RNA-seq分析流程分析:RNA-seq数据分析流程主要由序列对比,表达矩阵构建,差异基因鉴定等三大步骤构成【1-8】。目前,已有大量的软件被开发用于分析流程。而上游分析可能对下游分析或/与结果有实质上的影响。
文章主要探讨三种大步骤中各软件相互组合使用对结果所带来的差异【21】。文章选用的数据集为典型和非典型人类单核细胞数据集,且该数据集已被多个独立实验室进行了分析,以此数据集作为参考理论上具有高度的可信度。另,文章使用的测试数据集为真实数据,测试数据与参考数据集具有高度相关性。
研究显示不同的工作流的效果具有显著差异,只要体现在召回值和精确度两点,大体上呈现反义线性关系。提示我们在选取相应工作流时可首先考虑自己实验对这两项参数的需求。
代码:https://github.com/cckim47/kimlab/tree/master/rnaseq.
1、构建参考数据集
比对典型与非典型样本的测序数据参数,如测序质量、总读段数量、序列比对后读段数量,结果显示无显著差异。下载相关表达矩阵后,作者使用log2处理需要标准化的数据,并且使用了Significance Analysis of Microarrays (SAM)及limma进行差异基因分析并且进行比对后发现结果具有高度相关性,最终选取数据的交叉结果作为参考数据集
2、工作流组合
3、不同步骤对结果的影响
首先比对三大步骤内部差异基因数量,结果提示差异基因分析工具对显著基因数量产生的影响最大;另,各差异分析软件的结果稳定性也具有差异;
其次,作者比对了个工作流结果的召回值( recall ,intersecting significant genes divided by total number of significant reference genes)和精确度( precision ,intersecting significant genes divided by total number of significant genes identified by RNA-Seq),结果提示差异基因分析软件对结果影响最大。
4、工作流异质性
无论是转录本还是基因表达层面,召回值与差异基因数量呈线性关系,而精确度则呈反义线性关系。
基因表达水平的召回值,使用 SAMseq软件的工作流最高;转录水平的则以使用 baySeq及 NBPSeq软件的工作流最高;
基因表达水平的精确度,使用NOISeqBIO软件的工作流最高;转录水平中则具有多项,其中最常用的是Ballgown及NOISeqBIO。
值得注意的是,我们常用的TopHat2+cufflinks+cuffdiff工作流显示最高的精确度及第二低的差异基因数量。
5、工作流性能权衡
由于计算内部的关系,我们需要在召回值及精确度之间权衡,方能找到适合自己的工作流。研究中大部分工作流的召回值及精确度之间存在着反义线性关系。这一点在转录本及基因表达层面都是成立的。
Ballgown、DESeq2、 limma + voom、limma + vs及 and SAMseq最为接近该趋势,而baySeq和EBseq最为偏离。
SAM-seq(一种非参数方法)是一种高性能的软件【3,16】,尤其是在有大量重复数据可用时可使结果较为接近潜在的分布规律,但它趋向于牺牲精确度从而提高召回率;NOISeqBIO则倾向于在大规模的样本中鉴定更少的差异基因【3】并具有低召回值;baySeq和EBseq统计学方法最为接近,它们都以负二项模型(underlying negative binomial model)去估算每个基因差异表达的后验概率【46,48】,EBseq处理后的数据由于降低精确度而召回值未回升而偏离反义线性模型(EBseq在大样本中容易产生假阳性数据【16】)baySeq在处理基因层面数据时表现出与EBseq相似的倾向,可能是由于它们的计算模式相同;Ballgown是在limma的基础上发展的,它们三种软件性能良好并且趋向于反义线性模型。区别在于Ballgown倾向于更高的精确度,而limma+voom/vst更倾向于召回值。
序列对比和表达矩阵构建软件的选用通常对结果没有特殊影响,且差异基因分析所造成的影响远大于它们。除了以下两点:其一,BitSeq(表达矩阵构建软件)在与基于负二项模型的差异基因分析软件(BaySeq,DESeq2, edgeR, and NBPSeq)连用时会鉴定出大量的差异表达基因;其二,使用STAR进行读段对比会使一些高性能工作流的精确度和召回值达到平衡状态??,尤其是与Ballgown连用时。
工作流的权衡选择依据数据的下游分析和应用中对第一类错误和第二类错误的容忍度。( 当假设检验拒绝了实际上成立的零假设时,所犯的错误称为第一类错误,其概率用α表示;当假设检验接受实际上不成立的零假设时,所犯的错误称为第二类错误,其概率用β表示。),也即考虑召回值和精确度
附:其他研究者使用STAR+TPMCalculator+ DGA software进行测试分析,结果提示可能是TPMCalculator造成精确度的提高
https://ftp.ncbi.nlm.nih.gov/pub/RNASeqWF/notebooks/00%20-%20Project%20Notes.html