利用全基因組表達數據對比評價九種遺傳預測模型方法*

徐州醫科大學公共衛生學院流行病與衛生統計學系(221004) 余星皓 劉 兵 曾 平 黃水平

【提 要】 目的 研究稀疏模型(Lasso、ENET、ssLasso、貝葉斯變量選擇回歸模型(BVSR))與多基因模型(線性混合模型(LMM)、貝葉斯稀疏線性混合模型(BSLMM)、狄利克雷回歸模型(DPR))等九種遺傳預測方法在全基因組表達數據中對復雜疾病的遺傳預測表現。方法 通過模擬研究評價每種方法在不同的較大基因稀疏程度和不同的遺傳度下的預測精度,利用乳腺癌數據進行表型預測。結果 模擬結果顯示預測方法在滿足各自的模型假設時表現結果最好。在相同模擬假設情況下,隨著遺傳度的增高,模型的預測準確性也逐漸增高。BVSR運算速度和BSLMM運算速度相似,由于迭代次數的影響,BVSR與BSLMM的運算速度低于LMM。實際的乳腺癌數據顯示BSLMM和DPR的預測精度優于其他方法。結論 BSLMM和DPR在不同模擬情形下和真實數據中均表現出穩健的預測能力,值得在實際應用中推薦。

全基因組關聯研究(genome-wide association study,GWAS)已經發現了上萬個單核苷酸多態性(single-nucleotide polymorphism,SNP)位點與數百種復雜疾病存在關聯,為表型變異的遺傳基礎提供了前所未有的視角[1-7]。GWAS的成功也極大促進了利用遺傳信息(除環境和生活方式信息外)進行復雜疾病風險預測評估的研究[8]。在動物或植物中,具有遺傳標記的準確表型預測可以幫助選擇理想的育種個體,從而提高育種計劃的有效性[9];在人群中,利用遺傳標記的準確表型預測有利于復雜疾病的早期預防和干預、促進開發個性化藥物,從而制定治療方案和預測療效[10]。此外,在跨組學數據分析中,表型預測也被認為是將功能基因組測序研究與GWAS相結合的關鍵步驟,可以在GWAS中構建更高效和具有可解釋性的基因集測試,通過設置SNP權重,將其從預測模型中推斷出來,并用于表達數量性狀位點映射研究[11]。

研究者們發明了許多新的表型預測模型方法以利用全基因組遺傳位點信息提高預測精度。表型預測模型本質上是對基因效應大小分布的假設不同,例如Henderson等人[12]提出的線性混合模型(liner mixed model,LMM),Speed等人[13]提出的最佳無偏估計模型(multiple best linear unbiased prediction,MultiBLUP),Tang等人[14]提出的Bayes Lasso模型。上述方法多屬于參數模型,所假設的遺傳效應分布是高度參數化的、由少數幾個參數決定,其局限在于僅對特定滿足模型假設的疾病預測效果良好,對其它不滿足假設的疾病預測精度不佳。理論上,模型預測精度取決于遺傳效應先驗與真實分布的接近程度。然而,真實的遺傳結構往往未知,由于疾病的不同，在遺傳度、關聯位點個數、最小等位基因頻率和效應大小等方面存在諸多差別。另外,這些方法先前主要利用全基因組SNP位點信息進行復雜疾病的預測;而越來越多的研究表明基因表達在復雜疾病的發生、發展和轉歸中起著十分重要的作用,是SNP對復雜疾病起作用的重要分子中介。因此,通過全基因組表達數據預測復雜疾病具有重要的科學意義。通過基因表達預測復雜疾病與利用SNP預測復雜疾病的差別在于，基因表達數據的變量較少,效應可能更強,作為一個整體的功能單位,基因相對SNP更具有結果的可解釋性。那些通過全基因組SNP位點獲得復雜疾病預測結論未必適合基因表達預測。

本文主要對九種遺傳預測模型方法進行比較,這些方法大致可分為兩類:稀疏模型和多基因模型。前者包括Lasso、elastic net(ENET)、spike-and-slab Lasso GLMs(ssLasso)、貝葉斯變量選擇回歸模型(Bayesian variable selection regression,BVSR);后者包括線性混合模型(LMM)、貝葉斯稀疏線性混合模型(Bayesian sparse linear mixed model,BSLMM)、狄利克雷回歸模型(Dirichlet process regression,DPR)等。本文通過癌癥患者基因表達數據模擬各種復雜疾病可能的遺傳結構,比較上述方法在不同情況下對于復雜疾病的預測能力，進一步通過比較各方法在真實腫瘤數據中的預測能力評價這些方法在基因表達數據中的表型預測精度。

材料與方法

1.方法

比較九種模型(稀疏模型:BVSR[15]、Lasso[16]、ENET[17]和ssLasso,多基因模型:LMM、Linear-BSLMM、Probit-BSLMM[18]、DPR.VB和DPR.MCMC[19])在遺傳預測中的能力,這些模型方法被廣泛運用于人類復雜疾病的表型預測以及動植物繁殖的基因選擇中。具體而言,BVSR假設較少部分的基因變異對復雜疾病或性狀產生效應,單個基因變異對表型的效應較大,假設總體效應大小符合正態分布和零點分布的混合分布;Lasso是在模型的離均差平方和基礎上增加一個絕對值的懲罰項,使其總和小于等于一個常數,是一種壓縮估計方法;ENET是一種Lasso與嶺回歸組合后的回歸分析,在絕對值懲罰項基礎上同時施加平方和懲罰項;ssLasso假設基因的系數服從混合雙指數先驗,根據不同的基因型與表型數據,產生適合的收縮系數,并保留系數較大的基因;LMM假設所有的基因對表型均產生效應,且效應均較小,同時服從正態分布;BSLMM假設基因對表型的效應大小服從兩個正態分布組成的混合分布,兩個正態分布的混合程度決定于不同的混合比例,需要注意的是,針對病例對照研究(假設表型編碼為0-1),BSLMM可直接將0-1當作連續變量通過線性模型處理,稱為Linear-BSLMM,或將Probit鏈接函數通過廣義線性模型處理,稱為Probit-BSLMM,圖1中分別以BSLMM1和BSLMM2表示;DPR假設基因效應大小服從一個無窮混合正態分布,屬于非參數貝葉斯模型范疇,DPR的參數估計可通過傳統的馬爾科夫蒙特卡羅估計,稱為DPR.MCMC,也可通過變分貝葉斯方法估計,稱為DPR.VB。

LMM、BVSR、Linear-BSLMM和Probit-BSLMM方法采用GEMMA[18]軟件執行,DPR.VB、DPR.MCMC等方法采用DPR[19]軟件執行,Lasso、ENET通過glmnet軟件包執行,ssLasso通過BhGLM[20]軟件包執行。

2.模擬數據

為更加接近真實數據,使用916例乳腺癌患者的基因表達數據,模擬三種情境下的表型數據:(i)所有基因均對表型產生效應,分別服從正態分布、t分布或Laplace分布,設置遺傳度分別為0.2、0.5或0.8;(ii)僅有小部分基因(分別為10、50或100)對表型產生效應,服從正態分布,大部分剩余的基因對表型均不產生效應,遺傳度分別為0.2、0.5或0.8;(iii)所有基因對表型均產生效應,效應基因分為兩部分:小部分效應較大的基因和較大部分效應較小的基因,兩組基因的效應大小都服從正態分布,兩部分占遺傳度的比例分別為0.2和0.8;較大效應的基因數量分別為10、50或100;遺傳度分別為0.2、0.5或0.8。每種情景分別重復模擬20次,每一次重復中隨機將數據中90%的個體作為訓練集,剩余的10%作為測試集。通過訓練集數據估計每種方法的模型參數,并在測試集數據中比較不同方法的預測能力。通過表型與預測值的曲線下面積(area under the curve,AUC)評價預測效果。

3.真實數據

真實數據全部來源于癌癥和腫瘤基因圖譜(The Cancer Genome Atlas,TCGA),通過加利福尼亞大學基因組瀏覽器(UCSC Xena,https://xenabrowser.net/)下載乳腺癌(Breast Cancer,BRCA)數據。原始數據包括1215例患者臨床數據和1219例患者的20530個基因表達數據,首先對兩份數據進行合并,刪除重復的患者和男性患者,同時刪除2840個零表達值超過50%的基因,最終獲得916例患者的17690個基因表達數據。在數據分析前,首先將癌癥的基因表達數據分成兩部分:90%的患者為訓練集數據,剩余10%的患者為測試集數據,對數據集執行20次重復。通過AUC評價不同模型的預測能力。

結果

1.模擬數據結果

情景i中(圖1-i)當遺傳度為0.2時,效應大小的分布服從正態分布的情況下,各模型的預測能力接近。效應大小的分布服從t分布和Laplace分布的情況下,多基因模型的預測能力均高于稀疏模型,多基因模型的預測能力大小接近。當遺傳度為0.5或0.8時,三種情況下多基因模型中五種模型的預測準確性相對接近,優于稀疏模型,BVSR預測能力最低。當遺傳度為0.8時,各統計學模型的預測結果AUC均在0.7左右,高于遺傳度為0.2和0.5時的預測結果。

情景ii中(圖1-ii)當遺傳度為0.2或0.5時,稀疏模型的預測能力優于多基因模型,但是,隨著效應基因數量和遺傳度增加,稀疏模型與多基因模型之間的差異逐漸縮小;當遺傳度為0.8:效應基因數量為10的時候,BSLMM和DPR的預測能力較強,其余多基因模型精度稍低。隨著效應基因數量增加,稀疏模型與多基因模型之間的差異越小,當效應基因數量為100時,稀疏模型的預測能力要低于多基因模型。

情景iii中(圖1-iii)當遺傳度為0.2時,三種情況下多基因模型的預測效果均高于稀疏模型,其中BVSR的預測能力最差,linear-BSLMM、probit-BSLMM和DPR.MCMC的預測能力較強且較穩定。當遺傳度為0.5時,多基因模型的預測效果優于稀疏模型,BSLMM和DPR的預測能力較穩定,且精度較高。當較大效應基因數量增加時,Lasso和ENET的預測能力上升;當遺傳度為0.8時,較大效應基因的數量不同時,模型之間預測結果相似。多基因模型的預測結果相似。稀疏模型的預測結果相對不穩。在不同情況下,多基因模型的預測能力均稍高于稀疏模型。各統計學模型的預測結果AUC均在0.75左右,高于遺傳度為0.2或0.5時的預測結果。

圖1 遺傳度0.2(A)、0.5(B)、0.8(C)時三種模擬情景下Lasso(紅)、ENET(綠)、ssLasso(深藍)、BVSR(綠)、LMM(藍)、BSLMM1(Linear-BSLMM)(黃)、BSLMM2(Probit-BSLMM)(紫)、DPR.VB(淺藍)、DPR.MCMC(粉)之間AUC估計值的比較

2.真實數據結果

真實數據來源于上述的乳腺癌數據,將術后淋巴結陽性轉移作為預測的表型數據(編碼為0-1,視為連續變量應用線性模型處理)。經過計算,數據的真實遺傳度為0.15。各方法表型預測結果AUC均在0.61左右(圖2A),多基因模型的預測能力均高于稀疏模型。多基因模型中,五種方法的預測能力幾乎無差異(雖然DPR和BSLMM輕微略高)。稀疏模型中,ENET的預測能力最強,Lasso與BVSR的預測能力接近,ssLasso的預測能力最低。同時,采用Logistic回歸分析各基因的效應大小,并繪制相應的分布密度函數,與正態分布、Laplace分布和自由度為4的t分布相比較(圖2B)。

3.運算時間

最后對比各種方法的運算速度,以乳腺癌數據為例。這里需要說明的是,貝葉斯預測方法的運算時間與所使用的迭代次數有關系,為了節省時間,這里只進行了100次迭代。結果顯示(表1)，Lasso、ENET和ssLasso的運算速度最快,其次為LMM和DPR.VB;BVSR、BSLMM和DPR.MCMC運算速度相近,且時間均較長。

討論

本文所提及的多種統計學方法在遺傳統計學的多個方面都有較好的應用價值，例如估計遺傳度[18]和關聯映射[21]。本文主要致力于多種統計學方法在表型預測方面的研究。國內外已發表的文獻[13,18-19,22-24]主要研究表型預測方法在GWAS中的比較,已取得較大成果,本文通過模擬不同的疾病遺傳結構,比較了兩類統計學方法在基因表達數據中的預測能力。模擬數據分析結果顯示,各模型在表型預測上的表現不同。當模擬假設全部基因對表型均產生效應時,多基因模型的預測能力要高于稀疏模型;當模擬假設只有部分基因對表型產生效應時,稀疏模型的預測能力要高于多基因模型。當模擬假設小部分基因對表型產生較大的效應,其余大部分對表型產生較小的效應時,多基因模型中的BSLMM和DPR預測能力要高于其余模型。在相同模擬假設情況下,隨著遺傳度的增高,模型的預測準確性也逐漸增高。當假設全部基因對表型均產生效應、效應大小服從單個正態分布、遺傳度值較小時,多基因模型與稀疏模型的預測能力相似。當遺傳度較大時,多基因模型的準確性稍高于稀疏模型。當假設部分基因對表型產生效應,且效應大小服從單個正態分布時,隨著相應基因數量的增加,稀疏模型的預測能力逐漸由高于多基因模型逐漸變化為與多基因模型預測能力相似,甚至低于多基因模型。真實數據分析結果顯示不同統計學方法在表型預測的表現不同,各方法表型預測結果表型預測結果AUC均在0.61左右,多基因模型的預測能力均高于稀疏模型,其中BSLMM和DPR最優。由模擬研究和真實數據分析結果可以看出,模型預測的準確性與復雜疾病或性狀的真實遺傳結構有關,當復雜疾病的真實遺傳結構與模型假設相吻合時,模型的預測表現較好,而當復雜疾病的真實遺傳結構與模型假設差異較大時,模型的預測表現較差。這與以往研究的結論相同[13,18]。因此我們認為多基因的模型假設與復雜疾病的真實遺傳結構相符合,少數較大效應的基因和大多數小效應基因共同對復雜疾病產生影響。

圖2 BRCA九種方法預測能力的比較

表1 真實數據集中不同方法的運算時間(小時)

*:Lasso、ENET、ssLasso的調整參數均為通過100倍交叉驗證選擇的最佳值,BVSR、BSLMM、DPR.MCMC均為進行了100次MCMC迭代計算,括號中的值為標準差,均值和標準差都是基于模擬的20次重復計算。

本文研究的基因表達數據相對較小,三種方法的運算速度均較高,但是這些方法往往并不適用于較大數據集[22]。需要注意的是,Lasso、ENET和ssLasso的運算速度同時受到交叉驗證次數的影響,本研究中只通過100倍交叉驗證,并未進一步研究交叉驗證數量對于運算速度的影響。LMM、BVSR和BSLMM使用GEMMA軟件進行分析,同時適用于大規模數據集。BVSR和BSLMM的運算速度相似,由于迭代次數的影響,BVSR與BSLMM的運算速度低于LMM。有研究報道[18],當納入模型的自變量較大時,BSLMM的運算速度將遠高于BVSR。DPR通過DPR軟件進行分析,其中DPR.VB與LMM運算時間接近,而DPR.MCMC運算時間同樣受到迭代次數的影響,遠高于其它方法。此外,計算時間還受到工作環境、計算機語言、以及基因和表型之間的關系等因素的影響。雖然目前對于表型預測方法的研究層出不窮,但是尚未有一種方法對于不同遺傳結構的復雜疾病或性狀都具有較高的預測能力且運算時間較快。同時,遺傳度較低時,所有方法的預測表現均較差。此外,本文所納入的方法均為通過基因表達數據進行復雜疾病或性狀的預測,都未能考慮到在遺傳數據中非常重要的分組結構。例如:SNP可以通過基因或者功能注釋的形式分成不同的組;基因之間存在交互作用和共同作用,基因可以通過通路信息進行分組,功能相同的基因可以被劃分到一組[25]。這些被忽視的遺傳信息(分組結構、非線性效應和交互作用)可能會對預測結果造成一定的影響。已經有研究表明,充分利用多組學的信息可以提高遺傳預測的精度[26],這也將是接下來的研究重點之一。