基于多層感知機的綿羊限性性狀基因組選擇模擬研究

王萬年,陳思佳,郜金榮,溫中豪,袁夢嬌,張洪志,龐志旭,喬利英,劉文忠

(山西農業大學動物科學學院,太谷 030801)

基因組選擇(genomic selection, GS)是利用覆蓋于全基因組范圍內的遺傳標記,結合表型信息對個體育種值進行預測的方法[1]。隨著近年來基因分型成本的降低,基因組選擇在畜禽育種中得到了廣泛應用[2-3]。相對于最佳線性無偏預測(best linear unbiased prediction, BLUP)[4],基因組最佳線性無偏預測(genomic BLUP, GBLUP)模型[5]和一步法基因組最佳線性無偏預測(single-step GBLUP, SSGBLUP)[6]利用了基因型信息,提高了基因組估計育種值的準確性。除BLUP系列方法之外,由于深度學習(deep learning, DL)在許多其他科學領域的應用力量,人工神經網絡(artificial neural network, ANN)模型也開始用于GS[7],多層感知機(multilayer perceptron, MLP)是一種前饋神經網絡模型,利用了多層神經網絡擬合基因組信息和表型信息間的關系計算GEBV。在前人研究中[8],MLP在奶牛生長性狀中有更小的基因組預測誤差和更高的表型值預測精度[8]。

與生長性狀不同,限性性狀指在一個物種的一種性別中表達的性狀。這些性狀由性染色體或僅在一種性別中表達的常染色體基因控制,如單睪、產仔數、產蛋量等。相對于非限性性狀,限性性狀目前只能在單性別中測定表型,在育種中缺少一個性別的選種,選擇效率不高,利用基因組選擇可以實現雙性別都參與選種。GBLUP、SSGBLUP、MLP在非限性性狀基因組選擇中,由于可以使用全部的表型數據和與之對應的基因組數據,GS準確性可以達到較高水準,并且SSGBLUP和MLP優于GBLUP;而在限性性狀中,表型數據只來自于單性別,GBLUP通過偽表型,SSGBLUP結合了系譜信息,MLP通過神經元學習基因組和表型關系均同時使用了來自公畜和母畜基因組數據,這些不同導致估計方法的GEBV準確性可能會存在差異。許多研究表明,對于易于測量性狀和中高遺傳力性狀,BLUP取得了較好的應用效果[9-10],對于限性性狀,BLUP方法應用效果較差[11]。相比于傳統BLUP和GBLUP,使用SSGBLUP可以有效提高GS的準確性[12-14],而MLP法在限性性狀中的GS效果仍需進一步評估。由于在方法原理和使用條件中差異的存在,檢驗MLP在限性性狀中GS的效果是否優秀,對促進限性性狀的選擇具有重要意義。

為了估計育種值,需要對遺傳參數進行估計,3大遺傳參數指遺傳力(heritability,h2)、重復力和遺傳相關,通常使用估計遺傳力(遺傳參數)與當代群體真實遺傳力(遺傳參數)的接近程度衡量遺傳評估效果。近年來,遺傳參數估計的方法被不斷改進,從傳統的方差分析法[15]到適用于混合模型的Henderson方法[16]、最大似然法[17]和約束最大似然法(restricted maximum likelihood, REML)[18]。REML法避免了最大似然估計造成的偏低估計,是使用線性混合模型進行遺傳參數估計的常用方法。相對于REML法,ANN是非線性模型,沒有對SNP標記和表型之間的關系做出先驗假設,而是通過神經元的連接分析兩者之間的潛在關系,為適應輸入和輸出之間的復雜關聯提供了靈活性[19]。前人研究報道中[14],在計算得出EBV/GEBV后,必須還要對其準確性進行評估,通過EBV/GEBV準確性的高低評判各模型的選擇效果。

為了研究MLP在限性性狀GS中的有效性,檢驗MLP在限性性狀GS中的應用效果是否優秀,本研究基于多個綿羊模擬群體進行以下研究:1)比較不同情況下MLP與BLUP、GBLUP、SSGBLUP遺傳參數估計結果;2)預測不同情況下MLP的GEBV準確性,并對MLP在GS中的使用效果進行評估。

1 材料與方法

1.1 數據模擬

使用QMsim軟件[20]模擬3個水平(50、100、500)數量性狀基因座(quantitative trait loci, QTL)數目及3個水平(0.05、0.2、0.5)h2的限性性狀,性狀表型方差標準化為1,對整個模擬過程重復10次,以減少抽樣誤差。由于模擬群體經過幾代之后,當代群體h2與設置值有出入,所以需要計算真實的當代群體h2。為了得到與實際綿羊群體相似的連鎖不平衡(linkage disequilibrium, LD)和突變-漂變平衡[21],構建了2 000個世代的綿羊群體,起始規模3 000頭,群體規模大小不變,群體構建階段,公、母數目保持一致,采用隨機交配。從構建的群體最后一代隨機抽取50頭公羊,1 000頭母羊作為當代的初始群體,分別構建4個世代的當代群體Pop1和5個世代的當代群體Pop2。模擬中,公、母羊淘汰率均為20%,公、母羊增長率均為40%。公、母羊隨機交配,每頭母畜后裔數為1,公母后裔各半。選擇估計育種值高的個體留種,并按年齡進行淘汰。

本試驗模擬綿羊的基因組信息,設置26對常染色體,總長2 584 cM。為了模擬接近真實情況的LD,依據真實染色體長度及數目模擬基因組信息。分別模擬10 000(10K)和50 000(50K)個SNPs標記,假定SNP標記均勻分布在全基因組范圍內,且標記無效應。分別設置QTL數目為50、100和500個,QTL的最小等位基因頻率(minor allele frequency, MAF)為0.1。假定QTL在全基因組范圍內隨機分布,QTL效應值服從伽馬分布(形狀參數為0.4)。為了創建突變-漂變平衡,設置構建群體中QTL和SNP標記的回復突變率為10-5。

在Pop1中選擇第2～3世代作為參考群,第4世代作為候選群;在Pop2中選擇1～4世代作為參考群,第5世代作為候選群。保留Pop1和Pop2的完整系譜,和參考群中所有母畜的表型信息。選取參考群有后裔的公畜和所有母畜進行基因分型,隨機選取候選群的100個個體進行基因分型。本研究只對部分候選群體進行基因分型,更有實踐意義。

1.2 基因組預測模型

(1)

其中,

MME的求解通過DMU[22]軟件完成。

1.2.2 MLP MLP由輸入層(input layer, IL)、隱藏層(hidden layer, HL)和輸出層(output layer, OL)組成。IL接收標記信息,HL分別由n1, n2, …, nj神經元組成,OL由代表神經網絡輸出值的nol神經元組成,輸入經過IL被傳遞給HL1的神經元,然后每個隱藏的神經元產生一個輸出,用作HL2的每個神經元的輸入,以此類推最終得到輸出值yi。在這種結構中,每一層均產生權重w[i]矩陣,其中輸出由最后一個HL的線性組合產生。最終得到的預測值為:

(2)

其中,bi代表第i層的偏置,初始偏置設置為1,w[i]代表第i層權重,fi為激活函數,p代表IL神經元個數,t表示第t個神經元。在本試驗中,采用含3個HL的MLP架構,IL、HL1、HL2、HL3、OL神經元個數依次為64、32、32、16、1。為了防止模型過擬合,在第一隱藏層到第二隱藏層間增加了隨機失活層,訓練時設置為0.2,即保留80%的節點,其余的20%置為0。

本研究中使用修正線性單元(rectifier linear unit, ReLU)激活函數[23],增加了神經網絡各層之間的非線性關系。

(3)

采用均方誤差(mean square error, MSE)作為損失函數,平均絕對誤差(mean absolute error, MAE)作為評估函數,在反向傳播過程中,損失函數En對輸出層輸出yi的求偏導,得到輸出yi與真實值yt的差值δi=yi-yt,若設定zi=ReLU(x),則損失函數對隱藏層輸出ai的偏導為:

(4)

損失函數對網絡權重w[i]的梯度為:

(5)

為了使模型輸出逼近最優值,采用了均方根傳遞(root mean square prop, RMSProp)算法[24]來更新模型的網絡參數。RMSProp除了自適應調節學習率之外,對權重w和偏置b的梯度計算了微分平方加權平均數,用以修正擺動幅度,加快函數的收斂。

E[g2]t=αE[g2]t-1+(1-α)g2

(6)

(7)

(8)

1.3 遺傳參數估計方法

(9)

其中,

P=V-1-V-1X(X-1V-1X)-1XTV-1

(10)

(11)

IAE=λIEM+(1-λ)IAI

(12)

其中,IEM是EM算法的期望最大陣,IAI是AI算法的平均信息陣。λ∈[0,1],設置λ以一定步長增長。直到迭代后的方差組分估計值在參數空間內。

1.3.2 ANN 在ANN中,估計h2和SNP效應,需要獲得標記的相對重要性(relative importance, RI)。Olden等[29]提出了一種方法,即使神經網絡有多個HL,可以使用所有的連接權重來獲得RI。連接權值矩陣的優劣決定了SNP效應估計的準確性,為了計算所有標記的RI向量,將連接權重矩陣相乘。即

RI=w[1]*w[2]*…*w[j-1]

(13)

其中,將w[i=1]作為連接第j-1層到第j層的的估計權值矩陣,j是人工神經網絡的層數。使用如下線性近似[29],通過RI值估計加性SNP效應向量βα,即:

(14)

(15)

1.4 模型評價

1.4.1h2比較 10次重復下,各方法計算當代群體估計h2的均值和標準誤在相同條件下與真實h2進行比較并做顯著性檢驗,h2估值越接近真實h2越準確,h2估值越準確說明遺傳評估效果越好。各方法估計h2與真實h2差異越顯著,說明遺傳評估效果越差。

1.4.2 候選群體EBV/GEBV準確性比較 用EBV/GEBV準確性來比較不同模型的性能。準確性通過EBV/GEBV與QMsim軟件中模擬的TBV之間的相關系數衡量。即

從不同QTL數目、不同標記數與不同遺傳力3個角度,對Pop1和Pop2進行育種值估計,計算其候選群EBV/GEBV準確性。相同條件下,準確性越高說明模型性能越好。

2 結 果

2.1 參考群體和候選群體

2個獨立的當代群體Pop1和Pop2群體結構模擬結果見表1。其中Pop1的參考群個體為3 360個,基因分型個體為1 866個,候選群個體為2 744個,基因分型個體為100個;Pop2的參考群個體為7 154個,基因分型個體為3 351個,候選群個體為3 842,基因分型個體為100個。

表1 模擬數據的當代群體結構Table 1 Structure of recent population obtained from QMsim

2.2 MLP模型評估

如圖1所示,MSE最終達到了0.002 7,MAE最終達到了0.037 2,損失函數和評估函數在訓練600次之后基本均趨于平緩,模型擬合效果良好。

圖1 MLP模型損失函數與評估函數結果Fig.1 MLP model loss function and evaluation function results

2.3 不同情況下各模型的h2估計結果

設定初始群體特定限性性狀h2為0.05時,MLP、SSGBLUP、GBLUP和傳統BLUP的h2估計結果見表2。4種方法的估計值介于0.028～0.075,與當代群體遺傳力0.047～0.051相符。MLP的估計值介于0.044～0.057,相比SSGBLUP介于0.041～0.068范圍更小。在QTL數目和標記數相同的情況下,Pop1(Pop2)中當代群體遺傳力為0.049～0.051(0.047～0.050),與MLP的h2估值0.044～0.057(0.046～0.055)差異不顯著。在Pop2群體中標記數為10K且QTL數為100時,SSGBLUP與MLP差異不顯著,但兩種方法均顯著(P<0.05)優于GBLUP和BLUP;除此之外,MLP的h2估計結果均顯著(P<0.05)優于SSGBLUP、GBLUP和BLUP。

表2 設定初始群體特定限性性狀h2為0.05時各模型h2估計結果Table 2 Estimation results of h2 of each model when the prior population-specific limit trait h2 is set to 0.05

設定初始群體特定限性性狀h2為0.2時,各方法h2估計結果與當代群體h2的對比見表3。在QTL數和參考群大小相同的情況下,SSGBLUP和GBLUP法遺傳力估計結果在不同標記數下差別較大,MLP差別較小。在Pop1群體中標記數為10K且QTL數為50時,MLP的h2估值與當代群體h2差異顯著(P<0.05);除此之外,在QTL數目和標記數相同的情況下,Pop1(Pop2)中當代群體遺傳力為0.196～0.207(0.194～0.207),與MLP的h2估值0.196～0.216(0.185～0.203)差異不顯著。在Pop1中QTL數為500的兩個標記數下和Pop2中QTL數為10且標記數為50K兩種情況下,MLP的h2估值與SSGBLUP差異不顯著,但兩種方法均顯著(P<0.05)優于GBLUP和BLUP;除此之外,MLP的h2估計結果均顯著(P<0.05)優于SSGBLUP、GBLUP和BLUP。

表3 設定初始群體特定限性性狀h2為0.2時各模型h2估計結果Table 3 Estimation results of h2 of each model when the prior population-specific limit trait h2 is set to 0.2

設定初始群體特定限性性狀h2為0.5時,各方法h2估計結果見表4。在QTL數及標記數相同的情況下,4種方法在Pop2群體中的遺傳參數估計結果較Pop1更接近于真值。QTL數為100,Pop1中標記數10K下和Pop2中標記數50K下,MLP的h2估值與SSGBLUP差異不顯著,但兩種方法均顯著(P<0.05)優于GBLUP和BLUP。除此情況之外,MLP的h2估計結果均顯著(P<0.05)優于SSGBLUP、GBLUP和BLUP。

2.4 不同模型預測效果

表5為不同h2初值時,各方法的預測準確性范圍。在QTL數目和標記數相同的情況下,h2初值為0.05時,MLP的GEBV準確性較SSGBLUP和GBLUP明顯提高;h2初值為0.2時,MLP的GEBV預測準確性較SSGBLUP有略微提高,較GBLUP有明顯提高;h2初值為0.5時,MLP的GEBV準確性較SSGBLUP沒有明顯提高,較GBLUP和BLUP優勢明顯。

表5 不同h2初值時各方法預測準確性范圍Table 5 When different prior h2 values, the prediction accuracy range of each method

圖2為不同情況下MLP、SSGBLUP、GBLUP的GEBV準確性及BLUP的EBV準確性,在h2、QTL數和標記數相同的情況下,Pop1和Pop2群體中,MLP的預測準確性均為最高,BLUP均為最低;Pop2群體中各方法的預測準確性較Pop1群體均有提升。在QTL數、標記數和群體相同的情況下,隨著h2降低,MLP較GBLUP和SSGBLUP的GEBV準確性差值變大,優勢逐漸增加。在h2、QTL數和群體相同的情況下,MLP、SSGBLUP及GBLUP模型在標記數為50K時的準確性較標記數為10K時有提高。其中,在QTL數目相同及h2初值為0.05時,Pop1(Pop2)群體中MLP在標記數為10K時預測準確性介于0.269～0.278(0.287～0.300)略高于SSGBLUP在標記數為50K的預測準確性0.224～0.241(0.262～0.279)。在h2、標記數和群體相同的情況下,預測準確性隨QTL數改變并沒有明顯變化。

3 討 論

在綿羊育種中,限性性狀是影響經濟效益的重要性狀,這種復雜性狀的遺傳增益可以有效的提高綿羊的生產效率[30]。除產奶量外,乳房構型、產羔困難和乳腺炎等限性性狀對生產成本和動物福利影響同樣很大,也需要在育種目標中加以考慮[31]。在綿羊中實施基因組選擇最嚴峻的挑戰是品種問題,許多品種個體數少,導致參考群小,連鎖不平衡低,有時甚至缺乏表型記錄[32]。用MLP進行GS[8]已有報道,但是在綿羊限性性狀中還未應用,所以本研究在不同h2、QTL數目、標記數的兩個群體(Pop1和Pop2)中對MLP模型進行了全面的測試。研究發現,在QTL數、標記數和群體相同的情況下,在3個h2水平下,MLP較GBLUP和SSGBLUP優勢更大;隨著h2降低,MLP較GBLUP和SSGBLUP的優勢逐漸增加,表明DL方法的實施為改善綿羊限性性狀的遺傳增益提供了巨大潛力。

遺傳參數估計是家畜育種中的重要組成部分,群體遺傳參數估值的優劣對正確評定種畜育種值和育種規劃有非常重要的作用。MLP是通過神經元權值矩陣計算SNP效應向量從而進行遺傳參數估計;而REML法是對線性模型求解約束似然函數,通過迭代的方式求解其極值點計算方差組分。本研究中,4種方法在Pop2群體中的遺傳參數估計結果較Pop1更接近于真值,表明在估計遺傳參數時MLP與SSGBLUP、GBLUP和BLUP同樣受參考群體大小的影響,參考群體越大,遺傳參數估計結果與當代群體遺傳參數的差值越小[33]。除個別情況之外,MLP的遺傳力估計結果與其他3種方法差異顯著,且在相同h2、QTL和標記數下,對于Pop1和Pop2兩個群體MLP的遺傳參數估值與群體遺傳參數的差值較SSGBLUP、GBLUP和BLUP均更小,表明MLP在遺傳參數估計上優于其他3種方法,這可能是因為多層人工神經網絡采用非先驗假設很好的解釋了非加性效應,彌補了線性模型的局限性[34]。Gholizadeh等[35]在奶牛中應用神經網絡,得到了準確的產奶量遺傳力估計結果。Coelho de Sousa等[36]在咖啡中應用MLP進行遺傳力估計,得到的結果優于GBLUP,這是因為標記物之間的相互作用隱含在神經元多次加權的輸出中。

Pérez-Rodríguez 等[37]通過小麥數據集對神經網絡和幾種線性模型進行了評估,發現在植物中非線性神經網絡比線性模型有更好的預測精度。McDowell[38]使用擬南芥、玉米和小麥數據對常規基因組預測模型與MLP進行了比較,在3個性狀中,MLP的預測準確性均高于傳統的基因組預測模型。預測準確性結果顯示,相同h2、QTL和標記數下,相比于Pop1群體,具有更大參考群規模的Pop2群體中4種方法的GEBV準確性均有提高,由此推測低遺傳力性狀需要更多的表型數據來提高準確性。Luan等[39]使用挪威紅牛數據應用GBLUP及Bayes模型時也顯示出與本研究相似的結果,對于h2更高的性狀和更大參考群體,育種值估計更準確。本研究在不同QTL數下各模型獲得的準確性無明顯變化,可能是因為試驗中假設所有標記解釋相同的效應方差且無主效基因。在Daetwyler等[40]的模擬研究中,GBLUP在不同QTL數目時,準確性基本無變化,與本研究結果相符。隨著標記數從10K增加至50K,4種方法整體上基因組估計育種值準確性均呈上升趨勢,在張猛[41]的研究中也得到了相同的結論。在h2為0.05時,MLP在標記數為10K時的準確性略高于其他方法在50K時的準確性,這可能是由于標記數過少造成SSGBLUP法性能的下降。

相對于非限性性狀,在限性性狀中想要使遺傳評估和育種值估計結果達到更高水平難度更大,雖然使用GS利用雙性別的基因型數據彌補了一部分表型缺失的劣勢,但是育種值準確性的提高仍然有限[13]。在方法上,各模型對非限性性狀進行正常流程的遺傳評估和基因組選擇,而對于限性性狀,GBLUP必須使用EBV[42]或逆回歸育種值[43-44]作為偽表型值,通過偽表型值和基因組信息得到基因組估計育種值。SSGBLUP將系譜信息與基因組信息整合到模型[6]中,構建基因組關系矩陣(H陣),結合表型信息預測未進行基因分型個體的EBV及基因分型個體的GEBV。MLP是使用輸入數據為該網絡尋找最佳權重,它將輸入的多個數據集映射到單一的輸出數據集上,每個神經元的輸出表示為神經元全部輸入的加權和,通過線性和非線性激活函數[45-46]調節預測變量的特定權值從而學習基因型和表型之間的關系。在結果上,結合前人研究[47],非限性性狀的GEBV準確性較限性性狀更高,這是由于使用偽表型值或非線性學習模式均不能完全彌補表型的缺失。

在模擬研究中,非線性模型優于線性模型[48],而在真實數據中,非線性模型較線性模型只有輕微優勢或沒有優勢[49]。無論哪種模型都是簡化模型,與真實模型會有較大差異,可能更接近復雜生物現象的其他模型尚未找到,還需要繼續大量研究。本研究中MLP法的優勢和有效性仍然需要在真實綿羊群體中進行驗證。

4 結 論

對于綿羊限性性狀,MLP模型在遺傳參數估計方面較線性模型有更好的表現。MLP可以有效提高GEBV準確性,在綿羊限性性狀GS中有良好應用潛力。