科研 | NC：使用iDEA方法对单细胞转录组数据进行差异表达和基因富集分析
编译：夕夕，编辑：十九、江舜尧。
原创微文，欢迎转发转载。
论文ID
原名：Integrative differential expression and gene set enrichment analysis using summary statistics for scRNA-seq studies
译名：使用iDEA方法对单细胞转录组数据进行差异表达和基因富集分析
期刊：Nature communication
IF: 12.353
发表时间：2020.3.27
通讯作者：周翔
通讯作者单位: Center for Statistical Genetics, University of Michigan, Ann Arbor, MI 48109, USA.
DOI号： https://doi.org/10.1038/s41467-020-15298-6
实验设计
结果
1. 方法总结和仿真设计
IDEA方法的示意图见图1。简单来说，iDEA需要基因水平的总结统计数据作为输入文件。可以使用现有的差异表达分析的方法获得IDEA的输入文件。除了差异表达分析的总结统计信息外，iDEA还需要预编译的基因集。对于人类的数据，作者从七个现有的基因集和数据库中汇总了12033个基因集。对于小鼠的数据，作者从GO数据库中汇总了2851个基因集。对于这些输入文件，iDEA一次检查一个基因集，使用期望最大化算法进行推断并使用Louis方法计算校准p值，以检查该基因集是否包含差异表达基因。iDEA是一个开源的R包，可通过www.xzlab.org/software.html获得。
作者激进型模拟分析以评估iDEA在GSE分析和DE分析中的有效性。简单来说，作者通过一个真实的单细胞转录组的数据集推断出的参数构建零截尾负二项分布模拟了174个细胞上的10000个基因的零膨胀数据。模拟数据与真实的scRNA-seq数据具有相似的特征。在模拟基因中，其中有一定比例的基因属于一个基因集，作者将次百分比称为基因集覆盖率（CR）。在无效模拟的GSE分析中，将10000个基因中的每个基因都随机分配为DE基因的概率为exp(γ0)/(1+exp(γ0))。在GSE分析的替代模拟中，第j个基因被随机分配为DE基因的概率为exp(γ0+αjγ1)/(1+exp(γ0+αjγ1 ))。作者主要在γ0=-2和CR=10%的标准下进行模拟分析，并探索γ0，γ1和CR的不同组合以创建各种模拟情景。
作者比较了使用iDEA的GSE分析结果和常规的GSE方法的分析结果。作者发现，在不同的模拟情况下，iDEA会得到较好的p值（图2）。基因组控制因子（λgc），定义为经测试统计得到的经验中位数和期望中位数的比值。在其他方法中，例如fGSEA,PAGE和GSEA都产生了校准的p值，但是PAGE的p值稍微偏离对角线（图2d）。相反的，只有当CR值很低时，CAMERA得到的p值才会被校准。在CR较大时，CAMERA的p值不准确可能是因为CAMERA方法使用渐进法。作者注意到γ0值极低和CR极低时，所有方法得到的p值都偏离预期。在这些极端参数组合下，iDEA在Ⅰ类错误控制方面的次优性能可能是由于在你和罕见且不平衡的事件数据时遇到了潜在的参数可辨识性问题。其他方法的次优性能可能是由于这些方法中获取p值的渐近正态逼近不再准确。
图2 GSE分析中iDEA得到较好的p值
除了类型I错误控制之外，作者发现iDEA在一系列可选场景中比其他GSE方法更强大（图3a，c）。由于不同的方法具有不同的类型I错误控制，因此为了公平比较，作者以5%的FDR计算功效。在γ1=0.5和CR=10%的情况下，作者发现iDEA的功效为98%。相反的，fGSEA,GAMERA,PAGE和GSEA的功效为26%，0%，8%和26%（图3a）。iDEA和其他方法的功效都随着γ1的增加以及CR的增加而增加。在使用iDEA进行GSE分析时，通过ROC曲线观察到其AUC更高。iDEA优于现有的GSE分析方法可能是由于其他的GSE分析方法无法检测基因组中DE检验统计量与外部检验之间的差异。
图3 iDEA在GSE和DE分析中表现较好
作者发现，进行DE分析时，iDEA可以提高分析能力与其输入文件的来源无关（图3b和3d）。由于iDEA的模型需要输入文件中p值较好，因此，本文作者主要采用zingeR的分析结果作为输入文件。在分析中，作者发现iDEA得到的功能增益主要是由于DE和GSE分析的联合建模，而不是其对所有基因的联合建模。重要的是，iDEA可以在一系列模拟情况下产生合理的FDR值。ROC曲线也产生一致的结果，在进行DE分析时iDEA得到的AUC更大。除了直接检查DE分析能力外，作者还使用了Jaccard指数检查了不同DE分析方法的结果异质性。大概是由于iDEA的功能增益，作者发现iDEA可通过不同方法获得的DE基因列表来提高DE结果的一致性。
作者使用iDEA方法分析了3个公共的单细胞转录组数据集。第一个单细胞转录组数据集来自5种细胞类型的15280个基因组成。作者进行了DE和GSE分析。
作者首先应用IDEA方法和其他GSE方法在已构建好的11474个人类基因集数据库中检测显著富集的基因集（图4a）。作者通过对每个基因集的标签扫描10次构建empirical null p值分布。与模拟结果一致，作者发现iDEA分析，fGSEA分析，PAGE分析和GSEA分析的p值较好，但CAMERA的p值较差（图4b）。对于每种方法，作者都依赖于empirical null p值分布，基于FDR的经验值计算检测富集基因集的能力。与模拟结果一致，与其他GSE分析方法相比，iDEA鉴定出更丰富的基因集（图4c）。例如，经验FDR为5%，iDEA鉴定到2106个显著富集的基因，比fGSEA鉴定到的基因多20.9%。相反的，CAMERA,PAGE和GSEA分别鉴定到了537，1328和1079个基因。值得注意的是，除了统计结果外，iDEA鉴定到的基因还与胚胎发育密切相关。例如Wnt信号通路，TGF-beta受体信号通路和相关的GO条目，例如GO:0048514,GO:0001944和GO:0007492。为了量化不同GSE分析方法鉴定到的基因集的生物学意义，作者通过在PubMed中搜索相关文献，以一种无偏见的方法量化了基因集与胚胎细胞发育之间的相关性。实际上，在通过不同方法鉴定出的前50名富集基因集中，iDEA鉴定出的与胚胎细胞发育相关的基因最多，其次为fGSEA,CAMERA,PAGE和GSEA。iDEA检测到的与胚胎发育有关的基因数量最多，为iDEA可用于GSE分析提供了证据。
作者接下来使用iDEA进行DE分析，与模拟结果一致，iDEA鉴定到的DE基因比zingeR鉴定到的DE基因多。例如，设置经验FDR值1%，iDEA鉴定到2753个DE基因，比zingeR多64%（图4d）。iDEA鉴定到的50个DE基因明显划分为两种细胞类型：DECs和ECs（图4f）。重要的是，iDEA鉴定到内胚层细胞分化直接相关的1119个基因，而zingeR仅检测到706个。由于iDEA检测到的与内胚层细胞分化相关的DE基因数量最多，也为iDEA可用于DE分析提供了证据。iDEA分析鉴定到了涉及细胞分化过程的重要的DE基因，但zingeR没有鉴定到。包括SMAD3,GATA3,TGFBR1,WNT7B,HAND1,CCND1和HEY2。其中SMAD3对激活必要的转录调控网络以指导形成内胚层至关重要，GATA对于血管内皮小的信号传导通路必不可少，TGFBR1在激活SMAD2和SMAD3中起到重要重要，WNT7B对于配体-受体系统是必需的。最后，与模拟结果一样，iDEA改善了DE分析结果的一致性，在FDR设置为1%时，通过三种DE分析方法获得的DE基因的Jaccard指数仅为0.1，而使用iDEA分析后，Jaccard指数提高到0.25。
第二个单细胞转录组数据集是有622个小鼠神经元细胞的13598个基因组成，共划分为11个细胞类型。本文中作者主要研究NP1细胞类型和其余十种细胞类型的比较结果。
作者首先使用iDEA对之前构建的2851个小鼠基因集进行了GSE分析（图5a）。与模拟结果一致，使用iDEA,fGSEA,PAGE和GSEA进行GSE分析得到的p值较好，而CAMERA分析得到的p值较差（图5b）。与其他分析方法相比，iDEA鉴定到了更丰富的基因集（图5c）。例如，在设置FDR值为5%时，iDEA鉴定到1268个富集基因集，是GSEA分析结果的5倍。相反的，fGSEA,CAMERA和PAGE分别鉴定到了236，134和205个富集基因集。值得注意的是，由iDEA鉴定的基因具有生物学意义（图5e）。富集程度在前1%的GO条目大部分与神经系统，神经元反应和神经元功能有关，例如神经元投射（GO:0043005）和神经元部分（GO:0097458）。鉴定到的其他富集基因集包括轴突（GO:0030424）,突触（GO:0045202）和离子转运（GO:0006811）在神经元功能和活性方面也起到重要作用。fGSEA和CAMERA方法都没有检测到上述基因集，虽然PAGE和GSEA都检测到了上述基因集，排名较为靠后。
作者使用iDEA进行DE分析。同样的iDEA鉴定到的DE基因比zingeR的多（图5d）。设置FDR的值为1%，iDEA鉴定到1103个DE基因，其数量比zingeR高11%。作者选择了iDEA鉴定到的50个基因，其可明显划分为两种类型（图5f）。iDEA鉴定到的许多NP1神经元marker基因在zingeR中没有鉴定到。这些marker基因包括MRGPRB5, STX1B, FAM167A, KLK8 和STK32A。MRGPRB5是在初级伤害性感觉神经元中表达的与Mas相关的基因。KLK8介导伤害感受神经元中PAR1依赖性信号响应。重要的是，iDEA在已知的100个NP1的DE基因中检测到了79个，而zingeR检测到75个，再次说明了iDEA的强大功能。最后，与模拟结果一致，iDEA还提高了DE分析结果的一致性即通过三种DE方法计算得到的Jaccard指数为0.14，而iDEA的Jaccard指数为0.17。
结论
本文作者主要提出了一种新的可用于单细胞转录组分析中DE和GSE分析的计算方法。iDEA以现有的单细胞转录组差异表达分析的结果统计信息作为输入，为富集基因集生成校准的p值，并提供强大的DE和GSE分析。iDEA可以对汇总统计信息进行建模以避免对单个级别的 scRNA-seq数据进行显式建模的需要，从而可以将iDEA与现有的DE分析工具配对使用，以快速适应各种scRNA-seq数据类型。此外，使用汇总统计信息可以使iDEA合理的利用计算资源，其计算复杂度相对于基因数量呈线性增长。作者已经通过分析已发表的scRNA-seq数据集证明了iDEA的优点。
更多推荐
1 科研 | PNAS：转录组学揭示急性和慢性饮酒对肝脏昼夜新陈代谢有不同的影响