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

摘要

基因的進化史有助于預測其功能以及與表型特征的關系丧荐。盡管序列保守通常用于破譯基因功能和評估醫(yī)學相關性寄悯，但缺乏從比較表達數(shù)據(jù)進行功能推斷的方法赠堵。在這里银觅，我們使用來自17個哺乳動物物種的7個組織的RNA-seq來證明跨哺乳動物的表達進化是由 Ornstein-Uhlenbeck過程準確建模的，這是一種普遍提出的連續(xù)性狀進化模型扫沼。我們應用此模型來識別中性江掩、穩(wěn)定和定向選擇下的表達途徑。我們進一步論證了該模型的新應用核畴，以量化基因表達的穩(wěn)定選擇程度膝但，參數(shù)化每個基因的最佳表達水平分布，并檢測個體患者的表達數(shù)據(jù)中的有害表達水平谤草。我們的工作為解釋跨物種和疾病的表達數(shù)據(jù)提供了一個統(tǒng)計框架跟束。

比較基因組學已經(jīng)通過跨物種的進化模式識別和注釋了功能性遺傳元件（Rubin 等人 2000莺奸；Kellis 等人 2003；Siepel 等人 2005冀宴；Pollard 等人 2006灭贷；Lindblad-Toh 等人 2011）。當前的比較研究主要集中在基因組序列的分析上略贮，依賴于從觀察中性序列隨時間線性發(fā)散而發(fā)展起來的成熟理論框架（Harris 1966甚疟；Lewontin 和 Hubby 1966；Kimura 1968逃延；Jukes 和 King 1971）览妖。這些方法允許檢測在中性進化的無效模型下進化比預期更慢（例如，由于純化選擇）或更快（例如揽祥，由于正選擇或放松的選擇約束）的序列元件黄痪。

長期以來俄精，人們普遍認為基因調控的分歧，表現(xiàn)為基因表達的表型變化啤誊，在進化中也起著關鍵作用（King and Wilson 1975贩绕；Wang et al. 1996；Pierce and Crawford 1997斩松；Ferea et al. 1999；Fraser等人 2010）≌拘鳎基因表達的進化分析應該有助于解釋基因功能和進化過程，而這些方法不能僅通過序列來解決：在不同組織中穩(wěn)定選擇基因表達水平的程度可以揭示基因在其中發(fā)揮作用的程度丽柿。最重要的角色恢准；基因表達水平的進化約束強度有助于解釋臨床樣本中觀察到的表達水平；并鑒定其表達水平處于定向（正）選擇下的基因可以幫助評估譜系和物種特異性表型的基礎甫题。

多項研究使用各種啟發(fā)式方法分析了跨哺乳動物物種收集的表達數(shù)據(jù)馁筐，以確定保守和趨異的（基因）表達水平（Chan 等人 2009；Brawand 等人 2011坠非；Merkin 等人 2012敏沉；Perry 等人 2012）。然而炎码，目前對于解決與表達水平進化相關的功能問題的定量框架還沒有達成共識盟迟，部分原因是對于如何最好地模擬哺乳動物的表達進化缺乏共識。在果蠅潦闲，研究發(fā)現(xiàn)攒菠，與序列進化不同，基因表達水平的差異在進化時間上并不是連續(xù)線性的歉闰。相反辖众，由于穩(wěn)定選擇壓力而達到飽和卓起，需要比標準中性漂移模型更復雜的模型（Bedford 和 Hartl 2009；Kalinka 等人 2010）赵辕。相比之下既绩，哺乳動物的初始基因表達研究受到小數(shù)據(jù)集的阻礙，導致關于哺乳動物譜系中中性漂移和穩(wěn)定選擇的相對貢獻的報告不一致（Khaitovich 等人 2004还惠；Yanai 等人 2004饲握；Blekhman 等人. 2008 ; Brawand 等人 2011）。早期基于微陣列的研究觀察到靈長類動物的表達差異和發(fā)散時間之間存在線性關系蚕键，這表明進化是中性的（Enard 等人 2002救欧；Khaitovich 等人 2004 , 2005）。然而锣光，隨后的分析表明笆怠，這些觀察結果被僅包含人類 DNA 探針的微陣列所混淆（Gilad 等人，2006 年）誊爹，一旦考慮到這種情況蹬刷，靈長類動物表達水平幾乎沒有差異，突出了穩(wěn)定選擇作為表達進化的主要模式频丘。最近通過 RNA-seq 對九種哺乳動物的表達進化進行了大規(guī)模研究（Brawand et al. 2011)——減輕雜交技術的局限性——注意到更密切相關的物種確實具有更相似的表達水平（支持中性模型）办成，但表達進化的速度似乎很慢（支持純化選擇的主要作用）。表達進化的正確模型缺乏明確性導致純中性漂移模型 ( Perry et al. 2012 ) 和那些包含穩(wěn)定選擇的模型 ( Brawand et al. 2011 ) 使用沖突) 在跨哺乳動物的比較分析中搂漠，使不同的研究難以比較或解釋迂卢。此外，尚未充分探索如何使用此類模型（一旦擬合）得出基因功能的結論桐汤，而不是關于適應增益和選擇性效應的理論推論（Bedford 和 Hartl 2009而克；Nourmohammad 等人 2017）≌基因表達的進化模型在深入了解轉錄途徑及其與健康和疾病過程的關系方面的適用性尚未得到廣泛探索员萍。

在這里，我們使用來自 17 個哺乳動物物種和 7 種不同組織的綜合 RNA-seq 數(shù)據(jù)集來表征整個哺乳動物譜系的表達進化模式拣度。我們發(fā)現(xiàn)包含穩(wěn)定選擇的進化模型最適合描述哺乳動物的表達進化充活。我們基于先前提出的 Ornstein-Uhlenbeck (OU) 模型進一步開發(fā)了一個框架，以參數(shù)化進化最優(yōu)基因表達的分布蜡娶，我們使用該分布來量化穩(wěn)定選擇對基因表達的程度混卵，識別有害的表達水平在患者表達數(shù)據(jù)中，并檢測譜系特異性表達程序中的方向選擇窖张。

結果

哺乳動物物種之間的表達差異隨著進化時間的增加而飽和

為了系統(tǒng)地探索哺乳動物的表達進化幕随，我們編制了跨越 17 個物種和 7 種不同組織（腦、心臟宿接、肌肉赘淮、肺辕录、腎、肝梢卸、睪丸）的哺乳動物系統(tǒng)發(fā)育數(shù)據(jù)集走诞。(圖1A; 補充表 S1）。該數(shù)據(jù)集結合了 12 個物種的已發(fā)布數(shù)據(jù)（Harr 和 Turner 2010蛤高；Brawand 等人 2011蚣旱；Merkin 等人 2012；Pipes 等人 2013戴陡；Cortez 等人 2014塞绿；Wong 等人 2015）和另外 5 個物種的數(shù)據(jù)。我們在這里新收集的物種以提高系統(tǒng)發(fā)育覆蓋率（方法）恤批。我們專注于 10,899 個 Ensembl 注釋的哺乳動物一對一直系同源物 ( Aken et al. 2017）异吻。我們通過重新比對跨物種的轉錄組來確認基因注釋的質量，并發(fā)現(xiàn) Ensembl 的 95%–99% 的一對一直系同源物也被我們的程序確定為互惠最佳比對喜庞；此外诀浪，Ensembl 注釋的直系同源物與其人類對應物之間的平均序列同一性隨著進化時間線性下降（補充圖 S1）。此外延都，正如預期的那樣雷猪，表達譜首先按組織聚類，然后按物種聚類窄潭，它們的層次聚類與系統(tǒng)發(fā)育樹密切匹配（補充圖 S2、S3）酵颁。

圖1
跨哺乳動物譜系的表達進化由 OU 過程準確建模嫉你。(A)數(shù)據(jù)概覽。所有 17 種哺乳動物的系統(tǒng)發(fā)育樹（左圖）躏惋，以組織類型（彩色圓點）為標記幽污，其中包括（表達）譜。(*) 新生成的數(shù)據(jù)簿姨。( B ) 表達分歧不是線性的距误。顯示的是哺乳動物和人類之間肝臟樣本在進化時間內的成對均方表達距離（y軸），通過每 100 bp 的替換（x軸）估計扁位。（誤差線）重復平均值的標準偏差准潭；（實線）非線性（y = ax ^k）回歸擬合。( C )OU 模型域仇。描述 OU 模型的方程 (top ): (σ) 遺傳漂變率刑然；[dB(t)] 布朗運動；(θ) 最佳表達水平暇务；(α) 選擇強度泼掠。（左）在布朗運動（頂部）和 OU（底部）過程下怔软，在進化時間（x軸）上的模擬表達軌跡（y軸）。顯示了十個示例軌跡择镇。（右）從 1000 個模擬軌跡到初始值（y軸）在時間（x軸）上的均方距離挡逼。( D ) 最優(yōu)表達分布。（上）表達式（ y軸）的概率分布隨時間變化的圖示（x軸）在 OU 進程下腻豌。隨著時間接近無窮大家坎，分布趨于穩(wěn)定。（底部）具有低（ NRBP1）和高（APOA4）方差的兩個示例基因的所有肝臟樣本（x軸）的log₁₀TPM值（y軸）散點圖饲梭。（紅色實線和虛線）分別估計使用 OU 過程估計的每個基因表達值的漸近（最佳）分布的均值和方差乘盖。請注意，均值和方差是在對數(shù)空間中計算的憔涉。

平均而言订框，物種之間的成對表達差異（補充方法；補充圖 S4）在冪律關系中隨著進化時間而飽和（圖 1B)兜叨，與先前在果蠅中觀察到的進化趨勢一致( Bedford and Hartl 2009)穿扳。例如，當將每個物種的（表達）譜與相應的人類（表達）譜進行比較時国旷，差異最初會隨著進化距離的增加而發(fā)散矛物，但這種趨勢在靈長類譜系之外趨于穩(wěn)定。這種關系在我們擁有所有靈長類動物（大腦跪但、心臟履羞、腎臟、肝臟屡久、睪丸）的表達數(shù)據(jù)的五種組織中的每一種中都可以觀察到（補充圖 S5）忆首，并且不受不同數(shù)據(jù)源的批次效應或變異驅動每個物種可用的樣本數(shù)量（補充方法；補充圖 S6被环、S7）糙及。在我們擁有多個Mus物種表達數(shù)據(jù)的兩個組織中的每一個組織中，用Mus musculus作為參考物種時筛欢，我們觀察到同樣的關系浸锨。（補充圖 S8，S9）版姑。

表達進化可以被建模為一個Ornstein-Uhlenbeck過程

觀察到的表達差異模式對應于 Ornstein-Uhlenbeck (OU) 過程（圖 1 C,D)柱搜， Hansen (1997)最初提出作為一般連續(xù)表型進化模型的隨機過程，最近被建議作為專門用于果蠅基因表達水平進化的合適模型( Bedford and Hartl 2009 ) .

在表達式級別的上下文中剥险，OU 過程是對隨機游走的修改冯凹，通過 dX _t = σdB _t + α(θ – X _t)dt 來描述表達式 (dX _t ) 隨時間 (dt) 的變化，其中 dB _t表示布朗運動過程。該模型優(yōu)雅地量化了任何給定基因的漂移和選擇壓力的貢獻：（1）漂移由布朗運動建模宇姚，速率為 σ（圖匈庭。1C，頂部）浑劳，而（2）選擇壓力驅動表達回到最佳表達水平 θ 的強度由 α 參數(shù)化（圖阱持。1C，底部）魔熏。OU 過程結合了時間信息并充分考慮了系統(tǒng)發(fā)育關系衷咽，從而使我們能夠適應個體進化表達軌跡。在較長的時間尺度上蒜绽，漂移率 (σ) 和選擇強度 (α) 之間的相互作用達到平衡镶骗，并且隨著時間增加到無窮大，將表達式 X _t約束為穩(wěn)定的正態(tài)分布躲雅，具有平均 θ鼎姊，并且方差 σ ² /2α (圖。1D)相赁。

到目前為止相寇，OU 模型主要用于關于適應度增益和不斷變化的表達水平的選擇性影響的理論推論（Bedford 和 Hartl 2009；Kalinka 等人 2010钮科；Nourmohammad 等人 2017）唤衫。OU 模型在檢測較小哺乳動物系統(tǒng)發(fā)育和不完全基因注釋的表達選擇的應用有限（Brawand 等人 2011；Rohlfs 和 Nielsen 2015）绵脯。然而佳励，使用 OU 模型來表征基因表達的進化歷史以獲得生物學洞察力的完整功能尚未得到充分探索。

因此蛆挫，我們接下來開發(fā)了 OU 模型的應用程序赃承，以產(chǎn)生生物學上可解釋的結果，以解決有關基因表達水平璃吧、基因功能和疾病基因發(fā)現(xiàn)的進化問題楣导。首先废境，對于每個組織畜挨，我們根據(jù)我們的數(shù)據(jù)估計穩(wěn)定選擇下基因進化最優(yōu)表達的漸近分布。我們證明了這種分布的 OU 方差（我們稱之為“進化方差”）準確地描述了每個組織中基因表達水平的受限程度噩凹。其次巴元，我們將患者數(shù)據(jù)中觀察到的表達水平與進化模型估計的最佳表達分布進行比較，以檢測潛在有害的表達水平并指定致病基因驮宴。第三逮刨，我們對OU模型(Butler和King 2004年)進行了擴展，該模型解釋了系統(tǒng)發(fā)育中存在多種最優(yōu)表達分布，以確定可能與譜系特異性適應相關的遺傳路徑修己。我們將依次描述這些應用程序恢总。

大多數(shù)基因的表達在哺乳動物譜系內的穩(wěn)定選擇下進化

為了測試基因的表達是否在穩(wěn)定選擇下，我們使用了一個似然比測試來比較了沒有選擇的擬合（α= 0;棕色運動）（圖2A睬愤，頂部）片仿，與穩(wěn)定選擇的一個（情況）（α> 0，OU過程）（圖2A尤辱，底部）砂豌。（盡管我們直觀地展示了關于單個參考物種的表達進化模式[例如，圖 1B]光督，這兩種模型都考慮了整個系統(tǒng)發(fā)育樹的進化距離阳距，而不僅僅是與一個參考物種的相對距離。）因為表達式低表達基因的級別估計與高技術變化相關结借，并且它們的真實生物變化不太可能被精確推斷（Supplemental Fig. S10A; Silvestro et al. 2015）筐摘，在所有分析中，我們只關注每百萬分之一（TPM）表達超過 5 個轉錄本的基因映跟，導致分析的基因在 3428 到 5822 個之間蓄拣，具體取決于組織（Supplemental Methods; Supplemental Fig. S10B）。

平均而言，83% 的測試基因的表達（范圍：77%–90%挤巡；錯誤發(fā)現(xiàn)率 [FDR] < 0.05）更適合穩(wěn)定選擇模型（圖 2A剩彬、底部；補充圖S11）矿卑，盡管每個組織內數(shù)百個基因的表達似乎是中性進化的（圖 2A top; 補充圖S11）襟衰。57% (5669/8912) 基因的表達水平在所有表達它們的組織中處于穩(wěn)定選擇狀態(tài)，39% (2722) 僅在一些表達它們的組織中處于穩(wěn)定選擇狀態(tài)粪摘，只有6% （521）在我們研究的任何組織中都沒有處于穩(wěn)定選擇之下（圖 2B）瀑晒。
（上面這一段的寫法可以借鑒，多個數(shù)值的時候先說一個平均值/中位數(shù)/...徘意，然后說明范圍苔悦、相關qvalue等，一目了然數(shù)據(jù)分布椎咧，簡潔）

我們使用jackknifing程序評估了我們在表達穩(wěn)定選擇下檢測基因的敏感性和特異性玖详，我們進行二次抽樣以考慮3到16個物種的系統(tǒng)發(fā)育（補充方法）。正如預期的那樣勤讽，在穩(wěn)定選擇（即拒絕零假設）下調用的基因數(shù)量隨著包括更多物種而增加（補充圖 S12A）蟋座，但在 14 個物種時確實飽和（了）。重要的是脚牍，不一致率（相對于對完整數(shù)據(jù)集的分析）非常低：在子樣本系統(tǒng)發(fā)育中被發(fā)現(xiàn)處于選擇狀態(tài)的基因向臀，在完整系統(tǒng)發(fā)育中被發(fā)現(xiàn)為中性基因（即接受無效假設）的比例小于1%（補充圖S12B）。

基因表達水平的進化分布預測基因功能

OU 過程在最初被提議用于對果蠅的表達進化進行建模時被認為具有吸引力诸狭，因為它能夠區(qū)分中性選擇和穩(wěn)定選擇券膀。鑒于我們發(fā)現(xiàn)大多數(shù)哺乳動物基因處于穩(wěn)定選擇狀態(tài)，我們接下來探索了 OU 模型估計基因表達水平穩(wěn)定分布的能力驯遇，我們推斷這是對進化最優(yōu)表達分布的估計芹彬。因此，我們研究了使用 OU 模型的“進化方差”作為每個組織中基因表達的進化約束程度的定量測量叉庐。與上述相同的jackknifing程序（結果）表明舒帮，OU 模型的估計進化方差對子采樣具有高度魯棒性，這由在估計子采樣系統(tǒng)發(fā)育的方差時非常低的均方誤差 (MSE < 0.005) 確定（補充圖 S12C）陡叠。事實上玩郊，當使用來自少于六個物種的數(shù)據(jù)時，我們發(fā)現(xiàn)進化方差比樣本方差更加穩(wěn)健匾竿，這并沒有解釋物種之間的系統(tǒng)發(fā)育關系（補充圖 S12C）瓦宜。
（這個穩(wěn)轿低颉岭妖？個人理解是物種間都會有一個相對物種內更大的尺度，且相關指標比較穩(wěn)定）

我們首先檢查了跨組織的進化變異模式。為了控制每個組織中的樣本數(shù)量昵慌，我們在與跨組織相同數(shù)量樣本匹配的數(shù)據(jù)子集上重新擬合了 OU 進化參數(shù)（補充表 S1）假夺。我們發(fā)現(xiàn)大腦具有最多的低方差基因（約束最多），而睪丸最少斋攀，這與之前對這些組織表達進化速率的估計一致已卷。(圖 2C; 補充圖S13；陳等人淳蔼。2009 年侧蘸；布拉萬等人。2011 年）鹉梨』浒跨組織，方差合理相關（平均 Pearson 的r = 0.70）（補充圖 S14A）存皂。對于跨三個或更多組織表達的基因晌坤，表達水平和方差在體細胞組織中呈中度負相關(Pearson’s r =?0.25)，在27.2%的基因中旦袋，表達最高的組織與方差最低的組織相匹配(1263/4645)（補充圖 S14B骤菠，頂部）；在該分析中進一步包括睪丸導致表達水平和方差之間幾乎沒有相關性（中位數(shù) Pearson's r = -0.006）（補充圖 S14B疤孕，底部）商乎。

接下來，我們使用我們的完整數(shù)據(jù)集檢查了組織內的進化變異模式祭阀。為了避免數(shù)據(jù)源的多樣性引入的偏差截亦，我們沒有嘗試解釋方差的絕對值，而是專注于理解具有較低和較高方差的基因之間的相對關系柬讨。使用基于等級的基因本體論 (GO) 富集測試 ( Eden et al. 2009 )崩瓤，我們發(fā)現(xiàn)進化方差和功能密切相關，這與之前未使用基于系統(tǒng)發(fā)育方法估計方差的比較研究的結果一致 ( 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 )趁啸。當我們測試按表達水平排序的基因富集時强缘，我們沒有發(fā)現(xiàn)相同的富集類別。因此不傅，我們可以依賴進化方差的估計作為表達約束和基因功能的指標旅掂。

我們發(fā)現(xiàn)表達和序列約束之間只有適度的相關性（Pearson's r = -0.25）（圖 2E; 補充方法）。在表達和序列中保守的基因顯著富集了管家過程（FDR < 10 ^-4）（補充表 S3）访娶，并且兩者中不同的基因富集了免疫和炎癥反應（FDR < 10 ^-6）商虐。序列保守但表達不同的基因富含轉錄調節(jié)因子（FDR = 3.1 × 10 ^-5），尤其是那些參與胚胎形態(tài)發(fā)生的基因（FDR = 9.8 × 10 ^-8崖疤；例如秘车，IRX5、HAND2劫哼、NOTCH1）鲫尊。盡管表達水平的較高進化變異可能受環(huán)境、細胞類型組成的變化和遺傳差異的影響沦偎，但我們的分析支持這樣的假設疫向，即沒有蛋白質序列差異的基因調控差異可以解釋物種特異性表型。

使用基因表達的進化分布來預測有害水平

在分析罕見疾病時豪嚎，序列保守通常用于優(yōu)先考慮那些在發(fā)生突變時對罕見疾病更為重要且可能導致罕見疾病的基因突變（Alf?ldi 和 Lindblad-Toh 2013搔驼；Jordan 等人 2015；Richards 等人 2015）侈询。以此類推舌涨，我們假設表達保守性也應該可以預測基因的重要性。事實上扔字，基因的表達水平要么是培養(yǎng)所必需的 ( Hart et al. 2014 )囊嘉，要么是小鼠必需的 ( Georgi et al. 2013 )，要么是人類的單倍體不足 ( Rehm et al. 2015 )）在幾乎所有組織中的進化方差（更高的約束）顯著低于它們的非必要或單倍充足對應物（Wilcoxon 秩和檢驗P值 <0.01）（圖 3A），關系不是由表達水平驅動的（補充圖S15）。

圖 3
基因表達的進化分布有助于識別導致疾病的基因屈扎。（A）必需基因的進化變異性較低。箱線圖顯示了培養(yǎng)中必需基因（頂部）琢蛤、小鼠必需基因（中間）和人類單倍體不足（底部；深灰色）及其非必需或單倍體充足基因的 log（進化方差）（y軸）的分布（淺灰色）在七個組織中的每一個（x軸）：（***）P < 0.001抛虏；(**) P < 0.01博其。( B ) 疾病基因具有較低的進化變異性。箱線圖顯示對數(shù)的分布（進化方差）（y軸）與相關組織（大腦迂猴、心臟）中的高外顯率單基因自閉癥譜系障礙（上）慕淡、先天性心臟缺陷（中）和神經(jīng)肌肉疾病（下）相關（深灰色）和不相關（淺灰色）的基因和肌肉）沸毁。（左）該組織中表達受限的基因（在三個或更少的組織中>5 TPM）峰髓；（右）普遍表達的基因傻寂；(***) P < 0.001；(*) <0.05儿普。（C，D）使用進化分布或 GTEx RNA-seq 分布從肌營養(yǎng)不良患者的 RNA-seq 中識別異持谰螅基因表達的概述眉孩。( C）基于進化分布（左）或 GTEx RNA-seq 分布（右）的兩種評分方法。（D）表格顯示了使用進化數(shù)據(jù)（左）或 GTEx RNA-seq 數(shù)據(jù)（右）估計的分布時勒葱，所有肌營養(yǎng)不良患者的顯著異忱送簦基因數(shù)、-log₁₀FDR 評分和DMD顯著性等級凛虽。

我們接下來檢查了三種環(huán)境中疾病基因的差異：與非綜合征性自閉癥譜系障礙 (ASD)（大腦）（Banerjee-Basu 和 Packer 2010）死遭、先天性心臟缺陷（心臟）（Amberger 和 Hamosh ）直接相關的罕見單基因疾病基因2017 年；Blake 等人 2017 年）和神經(jīng)肌肉疾部（骨骼佳教丁）（補充方法；Cummings 等人 2017 年）至非。在每種情況下钠署，具有組織限制表達的疾病基因（在三個或更少的組織中TPM > 5）（補充方法；補充圖 S16）荒椭。在疾病相關組織中始終表現(xiàn)出比非疾病基因更低的方差（P值 <0.05） (圖 3B; 補充圖S17）谐鼎。僅在 ASD 相關基因中，我們還觀察到普遍表達的疾病基因與非疾病基因的差異顯著降低（圖 3B）趣惠。

接下來狸棍，我們假設每個基因的最佳 OU 分布的參數(shù)可以通過突出患者數(shù)據(jù)中的異常值、可能致病的基因表達水平來預測疾病基因味悄。這類似于通過使用保守性來識別全外顯子組測序中假定的致病序列突變來發(fā)現(xiàn)致病基因（Choi et al. 2009 ; O'Roak et al. 2011）草戈。為此，我們獲得了臨床診斷為神經(jīng)肌肉疾病的 93 名患者肌肉活檢的 RNA-seq（方法侍瑟；補充表 S4）猾瘸。對于每個患者樣本，我們計算每個基因的z值丢习，以評估它們如何偏離骨骼肌中該基因表達的(最佳)進化適合度牵触，并對多重假設檢驗進行校正（方法；補充圖 S18A）咐低。與來自 184 名健康人的 GTEx 肌肉樣本（The GTEx Consortium 2013）相比揽思，通過這種測量，患者的總體失調基因平均多 3.2 倍（Wilcoxon 秩和檢驗P值 = 2 × 10 ^-9）（補充圖S18B）见擦，表明OU模型擬合的進化參數(shù)可用于檢測更可能有害的異常值表達值钉汗。

然后羹令，我們測試了 OU 模型是否可用于識別罕見病分析中的致病基因。作為原則證明损痰，我們關注了肌肉疾病隊列中臨床診斷為 Becker 或 Duchenne 肌營養(yǎng)不良癥的 8 名患者的子集福侈，包括通過免疫印跡確認肌營養(yǎng)不良蛋白缺失或減少（Cummings 等人，2017 年）卢未。為了將我們的方法與標準差異表達分析進行比較肪凛，我們通過離群值表達對基因進行排序，Z分數(shù)基于（1）與從我們的進化數(shù)據(jù)估計的均值和方差進行比較辽社；或 (2) 與僅從健康 GTEx 人類數(shù)據(jù)估計的均值和方差進行比較（圖 3C）伟墙。根據(jù)我們的進化數(shù)據(jù)，在每位患者中被列為顯著異常值的基因較少（中位數(shù)：4滴铅，范圍：0-32）戳葵，并且DMD基因在 8 名患者中的 6 名中被列為最高或第二最異常表達的基因，每個都顯示出顯著表達不足 (FDR < 10 ^-3 ) (圖 3D)汉匙。相比之下拱烁，參考 GTEx 表達數(shù)據(jù)的評分并未產(chǎn)生如此具體的結果：異常值的基因數(shù)中位值為14.5（范圍：0-250），8 名患者中只有 4 名被稱為DMD顯著低表達（FDR < 10 ^-3） , 及其在這些患者中的重要性排名在 1 到 50 之間噩翠。這種特異性差異可能反映了在使用進化數(shù)據(jù)與人類數(shù)據(jù)時對健康（可容忍）方差的更準確估計：盡管對平均表達的估計在兩種方法之間高度一致（補充圖. S19邻梆，左），當使用進化數(shù)據(jù)集估計時绎秒，表達方差幾乎總是更大（補充圖 S19浦妄，右），反映了在較長時間周期上见芹，每個基因必須充分探索生理上可接受的表達水平空間剂娄。因此，使用人類 GTEx 分布導致更多的假陽性基因似乎異常表達玄呛。相反阅懦，使用 OU 模型對最佳基因表達的進化均值和方差的估計有助于檢測實際疾病基因的失調，并有助于發(fā)現(xiàn)新的疾病基因徘铝。重要的是耳胎，與患者和健康對照之間差異表達的方法相比，這種方法不需要對照群體惕它，并且可以針對個體患者樣本進行怕午。

多變量 OU 模型可用于檢測譜系特異性表達變化

最后，我們探索了使用 OU 模型來檢測基因表達中的方向選擇淹魄。我們使用了模型的擴展郁惜，該模型通過將表達水平的分布建模為多元正態(tài)分布來解釋單個系統(tǒng)發(fā)育中的多種選擇方案，其均值和方差是針對每個（預定義的）子進化枝估計的甲锡。(圖 4A; 巴特勒和金 2004兆蕉；羅爾夫斯等人羽戒。2014 年）。該擴展 OU 模型的先前應用確定了哺乳動物系統(tǒng)發(fā)育中超過 9000 種表達變化（Brawand 等人虎韵，2011 年）易稠，但該分析依賴于較小的系統(tǒng)發(fā)育，因此專注于識別基因表達的物種特異性變化包蓝，不幸的是驶社，這種變化可能是容易被環(huán)境原因或技術影響所混淆。

圖 4
多變量 OU 過程可以檢測譜系特定的表達變化养晋。( A ) 多元 OU 過程衬吆。在多變量 OU 過程下梁钾，隨著時間（x軸）的模擬表達軌跡（y軸）绳泉。灰色軌跡是從同一分布 (N ₀ ) 跨時間采樣的，而橙色軌跡從相同的祖先分布 (N ₀ ) 開始姆泻，但在物種形成事件后在新分布 (N _{1 ) 下演變零酪。}( B ) 表達進化的三個測試假設：從左到右：單變量OU _all模型，其中基因表達在整個系統(tǒng)發(fā)育的單一穩(wěn)定機制下進化（黑色）拇勃，以及兩個多變量 OU 模型四苇，OU_{靈長類動物}和 OU_嚙齒動物，其中基因表達在祖先機制下進化（黑色）和指定的新機制子進化枝（橙色）方咆。( C ) 肝臟中的譜系特異性表達月腋。肝臟樣本中的成對均方表達距離 ( y軸）在參考物種（標記為黑點）和其他哺乳動物之間，用于分配給三個測試 OU 模型中的每一個的基因瓣赂。（黑點）在祖先分布下進化的物種榆骚；（標記為橙色點）譜系分裂后在新制度下進化的物種；（實線）非線性回歸擬合在祖先分布下進化的物種煌集。（D）為譜系特異性表達豐富的示例過程妓肢。熱圖顯示來自基因（列）的列歸一化表達（紅色：高；藍色：低）苫纤，在三個富集的 GO 類別（FDR < 0.05）中具有譜系特異性表達模式：肝臟中的脂質轉運（左）碉钠，肝臟中的免疫調節(jié)（中）和睪丸中的微管運動（右）。

我們利用我們更全面的系統(tǒng)發(fā)育覆蓋范圍卷拘，專注于檢測三個或更多物種的整個亞進化枝的方向和幅度一致的表達變化喊废，它們的樣本被收集并跨多個來源測序，以減輕非遺傳混雜因素栗弟。我們根據(jù)Butler 和 King (2004)建議的方法確定了整個樹的“差異基因表達”（方法）：對于每個基因操禀，我們應用了標準的單變量 OU 模型，該模型對所有物種使用單一最優(yōu)值横腿，以及擴展 OU 模型颓屑，該模型使用兩個最優(yōu)值——一個用于祖先分布斤寂，一個用于內部分布感興趣的進化枝——并使用擬合優(yōu)度測試分配了最佳 OU 模型。作為保守措施揪惦，我們只保留了那些在亞進化枝之間也發(fā)生至少兩倍變化并且在其中一個亞進化枝中平均表達水平至少為 1 TPM 的基因遍搞。

我們首先評估了這種方法在增加系統(tǒng)發(fā)育距離時檢測譜系特異性表達的能力，方法是測試所有靈長類動物（分支長度 = 0.121）器腋、嚙齒動物（分支長度 = 0.177）溪猿、勞亞動物（分支長度= 0.407）或兔形目動物（分支長度 = 0.575）。我們構建我們的數(shù)據(jù)集纫塌，以便我們針對八個比較物種測試感興趣的進化枝內三個物種之間的共享差異表達變化诊县，并且為了避免批次效應的混雜因素，我們選擇三個感興趣的物種措左，以便獲得每個物種的數(shù)據(jù)來自不同的來源（補充方法依痊；補充表 S1）。正如預期的那樣怎披，我們發(fā)現(xiàn)檢測到的差異表達基因的數(shù)量隨著感興趣的進化枝內距離的增加而單調減少胸嘁，范圍從靈長類進化枝中的 470 個基因到兔形進化枝中的 327 個（補充圖 S20A）。不幸的是凉逛，我們通過改組物種表達數(shù)據(jù)（方法）估計的這種分析的錯誤發(fā)現(xiàn)率從 54% 到 78% 不等性宏。

為了提高我們檢測差異表達的能力，我們轉向我們的完整數(shù)據(jù)集并測試靈長類動物進化枝（OU_{靈長類動物}状飞，5-7 靈長類動物與 8-10 個比較物種）和嚙齒動物進化枝（OU_嚙齒動物毫胜，3-5嚙齒動物與 10-12 個比較物種）（圖 4B）。即使有了這個更大的系統(tǒng)發(fā)育诬辈，我們的 FDR 范圍從 18% 到 52%酵使，這表明未來對差異表達的研究應該依賴于更大的系統(tǒng)發(fā)育。不同組織的不同富集大小可能很大程度上是由樣本大小的差異驅動的自晰。當使用與跨組織樣本大小匹配的數(shù)據(jù)集時凝化，差異表達的表達基因的百分比在組織之間是相當一致的（補充圖 S20B）。

盡管我們的分析能力有限酬荞，但我們能夠在肝臟（靈長類動物：FDR = 0.18搓劫；嚙齒動物：FDR = 0.27）和睪丸（靈長類動物：FDR = 0.29；嚙齒動物：FDR = 0.29）中實現(xiàn) FDR < 30%與肺（FDR = 0.26）和腦（FDR = 0.18）的靈長類進化枝一樣（補充圖S21）混巧，我們進一步檢查了這些差異表達基因組枪向。例如，在肝臟中咧党，我們分別在靈長類動物和嚙齒動物中鑒定了 640 個和 794 個具有譜系特異性表達變化的基因秘蛔，突出了每個進化枝中調節(jié)不同的特定代謝過程。僅當存在進化枝特異性選擇時，譜系特異性基因的表達水平才顯著偏離預期（圖 4C）深员。由于與以前的應用相比负蠕，差異表達基因的集合更大，我們可以識別譜系特異性基因之間的功能富集（補充表 S5）倦畅。我們發(fā)現(xiàn)與肝臟中許多脂質代謝過程相關的基因的靈長類動物特異性下調（FDR = 1.88 × 10 ^-11）遮糖。這些過程包括過氧化物酶體功能（FDR = 2.45 × 10 ^-8）、脂肪酸代謝（FDR = 1.52 × 10 ^-8）和脂質轉運（FDR = 3.36 × 10 ^-3）（圖 4D)叠赐，并含有已知的脂質代謝調節(jié)劑欲账，例如 LDL 受體 ( LDLR ) ( Brown and Goldstein 1979 )、肝脂肪酶 ( LIPC ) ( Guerra et al. 1997 ) 和轉錄因子PPARA ( Kersten 2014 )芭概。因此赛不，多種途徑的表達可能在祖先靈長類動物分支上有所不同，這與沒有進一步基因改造的小鼠無法很好地模擬人類脂質血癥的觀察結果一致（von Scheidt et al. 2017）罢洲。在其他例子中踢故，參與調節(jié)免疫反應的基因在嚙齒動物肝臟（FDR = 6.97 × 10 ^-4）和基于微管的運動基因（FDR = 2.82 × 10^-3 ) 和精子發(fā)生 (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 ) 的快速進化奏路，分別畴椰。

討論

在這里臊诊，我們將跨哺乳動物的比較基因表達譜的大型數(shù)據(jù)集與系統(tǒng)分析相結合鸽粉，并表明一對一哺乳動物直系同源物的基因表達差異在進化時間中飽和，這可以通過 OU 過程很好地建模. 我們進一步展示了如何使用 OU 模型來查詢基因功能抓艳，評估候選疾病基因的有害表達水平触机，以及識別譜系特異性的表達進化。

與任何比較物種分析一樣玷或，批次效應引入的偽影儡首，或直系同源物或轉錄本注釋中的錯誤可能會使我們的數(shù)據(jù)產(chǎn)生偏差。然而偏友，我們觀察到的進化模式不僅僅來自單個批次蔬胯，因為每個數(shù)據(jù)源都包含跨越整個系統(tǒng)發(fā)育樹的物種。此外位他，我們僅使用一對一的哺乳動物直系同源物有助于減輕直系分配中的錯誤氛濒，并且我們確認我們的注釋轉錄本的序列身份在進化時間內呈線性發(fā)散，因此不會驅動觀察到的表達進化模式鹅髓。

表達進化的非線性模式由先前提出的 OU 過程 ( Hansen 1997 ) 準確建模舞竿，該模型結合了中性漂移和穩(wěn)定選擇。盡管我們發(fā)現(xiàn)穩(wěn)定選擇在哺乳動物譜系中的表達進化中起主導作用窿冯，但我們注意到適當?shù)哪Ｐ腿Q于進化距離：在靈長類動物譜系中骗奖，我們確實發(fā)現(xiàn)表達差異接近線性，證實了原始研究提出了表達進化的中性模型（Enard et al. 2002 ; Khaitovich et al. 2004）；但在更大的進化距離上执桌，穩(wěn)定選擇的影響越來越大鄙皇，正如最近在哺乳動物中進行的 RNA-seq 研究中所指出的那樣。布拉萬等人仰挣。2011 年）和果蠅（Bedford 和 Hartl 2009 年育苟；Nourmohammad 等人 2017 年）。

重要的是椎木，盡管 OU 過程準確地模擬了我們的數(shù)據(jù)违柏，但還有其他因素可能會限制基因的表達。這些因素包括 (1) 定義的基因表達下限（即 0 TPM）香椎；(2)基因表達的上限漱竖；(3) 在一個組織中對一個基因的選擇壓力對同一基因在其他組織中的表達有限制作用；(4) 對一個形成反式的基因的選擇壓力對其他基因的表達限制的影響畜伐。盡管案例 (1) 和 (2) 代表了我們研究中測試的所有基因的一小部分馍惹，尤其是使用我們的過濾器過濾表達 > 5 TPM 的基因，我們無法分離由于間接選擇力而表達受限的基因玛界，如案例（3）和（4）万矾，來自使用我們當前數(shù)據(jù)的直接選擇壓力下的表達水平。盡管如此慎框，OU 模型仍然是描述基因表達進化歷史的重要定量工具良狈，并為進一步探究表達進化的機制奠定了基礎。

盡管以前使用 OU 模型和基因表達數(shù)據(jù)的研究集中在表達進化的理論方面笨枯，但我們現(xiàn)在展示如何使用 OU 模型來估計最佳基因表達水平的分布薪丁，并回答有關基因功能的生理和臨床相關問題，包括從個體患者數(shù)據(jù)中檢測有害表達水平馅精。因為我們的數(shù)據(jù)來自可能影響我們對進化分布估計的準確性的各種來源严嗜，所以我們小心地只分析相對進化方差（例如，基于等級的 GO 富集測試）或構建我們數(shù)據(jù)集的子集進行分析避免批次效應（例如洲敢，測試多個物種的表達變化漫玄，每個物種都從不同的來源收集）。然而压彭，當直接使用進化均值和方差估計進行疾病表達分析時睦优，我們發(fā)現(xiàn)我們的估計在明確識別疾病基因的異常表達方面優(yōu)于健康人類數(shù)據(jù)。這表明進化估計對技術差異具有穩(wěn)健性哮塞，并且最終可以跨多種組織類型提供刨秆，以幫助科學和臨床發(fā)現(xiàn)。

我們最終應用了一個多元 OU 模型（Butler and King 2004忆畅；Rohlfs et al. 2014) 以識別我們數(shù)據(jù)中跨物種進化枝的譜系特異性表達變化衡未。盡管我們改進了以前依賴較小數(shù)據(jù)集的研究尸执，即使有 17 個物種，我們也無法達到低于 18% 的 FDR缓醋，并且我們發(fā)現(xiàn)在更遠的進化枝之間共享的差異表達變化更少如失。這表明未來研究譜系特異性表達的研究將需要從更大的系統(tǒng)發(fā)育中取樣的數(shù)據(jù)，并且在檢測密切相關物種之間共享的譜系特異性表達變化時表現(xiàn)更好送粱。我們還注意到褪贵，這種方法的一個缺點是必須首先構建經(jīng)過檢驗的假設（例如，靈長類動物與非靈長類動物抗俄；嚙齒類動物與非嚙齒類動物）脆丁，并且最佳擬合模型可能無法真正反映潛在的進化歷史。

展望未來动雹，我們預計 OU 模型可以進一步開發(fā)用于其他生物學查詢槽卫，例如，測試跨基因或旁系同源家族途徑的穩(wěn)定選擇胰蝠，估計祖先表達狀態(tài)歼培，或應用于其他轉錄區(qū)域，如短或長非編碼 RNA茸塞。正如我們的分析所示躲庄，在不同的發(fā)育和環(huán)境背景下，其他組織類型和物種的表達特征將為基因表達的進化以及基因型和表型之間的關系提供更大的動力和進一步的洞察力钾虐。

方法

用于進化數(shù)據(jù)集的 RNA-seq

來自狗和兔組織的 RNA 樣品購自 Zyagen噪窘。來自雪貂組織的 RNA 樣本是 John Englehardt（愛荷華大學）的禮物。來自負鼠組織的 RNA 樣本是來自 Paul Samollow (Texas A&M) 的禮物禾唁。來自犰狳組織的 RNA 樣本是 Jason Merkin 和 Christopher Burge (MIT) 的禮物效览。所有組織收集均經(jīng) IACUC 批準无切，并按照各自的機構指南進行荡短。RNA-seq 文庫的制備如Levin 等人所述。（2010）（補充方法）哆键。樣本在 Illumina HiSeq 2000 測序儀上進行測序掘托，最小深度為 3500 萬次讀數(shù)。

比對和表達量化

RSEM v1.2.12 ( Li and Dewey 2011 ) 與默認參數(shù)一起使用籍嘹，以將讀數(shù)與每個物種的轉錄組對齊并量化每個基因的 TPM闪盔。

基因表達值的標準化

使用來自 Bioconductor package edgeR ( Robinson et al. 2010 ) 的 TMM 歸一化 ( Robinson and Oshlack 2010 ) 對基因表達值 (log ₁₀ TPM) 進行歸一化。簡而言之辱士，TMM 標準化假設大多數(shù)基因在樣本之間沒有差異表達 (DE)泪掀，并估計一對樣本之間的比例因子，使得對數(shù)表達比的修剪平均值（M 值的修剪平均值 [TMM]）相等到 1. 假設物種對之間的大多數(shù)基因不是 DE 是合理的颂碘，因為即使在遙遠的哺乳動物（例如人類和負鼠）之間异赫，給定組織中表達水平的 Pearson 相關性也大于 0.75。

擬合 OU 工藝參數(shù)

BM 和 OU 模型使用帶有默認參數(shù)的 R 包ouch ( Butler and King 2004 )擬合歸一化表達式值。使用比較 OU（替代假設）與 BM（零假設）模型的似然比檢驗計算每個基因的*P值塔拳，然后使用 Benjamini-Hochberg FDR 程序（ *Benjamini 和 Hochberg 1995）對多假設檢驗進行校正鼠证。

神經(jīng)肌肉疾病數(shù)據(jù)集的樣本

本研究中描述的神經(jīng)肌肉疾病患者 RNA-seq 隊列是Cummings 等人描述的隊列的超集。（2017 年）（dbGaP 加入 phs000655.v3.p1）和另外 30 名患者靠抑。組織是根據(jù)機構審查委員會 (IRB) 批準的協(xié)議在國家神經(jīng)疾病和中風研究所 (協(xié)議 #12-N-0095)量九、紐卡斯爾大學 (CF01.2011)、波士頓兒童醫(yī)院 (03-12-205R) 采購的颂碧、倫敦大學學院 (08ND17)荠列、加州大學洛杉磯分校 (15-001919) 和圣裘德兒童研究醫(yī)院 (10/CHW/45)≡爻牵患者在活檢前就診時同意這些方案弯予。如Cummings 等人所述收集患者肌肉活檢。（2017）并使用與 GTEx 項目相同的協(xié)議進行測序（補充方法; GTEx 聯(lián)盟 2013 年）个曙。

人體肌肉數(shù)據(jù)的對齊和表達量化

GTEx BAM 文件從登錄 ID phs000424.v6.p1 下的 dbGaP 下載锈嫩，并在使用 Picard SamToFastq 轉換為 FASTQ 文件后重新對齊。使用 hg19 作為基因組參考垦搬，使用 STAR 2-Pass v.2.4.2a (Dobin et al. 2013 )比對患者和 GTEx 讀數(shù)呼寸。使用 GENCODE v19 注釋使用RNA-SeQC v1.1.8 ( DeLuca et al. 2012 ) 對表達進行量化。

檢測患者樣本中的異常值表達

首先通過 TMM 標準化（Robinson 和 Oshlack 2010 ）將基因表達值（log ₁₀ TPM）標準化為來自進化數(shù)據(jù)集的人類骨骼肌表達值猴贰。對于每個患者樣本中的每個基因对雪，使用從進化數(shù)據(jù)估計的漸近平均值和方差計算Z分數(shù)。Z分數(shù)僅針對被評估為在 OU 而不是 BM 模型下更適合的基因計算米绕，并且其漸近平均值估計為 5 TPM 或更高瑟捣。Z分數(shù)轉換為P-值，然后使用 Benjamini-Hochberg FDR 程序對多個假設檢驗進行校正栅干。我們使用 0.01 的 FDR 閾值來初步定義顯著性迈套。其中，我們刪除了另外 330 個基因碱鳞，這些基因在超過 25% 的 GTEx 樣本中得分為顯著異常值桑李。作為比較器，還使用從健康人類 GTEx 樣本估計的樣本均值和方差計算Z分數(shù)窿给。為了確保這兩種方法之間的可比性贵白，我們只計算了在上述進化方法中的任何步驟都沒有被過濾掉的基因的Z分數(shù)。

檢測譜系特異性表達程序

在所有譜系特異性差異表達分析中崩泡，使用將每個 OU 模型與 BM 模型進行比較的似然比檢驗計算P值禁荒，然后使用 Benjamini-Hochberg FDR 程序針對多重假設檢驗進行調整。對于每個基因角撞，在所有模型上計算 Akaike 和貝葉斯信息標準（AIC 和 BIC）得分呛伴，這些模型對 null 具有顯著性寥掐，以確定最佳擬合模型。在所有情況下磷蜀，這兩個分數(shù)都同意最佳擬合模型召耘。為了估計 FDR，我們使用改組的表達數(shù)據(jù)在每個組織中執(zhí)行相同的程序褐隆，其中來自一個物種的表達數(shù)據(jù)被隨機重新分配給不同的物種污它。（在相關方法中 [ Brawand et al. 2011 ; Rohlfs and Nielsen 2015 ]，Q 值是通過針對零假設 θ _subclade = θ _ancestral直接測試備擇假設 θ _subclade ≠ θ _ancestral得出的庶弃，并針對多個假設檢驗進行調整衫贬。然而，我們發(fā)現(xiàn)這種嚴格的方法導致多重測試的負擔更大歇攻。）我們分別對每組上調和下調的基因進行了 GO 富集分析固惯，使用平均表達至少為 1 的背景基因組適當組織中所有物種的 TPM。

數(shù)據(jù)訪問

本研究的 RNA-seq 原始數(shù)據(jù)（來自兔子缴守、狗葬毫、雪貂和負鼠）已提交給 NCBI 基因表達綜合系統(tǒng)（GEO；http: //www.ncbi.nlm.nih.gov/geo/ ）編號GSE106077屡穗。每個組織環(huán)境中所有一對一哺乳動物直系同源物的處理表達數(shù)據(jù)和進化表達分布可在https://portals.broadinstitute.org/evee/獲得贴捡。

end