A survey of best practices for RNA-seq data analysis

image.png

数据质控及过滤

fastqc , multiqc, trimmomatic, cutadapt ,trim-galore

FastQC数据质控

cutadapt数据过滤

去文库构建时加入的接头、reads两端低质量碱基、含N多的碱基、去短的reads。

参考基因组比对

star, hisat2, bowtie2, tophat, bwa, subread
利用转录组比对软件，将上一步得到的clean reads与基因组比对。

「比对的作用」：1、评估文库构建的质量（文库构建时随机打断，如果文库优质，打断随机性好，比对时会发现reads均匀分布在基因组上）；2、评估mRNA测序是否有意义（理论上reads应主要比对到CDS_Exons ，CDS_Exons占据mRNA绝大部分）。

bowtie2作用

比对之前需要先对参考基因组建index
推荐使用 Hisat2或 STAR。STAR需要更大存，运行时间也更长。准确性相差不大。
因为可变剪切，RNA-Seq reads由于不包含内含子，所以来自外显子边界处的 reads被重新 mapping回基因组时，会被中间的内含子分开，这种情况叫做splice-alignment。

转录本组装

bam文件只有序列比对的情况，而注释信息会标注哪里是基因

「组装原因」，由于二代测序读长限制，必须把mRNA打断成小片段进行测序，组装的目地就是利用生信方法重新拼接全长转录本。「主流组装软件」为「Stingtie」和「Cufflinks」，两种软件结合注释文件gff/gtf，将转录本构建出来

cufflinks作用

基因表达量分析

基因表达量即在当前研究条件下，一个细胞中、或者是一定摩尔量的 RNA中某个特定基因比对上了多少条转录本（绝对定量，需要知道细胞个数）或者比对上转录本的比例（相对定量），更直接一点讲就是比对上reads的count数；

Manuel Garber et al., Nature Methods, 2011

基因表达量与「基因长度」，「测序深度」成正相关：一个样本中，A基因越长，建库时被随机打断的片段越多，被测序的概率越大，比对到A基因的reads就越多；不同样本中，样本的测序深度越高，A基因被测到的次数越多，比对到A基因的 reads 就越多。所以，直接数reads来算表达量的方法是有问题的，有TPM/FPKM/RPKM三种标准化方法，将reads数除以基因长度，测序深度来校正二者对表达量影响。现在比较常用的是TPM。

软件

虽然认为TPM更准确，但是由于三者可相互转换，所以都在用。根据关注点不同，可以使用不同的软件组合：

关注已知和新转录本定量，可用Cufflinks或StringTie；
关注转录本水平定量，可用RSEM或eXpress直接将reads比对到参考转录本；
不经过比对的定量，节省计算资源，可用Sailfish或Salmon或quasi-mapping或kallisto；

基因差异表达分析

「转录组分析的重要目标」就是找case组和control组样本之间差异表达的基因；差异表达分析依赖上一步所获的各基因表达量；分析工具有很多，根据「依赖技术」可划分为：

count-based 方法，可用DESeq、limma和edgeR；
assembly-based方法，可用Cuffdiff和Ballgown；
alignment-free方法，sleuth；

根据「有无生物学重复」，无生物学重复可用DESeq，有生物学重复可用DESeq2；

差异基因富集分析

评估差异基因主要影响的生物学功能和通路。

KEGG

找差异表达基因主要显著影响了哪些「生化代谢途径和信号转导途径」。

GO（GENE ONTOLOGY）

找差异表达基因主要富集在哪些GO term（分子到生物过程，分三类：molecular function、cellular component、biological process）中，评估「差异表达基因与哪些生物学功能显著相关」，对生物学功能起上调还是下调作用。

基因融合分析

可变剪切分析

引用自直观理解RNA-seq分析

以前 bowtie2+samtools+cufflinks+deseq2，
现在 hisat2+stringtie+ballgown
其本质的工作流程其实是一样的，只不过使用的算法不同

实验验证

用qPCR进行实验验证RNA-seq的可靠性，目前还没有人会把RNA-seq作为最可靠的证据，一般都是要进行qPCR实验验证的

转录组分析流程