转录组及其他RNA测序最低需要几个生物学重复？差异分析软件到底哪家强？ | 转录调控专题
RNA-seq现已取代表达谱芯片成为全基因组差异基因表达（DGE）实验的首选技术。在任何实验设计中，选择适当数量的生物学重复（指经过相同方式处理的相同样品）都是成本和精度之间的权衡。RNA-seq技术的兴起后随之而来便是诞生了多款用于差异分析的软件。许多已发表的高质量RNA-seq研究中，每组或每个条件下至少有不低于5次生物学重复，但由于每组设置的生物学重复数量过少会干扰到DGE软件精准识别差异表达基因的能力。例如联川生物用户近期发表的部分转录组大样本研究中，生物学重复普遍不低于5个。（详情点击：大样本用户案例Cell Rep：600例转录组揭示禁食激活蛋白酶体开关| 转录调控专题）。其次，目前还没有文献证实，差异基因的变化倍数越大越能说明其发挥更重要的生物学功能。相反，有一些文献证实部分差异变化倍数不大的显著差异基因也能发挥重要的生物学功能。因此，如果研究人员比较关注差异倍数变化不大，但行使重要生物学功能的差异基因时，建议每组最低生物学重复样本数≥12，若受到经费限制或样本来源限制，最低生物学重复数也要≥6。
针对上述问题，2016年英国邓迪大学的Geoffrey J. Barton教授在冷泉港出版的专业学术期刊《RNA》上发表了一篇文章专门进行了评估。下面就让我们一起来详细了解一下：
论文标题：How many biological replicates are needed in an RNA-seq experiment and which differential expression tool should you use?
刊登日期：2016年06月
发表期刊：Journal of Extracellular Vesicles
影响因子：5.636
研究机构：英国邓迪大学生命科学学院计算生物学系
技术手段：RNA-seq、11款差异基因表达鉴定软件（baySeq、Cuffdiff、DEGseq、DESeq、DESeq2、EBSeq、edgeR、limma、NOISeq、PoissonSeq、SAMSeq）
原文链接：https://rnajournal.cshlp.org/content/22/6/839
这篇文章针对生物学重复数和正确识别差异表达基因的能力之间的确切关系进行了探索。作者分别对野生型和snf2突变型酵母样品的48个生物学重复进行RNA-seq检测。质控后，野生型样品保留42个生物学重复，突变型样品保留44个生物学重复。在控制假阳性率为0.05的标准下，用所有可用的生物学重复计算出的差异基因被定义为该DGE工具的“金标准”差异基因集。作者通过高度重复的RNA-seq实验评估了11个DGE工具的性能，最后得出结论：性能最好的是edgeR [exact], DESeq2和limma。该研究提供了以下两方面的建议：
1. 未来的RNA-seq实验需要多少次重复才能最大限度地提高DGE识别和量化的敏感性和准确性？
2. 在给定重复数的RNA-seq实验中，使用下面列表中哪种DGE工具来检测差异表达基因最合适？
表1 RNA-seq差异基因表达工具及统计检验方法
主要研究内容：
01
特定工具的“金标准”
本文主要通过比较来自重复子集的DGE结果与11个工具特定的“金标准”DGE集，来评估每个DGE工具的性能，并且该性能可作为重复数和表达差异倍数的函数。这些工具特定的“金标准”是通过在分析工具中运行全部clean data，然后对每个基因的读取计数的测量值进行计算得到的。此外还将P值或FDR≤0.05的差异表达基因标记为“显著差异表达”（SDE）。这些“金标准”运行后，来自酿酒酵母（Saccharomyces cerevisiae）被注释的7126个基因中通常有60%-75%被鉴定为SDE。对各差异基因鉴定工具特定的“金标准”进行定义后，从clean重复样本中抽取的重复子集上利用DGE算法进行迭代运行。对于这11种工具，在每个条件下以i=100次迭代和nr=2, …, 40次迭代(nr: number of biological replicates)进行自举运行（cuffdiff的运行速度明显比其他工具慢，所以这个工具的迭代次数减少到i=30）。对于一个给定的nr值，所有自举迭代中计算出的平均对数转换成了倍数变化[log2(FC)]、P值或FDR，并代表每个基因的测量值。此外，作者通过比较每个自举的SDE基因和各DGE工具“金标准”的SDE基因，将真阳性、真阴性、假阳性和假阴性率（TPR、TNR、FPR、FNR）计算为四个任意倍变阈值（|log2（FC）| = T∈ {0, 0.3, 1, 2}）。02
工具性能
图1
图1展示了利用edgeR运行数据后的一个示例。图1A展示了在控制相同的假阳性率水平下，不同的生物学重复鉴定出的差异基因数目。作者从所有生物学重复中随机抽取了2组、3组…40组生物学重复，分别计算差异基因，结果发现：1、差异基因的数目与生物学重复数量正相关；
2、差异基因数目的稳定性与生物学重复数量正相关；
3、生物重复较少时，不同的抽样导致的差异基因数目波动较大；
4、生物重复较多时，检测出的差异基因数目受抽样影响较小，体现在箱线图数据分布更集中。
图1B和1C展示了不同生物学重复与鉴定的差异基因真阳性率的关系。不同实线代表不同的差异基因倍数变化（T=|Log2(FC)|）条件下和不同的生物学重复条件下的真阳性率。虚线代表假阳性率。而虚线近乎一条直线，说明edgeR对假阳性率的控制还是比较好的，且不受生物重复数影响。如果倍变阈值比较高，比如4倍差异（T=2）时，较低的重复数即可获得较高的真阳性率。当倍变阈值较低（T=0）时，真阳性率受生物学重复影响较大，表现为生物学重复越少，真阳性率越低。值得注意的是：常规倍变筛选标准2倍差异(T=1)时，需要20个生物重复才能达到与4倍差异相同的真阳性率。生物重复越多，真阳性率越高，并受倍变阈值影响越小。
图1D也展示了真阳性、真阴性 (非“金标准”差异基因定义为真阴性（无差异）基因)、假阳性和假阴性的基因数目随生物学重复数的变化而发生变化。总的来说，生物学重复越多，漏掉的差异基因（假阴性基因）就越少（图1D）。
图2
图2总结了该研究中所有11种工具的性能与生物学重复数（nr）和倍变阈值（T）的关系。除NOISeq和DEGSeq外，所有工具的3次重复和0倍变阈值（nr = 3 , T = 0）自举子集的阳性率仅为20%～40%，表明无论用什么工具分析数据，在这么少的重复中，这些实验都无法识别大部分的差异表达基因（图2A）。
DEGSeq和NOISeq都显示出较高的真阳性率，但同时也具有较高的假阳性率（DEGSeq：FPR=17％，NOISeq：FPR=9％）。除去DEGSeq，其余所有工具的真阳性率都是较低波动变化（图2C）。
T=2时，这些工具显示的真阳性率均≥85%，除cuffdiff外，FPR基本为0（图2E）。高倍数变化（T=2）的SDE基因队列中，真阳性率对重复数基本不敏感。在T=2的情况下，无论哪种工具，生物学重复的数目对差异基因的影响均较小（图1B，2F）。
而重复数的增加对具有较小倍数变化的基因的检测率有很大的影响。除了DEGseq，倍变阈值的降低都会使其它所有工具的真阳性率降低，并与重复数无关（图2A-D）。而降低倍变阈值导致的真阳性率降低可以通过增加重复数来恢复。
通过上述结果我们可以得出，性能最好的工具为DESeq、DESeq2、EBSeq、edgeR [exact]和limma，它们成功地控制了假阳性率。无论倍变阈值或重复数如何，都能保持在接近或低于5%的水平，再次强调了增加重复数的主要作用是提高这些工具的敏感性，将假阴性转化为真阳性。
接下来，作者介绍了一种更可靠的方法来探测DGE工具的FPR性能。
03
工具的一致性与高重复数据
图3
作者对DGE工具的性能进行进一步测试，并假设：在生物学重复足够时，这些工具收敛于数据中真正的潜在差异表达信号。这一假设是通过对每个工具的 "金标准 "的DGE测量值和应用于相同数据的另外5个简单统计测试（t-test、long t-test、Mann-Whitney、Permutation、Bootstrap）的结果进行聚类来检验的（关于统计测试的详细描述，见论文中的材料与方法）。对于每个工具或测试，作者构建了一个由1和0组成的含有7126个元素的长向量，代表注释中的每个基因是否被该工具鉴定为SDE（调整后的P值或FDR阈值≤0.05）。每个工具或测试的向量按基因id排序，然后用R软件包pvclust按完全连锁的相关距离进行分层聚类。其中，为聚类中每个分支计算出的近似无偏的P值百分比（AU%-图3，括号内的数值）表明了每个分支对抽样误差的稳健程度。用标准的统计学检验后，有三个广泛使用的工具（DESeq2，edgeR [exact] 和 limma）聚集在一起（图3，群组3），cuffdiff，DESeq和EBSeq也紧紧地聚在一起（图3，群组4），此外，NOISeq、DEGSeq、baySeq和edgeR [GLM] 也形成了一个独特的独立群组（图3，群组2），表明这些工具得出的结果与群组1中的工具大不相同。
测试工具误报率
图4
RNA-seq差异表达工具最重要的性能指标可能还是其错误检测率。该研究中，大量的重复足以允许对每个工具的假阳性率进行简单测试。
作者从WT条件下随机选择两组nr重复，并假设两组之间没有表达差异，根据定义，每一个被鉴定为SDE的基因都是一个假阳性。此外，对于工具的每一次自举运行，都要计算总基因组中被鉴定为SDE的这一部分。图4显示了每个DGE工具的假阳性比例与重复数的关系。结果表明，即使有大量的重复，DEGSeq、NOISeq和SAMSeq的表现也很差。特别是DEGSeq，其假阳性率的表现更差，每一次自举迭代都会将>5%的基因识别为假阳性，无论重复多少次，假阳性率的中位数都在50%以下。此外，大约10%的cuffdiff、PoissonSeq和40%的SAMSeq自举迭代结果中将>5%的基因识别为了假阳性，这表明这些工具也没有很好地控制其假阳性率。BaySeq、DESeq和EBSeq在这项测试中表现非常好，edgeR、DESeq2和limma也表现良好。
相关阅读
植物RNA多组学研究的正确打开方式全部在这里了，套餐组合任意挑选 | 转录调控专题用户文章-PBJ：m6A测序揭示阅读蛋白MhYTP2调控苹果白粉病抗性的分子机制-转录调控专题
干货：5千字解密转录组发高分常见应用场景（含推荐样本数） | 转录调控专题
所见即所得，绘图高规格联川云平台，让科研更自由