从0开始的RNA-seq教程

从0开始的RNA-seq教程

五香可达鸭

2020.10.26 10:35:56字数 1,194阅读 135

本次RNA-seq分析目标明确,得到基因的表达矩阵即可,不涉及其他分析。

样本:植物的叶和茎的转录组;测序方法:双端测序;有参考基因组文件

1.数据质控
拿到测序数据,第一步进行md5检测,确定文件没有损坏。然后使用fastqc对测序数据进行质控,multiqc可以聚合多个qc结果进行展示。双端测序文件一般是两个,命名时一般会在末尾加上1,2加以区分。
从 SRA下载了raw data,质控很多人用fastqc,这里为了方便,我使用了fastp,直接生成过滤后的文件,省去了过滤的步骤,并且速度很快。

#SRA数据提取fastq文件,详细参数自行搜索
fastq-dump --gzip --split-e SRR_ID

#fastqc 命令
fastqc -t 8 -o out_path sample1_1.fq sample1_2.fq

#fastp命令  输入文件、输出文件、输入文件、输出文件、线程数
fastp -i  A1.fq.gz  -o fastp_A1.fq.gz -I A2.fq.gz -O fastp_A2.fq.gz -w 4
12345678

2.参考基因组和基因注释文件
我的样本特殊,有参考基因组和注释文件,但是注释文件不完善,只进行了基因的结构预测,并没有给gene_Id,所以要在后来的处理中花费额外的时间去做(第一个坑)。
从ncbi下载组装注释好的基因组。

3.序列比对

目的:这一步的目的是把测序的reads比对到参考基因组上。

RNA-Seq数据分析分为很多种,比如说找差异表达基因或寻找新的可变剪切。如果找差异表达基因单纯只需要确定不同的read计数就行的话,我们可以用bowtie, bwa这类比对工具,或者是salmon这类align-free工具,并且后者的速度更快。

但是如果你需要找到新的isoform,或者RNA的可变剪切,看看外显子使用差异的话,你就需要TopHat, HISAT2或者是STAR这类工具用于找到剪切位点。因为RNA-Seq不同于DNA-Seq,DNA在转录成mRNA的时候会把内含子部分去掉。所以mRNA反转的cDNA如果比对不到参考序列,会被分开,重新比对一次,判断中间是否有内含子。
链接:https://www.jianshu.com/p/681e02e7f9af

就唯一比对而言,STAR是三者最佳的,主要是因为它不会像TopHat和HISAT2一样在PE比对不上的情况还强行把SE也比对到基因组上。而且在处理较长的read和较短read的不同情况,STAR的稳定性也是最佳的。
就速度而言,HISAT2比STAR和TopHat2平均快上2.5~100倍
链接:https://www.jianshu.com/p/681e02e7f9af

工具:目前比较推荐的是hisat2,hisat2正确率高,当然总数量会降低,二类错误率低了,一类错误率就会高。STAR好像也行,但是我还是用了使用比较多的hisat2

正式开始:
①索引,为了提高比对效率,通过BWT算法对基因组建立索引去进行比对。有现成的就下载现成的,没有现成的就自己建一个索引(我用的服务器比较富裕,所以时间也没有很久,十几分钟吧)。hisat2快速上手教程

 # 其实hisat2-buld在运行的时候也会自己寻找exons和splice_sites,但是先做的目的是为了提高运行效率
 extract_exons.py gencode.v26lift37.annotation.sorted.gtf > hg19.exons.gtf 
 extract_splice_sites.py gencode.v26lift37.annotation.gtf > hg19.splice_sites.gtf 
 # 建立index, 必须选项是基因组所在文件路径和输出的前缀
 hisat2-build --ss hg19.splice_sites.gtf --exon hg19.exons.gtf genome/hg19/hg19.fa hg19
12345

这里请注意一点,hisat2使用gtf文件生成exon和ss文件,gff文件不可以,所以要先把gff文件使用gffread转化成gtf格式。

②开始比对
③sam文件转换为bam文件并进行排序,建立索引
④bam文件质控(因为后续输入文件可以用sam文件,所以省去了转换格式那一步)


hisat2 -p 4 -x ../index/Rosa -1 A_1.fastq.gz -2 A_2.fastq.gz –S A.sam

123

4.reads计数

featurecounts 计数。


1

这里要说一下,目前FPKM和RPKM不被用来做差异分析,TPM和TMM标准化用来做差异分析的比较多,且DEseq的输入文件是标准化之前的,所以推荐featurecounts计数。

5.基因差异表达分析
6.富集分析

写在最后:
推荐的两篇综述文献(虽然很难读):
①Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis
②A survey of best practices for RNA-seq data analysis

1人点赞

随笔

“小礼物走一走,来简书关注我”

还没有人赞赏,支持一下

五香可达鸭

总资产0.249 (约0.02元)共写了1194字获得1个赞共3个粉丝

点击数:0

发表评论

邮箱地址不会被公开。 必填项已用*标注