基因組預測中先驗分布的兩種超參數設定策略比較

董林松,方 銘,王志勇

(集美大學農業部東海海水健康養殖重點實驗室,福建 廈門 361021)

基因組選擇(GS)是一種通過全基因組范圍內的標記進行基因組估計育種值(GEBV)并加以選擇的方法[1].傳統的最佳線性無偏預測(BLUP)方法可以結合系譜和表型信息來預測個體的育種值[2];相比簡單的表型選擇方法,BLUP方法可以提高育種值估計的準確性.然而,系譜記錄無法準確捕獲等位基因在個體間傳遞的具體信息.換言之,傳統方法無法很好地估計出個體育種值中所包含的孟德爾抽樣誤差項[3].與傳統BLUP方法相比,基因組選擇方法捕獲了個體的標記基因型信息,可以更加準確地了解到個體之間的親緣關系,從而提高GEBV的準確性[4-6].基于其較高的GEBV準確性,并且可以進行早期選種[7],目前基因組選擇已經成為奶牛育種的標準方法,并在其他動植物育種中被廣泛應用[8].但由于基因組的標記數量往往多于參考群個體數量,所以基因組預測在算法上存在挑戰,即所謂的“大P小N”的問題[9].

為解決該問題,近年來出現大量的關于基因組預測的方法,其中較為典型的是Meuwissen等[1]在2001年提出的嶺回歸BLUP(RRBLUP)、BayesA和BayesB方法,在此基礎上又有學者提出基因組BLUP(GBLUP)[10]、GBLUP-最小絕對值收縮與選擇操作(LASSO)[11]、性狀特異的關系矩陣BLUP(TABLUP)[12]、BayesCπ[13]、BayesLASSO[14-15]、Bayes-隨機搜索變量選擇(SSVS)[16]、快速BayesB(FBayesB)[17]、期望最大化BayesB(emBayesB)[18]和基于帕雷托法則的混合正態分布(MixP)[19-20]等方法.這些方法都是在一些先驗假設下進行推導的,RRBLUP和GBLUP假定基因組內所有標記效應值的方差相等,而Bayes方法則假定標記的方差(或效應)服從某一分布,如逆卡方分布[1,13]或拉普拉斯(雙指數)分布[15,17]等.但不論何種分布,都需要在計算之前對超參數進行設置.朱波等[21]已經證明當BayesA的先驗分布超參數取不同數值時,GEBV的預測準確性會產生變化,說明設定合理的超參數值是非常重要的.同時,由于部分研究者對先驗超參數的設定產生一些誤解,因此有必要對基因組選擇算法中先驗分布的超參數設定策略加以探討.

本研究提出基因組選擇算法的2種設定先驗超參數的策略,并利用QTLMAS2012的模擬數據,采用7種基因組選擇算法來探討和驗證超參數設定策略的可行性.

1 計算與實驗方法

1.1 計算模型

基因組預測模型分為兩個水平,第一個水平是對觀測值y的剖分,基因組預測模型的第二個水平是對單核苷酸多態性(SNP)效應或其方差進行解釋.對于表型值的剖分可以用如下混合模型表示:

y=Xb+Za+e.

其中，π是SNP效應值方差為0的概率，當π取值為0時即是BayesA方法[1].BayesCπ的先驗分布與BayesB相似，區別在于所有非零效應SNP的方差被認為是相等的，并且π值可以通過程序進行更新，而不是一個固定值[13].FBayesB是2009年由Meuwissen等[17]提出的一種基于迭代條件期望(ICE)算法的快速Bayes方法，該方法同樣認為只有一部分SNP存在效應，并且這部分SNP效應值的方差相等，但SNP效應服從拉普拉斯分布：

1.2 加性遺傳方差與σ2的關系

這里先對單個SNP位點進行研究，然后將問題擴展到所有SNP位點.假設基因組中的某一SNP位點存在2種等位基因M和m，就有可能出現3種基因型，即MM、Mm和mm.為便于計算，這里將3種基因型轉化為0，1，2，其中Mm為1，MM和mm分別為0和2.假定該位點基因型符合哈代-溫伯格平衡(HWE)；等位基因M的頻率為pi，m的頻率為qi=1-pi，則該位點各種基因型期望值E(w)和方差V(w)分別為[23]：

E(w)=2qi，

(1)

V(w)=E(w2)-E2(w)=2piqi.

(2)

(3)

根據式(3)可以看出一個位點的期望加性遺傳方差與該位點的等位基因頻率和SNP效應值的方差有關.

現將該問題擴展到所有SNP位點，設所有位點總的期望加性遺傳方差為VA，且每個位點效應值方差都等于σ2，則[24]

故：

(4)

1.3 先驗分布超參數的設定

RRBLUP可通過求解下列混合模型方程組來獲得SNP的效應值[1]：

(5)

式(5)中有兩個未知數，即v和s2.為獲得這兩個參數值，通常需對v任意設定一個參數值，例如設定v等于4.2[13]或5.0[32]等，這樣s2就可以通過式(5)求出.

在BayesB[1]和BayesCπ[13]中，先驗分布假設只有一部分的標記存在效應，而其余標記效應的方差為0.此時s2可以通過如下式獲得[13]：

其中π為一個位點效應值為0的概率.

1.4 先驗超參數的另一種設定方法

鑒于不同研究者對SNP基因型的編碼習慣不同，為統一超參數設定公式，可以在計算參數前對標記基因型進行優化.優化的方法就是將基因型進行如下轉換[17]：

1.5 模擬驗證

本研究采用QTLMAS2012的模擬數據[34]對比2種先驗分布參數設定方案的結果.該群體參考群個體數為3 000，估計群為1 020.共模擬了5條染色體，每條染色體的長度均為100 Mb.基因組共模擬產生10 000個SNP位點，剔除最小等位基因頻率(MAF)小于0.05的位點之后，共保留9 042個有效的SNP標記.基因組中共模擬50個數量性狀基因座(QTL)，每個QTL的效應從gamma(0.42，5.4)抽樣而來，其中0.42為分布的形狀參數，5.4為尺度參數.本研究選取其模擬的第一個性狀(產奶量)進行分析，該模擬性狀的遺傳力為0.35.

利用上述提到的7個預測方法(RRBLUP、BayesA、BayesB、BayesCπ、FBayesB、FMixP、MMixP)估計SNP的效應值.在BayesB，FBayesB和FMixP模型中，π值設定為0.99(假定一個QTL大概與周圍兩個SNP發生緊密連鎖).

本研究對SNP基因型采用進行標準化和不進行標準化2種策略.在基于MCMC算法的Bayes方法中，逆卡方分布的自由度v取值為5.0，相應的尺度參數s2的計算和其他Bayes方法的超參數計算公式如表1所示.RRBLUP采用高斯-賽德爾(Gauss-Seidel)迭代法求解，FBayesB和FMixP采用ICE算法獲得SNP的效應值，判斷迭代收斂的標準為：(at-at-1)T(at-at-1)/((at)Tat)<10-8，其中t表示第t次迭代.BayesA、BayesB、BayesCπ和MMixP中，運行Gibbs抽樣20 000次循環，并舍棄前5 000次循環，取后面15 000次的平均值作為SNP的效應值.在BayesB的每個Gibbs抽樣循環中，對SNP效應值的方差通過Metropolis-Hastings算法抽樣100次.估計群體的GEBV通過各個位點的SNP效應值的累加獲得，估計準確性定義為GEBV與真實育種值的相關系數，并作為2種先驗超參數設定結果的比較依據.所有的計算程序均使用Fortran90代碼來運行.

表1 2種SNP基因型編碼策略下的各種預測模型的先驗超參數的計算方法

Tab. 1 Calculations of prior hyper-parameters in various prediction models in the two strategies of SNP genotypes coding

預測模型先驗分布類型策略1策略2BayesAσ2ai～χ-2(v,s2)s2=(v-2)VAvks2=(v-2)VAv∑ki=12piqiBayesBσ2ai～χ-2(v,s2)s2=(v-2)VAv(1-π)ks2=(v-2)VAv(1-π)∑ki=12piqiBayesCπσ2ai～χ-2(v,s2)s2=(v-2)VAv(1-π)ks2=(v-2)VAv(1-π)∑ki=12piqiFBayesBai～Laplace(0,λ)λ=2(1-π)kVAλ=2(1-π)∑ki=12piqiVAFMixPai～N(0,σ21)或～N(0,σ22)σ21=(1-π)VAπkσ22=πVA(1-π)kσ21=(1-π)VAπ∑ki=12piqiσ22=πVA(1-π)∑ki=12piqiMMixPσ2ai～χ-2(v,s21)或～χ-2(v,s22)s21=(v-2)(1-π)VAvπks22=(v-2)πVAv(1-π)ks21=(v-2)(1-π)VAvπ∑ki=12piqis22=(v-2)πVAv(1-π)∑ki=12piqi

注：策略1為SNP基因型進行標準化，策略2為基因型不進行標準化.

2 結 果

通過各種預測模型獲得的2種基因型編碼策略下的預測結果如表2所示，從結果中可以看出，不同方法間的GEBV預測準確性相差較大，其中預測準確性較高的為4種基于MCMC算法的Bayes方法(其中MMixP和BayesCπ預測的準確性相對較高，而BayesB相對較低)，而2種快速Bayes方法的預測結果與RRBLUP相似.在BayesCπ結果中，π值(無效應的SNP比例)在SNP基因型進行標準化和不進行標準化下的預測值非常接近，分別為98.22%和98.24%.基于該結果，本試驗中重新將BayesB中的π值設定為98.2%，再次運行BayesB程序，在SNP基因型標準化和非標準化下獲得的GEBV準確性結果分別為0.783和0.780，較之前有所提高.Usai等[34]在使用相同的模擬數據下，通過RRBLUP、BayesA、BayesB和BayesC方法獲得的GEBV預測準確性分別為0.74，0.79，0.79和0.79，這個結果與本試驗的結果比較吻合，出現稍許差異的原因可能是因為本研究剔除了MAF小于0.05的SNP位點，導致參考群中可用于分析的SNP數量減少.但是在任意一種計算方法下，2種參數設定策略獲得的GEBV預測準確性差異都非常小，幾乎可以忽略，說明這2種參數設定策略都是可行的.

表2 兩種SNP基因型編碼策略下的GEBV預測準確性

Tab. 2 GEBV predictive accuracies in the two strategies of SNP genotypes coding

編碼策略RRBLUPBayesABayesBBayesCπFBayesBFMixPMMixP策略10.7340.7800.7710.7870.7360.7360.795策略20.7320.7770.7680.7890.7350.7340.790

注：策略1為SNP基因型進行標準化，策略2為基因型不進行標準化.

3 討 論

本研究使用了7種基因組預測方法來觀察GEBV的預測準確性，總體上，基于MCMC算法的Bayes方法的預測結果要優于基于ICE算法的Bayes方法和RRBLUP方法.在4種基于MCMC算法的Bayes方法中，BayesCπ和MMixP的結果較為相似并且最好，原因可能是因為這兩種方法在計算過程中可以更新π值，從而獲得一個最優的解.從BayesCπ方法預測的π值也可以看出，無效應的SNP比例接近98.2%，而本研究在BayesB中將該比例設定為99%，這可能是導致BayesB準確性略低于BayesCπ和MMixP的原因.重新設定π值之后，再次運行BayesB就提高了GEBV預測準確性，驗證了上述猜測.快速Bayes方法的預測結果低于基于MCMC方法的預測結果，這與Meuwissen等[17]的研究結果吻合，原因可能是因為ICE本身是一種近似算法[17]，此外π在計算時取值不夠合理，因此更加降低了這些方法的預測效果.雖然BayesA和RRBLUP都認為所有位點存在效應，但BayesA認為不同SNP位點效應值的方差是有差異的，因此能夠將無效SNP位點的效應值壓縮至接近0，而RRBLUP受限于其先驗分布的問題，壓縮能力明顯不如BayesA.換言之，BayesA的先驗假設更接近模擬的QTL的分布情況，因此預測準確性高于RRBLUP.

4 結 論

本研究對比了基因組選擇算法中先驗分布的兩種超參數設定策略的結果，闡述了標記效應方差與位點加性遺傳方差的區別，并給出了相應的計算公式.從模擬結果來看，在計算SNP效應之前，基因型是否進行標準化對最終的結果影響很小.但需要注意的是，如果SNP基因型不進行標準化，在計算標記效應的期望方差時，分母是所有SNP位點的基因型方差之和；如果進行標準化，該項變為k，即總的標記位點數.由于基因型標準化可以突出效應較大的SNP位點，故建議先將SNP基因型標準化，然后再設定先驗分布的超參數值.