从0开始的RNA-seq教程

本次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

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

这里请注意一点，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

4.reads计数

featurecounts 计数。

featureCounts -g transcript_id -Q 20 -p -B -C -T 4 -a corrected.gtf -o feature_counts T9_sorted.bam T10_sorted.bam T11_sorted.bam T12_sorted.bam  > feature_count.log   ###Q 进行质量过滤（hisat2的比对文件不要选择此选项）-p 双端测序 -B -C 只计算两条序列都比对上同一位置的序列  -T 线程

这里要说一下，目前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