2023-03-11

#RNA-seq实战代码

#正常情况下公司给我们的数据就直接是fastq格式文件

#使用fastqc查看测序数据质量,就会得到*.html文件

fastqc *.fastq.gz

#在html文件所在目录下将*.html文件拷入一个单独的文件夹:cp *.html 新建的单独的文件夹所在的目录

#将linux中的文件复制到桌面

cp -r 要复制的文件所在目录 /mnt/c/Users/happy/Desktop/

#将所有的质量报告合并:对html文件所在目录跑一下multiqc:multiqc 新建的单独的文件夹所在的目录

multiqc 新建的单独的文件夹所在的目录   (要在rnaseq3.7环境中跑这个代码)

#使用trim_galore 过滤fastq文件中质量不好的reads和remove接头序列

ls /home/happy/project/airway/raw/*_1.fastq.gz >1

ls /home/happy/project/airway/raw/*_2.fastq.gz >2

#生成只包含*—1.fastq.gz和*—1.fastq.gz文件的目录

paste 1 2 >config

#过滤文件打算存储的位置

dir='/home/happy/project/airway/clean/'

#循环如下:

cat config |while read id;

do

arr=($id)

fq1=${arr[0]}

fq2=${arr[1]}

nohup trim_galore -q 25 --phred33 --stringency 3 --length 36 -e 0.1 --paired $fq1 $fq2 --gzip -0 $dir &

done

#3、再次使用fastqc查看过滤后的文件

#比对,转录组比对主要是使用hisat2、subjunc、star;基因组主要是使用bwa和bowtie2

#比对首先要下载参考基因组构建索引index,也可以直接在官网下载index,人类和小鼠的index有现成的

#方法一:直接在hisat2官网下载index,选择合适的物种

#hisat2官网:http://daehwankimlab.github.io/hisat2/download/

#构建存储索引的文件夹

mkdir index

cd index

wget 复制链接

#解压index文件

tar -zxvf *.tar.gz

#删除压缩包

rm -rf *.tar.gz

#方法二:自己下载参考基因组和注释文件构建index

#下载基因组文件

wget 参考基因组链接

tar -zxvf  *.tar.gz

#下载注释文件

wget 注释文件链接

gunzip *.gtf.gz

extract_exons.py hg19.ncbiRefSeq.gtf > genome.exon        其中hg19.ncbiRefSeq.gtf是解压后的注释文件

extract_splice_sites.py hg19.ncbiRefSeq.gtf > genome.ss

#建立index

hisat2-build hg19.fa --ss genome.ss --exon genome.exon genome_tran

#4.比对

hisat2 --dta -t -x /data/RNAseq/mm10/genome -1 /data/RNAseq/cleandata/trim_galoredata/CK-4_1_val_1.fq.gz -2 /data/RNAseq/cleandata/trim_galoredata/CK-4_2_val_2.fq.gz -S cleandata/hisat2_mm10data/CK4.sam

#/data/RNAseq/mm10/ 是参考基因组文件所在目录,genome这个不是文件夹,是index文件的前缀,mm10文件下并没有这个文件,如果不加genome就会发生报错

#-1和-2分别表示双端测序的1个文件,后面跟的是文件路径

#-S 后面跟的是输出数据所在文件夹

#-x 后面接的是索引文件

#命令没有问题的话,出现下面提示表示程序正在执行,就去干其他的,等执行结束再往下:

Time loading forward index: 00:02:51

Time loading reference: 00:00:16

#写循环批量处理(for  in 循环和while循环各写一个)

#for in 循环,创建脚本文件

vim hisat2_mm10_batch.sh

#输入下面内容:

for i in CK-7 CK-8 HGJ-10 HGJ-6 HGJ-9

do

hisat2 --dta -t -p 8 -x /data/RNAseq/mm10/genome

-1 /data/RNAseq/cleandata/trim_galoredata/"$i"_1_val_1.fq.gz

-2 /data/RNAseq/cleandata/trim_galoredata/"$i"_2_val_2.fq.gz

-S /data/RNAseq/cleandata/hisat2_mm10data/"$i".sam;done

#保存脚本,再运行脚本

bash hisat2_mm10_batch.sh

#while 循环

#先创建存储了SRR号的文件。例如sra.txt

less sra.txt|while read id ;

do

hisat2 --dta -t -p 8 -x /data/RNAseq/mm10/genome

-1 /data/RNAseq/cleandata/trim_galoredata/"$i"_1_val_1.fq.gz

-2 /data/RNAseq/cleandata/trim_galoredata/"$i"_2_val_2.fq.gz

-S /data/RNAseq/cleandata/hisat2_mm10data/"$i".sam;done

#subjunc比对会直接输出bam文件,就不用后续的sam转bam文件

subjunc -T 5 -i /home/happy/project/airway/reference/mm10/genome -r ${id}_1_val_1.fq.gz -R ${id}_2_val_2.fa.gz -o ${id}.subjunc.bam

#5.samtools将sam文件转bam文件,samtools可以进行格式转换、排序、索引

#第一步:转换

#SRR.txt为存储了SRR号的文件

单样本处理:samtools view -@8 -b SRR1039509.sam > SRR1039509.bam

多样本处理:less /home/happy/project/airway/hisat2/SRR.txt|while read id ;do samtools view -@8 -b {$id}.sam > {$id}.bam done

#第二步:排序 samtools sort

#其中 -l是(-L的小写)

单样本处理:samtool sort -@ 8 -l -o SRR1039508.bam.sort SRR1039508.bam

多样本批量处理:ls /home/happy/project/airway/SRR.txt |while read id ;do samtools sort -@ 8 -l {$id}.bam -o {$id}.bam.sort;done

#第三步:建立索引(index命令)在SRR1039508.bam.sort所在文件夹下操作

单样本处理:samtools index  SRR1039508.bam.sort SRR1039508.bam.index

查看sam文件直接用less,查看bam文件用samtools view SRR1039508.bam |less -S(分页查看) 可以看出来是按pos排序还是按name排序,这里是按pos(染色体位置)

#第四步:查看reads比对情况(flagstat命令),输出.falgstat文件

samtools flagstat -@ 8 SRR1039508.bam.sort > SRR1039508.bam.flagstat

#第五步:igv可视化比对结果

先将bam文件转为bw文件,要先构建索引 samtools index SRR1039508.bam.sort

再使用bamCoverage转换格式:bamCoverage -b SRR1039508.bam.sort -o SRR1039508.sort.bw

再将bw文件复制到桌面:cp *.bw /mnt/c/Users/happy/Desktop/    (bw文件要指定文件所在目录)

igv里面选择file 载入bw文件

#6.计数(htseq、featureCounts )通常用featureCounts计数

htseq计数代码(很慢)

htseq-count -f bam -r pos /home/happy/project/airway/hisat2data/SRR1039508.bam.sort /home/happy/project/airway/reference/gencode.v10.annotation.gtf.gz

#-f 指定bam文件格式

#-r 指定排序方式,这里是按位置position

#/home/happy/project/airway/hisat2data/SRR1039508.bam.sort  进行计数的bam文件

#/home/happy/project/airway/reference/gencode.v10.annotation.gtf.gz  注释文件所在目录

featureCounts计数(很快):参考基因组和GTF/GFF/SAF注释文件来自同一个网站,同一个版本

featureCounts 需要两个输入文件:比对产生的BAM/ SAM文件(一般会用bam文件,因为所占空间小);区间注释文件(GTF格式, SAF格式)

代码1:featureCounts -T 10 -a ../reference/gencode.v10.annotation.gtf.gz -o count.txt -p -B -C -f -t exon -g gene_id /home/happy/project/airway/hisat2data/*sort

代码2:featureCounts -T 10 -a ../reference/gencode.v10.annotation.gtf.gz -o count1.txt -p -B -C -t exon -g gene_name /home/happy/project/airway/hisat2data/*sort

-g  后面也可以接gene_name

#如果代码1比对率不高,说明不是全外显子测序,则用代码2

-T 指定线程数    -f 如果-f被设置,那将会统计feature层面的数据,如exon-level,否则会统计meta-feature层面的数据,如gene-levels(feature数目指的是exon数目   meta-feature 指的是基因数目)

-a 参考GTF文件名     -F 参考文件的格式,一般为GTF/SAF ,默认为GTF    -p 双端测序    -B  在-p选择的条件下,只有两端reads都被比对上的fragment才会被统计

-t 设置feature-type,-t指定的必须是gtf中有的feature,同时read只有落到这些feature上才会被统计到,默认是“exon”

-g 默认为gene_id    -o 输出文件的名字,可输出被统计数据的txt文本及summary文本

/home/happy/project/airway/hisat2data/*sort  被统计的文件所在目录及名称

#featureCounts查看及结果解析:less -SN  count.txt  (-SN可以使查看是排版更好)

Geneid:基因的ensemble基因号;   Chr:多个feature所在的染色体编号;  Start:多个feature起始位点,与前面一一对应;  End:多个feature终止位点,与前面一一对应

Strand:正负链   Length:基因长度  sampleID:一列代表一个样本,数值表示比对到该基因上的read数目

feature数目指的是exon数目   meta-feature 指的是基因数目

#7.表达矩阵探索

multiqc 要在3.7版本的python才不会报错,所以要创建pyhton3.7环境下下载multiqc :conda create --name rnaseq3.7 python=3.7

conda activate rnaseq3.7

conda install -c bioconda -c conda-forge multiqc

multiqc count1.txt.summary    就会得到一个html文件  复制到电脑桌面看

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 215,723评论 6 498
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 92,003评论 3 391
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 161,512评论 0 351
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 57,825评论 1 290
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 66,874评论 6 388
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 50,841评论 1 295
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,812评论 3 416
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,582评论 0 271
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,033评论 1 308
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,309评论 2 331
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,450评论 1 345
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,158评论 5 341
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,789评论 3 325
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,409评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,609评论 1 268
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,440评论 2 368
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,357评论 2 352

推荐阅读更多精彩内容

  • 生信学习笔记 linux部分功能 查看文件夹 工具 选项 可以设置鼠标功能 可以设置右键粘贴 双击这个窗口可以再打...
    Vikenn阅读 1,150评论 1 4
  • 一、测序数据的准备 新建工作目录,RNASeq-analysis mkdirRNASeq-analysiscdRN...
    胡小兵plus阅读 29,329评论 1 36
  • control+左右键:光标在单词之间跳转control+a:光标跳到首control+e:光标跳到尾more 查...
    王琪_738c阅读 660评论 0 0
  • 一、文章数据下载 安装miniconda sudo apt-get install wgetwget https:...
    黑白配_b74c阅读 158评论 0 1
  • RNAseq实际操作(实战) 首先声明,虽然是实战,但是其实是学习笔记而已,初学,参考了大量大神的博客和帖子,还有...
    zd200572阅读 4,397评论 0 29