中間偃麥草全基因組的選擇

郭 鵬， 曹 晟

(天津農學院計算機與信息工程學院，天津 300384)

中間偃麥草[(Host) Nevski]原產于東歐，是一種根系發(fā)達、抗寒性強、耐旱、耐鹽、適應性廣、再生性好的優(yōu)良牧草。作為多年生牧草，中間偃麥草產量高，已經在畜牧業(yè)發(fā)達的國家得到了廣泛應用。中間偃麥草的果實可食用，其谷物制成的食品在美國某些市場上已有銷售。此外，中間偃麥草對小麥銹病、白粉病等病害免疫，對黑穗病、葉枯病和根腐病等病害具有高抗性，能與小麥雜交，已經成為小麥遺傳改良中具有重要利用價值的野生親本之一。李振聲院士獲得了偃麥草與小麥的雜交種子，并將其廣泛用作與小麥遠緣雜交的野生親本材料，通過雜交將其優(yōu)良基因轉移進小麥，已經成為小麥遺傳改良的有效途徑。

黑龍江省分子細胞遺傳與遺傳育種重點實驗室研究團隊對中間偃麥草種質資源在直接馴化改良、遠緣雜交改良等方面進行了研究，并對利用基因組學育種技術進行中間偃麥草遺傳育種研究的重點和選育策略進行了展望。全基因組選擇使用單核苷酸多態(tài)性(single nucleotide polymorphism，SNP)標記信息和表型特征值估計候選個體育種值，根據育種值的大小對候選個體進行篩選。作為育種領域的新興技術，與傳統的僅使用個體表型的估計育種值相比，全基因組選擇具有準確度更高、世代間隔短、育種效率高等優(yōu)點，目前已經被越來越多地應用在農作物的育種領域。

在國外，Zhang等利用基因組選擇技術構建了基因標記數為3 883個、大小為1 126個的育種群體，并用不同的基因組選擇方法進行了研究，證明了使用基因組選擇技術進行中間偃麥草遺傳育種的可行性。Bajgain等使用全基因組關聯分析尋找中間偃麥草表型發(fā)育相關的基因座，利用顯著位點作為固定效應預測育種值，使用rrBLUP在SNP標記數為8 899個的數據集中進行育種值的估計，結果表明，該方法可使預測能力提高14%?？紤]到環(huán)境因素對表型發(fā)育的影響，Bajgain等此后使用GBLUP、貝葉斯方法研究加性、顯性效應以及環(huán)境效應對中間偃麥草性狀發(fā)育的影響，在標記數為 8 899 個、群體大小為451個的基因數據集中進行了穗質量、穗長、每個花序小穗數、產量等9種性狀的全基因組選擇研究。Crain等研究了中間偃麥草 3 658 個個體的46個性狀，對18 357個SNP標記的基因組數據使用GBLUP方法進行了研究，發(fā)現預測能力隨著訓練群體規(guī)模、標記數量的增加而提高。

本研究對標記數為23 495個、研究群體大小為5 521個的中間偃麥草數據采用隨機方式抽取個體構建訓練群體和校驗群體，使用交叉驗證的方式計算候選育種個體全基因組估計育種值和育種準確度，對比GBLUP、BayesA、BayesB、BayesCπ法的基因組育種值估計結果，以期找到適合中間偃麥草性狀育種分析的全基因組選擇方法來提高育種效率，為中間偃麥草全基因組育種提供參考。

1 材料與方法

本研究所用數據是從美國國立生物技術信息中心(NCBI)網站(https：//www.ncbi.nlm.nih.gov/bioproject/)上下載的美國堪薩斯州薩利納市土地研究所的中間偃麥草的公共數據。在試驗中選擇中間偃麥草的自由脫粒率、穗產量、株高、種子質量、每個花序小穗數、穗長、落粒性等7種性狀進行研究，研究群體包含5 521個個體，基因組數據包括23 495個SNP標記位點。

生物育種中的遺傳力也稱遺傳率，表示遺傳因素對表型發(fā)育所起作用的大小。由于表型數據中存在表型值缺失或無效的情況，因此需要對數據進行篩選，經過篩選后具有有效表型數據的7種性狀數據值的平均值、標準差、最大值和最小值的統計結果及對應的遺傳力如表1所示，其中自由脫粒、穗產量、株高、種子質量、每個花序小穗數、穗長、落粒性具有有效表型數據和SNP標記數據的群體中的個體數分別為5 019、5 230、5 488、4 703、5 013、4 512、4 694個。

表1 表型數據及遺傳力

1.1 統計模型

本研究選用的統計模型如下所示，

式中：為個體的性狀表型值；為標記的數量；為平均表型值；為第個分子標記的效應值；為個體第個標記的基因型編碼(編碼符號為0、1或2)；為個體的隨機殘差效應值。

1.2 BayesA

1.3 BayesB

BayesB假設只有部分位點效應能夠影響到遺傳效果，在所有位點的遺傳方差分布中，只有少數位點具有遺傳方差。因此，在BayesB中使用參數控制位點是否具有標記效應。

1.4 BayesCπ

BayesCπ假設所有位點具有1個公共方差，而不是每個位點都有自己的方差，并且假設控制標記效應為0的參數()未知。公共方差符合自由度=42、縮放因子為的尺度逆卡方分布，其中的推導過程與BayesB中縮放因子()的推導過程相同。概率()未知，并且假設符合(0，1)間的均勻分布，符合概率為(1-)的SNP效應由多元分布的混合模型計算得到，標記效應值也是使用正態(tài)分布取樣的方式進行估計的。

春秋末期，宗法制逐步被破壞，導致天子與各諸侯國的關系也發(fā)生了重大變化。天子作為天下共主的地位和權威進一步喪失?！按呵飼r猶嚴祭祀，重聘享，而七國則無其事矣；春秋時，猶論宗姓氏族，而七國則無一言及之矣；春秋時，猶宴會賦詩，而七國則不聞矣；春秋時猶赴告策書，而七國則無有矣。邦無定交，土無定主?！盵23]P715宗法、會盟、祭祀、等諸多方面都發(fā)生了重大變化。

1.5 GBLUP

在GBLUP中，基于系譜的親緣關系矩陣由遺傳關系矩陣(G陣)代替，然后使用最佳線性無偏估計計算育種值。G矩陣的公式如下：

式中：為標記的數量；為標記的等位基因頻率；為SNP標記效應的中心關聯矩陣。

1.6 GEBV計算

全基因組估計育種值(genomic estimated breeding value，Gebv)使用如下公式計算：

式中：為個體的基因組估計育種值。

1.7 交叉驗證

本研究選用隨機掩蔽交叉驗證方法將中間偃麥草的表型數據、基因型數據分成校驗集合和訓練集合，使用5倍交叉驗證方法對群體進行分析。每次隨機抽取約1/5個體構成校驗集合，剩余的個體構成訓練集合，具體分組見表2。對于每種性狀，重復基因組育種值估計10次，以10次估計結果的平均準確度作為全基因組估計的準確度。

表2 訓練群體和校驗群體分組情況

數據分組后，使用訓練群體中的表型數據、基因型數據估計位點效應值，以BayesA為例，交叉驗證全基因組估計實現流程如圖1所示。驗證階段使用估計的效應值，校驗群體中的基因型編碼計算校驗群體中個體的育種值。

1.8 估計準確性標準

將全基因組估計所得育種值與表型數據求相關系數，作為評判標準，計算公式：

2 結果與分析

2.1 試驗環(huán)境與試驗相關參數

本研究的環(huán)境：AMD RyZen 5 1600六核(12線程)中央處理器(CPU)，32 G DDR4 2 666 MHz內存，希捷2 TB(1萬轉)硬盤；中文版64位Windows 7操作系統。

在試驗中，BayesA、BayesB、BayesCπ中的burnin和總取樣運行次數分別設置成400 0、20 000。=4012、=0.002 0作為BayesA中的先驗分布參數；=4.2作為BayesCπ中的先驗分布參數，=4234、=0.042 9作為BayesB中的先驗分布參數，在BayesB、BayesCπ算法中，有效基因位占比()=0.01。

2.2 估計準確度

本試驗采用5倍交叉驗證，經過10次運行后，計算每次運行產生的估計準確度。由于訓練群體和校驗群體的構建采用隨機抽取的方式選擇個體，不同個體的SNP標記數據和表型數據存在差異，使得每次最終估計的育種值準確度存在差異，因此選用育種值估計準確度的“均值±標準偏差”作為最終估計結果。由7種中間偃麥草全基因組育種值估計的準確度(表3)可以看出，GBLUP預測的準確度最低，標準偏差也最低，說明GBLUP的性能最穩(wěn)定。在用BayesB方法得出的結果中，4個最優(yōu)，1個次優(yōu)；在用BayesCπ方法得出的結果中，2個最優(yōu)，2個次優(yōu)；在用BayesA方法得出的結果中，1個最優(yōu)，3個次優(yōu)；在用GBLUP方法得出的結果中，1個次優(yōu)?？傮w而言，BayesB方法的估計準確度表現最優(yōu)異。

表3 全基因組育種值估計的準確度

2.3 運行時間

在10次交叉驗證運行后，產生的10次運行時間如表4所示?？梢钥闯?，同一性狀的3種貝葉斯方法的運行時間均高于GBLUP法的運行時間，而且3種貝葉斯方法的運行時間也有差異，排序為BayesB>BayesA>BayesCπ。

表4 不同方法的運行時間

在分析單次貝葉斯全基因組估計的運行時間時，選擇7種性狀的第1次運行時間進行比較。如圖2所示，BayesA法的運行時間大約是BayesB法的75%，BayesCπ法的運行時間大約是BayesB法的60%。

2.4 中間偃麥草性狀全基因組估計育種值

由于不同全基因組選擇方法估計的中間偃麥草育種值不同，因此選擇10個不同個體，使用對各性狀估計精度最高的方法計算全基因組估計育種值，結果見表5。在育種實踐中，不同性狀的表型值間存在較大差異，使得估計所得育種值在數量級上的差異較大。因此在育種過程中可以針對某一種性狀進行個體的選擇，也可以針對所有性狀對候選個體進行取舍。在對所有研究性狀的估計育種值進行選擇時，考慮到不同育種值間數量級的差異，通常情況下對不同性狀賦予不同的權重，采用加權求和的方式或者根據綜合公式計算綜合育種值，再根據綜合育種值篩選候選個體。

表5 中間偃麥草不同性狀全基因組的估計育種值

3 討論

育種過程的關鍵環(huán)節(jié)是估計育種值，作為新興的育種技術，全基因組選擇在育種值估計準確度、縮短育種世代間隔、遺傳力低的性狀的選育中表現出了優(yōu)異性能。中間偃麥草作為一種優(yōu)良的牧草，對小麥病害具有高抗性，甚至具有免疫性，是雜交小麥重要的野生親本，進行中間偃麥草的全基因組育種研究既能促進優(yōu)良牧草的改良，又能為雜交小麥育種提供重要的遺傳親本，具有重要的實用價值。

本試驗對中間偃麥草自由脫粒率、穗產量、株高、種子質量、每個花序小穗數、穗長、落粒性7種性狀的全基因組進行了選擇性分析。結果顯示，從育種值估計準確度上看，BayesB法對自由脫粒率、穗產量、株高、落粒性4種性狀的估計準確度最高；BayesCπ法對種子質量、每個花序小穗數的估計準確度最高；BayesA法對穗長性狀的估計準確度最高。究其原因，全基因組選擇育種值估計的準確度受到群體大小、全基因組選擇方法、遺傳力、種群數據等因素的影響，很少有某種方法能夠保證在所有性狀育種值的估計中都能獲得最優(yōu)的準確度。

GBLUP使用遺傳關系矩陣代替親緣關系矩陣，利用最佳線性無偏估計方法估計育種值，相校于傳統的基于表型的BLUP方法，更能反映個體發(fā)育過程中性狀的遺傳性能，在保持BLUP較短運行時間的情況下，準確性更高，是當前動植物育種領域全基因組選擇的重要工具。GBLUP假設所有位點都具有遺傳效應，并未考慮單個SNP標記效應，在估計過程中不會出現由于隨機取樣產生的偶然性因素，因此性能穩(wěn)定，中間偃麥草7種性狀的全基因組估計結果的偏差最小。

貝葉斯方法在基因組選擇過程中求解每個SNP位點的效應值，更能反映基因位點對遺傳影響的真實情況。BayesA假設所有位點對性狀發(fā)育產生遺傳效應，BayesB、BayesCπ假設只有少量基因位點可以影響性狀的發(fā)育，其余的位點效應可以忽略不計，從目前全基因組關聯分析的結果可以看出，少量顯著位點效應影響性狀的發(fā)育更能反映遺傳效果。BayesB假設每個位點都具有自己的方差，BayesCπ假設所有位點都有1個公共的方差，BayesB更能反映SNP效應的真實情況，因此在中間偃麥草7種性狀全基因組育種值的估計精度中，BayesB法得出的最優(yōu)結果最多。

就運行時間而言，使用SNP標記數據計算遺傳關系矩陣占據了GBLUP絕大部分運行時間，在計算機運行環(huán)境相同的情況下，全基因組選擇群體大小、個體SNP標記數量決定了基因型數據，進而主導著GBLUP的運行時間。BayesA、BayesB、BayesCπ全基因組選擇方法假設SNP標記效應值符合正態(tài)分布，該正態(tài)分布的參數符合尺度逆卡方分布，通過隨機取樣的方式從尺度逆卡方分布中產生正態(tài)分布函數的參數，然后使用馬爾科夫鏈蒙特卡羅(Markov chain Monte Carlo，MCMC)取樣技術(BayesA使用Gibbs取樣，BayesB、BayesCπ使用Metropolis Hastings取樣)，利用正態(tài)分布取樣函數產生標記效應值。MCMC的運行過程包括burn-in訓練階段和sampling取樣階段，為了消除取樣偶然性的影響，burn-in階段、sampling階段的迭代次數通常設成較大的值，貝葉斯全基因組選擇方法屬于高度密集型計算任務，往往需要較長的運行時間。BayesCπ假設所有位點具有1個公共的方差，BayesB、BayesA假設每個位點具有自己的方差，需要為每個位點估計正態(tài)分布方差，而且解決相同的問題，Metropolis Hastings取樣運行時間比Gibbs 取樣運行時間更長，因此在中間偃麥草群基因組選擇中，GBLUP的運行時間最短。在3種貝葉斯方法中，BayesB的運行時間最長，BayesCπ的運行時間最短。貝葉斯方法的全基因組選擇準確度高，但漫長的運行時間是育種實踐應用過程中的不利因素，使用更高主頻的高性能處理器結合并行計算技術縮短運行時間，將有助于提高貝葉斯方法全基因組選擇技術的實用性。

4 結論

本研究使用GBLUP、BayesA、BayesB和BayesCπ進行中間燕麥草的自由脫粒率、穗產量、株高、種子質量、每個花序小穗數、穗長、落粒性7種性狀的全基因組選擇，經過5倍的交叉驗證發(fā)現，在4種方法育種值的估計結果中，BayesB的準確度最優(yōu)，運行時間最長。進行中間偃麥草全基因組的選擇需要針對不同性狀選擇最優(yōu)的全基因組選擇方法。在育種實踐中，綜合考慮不同方法育種性狀估計準確度和算法運行時間才能在中間偃麥草育種中高效地利用全基因組選擇技術，提高育種效率。