Chen J, Swofford R, Johnson J, et al. A quantitative framework for characterizing the evolutionary history of mammalian gene expression. Genome Res. 2019;29(1):53-63. doi:10.1101/gr.237636.118

摘要

基因的进化史有助于预测其功能以及与表型特征的关系。尽管序列保守通常用于破译基因功能和评估医学相关性，但缺乏从比较表达数据进行功能推断的方法。在这里，我们使用来自17个哺乳动物物种的7个组织的RNA-seq来证明跨哺乳动物的表达进化是由 Ornstein-Uhlenbeck过程准确建模的，这是一种普遍提出的连续性状进化模型。我们应用此模型来识别中性、稳定和定向选择下的表达途径。我们进一步论证了该模型的新应用，以量化基因表达的稳定选择程度，参数化每个基因的最佳表达水平分布，并检测个体患者的表达数据中的有害表达水平。我们的工作为解释跨物种和疾病的表达数据提供了一个统计框架。

比较基因组学已经通过跨物种的进化模式识别和注释了功能性遗传元件（Rubin 等人 2000；Kellis 等人 2003；Siepel 等人 2005；Pollard 等人 2006；Lindblad-Toh 等人 2011）。当前的比较研究主要集中在基因组序列的分析上，依赖于从观察中性序列随时间线性发散而发展起来的成熟理论框架（Harris 1966；Lewontin 和 Hubby 1966；Kimura 1968；Jukes 和 King 1971）。这些方法允许检测在中性进化的无效模型下进化比预期更慢（例如，由于纯化选择）或更快（例如，由于正选择或放松的选择约束）的序列元件。

长期以来，人们普遍认为基因调控的分歧，表现为基因表达的表型变化，在进化中也起着关键作用（King and Wilson 1975；Wang et al. 1996；Pierce and Crawford 1997；Ferea et al. 1999；Fraser等人 2010）。基因表达的进化分析应该有助于解释基因功能和进化过程，而这些方法不能仅通过序列来解决：在不同组织中稳定选择基因表达水平的程度可以揭示基因在其中发挥作用的程度。最重要的角色；基因表达水平的进化约束强度有助于解释临床样本中观察到的表达水平；并鉴定其表达水平处于定向（正）选择下的基因可以帮助评估谱系和物种特异性表型的基础。

多项研究使用各种启发式方法分析了跨哺乳动物物种收集的表达数据，以确定保守和趋异的（基因）表达水平（Chan 等人 2009；Brawand 等人 2011；Merkin 等人 2012；Perry 等人 2012）。然而，目前对于解决与表达水平进化相关的功能问题的定量框架还没有达成共识，部分原因是对于如何最好地模拟哺乳动物的表达进化缺乏共识。在果蝇，研究发现，与序列进化不同，基因表达水平的差异在进化时间上并不是连续线性的。相反，由于稳定选择压力而达到饱和，需要比标准中性漂移模型更复杂的模型（Bedford 和 Hartl 2009；Kalinka 等人 2010）。相比之下，哺乳动物的初始基因表达研究受到小数据集的阻碍，导致关于哺乳动物谱系中中性漂移和稳定选择的相对贡献的报告不一致（Khaitovich 等人 2004；Yanai 等人 2004；Blekhman 等人. 2008 ; Brawand 等人 2011）。早期基于微阵列的研究观察到灵长类动物的表达差异和发散时间之间存在线性关系，这表明进化是中性的（Enard 等人 2002；Khaitovich 等人 2004 , 2005）。然而，随后的分析表明，这些观察结果被仅包含人类 DNA 探针的微阵列所混淆（Gilad 等人，2006 年），一旦考虑到这种情况，灵长类动物表达水平几乎没有差异，突出了稳定选择作为表达进化的主要模式。最近通过 RNA-seq 对九种哺乳动物的表达进化进行了大规模研究（Brawand et al. 2011)——减轻杂交技术的局限性——注意到更密切相关的物种确实具有更相似的表达水平（支持中性模型），但表达进化的速度似乎很慢（支持纯化选择的主要作用）。表达进化的正确模型缺乏明确性导致纯中性漂移模型 ( Perry et al. 2012 ) 和那些包含稳定选择的模型 ( Brawand et al. 2011 ) 使用冲突) 在跨哺乳动物的比较分析中，使不同的研究难以比较或解释。此外，尚未充分探索如何使用此类模型（一旦拟合）得出基因功能的结论，而不是关于适应增益和选择性效应的理论推论（Bedford 和 Hartl 2009；Nourmohammad 等人 2017）。基因表达的进化模型在深入了解转录途径及其与健康和疾病过程的关系方面的适用性尚未得到广泛探索。

在这里，我们使用来自 17 个哺乳动物物种和 7 种不同组织的综合 RNA-seq 数据集来表征整个哺乳动物谱系的表达进化模式。我们发现包含稳定选择的进化模型最适合描述哺乳动物的表达进化。我们基于先前提出的 Ornstein-Uhlenbeck (OU) 模型进一步开发了一个框架，以参数化进化最优基因表达的分布，我们使用该分布来量化稳定选择对基因表达的程度，识别有害的表达水平在患者表达数据中，并检测谱系特异性表达程序中的方向选择。

结果

哺乳动物物种之间的表达差异随着进化时间的增加而饱和

为了系统地探索哺乳动物的表达进化，我们编制了跨越 17 个物种和 7 种不同组织（脑、心脏、肌肉、肺、肾、肝、睾丸）的哺乳动物系统发育数据集。(图1A; 补充表 S1）。该数据集结合了 12 个物种的已发布数据（Harr 和 Turner 2010；Brawand 等人 2011；Merkin 等人 2012；Pipes 等人 2013；Cortez 等人 2014；Wong 等人 2015）和另外 5 个物种的数据。我们在这里新收集的物种以提高系统发育覆盖率（方法）。我们专注于 10,899 个 Ensembl 注释的哺乳动物一对一直系同源物 ( Aken et al. 2017）。我们通过重新比对跨物种的转录组来确认基因注释的质量，并发现 Ensembl 的 95%–99% 的一对一直系同源物也被我们的程序确定为互惠最佳比对；此外，Ensembl 注释的直系同源物与其人类对应物之间的平均序列同一性随着进化时间线性下降（补充图 S1）。此外，正如预期的那样，表达谱首先按组织聚类，然后按物种聚类，它们的层次聚类与系统发育树密切匹配（补充图 S2、S3）。

图1
跨哺乳动物谱系的表达进化由 OU 过程准确建模。(A)数据概览。所有 17 种哺乳动物的系统发育树（左图），以组织类型（彩色圆点）为标记，其中包括（表达）谱。(*) 新生成的数据。( B ) 表达分歧不是线性的。显示的是哺乳动物和人类之间肝脏样本在进化时间内的成对均方表达距离（y轴），通过每 100 bp 的替换（x轴）估计。（误差线）重复平均值的标准偏差；（实线）非线性（y = ax ^k）回归拟合。( C )OU 模型。描述 OU 模型的方程 (top ): (σ) 遗传漂变率；[dB(t)] 布朗运动；(θ) 最佳表达水平；(α) 选择强度。（左）在布朗运动（顶部）和 OU（底部）过程下，在进化时间（x轴）上的模拟表达轨迹（y轴）。显示了十个示例轨迹。（右）从 1000 个模拟轨迹到初始值（y轴）在时间（x轴）上的均方距离。( D ) 最优表达分布。（上）表达式（ y轴）的概率分布随时间变化的图示（x轴）在 OU 进程下。随着时间接近无穷大，分布趋于稳定。（底部）具有低（ NRBP1）和高（APOA4）方差的两个示例基因的所有肝脏样本（x轴）的log₁₀TPM值（y轴）散点图。（红色实线和虚线）分别估计使用 OU 过程估计的每个基因表达值的渐近（最佳）分布的均值和方差。请注意，均值和方差是在对数空间中计算的。

平均而言，物种之间的成对表达差异（补充方法；补充图 S4）在幂律关系中随着进化时间而饱和（图 1B)，与先前在果蝇中观察到的进化趋势一致( Bedford and Hartl 2009)。例如，当将每个物种的（表达）谱与相应的人类（表达）谱进行比较时，差异最初会随着进化距离的增加而发散，但这种趋势在灵长类谱系之外趋于稳定。这种关系在我们拥有所有灵长类动物（大脑、心脏、肾脏、肝脏、睾丸）的表达数据的五种组织中的每一种中都可以观察到（补充图 S5），并且不受不同数据源的批次效应或变异驱动每个物种可用的样本数量（补充方法；补充图 S6、S7）。在我们拥有多个Mus物种表达数据的两个组织中的每一个组织中，用Mus musculus作为参考物种时，我们观察到同样的关系。（补充图 S8，S9）。

表达进化可以被建模为一个Ornstein-Uhlenbeck过程

观察到的表达差异模式对应于 Ornstein-Uhlenbeck (OU) 过程（图 1 C,D)， Hansen (1997)最初提出作为一般连续表型进化模型的随机过程，最近被建议作为专门用于果蝇基因表达水平进化的合适模型( Bedford and Hartl 2009 ) .

在表达式级别的上下文中，OU 过程是对随机游走的修改，通过 dX _t = σdB _t + α(θ – X _t)dt 来描述表达式 (dX _t ) 随时间 (dt) 的变化，其中 dB _t表示布朗运动过程。该模型优雅地量化了任何给定基因的漂移和选择压力的贡献：（1）漂移由布朗运动建模，速率为 σ（图。1C，顶部），而（2）选择压力驱动表达回到最佳表达水平 θ 的强度由 α 参数化（图。1C，底部）。OU 过程结合了时间信息并充分考虑了系统发育关系，从而使我们能够适应个体进化表达轨迹。在较长的时间尺度上，漂移率 (σ) 和选择强度 (α) 之间的相互作用达到平衡，并且随着时间增加到无穷大，将表达式 X _t约束为稳定的正态分布，具有平均 θ，并且方差 σ ² /2α (图。1D)。

到目前为止，OU 模型主要用于关于适应度增益和不断变化的表达水平的选择性影响的理论推论（Bedford 和 Hartl 2009；Kalinka 等人 2010；Nourmohammad 等人 2017）。OU 模型在检测较小哺乳动物系统发育和不完全基因注释的表达选择的应用有限（Brawand 等人 2011；Rohlfs 和 Nielsen 2015）。然而，使用 OU 模型来表征基因表达的进化历史以获得生物学洞察力的完整功能尚未得到充分探索。

因此，我们接下来开发了 OU 模型的应用程序，以产生生物学上可解释的结果，以解决有关基因表达水平、基因功能和疾病基因发现的进化问题。首先，对于每个组织，我们根据我们的数据估计稳定选择下基因进化最优表达的渐近分布。我们证明了这种分布的 OU 方差（我们称之为“进化方差”）准确地描述了每个组织中基因表达水平的受限程度。其次，我们将患者数据中观察到的表达水平与进化模型估计的最佳表达分布进行比较，以检测潜在有害的表达水平并指定致病基因。第三，我们对OU模型(Butler和King 2004年)进行了扩展，该模型解释了系统发育中存在多种最优表达分布，以确定可能与谱系特异性适应相关的遗传路径。我们将依次描述这些应用程序。

大多数基因的表达在哺乳动物谱系内的稳定选择下进化

为了测试基因的表达是否在稳定选择下，我们使用了一个似然比测试来比较了没有选择的拟合（α= 0;棕色运动）（图2A，顶部），与稳定选择的一个（情况）（α> 0，OU过程）（图2A，底部）。（尽管我们直观地展示了关于单个参考物种的表达进化模式[例如，图 1B]，这两种模型都考虑了整个系统发育树的进化距离，而不仅仅是与一个参考物种的相对距离。）因为表达式低表达基因的级别估计与高技术变化相关，并且它们的真实生物变化不太可能被精确推断（Supplemental Fig. S10A; Silvestro et al. 2015），在所有分析中，我们只关注每百万分之一（TPM）表达超过 5 个转录本的基因，导致分析的基因在 3428 到 5822 个之间，具体取决于组织（Supplemental Methods; Supplemental Fig. S10B）。

平均而言，83% 的测试基因的表达（范围：77%–90%；错误发现率 [FDR] < 0.05）更适合稳定选择模型（图 2A、底部；补充图S11），尽管每个组织内数百个基因的表达似乎是中性进化的（图 2A top; 补充图S11）。57% (5669/8912) 基因的表达水平在所有表达它们的组织中处于稳定选择状态，39% (2722) 仅在一些表达它们的组织中处于稳定选择状态，只有6% （521）在我们研究的任何组织中都没有处于稳定选择之下（图 2B）。
（上面这一段的写法可以借鉴，多个数值的时候先说一个平均值/中位数/...，然后说明范围、相关qvalue等，一目了然数据分布，简洁）

我们使用jackknifing程序评估了我们在表达稳定选择下检测基因的敏感性和特异性，我们进行二次抽样以考虑3到16个物种的系统发育（补充方法）。正如预期的那样，在稳定选择（即拒绝零假设）下调用的基因数量随着包括更多物种而增加（补充图 S12A），但在 14 个物种时确实饱和（了）。重要的是，不一致率（相对于对完整数据集的分析）非常低：在子样本系统发育中被发现处于选择状态的基因，在完整系统发育中被发现为中性基因（即接受无效假设）的比例小于1%（补充图S12B）。

基因表达水平的进化分布预测基因功能

OU 过程在最初被提议用于对果蝇的表达进化进行建模时被认为具有吸引力，因为它能够区分中性选择和稳定选择。鉴于我们发现大多数哺乳动物基因处于稳定选择状态，我们接下来探索了 OU 模型估计基因表达水平稳定分布的能力，我们推断这是对进化最优表达分布的估计。因此，我们研究了使用 OU 模型的“进化方差”作为每个组织中基因表达的进化约束程度的定量测量。与上述相同的jackknifing程序（结果）表明，OU 模型的估计进化方差对子采样具有高度鲁棒性，这由在估计子采样系统发育的方差时非常低的均方误差 (MSE < 0.005) 确定（补充图 S12C）。事实上，当使用来自少于六个物种的数据时，我们发现进化方差比样本方差更加稳健，这并没有解释物种之间的系统发育关系（补充图 S12C）。
（这个稳健？个人理解是物种间都会有一个相对物种内更大的尺度，且相关指标比较稳定）

我们首先检查了跨组织的进化变异模式。为了控制每个组织中的样本数量，我们在与跨组织相同数量样本匹配的数据子集上重新拟合了 OU 进化参数（补充表 S1）。我们发现大脑具有最多的低方差基因（约束最多），而睾丸最少，这与之前对这些组织表达进化速率的估计一致。(图 2C; 补充图S13；陈等人。2009 年；布拉万等人。2011 年）。跨组织，方差合理相关（平均 Pearson 的r = 0.70）（补充图 S14A）。对于跨三个或更多组织表达的基因，表达水平和方差在体细胞组织中呈中度负相关(Pearson’s r =−0.25)，在27.2%的基因中，表达最高的组织与方差最低的组织相匹配(1263/4645)（补充图 S14B，顶部）；在该分析中进一步包括睾丸导致表达水平和方差之间几乎没有相关性（中位数 Pearson's r = -0.006）（补充图 S14B，底部）。

接下来，我们使用我们的完整数据集检查了组织内的进化变异模式。为了避免数据源的多样性引入的偏差，我们没有尝试解释方差的绝对值，而是专注于理解具有较低和较高方差的基因之间的相对关系。使用基于等级的基因本体论 (GO) 富集测试 ( Eden et al. 2009 )，我们发现进化方差和功能密切相关，这与之前未使用基于系统发育方法估计方差的比较研究的结果一致 ( Chan 等人 2009 年；Yue 等人 2014 年)：在所有组织中，变异低的基因富含管家功能（例如，RNA 结合和剪接、染色质组织、细胞周期），而变异高的基因富含细胞外蛋白（FDR < 0.001）。

一些过程仅在特定组织中富含具有低或高变异的基因（图 2D; 补充表 S2）：在具有组织特异性保护（低方差）的过程中，大脑中的突触蛋白（FDR = 0.011）和睾丸中的 Wnt 信号传导（FDR = 0.014）；具有高方差的过程包括心脏中的收缩纤维部分 (FDR = 0.005)、肾脏中的氧化还原酶活性 (FDR = 6.10 × 10 ^-6 ) 和肝脏中的脂质代谢 (FDR = 2.31 × 10 ^-9 )。当我们测试按表达水平排序的基因富集时，我们没有发现相同的富集类别。因此，我们可以依赖进化方差的估计作为表达约束和基因功能的指标。

我们发现表达和序列约束之间只有适度的相关性（Pearson's r = -0.25）（图 2E; 补充方法）。在表达和序列中保守的基因显着富集了管家过程（FDR < 10 ^-4）（补充表 S3），并且两者中不同的基因富集了免疫和炎症反应（FDR < 10 ^-6）。序列保守但表达不同的基因富含转录调节因子（FDR = 3.1 × 10 ^-5），尤其是那些参与胚胎形态发生的基因（FDR = 9.8 × 10 ^-8；例如，IRX5、HAND2、NOTCH1）。尽管表达水平的较高进化变异可能受环境、细胞类型组成的变化和遗传差异的影响，但我们的分析支持这样的假设，即没有蛋白质序列差异的基因调控差异可以解释物种特异性表型。

使用基因表达的进化分布来预测有害水平

在分析罕见疾病时，序列保守通常用于优先考虑那些在发生突变时对罕见疾病更为重要且可能导致罕见疾病的基因突变（Alföldi 和 Lindblad-Toh 2013；Jordan 等人 2015；Richards 等人 2015）。以此类推，我们假设表达保守性也应该可以预测基因的重要性。事实上，基因的表达水平要么是培养所必需的 ( Hart et al. 2014 )，要么是小鼠必需的 ( Georgi et al. 2013 )，要么是人类的单倍体不足 ( Rehm et al. 2015 )）在几乎所有组织中的进化方差（更高的约束）显着低于它们的非必要或单倍充足对应物（Wilcoxon 秩和检验P值 <0.01）（图 3A），关系不是由表达水平驱动的（补充图S15）。

图 3
基因表达的进化分布有助于识别导致疾病的基因。（A）必需基因的进化变异性较低。箱线图显示了培养中必需基因（顶部）、小鼠必需基因（中间）和人类单倍体不足（底部；深灰色）及其非必需或单倍体充足基因的 log（进化方差）（y轴）的分布（浅灰色）在七个组织中的每一个（x轴）：（***）P < 0.001；(**) P < 0.01。( B ) 疾病基因具有较低的进化变异性。箱线图显示对数的分布（进化方差）（y轴）与相关组织（大脑、心脏）中的高外显率单基因自闭症谱系障碍（上）、先天性心脏缺陷（中）和神经肌肉疾病（下）相关（深灰色）和不相关（浅灰色）的基因和肌肉）。（左）该组织中表达受限的基因（在三个或更少的组织中>5 TPM）；（右）普遍表达的基因；(***) P < 0.001；(*) <0.05。（C，D）使用进化分布或 GTEx RNA-seq 分布从肌营养不良患者的 RNA-seq 中识别异常基因表达的概述。( C）基于进化分布（左）或 GTEx RNA-seq 分布（右）的两种评分方法。（D）表格显示了使用进化数据（左）或 GTEx RNA-seq 数据（右）估计的分布时，所有肌营养不良患者的显着异常基因数、-log₁₀FDR 评分和DMD显着性等级。

我们接下来检查了三种环境中疾病基因的差异：与非综合征性自闭症谱系障碍 (ASD)（大脑）（Banerjee-Basu 和 Packer 2010）、先天性心脏缺陷（心脏）（Amberger 和 Hamosh ）直接相关的罕见单基因疾病基因2017 年；Blake 等人 2017 年）和神经肌肉疾病（骨骼肌）（补充方法；Cummings 等人 2017 年）。在每种情况下，具有组织限制表达的疾病基因（在三个或更少的组织中TPM > 5）（补充方法；补充图 S16）。在疾病相关组织中始终表现出比非疾病基因更低的方差（P值 <0.05） (图 3B; 补充图S17）。仅在 ASD 相关基因中，我们还观察到普遍表达的疾病基因与非疾病基因的差异显着降低（图 3B）。

接下来，我们假设每个基因的最佳 OU 分布的参数可以通过突出患者数据中的异常值、可能致病的基因表达水平来预测疾病基因。这类似于通过使用保守性来识别全外显子组测序中假定的致病序列突变来发现致病基因（Choi et al. 2009 ; O'Roak et al. 2011）。为此，我们获得了临床诊断为神经肌肉疾病的 93 名患者肌肉活检的 RNA-seq（方法；补充表 S4）。对于每个患者样本，我们计算每个基因的z值，以评估它们如何偏离骨骼肌中该基因表达的(最佳)进化适合度，并对多重假设检验进行校正（方法；补充图 S18A）。与来自 184 名健康人的 GTEx 肌肉样本（The GTEx Consortium 2013）相比，通过这种测量，患者的总体失调基因平均多 3.2 倍（Wilcoxon 秩和检验P值 = 2 × 10 ^-9）（补充图S18B），表明OU模型拟合的进化参数可用于检测更可能有害的异常值表达值。

然后，我们测试了 OU 模型是否可用于识别罕见病分析中的致病基因。作为原则证明，我们关注了肌肉疾病队列中临床诊断为 Becker 或 Duchenne 肌营养不良症的 8 名患者的子集，包括通过免疫印迹确认肌营养不良蛋白缺失或减少（Cummings 等人，2017 年）。为了将我们的方法与标准差异表达分析进行比较，我们通过离群值表达对基因进行排序，Z分数基于（1）与从我们的进化数据估计的均值和方差进行比较；或 (2) 与仅从健康 GTEx 人类数据估计的均值和方差进行比较（图 3C）。根据我们的进化数据，在每位患者中被列为显着异常值的基因较少（中位数：4，范围：0-32），并且DMD基因在 8 名患者中的 6 名中被列为最高或第二最异常表达的基因，每个都显示出显着表达不足 (FDR < 10 ^-3 ) (图 3D)。相比之下，参考 GTEx 表达数据的评分并未产生如此具体的结果：异常值的基因数中位值为14.5（范围：0-250），8 名患者中只有 4 名被称为DMD显着低表达（FDR < 10 ^-3） , 及其在这些患者中的重要性排名在 1 到 50 之间。这种特异性差异可能反映了在使用进化数据与人类数据时对健康（可容忍）方差的更准确估计：尽管对平均表达的估计在两种方法之间高度一致（补充图. S19，左），当使用进化数据集估计时，表达方差几乎总是更大（补充图 S19，右），反映了在较长时间周期上，每个基因必须充分探索生理上可接受的表达水平空间。因此，使用人类 GTEx 分布导致更多的假阳性基因似乎异常表达。相反，使用 OU 模型对最佳基因表达的进化均值和方差的估计有助于检测实际疾病基因的失调，并有助于发现新的疾病基因。重要的是，与患者和健康对照之间差异表达的方法相比，这种方法不需要对照群体，并且可以针对个体患者样本进行。

多变量 OU 模型可用于检测谱系特异性表达变化

最后，我们探索了使用 OU 模型来检测基因表达中的方向选择。我们使用了模型的扩展，该模型通过将表达水平的分布建模为多元正态分布来解释单个系统发育中的多种选择方案，其均值和方差是针对每个（预定义的）子进化枝估计的。(图 4A; 巴特勒和金 2004；罗尔夫斯等人。2014 年）。该扩展 OU 模型的先前应用确定了哺乳动物系统发育中超过 9000 种表达变化（Brawand 等人，2011 年），但该分析依赖于较小的系统发育，因此专注于识别基因表达的物种特异性变化，不幸的是，这种变化可能是容易被环境原因或技术影响所混淆。

图 4
多变量 OU 过程可以检测谱系特定的表达变化。( A ) 多元 OU 过程。在多变量 OU 过程下，随着时间（x轴）的模拟表达轨迹（y轴）。灰色轨迹是从同一分布 (N ₀ ) 跨时间采样的，而橙色轨迹从相同的祖先分布 (N ₀ ) 开始，但在物种形成事件后在新分布 (N _{1 ) 下演变。}( B ) 表达进化的三个测试假设：从左到右：单变量OU _all模型，其中基因表达在整个系统发育的单一稳定机制下进化（黑色），以及两个多变量 OU 模型，OU_{灵长类动物}和 OU_啮齿动物，其中基因表达在祖先机制下进化（黑色）和指定的新机制子进化枝（橙色）。( C ) 肝脏中的谱系特异性表达。肝脏样本中的成对均方表达距离 ( y轴）在参考物种（标记为黑点）和其他哺乳动物之间，用于分配给三个测试 OU 模型中的每一个的基因。（黑点）在祖先分布下进化的物种；（标记为橙色点）谱系分裂后在新制度下进化的物种；（实线）非线性回归拟合在祖先分布下进化的物种。（D）为谱系特异性表达丰富的示例过程。热图显示来自基因（列）的列归一化表达（红色：高；蓝色：低），在三个富集的 GO 类别（FDR < 0.05）中具有谱系特异性表达模式：肝脏中的脂质转运（左），肝脏中的免疫调节（中）和睾丸中的微管运动（右）。

我们利用我们更全面的系统发育覆盖范围，专注于检测三个或更多物种的整个亚进化枝的方向和幅度一致的表达变化，它们的样本被收集并跨多个来源测序，以减轻非遗传混杂因素。我们根据Butler 和 King (2004)建议的方法确定了整个树的“差异基因表达”（方法）：对于每个基因，我们应用了标准的单变量 OU 模型，该模型对所有物种使用单一最优值，以及扩展 OU 模型，该模型使用两个最优值——一个用于祖先分布，一个用于内部分布感兴趣的进化枝——并使用拟合优度测试分配了最佳 OU 模型。作为保守措施，我们只保留了那些在亚进化枝之间也发生至少两倍变化并且在其中一个亚进化枝中平均表达水平至少为 1 TPM 的基因。

我们首先评估了这种方法在增加系统发育距离时检测谱系特异性表达的能力，方法是测试所有灵长类动物（分支长度 = 0.121）、啮齿动物（分支长度 = 0.177）、劳亚动物（分支长度= 0.407）或兔形目动物（分支长度 = 0.575）。我们构建我们的数据集，以便我们针对八个比较物种测试感兴趣的进化枝内三个物种之间的共享差异表达变化，并且为了避免批次效应的混杂因素，我们选择三个感兴趣的物种，以便获得每个物种的数据来自不同的来源（补充方法；补充表 S1）。正如预期的那样，我们发现检测到的差异表达基因的数量随着感兴趣的进化枝内距离的增加而单调减少，范围从灵长类进化枝中的 470 个基因到兔形进化枝中的 327 个（补充图 S20A）。不幸的是，我们通过改组物种表达数据（方法）估计的这种分析的错误发现率从 54% 到 78% 不等。

为了提高我们检测差异表达的能力，我们转向我们的完整数据集并测试灵长类动物进化枝（OU_{灵长类动物}，5-7 灵长类动物与 8-10 个比较物种）和啮齿动物进化枝（OU_啮齿动物，3-5啮齿动物与 10-12 个比较物种）（图 4B）。即使有了这个更大的系统发育，我们的 FDR 范围从 18% 到 52%，这表明未来对差异表达的研究应该依赖于更大的系统发育。不同组织的不同富集大小可能很大程度上是由样本大小的差异驱动的。当使用与跨组织样本大小匹配的数据集时，差异表达的表达基因的百分比在组织之间是相当一致的（补充图 S20B）。

尽管我们的分析能力有限，但我们能够在肝脏（灵长类动物：FDR = 0.18；啮齿动物：FDR = 0.27）和睾丸（灵长类动物：FDR = 0.29；啮齿动物：FDR = 0.29）中实现 FDR < 30%与肺（FDR = 0.26）和脑（FDR = 0.18）的灵长类进化枝一样（补充图S21），我们进一步检查了这些差异表达基因组。例如，在肝脏中，我们分别在灵长类动物和啮齿动物中鉴定了 640 个和 794 个具有谱系特异性表达变化的基因，突出了每个进化枝中调节不同的特定代谢过程。仅当存在进化枝特异性选择时，谱系特异性基因的表达水平才显着偏离预期（图 4C）。由于与以前的应用相比，差异表达基因的集合更大，我们可以识别谱系特异性基因之间的功能富集（补充表 S5）。我们发现与肝脏中许多脂质代谢过程相关的基因的灵长类动物特异性下调（FDR = 1.88 × 10 ^-11）。这些过程包括过氧化物酶体功能（FDR = 2.45 × 10 ^-8）、脂肪酸代谢（FDR = 1.52 × 10 ^-8）和脂质转运（FDR = 3.36 × 10 ^-3）（图 4D)，并含有已知的脂质代谢调节剂，例如 LDL 受体 ( LDLR ) ( Brown and Goldstein 1979 )、肝脂肪酶 ( LIPC ) ( Guerra et al. 1997 ) 和转录因子PPARA ( Kersten 2014 )。因此，多种途径的表达可能在祖先灵长类动物分支上有所不同，这与没有进一步基因改造的小鼠无法很好地模拟人类脂质血症的观察结果一致（von Scheidt et al. 2017）。在其他例子中，参与调节免疫反应的基因在啮齿动物肝脏（FDR = 6.97 × 10 ^-4）和基于微管的运动基因（FDR = 2.82 × 10^-3 ) 和精子发生 (FDR = 2.82 × 10 ^-2 ) 在灵长类动物的睾丸中下调 (图 4D)，反映了已知的免疫相关基因 ( Kosiol et al. 2008 ; Areal et al. 2011 ; Yue et al. 2014 ) 和生殖相关基因 ( Swanson et al. 2001 ; Torgerson et al. 2002 ) 的快速进化，分别。

讨论

在这里，我们将跨哺乳动物的比较基因表达谱的大型数据集与系统分析相结合，并表明一对一哺乳动物直系同源物的基因表达差异在进化时间中饱和，这可以通过 OU 过程很好地建模. 我们进一步展示了如何使用 OU 模型来查询基因功能，评估候选疾病基因的有害表达水平，以及识别谱系特异性的表达进化。

与任何比较物种分析一样，批次效应引入的伪影，或直系同源物或转录本注释中的错误可能会使我们的数据产生偏差。然而，我们观察到的进化模式不仅仅来自单个批次，因为每个数据源都包含跨越整个系统发育树的物种。此外，我们仅使用一对一的哺乳动物直系同源物有助于减轻直系分配中的错误，并且我们确认我们的注释转录本的序列身份在进化时间内呈线性发散，因此不会驱动观察到的表达进化模式。

表达进化的非线性模式由先前提出的 OU 过程 ( Hansen 1997 ) 准确建模，该模型结合了中性漂移和稳定选择。尽管我们发现稳定选择在哺乳动物谱系中的表达进化中起主导作用，但我们注意到适当的模型取决于进化距离：在灵长类动物谱系中，我们确实发现表达差异接近线性，证实了原始研究提出了表达进化的中性模型（Enard et al. 2002 ; Khaitovich et al. 2004）；但在更大的进化距离上，稳定选择的影响越来越大，正如最近在哺乳动物中进行的 RNA-seq 研究中所指出的那样。布拉万等人。2011 年）和果蝇（Bedford 和 Hartl 2009 年；Nourmohammad 等人 2017 年）。

重要的是，尽管 OU 过程准确地模拟了我们的数据，但还有其他因素可能会限制基因的表达。这些因素包括 (1) 定义的基因表达下限（即 0 TPM）；(2)基因表达的上限；(3) 在一个组织中对一个基因的选择压力对同一基因在其他组织中的表达有限制作用；(4) 对一个形成反式的基因的选择压力对其他基因的表达限制的影响。尽管案例 (1) 和 (2) 代表了我们研究中测试的所有基因的一小部分，尤其是使用我们的过滤器过滤表达 > 5 TPM 的基因，我们无法分离由于间接选择力而表达受限的基因，如案例（3）和（4），来自使用我们当前数据的直接选择压力下的表达水平。尽管如此，OU 模型仍然是描述基因表达进化历史的重要定量工具，并为进一步探究表达进化的机制奠定了基础。

尽管以前使用 OU 模型和基因表达数据的研究集中在表达进化的理论方面，但我们现在展示如何使用 OU 模型来估计最佳基因表达水平的分布，并回答有关基因功能的生理和临床相关问题，包括从个体患者数据中检测有害表达水平。因为我们的数据来自可能影响我们对进化分布估计的准确性的各种来源，所以我们小心地只分析相对进化方差（例如，基于等级的 GO 富集测试）或构建我们数据集的子集进行分析避免批次效应（例如，测试多个物种的表达变化，每个物种都从不同的来源收集）。然而，当直接使用进化均值和方差估计进行疾病表达分析时，我们发现我们的估计在明确识别疾病基因的异常表达方面优于健康人类数据。这表明进化估计对技术差异具有稳健性，并且最终可以跨多种组织类型提供，以帮助科学和临床发现。

我们最终应用了一个多元 OU 模型（Butler and King 2004；Rohlfs et al. 2014) 以识别我们数据中跨物种进化枝的谱系特异性表达变化。尽管我们改进了以前依赖较小数据集的研究，即使有 17 个物种，我们也无法达到低于 18% 的 FDR，并且我们发现在更远的进化枝之间共享的差异表达变化更少。这表明未来研究谱系特异性表达的研究将需要从更大的系统发育中取样的数据，并且在检测密切相关物种之间共享的谱系特异性表达变化时表现更好。我们还注意到，这种方法的一个缺点是必须首先构建经过检验的假设（例如，灵长类动物与非灵长类动物；啮齿类动物与非啮齿类动物），并且最佳拟合模型可能无法真正反映潜在的进化历史。

展望未来，我们预计 OU 模型可以进一步开发用于其他生物学查询，例如，测试跨基因或旁系同源家族途径的稳定选择，估计祖先表达状态，或应用于其他转录区域，如短或长非编码 RNA。正如我们的分析所示，在不同的发育和环境背景下，其他组织类型和物种的表达特征将为基因表达的进化以及基因型和表型之间的关系提供更大的动力和进一步的洞察力。

方法

用于进化数据集的 RNA-seq

来自狗和兔组织的 RNA 样品购自 Zyagen。来自雪貂组织的 RNA 样本是 John Englehardt（爱荷华大学）的礼物。来自负鼠组织的 RNA 样本是来自 Paul Samollow (Texas A&M) 的礼物。来自犰狳组织的 RNA 样本是 Jason Merkin 和 Christopher Burge (MIT) 的礼物。所有组织收集均经 IACUC 批准，并按照各自的机构指南进行。RNA-seq 文库的制备如Levin 等人所述。（2010）（补充方法）。样本在 Illumina HiSeq 2000 测序仪上进行测序，最小深度为 3500 万次读数。

比对和表达量化

RSEM v1.2.12 ( Li and Dewey 2011 ) 与默认参数一起使用，以将读数与每个物种的转录组对齐并量化每个基因的 TPM。

基因表达值的标准化

使用来自 Bioconductor package edgeR ( Robinson et al. 2010 ) 的 TMM 归一化 ( Robinson and Oshlack 2010 ) 对基因表达值 (log ₁₀ TPM) 进行归一化。简而言之，TMM 标准化假设大多数基因在样本之间没有差异表达 (DE)，并估计一对样本之间的比例因子，使得对数表达比的修剪平均值（M 值的修剪平均值 [TMM]）相等到 1. 假设物种对之间的大多数基因不是 DE 是合理的，因为即使在遥远的哺乳动物（例如人类和负鼠）之间，给定组织中表达水平的 Pearson 相关性也大于 0.75。

拟合 OU 工艺参数

BM 和 OU 模型使用带有默认参数的 R 包ouch ( Butler and King 2004 )拟合归一化表达式值。使用比较 OU（替代假设）与 BM（零假设）模型的似然比检验计算每个基因的*P值，然后使用 Benjamini-Hochberg FDR 程序（ *Benjamini 和 Hochberg 1995）对多假设检验进行校正。

神经肌肉疾病数据集的样本

本研究中描述的神经肌肉疾病患者 RNA-seq 队列是Cummings 等人描述的队列的超集。（2017 年）（dbGaP 加入 phs000655.v3.p1）和另外 30 名患者。组织是根据机构审查委员会 (IRB) 批准的协议在国家神经疾病和中风研究所 (协议 #12-N-0095)、纽卡斯尔大学 (CF01.2011)、波士顿儿童医院 (03-12-205R) 采购的、伦敦大学学院 (08ND17)、加州大学洛杉矶分校 (15-001919) 和圣裘德儿童研究医院 (10/CHW/45)。患者在活检前就诊时同意这些方案。如Cummings 等人所述收集患者肌肉活检。（2017）并使用与 GTEx 项目相同的协议进行测序（补充方法; GTEx 联盟 2013 年）。

人体肌肉数据的对齐和表达量化

GTEx BAM 文件从登录 ID phs000424.v6.p1 下的 dbGaP 下载，并在使用 Picard SamToFastq 转换为 FASTQ 文件后重新对齐。使用 hg19 作为基因组参考，使用 STAR 2-Pass v.2.4.2a (Dobin et al. 2013 )比对患者和 GTEx 读数。使用 GENCODE v19 注释使用RNA-SeQC v1.1.8 ( DeLuca et al. 2012 ) 对表达进行量化。

检测患者样本中的异常值表达

首先通过 TMM 标准化（Robinson 和 Oshlack 2010 ）将基因表达值（log ₁₀ TPM）标准化为来自进化数据集的人类骨骼肌表达值。对于每个患者样本中的每个基因，使用从进化数据估计的渐近平均值和方差计算Z分数。Z分数仅针对被评估为在 OU 而不是 BM 模型下更适合的基因计算，并且其渐近平均值估计为 5 TPM 或更高。Z分数转换为P-值，然后使用 Benjamini-Hochberg FDR 程序对多个假设检验进行校正。我们使用 0.01 的 FDR 阈值来初步定义显着性。其中，我们删除了另外 330 个基因，这些基因在超过 25% 的 GTEx 样本中得分为显着异常值。作为比较器，还使用从健康人类 GTEx 样本估计的样本均值和方差计算Z分数。为了确保这两种方法之间的可比性，我们只计算了在上述进化方法中的任何步骤都没有被过滤掉的基因的Z分数。

检测谱系特异性表达程序

在所有谱系特异性差异表达分析中，使用将每个 OU 模型与 BM 模型进行比较的似然比检验计算P值，然后使用 Benjamini-Hochberg FDR 程序针对多重假设检验进行调整。对于每个基因，在所有模型上计算 Akaike 和贝叶斯信息标准（AIC 和 BIC）得分，这些模型对 null 具有显着性，以确定最佳拟合模型。在所有情况下，这两个分数都同意最佳拟合模型。为了估计 FDR，我们使用改组的表达数据在每个组织中执行相同的程序，其中来自一个物种的表达数据被随机重新分配给不同的物种。（在相关方法中 [ Brawand et al. 2011 ; Rohlfs and Nielsen 2015 ]，Q 值是通过针对零假设 θ _subclade = θ _ancestral直接测试备择假设 θ _subclade ≠ θ _ancestral得出的，并针对多个假设检验进行调整。然而，我们发现这种严格的方法导致多重测试的负担更大。）我们分别对每组上调和下调的基因进行了 GO 富集分析，使用平均表达至少为 1 的背景基因组适当组织中所有物种的 TPM。

数据访问

本研究的 RNA-seq 原始数据（来自兔子、狗、雪貂和负鼠）已提交给 NCBI 基因表达综合系统（GEO；http: //www.ncbi.nlm.nih.gov/geo/ ）编号GSE106077。每个组织环境中所有一对一哺乳动物直系同源物的处理表达数据和进化表达分布可在https://portals.broadinstitute.org/evee/获得。

end