《PLOS Computational Biology》:Not every gene is special: Modelling scale controls the false discovery rate when analysing high-throughput sequencing data
编辑推荐:
差异表达/丰度分析在使用高通量测序(high-throughput sequencing,HTS)的研究中十分常见;然而,当不同工具应用于同一数据集时,往往无法返回可比的结果。大多数工具采用标准化方法尝试校正计数数据中的技术变异。此前,研究人员证明了这些标准化
差异表达/丰度分析在使用高通量测序(high-throughput sequencing,HTS)的研究中十分常见;然而,当不同工具应用于同一数据集时,往往无法返回可比的结果。大多数工具采用标准化方法尝试校正计数数据中的技术变异。此前,研究人员证明了这些标准化方法常因对生物系统总体尺度(即大小)的错误假设而不适用。在本研究中,研究人员结合二项式稀疏化(binomial thinning)和样本分组置换,对11个RNA-seq及其他HTS数据集产生了100次分析迭代,其中预期5%的特征在组间存在显著差异。这使研究人员能够计算各次迭代的错误发现率(false discovery rate,FDR)和灵敏度。模拟表明,尺度误设导致若干常用工具对FDR的控制不佳,并且与直觉相反的是,随着组间模型化差异的增大,FDR反而升高。在ALDEx2或ALDEx3中实施尺度模型改善了高得不可接受的FDR;然而,在满意的FDR控制和高灵敏度之间存在固有的权衡——没有工具能同时提供两者。研究人员确定,增加尺度不确定性也增加了特征被报告为差异表达所需的组间最小差异。在不同类型的HTS数据中一致地观察到这一现象,且高度一致。关键的是,研究人员利用一个“真实世界”的非置换RNA-seq数据集分析,证明了后一种效应并非由稀疏化/置换方法导致。总体而言,该工作突出了所有研究人员在分析测序数据时可能无意中做出的灵敏度与FDR控制之间的选择,并为HTS数据分析中选择合适的尺度不确定性量提供了指导。
高通量测序(HTS)已广泛应用于生物医学研究,差异表达/丰度分析是常见目标。然而,不同分析工具应用于同一数据集时结果常不一致,原因在于大多数工具依赖标准化来校正技术变异,但这些标准化往往对生物系统总体尺度(即系统规模)作出错误假设。此前研究已证明,尺度误设会导致错误发现率(FDR)失控,但相关研究多局限于扩增子或宏基因组数据,或仅使用单一转录组数据集。因此,本研究旨在跨多种HTS数据类型系统评估尺度模型(scale model)对FDR和灵敏度的影响,并为选择合适的尺度不确定性提供指导。
研究人员使用了11个公开HTS数据集,涵盖bulk转录组(实验室来源与临床来源)、单细胞转录组、宏转录组和16S rRNA扩增子测序数据。通过seqgendiff包进行样本分组置换和基于二项式稀疏化的信号注入,使5%的特征被模拟为组间差异,从而在100次置换迭代中计算各工具的FDR和灵敏度。比较了ALDEx2(v1.4.1)、ALDEx3(v0.3.0)与DESeq2(v1.49.1)、edgeR(v4.7.2)、limma-voom(v3.65.1)的性能。尺度不确定性通过ALDEx的gamma(γ)参数从0.001到0.5递增建模。真实数据验证使用了未置换的乳腺癌(BRCA)数据集。
### 高错误发现率在尺度na?ve工具中普遍存在
通过置换和稀疏化模拟,发现DESeq2、edgeR和limma-voom在所有测试数据集中的FDR最高,常超过名义水平的两倍以上,而ALDEx2和ALDEx3即使不建模尺度也表现出更好的FDR控制。增加尺度不确定性(γ值)进一步降低FDR,当γ=0.1时,在模拟log2倍数变化(L2FC)达到约0.5之前几乎没有假阳性;当γ=0.4时,FDR低于5%甚至接近零。值得注意的是,FDR随着L2FC阈值的升高而增加,在最大阈值处可超过80%。
### 常用工具在灵敏度和FDR控制之间存在权衡
灵敏度分析显示,DESeq2、edgeR和limma在所有数据集中的灵敏度均高于两个ALDEx版本。增加尺度不确定性会显著降低灵敏度,γ≥0.4时ALDEx2在某些数据集中仅检出不到50%的真实阳性。ALDEx3在所有γ值下均优于ALDEx2,在γ=0.2时与DESeq2等工具相当,尤其是在bulk转录组数据中。单细胞和阴道宏转录组数据中ALDEx工具的灵敏度普遍偏低,且所有工具的灵敏度在L2FC≥0.5时趋于平台。
### 增加尺度不确定性提高了显著所需的最小倍数变化
以PD1免疫治疗数据集为例,γ值增加导致ALDEx2和ALDEx3报告的差异表达基因数减少。火山图和效应图显示,尺度模型通过动态截止值排除了两类特征:一类是L2FC大但校正后P值边缘的特征,另一类是L2FC极小但P值极小的特征。跨数据集分析表明,不同数据集的初始L2FC阈值差异与组内离散度相关;对于同一工具,最小L2FC与γ之间的关系近似恒定,线性回归斜率在ALDEx2平均为2.34,ALDEx3平均为1.59。
### 将尺度模型应用于真实数据集凸显尺度na?ve分析中的虚假发现
在未置换的BRCA真实数据中,尺度na?ve分析报告了约15,837(ALDEx2)和15,825(ALDEx3)个差异表达基因,而DESeq2、edgeR和limma报告更多。增加尺度不确定性后,γ=0.5时ALDEx2和ALDEx3分别仅报告2,777和4,437个基因。被排除的基因超过99%不属于实验验证的乳腺癌驱动基因。典型例子是肿瘤抑制基因PTEN,其L2FC仅为-0.0594,校正后P值为0.034,组内离散度较高,其显著差异由少数离群样本驱动,加入少量尺度不确定性后即不再显著。然而,过高的γ(如0.5)会使ERBB2和KRAS等已知癌症基因也丢失,说明γ值不宜过高。
讨论部分强调,常用工具(如DESeq2、edgeR)在灵敏度和FDR控制之间严重偏向灵敏度,导致假阳性过多。尺度模型能够在几乎不损失灵敏度的情况下显著控制FDR,且通过γ值可灵活调整。研究人员建议,当分析者希望为L2FC提供约0.5的缓冲时,ALDEx2使用γ=0.2,ALDEx3使用γ=0.3;若仅关注大于1的L2FC,则两种工具均可使用γ=0.4。总体而言,尺度模型为HTS数据分析提供了一种比传统双重P值/L2FC截止更可靠的FDR控制策略。
结论:研究人员总体证明了ALDEx2和ALDEx3所实现的尺度模型在高通量测序数据分析中的实用性及其出色的FDR控制能力。与其他研究一致,该工作强调常用RNA-seq分析工具返回的假阳性率高得不可接受,并指出所有研究人员在考虑如何分析数据时都需要意识到灵敏度与FDR控制之间固有的权衡。