楊靜蕾 吳冰杰 王安洲 肖英杰,2,*

1 作物遺傳改良全國重點實驗室 / 華中農(nóng)業(yè)大學(xué), 湖北武漢 430070; 2 湖北洪山實驗室, 湖北武漢 430070

雜種優(yōu)勢利用和基于表型選擇的傳統(tǒng)育種為糧食產(chǎn)量的提升作出了巨大貢獻。近年來, 氣候變化和人口增加給全球糧食安全帶來了嚴(yán)峻挑戰(zhàn), 預(yù)計到2050 年, 世界人口將增至95 億, 屆時需要多生產(chǎn)70%的糧食才能滿足人口需求[1-2], 因此亟需新的育種技術(shù)和方法來提高作物產(chǎn)量[3-4]。高通量基因分型平臺和表型平臺、基因組選擇、機器學(xué)習(xí)等新技術(shù),為育種效率的提升和作物新品種的培育提供了新的契機[5]。

Meuwissen 在 2001 年首次提出基因組選擇(genomic selection, GS)或基因組預(yù)測的概念[6], 其核心是利用參考群體的全基因組分子標(biāo)記基因型數(shù)據(jù)與表型數(shù)據(jù)建立統(tǒng)計模型, 在表型未知而基因型已知的候選群體中, 利用該模型估計每個材料的基因組育種值, 進而在早期實施選擇, 大幅度提高了遺傳增益。GS 在理論上類似于MAS (marker assisted selection), 但MAS 僅利用少量顯著性標(biāo)記對主效基因進行選擇。動植物重要性狀一般為數(shù)量性狀, 由大量微效基因控制, MAS 無法捕捉到微效基因的貢獻, 因此對數(shù)量性狀改良的作用有限。而基于覆蓋全基因組的高密度標(biāo)記, 直接估計全基因組中的所有標(biāo)記效應(yīng)的基因組選擇, 可以捕獲具有微小影響的遺傳位點[7], 能夠更好地解釋表型變異, 提高復(fù)雜農(nóng)藝性狀地預(yù)測精度和選擇效率。

隨著下一代測序技術(shù)和高密度單核苷酸多態(tài)性(single nucleotide polymorphism, SNP)基因分型技術(shù)的快速發(fā)展, 目前, 基因組預(yù)測已經(jīng)成為革命性的育種手段?；蚪M選擇的首次利用是在奶牛中, 極大縮短了奶牛選育的世代間隔, 后來又成功用于豬、羊等主要動物的育種[8-9]。據(jù)報道[10], 動物中基因組預(yù)測的精度可達0.8, 這為基于預(yù)測結(jié)果對動物幼崽進行早期篩選, 加快種畜選育提供了可能。在玉米、水稻、小麥等作物中[11-14], 也陸續(xù)開展了該類研究, 比如, Wang 等[13]在水稻雜交育種中開發(fā)的MV-ADV 模型有助于對低遺傳力性狀(如產(chǎn)量)進行基因組預(yù)測; Ma 等[14]曾使用不同遺傳結(jié)構(gòu)的群體構(gòu)建訓(xùn)練集對大豆百粒重進行全基因組預(yù)測, 預(yù)測精度最高可達0.75; Charmet 等[15]評估了基因組選擇在小麥3 個育種群體的3 個性狀的表現(xiàn), 具有較高遺傳力的抽穗期, 其預(yù)測精度可高達0.7。

在玉米中, 有研究發(fā)現(xiàn), 基因組選擇優(yōu)于分子標(biāo)記輔助選擇和基于系譜的傳統(tǒng)表型選擇, 可顯著提高低遺傳力性狀的預(yù)測精度和選擇效率[16-17]。Cao等[18]利用玉米的1 個關(guān)聯(lián)群體和3 個雙單倍體群體,對其焦油斑點病抗性性狀進行GS 研究, 發(fā)現(xiàn)預(yù)測精度受訓(xùn)練群體大小和標(biāo)記密度影響, 但總體而言預(yù)測精度較高, 說明GS 對玉米抗病性篩選具有較強的應(yīng)用潛力?；诰哂袩釒а壍?2 個玉米雙親分離群體的基因型和干旱和正常情況的表型數(shù)據(jù),Zhang 等[19]對產(chǎn)量、株高和開花期3 個性狀進行基因組預(yù)測, 發(fā)現(xiàn)使用200 個SNPs 標(biāo)記及50%群體作訓(xùn)練集時, 3 個性狀的預(yù)測精度分別是0.28、0.32 和0.29, 且隨著訓(xùn)練群體大小和標(biāo)記密度的增加, 預(yù)測精度也會隨之提升。目前, GS 已在玉米的多種遺傳和育種群體中得到廣泛應(yīng)用[20-22]。

統(tǒng)計模型是基因組預(yù)測的核心, 影響基因組預(yù)測精度。常用的 GS 模型有 LSE (least-squares estimation)、BLUP (best linear unbiased prediction)、Bayes (bayesian analysis)、LASSO (least absolute shrinkage prediction operator)等[23]。其中 rrBLUP(ridge regression best linear unbiased prediction)的利用最為廣泛, 其是利用訓(xùn)練群體估計標(biāo)記效應(yīng), 在預(yù)測群體中將標(biāo)記效應(yīng)累加, 進而預(yù)測未知表型個體的基因組育種值。據(jù)報道, 基于該方法對玉米的開花時間進行GS, 預(yù)測精度可達0.64[24]。Bhering等運用模擬數(shù)據(jù)集評估rrBLUP、GBLUP (genomic best linear unbiased prediction)、Bayesian LASSO 三種統(tǒng)計方法的基因組預(yù)測性能, 發(fā)現(xiàn)該方法在分析效率和預(yù)測精度上表現(xiàn)均最好[25]。Yan 等[26]基于rrBLUP、Bayes 和RF (random forest) 3 種統(tǒng)計方法,對玉米籽粒中鎘元素含量進行預(yù)測, 預(yù)測精度分別達到0.89、0.83 和0.75, 其中rrBLUP 方法的預(yù)測精度最高。LASSO 由Tibshirani 于1996 年首次提出[27],該方法通過最小化殘差平方和的約束, 實現(xiàn)對高維數(shù)據(jù)的變量選擇, 從而在基因表達分析中被廣泛應(yīng)用, 具有較高的預(yù)測準(zhǔn)確性、良好的解釋能力和穩(wěn)健性。在利用衛(wèi)星遙感數(shù)據(jù)預(yù)測玉米年產(chǎn)量變化的研究中, LASSO 的預(yù)測精度達0.78[28]。在陸地棉纖維品質(zhì)的預(yù)測研究中, Islam 等[29]發(fā)現(xiàn) Bayesian LASSO 比 GBLUP、rrBLUP、Bayes、PKHS(reproducing kernel space) 4 種統(tǒng)計模型的預(yù)測精度要更高。Tsai 等[30]使用春大麥和冬小麥2 套數(shù)據(jù)驗證rrBLUP 和Bayesian LASSO 兩種模型的預(yù)測性能,發(fā)現(xiàn)在春大麥中, Bayesian LASSO 對白粉病和產(chǎn)量的預(yù)測精度高于rrBLUP, 而在冬小麥中, 2 個模型對產(chǎn)量的預(yù)測精度相似。以上研究表明, 對于不同性狀、不同物種和群體, 基因組預(yù)測方法表現(xiàn)可能存在差異, 但是rrBLUP 與LASSO 模型在眾多統(tǒng)計模型中表現(xiàn)出更好的性能, 于是本研究主要基于以上2 種統(tǒng)計模型進行GS 研究。但前人研究往往基于少量性狀和同種類型數(shù)據(jù), 對農(nóng)藝和品質(zhì)性狀的預(yù)測能力系統(tǒng)評估還較缺乏。

近年來, 隨著多維組學(xué)技術(shù)的迅速發(fā)展, 植物在轉(zhuǎn)錄、翻譯和代謝水平上的變化都可以進行定量檢測。Azodi 等[24]基于轉(zhuǎn)錄組數(shù)據(jù)和基因組數(shù)據(jù)來構(gòu)建開花時間的預(yù)測模型, 通過評估模型特征的重要性, 發(fā)現(xiàn)模型中最重要的2個特征都為轉(zhuǎn)錄組特征, 驗證了轉(zhuǎn)錄組數(shù)據(jù)相比基因組數(shù)據(jù)對預(yù)測模型的重要性, 這也表明轉(zhuǎn)錄組提升基因組預(yù)測方面的巨大潛力。此外, Zhang 等[31]對由385 個自交系組成的玉米關(guān)聯(lián)群體進行全基因組代謝物分析, 共檢測到1035 種顯著變化的代謝物, 發(fā)現(xiàn)其中15 種代謝物性狀可解釋超過60%的玉米苗期干旱處理后存活率表型變異。Qin 等[32]曾通過代謝物分析發(fā)現(xiàn)硼元素通過增加抗氧化酶的活性和改變代謝產(chǎn)物, 減輕了鎘毒害對小麥的影響。Hu 等[33]發(fā)現(xiàn)使用轉(zhuǎn)錄組、代謝組數(shù)據(jù)對水稻產(chǎn)量的預(yù)測精度可達 0.4869 和0.4593, 遠高于基因組預(yù)測的精度, 這說明對于受微效多基因控制的產(chǎn)量性狀, 轉(zhuǎn)錄組和代謝組數(shù)據(jù)可能會提供比基因組更豐富的信息進行預(yù)測。

目前, 在玉米中, 利用多組學(xué)數(shù)據(jù)對玉米重要性狀的預(yù)測研究報道較少, 在不同模型和不同性狀上的系統(tǒng)評估也較為缺乏。為系統(tǒng)研究組學(xué)數(shù)據(jù)對玉米性狀預(yù)測的效果, 本研究以具有廣泛多樣性的368 個玉米自交系的基因組、轉(zhuǎn)錄組和代謝組數(shù)據(jù),基于rrBLUP 與LASSO 模型, 對55 個農(nóng)藝和品質(zhì)性狀進行預(yù)測分析, 系統(tǒng)評估了各組學(xué)數(shù)據(jù)和統(tǒng)計模型, 對農(nóng)藝性狀和品質(zhì)性狀預(yù)測能力的差異, 為后續(xù)玉米重要性狀的基因組育種提供了理論依據(jù)。

1 材料與方法

1.1 試驗材料

1.1.1 群體及表型數(shù)據(jù) 本研究使用來自368 個不同玉米自交系組成關(guān)聯(lián)群體作為試驗材料, 分別在海南(三亞, 18°25′N, 109°51′E, 2010 年)、云南(昆明, 24°25′N, 102°30′E, 2011 年)和重慶(29°25′N,106°50′E, 2011 年) 3 個地點進行田間試驗。根據(jù)系譜信息, 將所有自交系分為2 組, 即溫帶(temperate,TEM)和熱帶/亞熱帶(tropic/sub-tropic, TST)。所有自交系根據(jù)不完全隨機區(qū)組試驗設(shè)計, 單行種植在試驗地中。對收獲的玉米植株進行表型鑒定, 獲取包括株高、穗行數(shù)在內(nèi)的20 個農(nóng)藝性狀和包括籽粒生育酚、脂肪酸含量在內(nèi)的35 個品質(zhì)性狀。本研究使用的表型數(shù)據(jù)均來自于已發(fā)表文章[32-34]。對表型數(shù)據(jù)進行預(yù)處理, 刪除各表型的異常值并計算55 個性狀的描述性統(tǒng)計。

1.1.2 基因組和轉(zhuǎn)錄組數(shù)據(jù) 本研究所用群體已使用Illumina Maize SNP50 芯片和轉(zhuǎn)錄組進行測序,并鑒定到103 萬個覆蓋全基因組的高質(zhì)量SNPs 和28,769 個基因的表達數(shù)據(jù)[34-36]。本研究從中隨機挑選5 萬個SNPs 作為基因組數(shù)據(jù); 將基因表達數(shù)據(jù)進行標(biāo)準(zhǔn)化處理后作為轉(zhuǎn)錄組數(shù)據(jù)?；蚪M和轉(zhuǎn)錄組數(shù)據(jù)將用于后續(xù)的多組學(xué)預(yù)測研究。

1.1.3 代謝組數(shù)據(jù) 本研究使用群體的成熟籽粒已經(jīng)進行了靶向和非靶向代謝組分析。其中靶向代謝組是驗證測試樣品中是否存在目標(biāo)代謝物的檢測方式, 本研究使用的靶向代謝組是指以關(guān)聯(lián)群體在云南和重慶2 個環(huán)境的成熟玉米籽粒為測試樣品,檢測其中部分氨基酸含量所形成的數(shù)據(jù)集: 包括17種重要的氨基酸、47 個氨基酸的衍生性狀以及它們各種氨基酸的總和, 共計130 種靶向代謝物[37]。非靶向代謝組是通過對生物體內(nèi)代謝物進行全面分析,找出差異代謝物的研究方式。本研究使用的非靶向代謝組是關(guān)聯(lián)群體在海南、云南、重慶3 個環(huán)境下種植, 對成熟籽粒進行基于氣相色譜或液相色譜結(jié)合質(zhì)譜的代謝物分析獲取的數(shù)據(jù), 在3 個環(huán)境下共鑒定到2031 種代謝物[38]。以上代謝物數(shù)據(jù)將用于本研究的多組學(xué)預(yù)測研究。研究中進行了靶向代謝物與非靶向代謝物的整合, 通過逐步整合得到5 個數(shù)據(jù)集, 包括海南非靶向代謝物的數(shù)據(jù)集E1; 海南和云南非靶向代謝物的數(shù)據(jù)集E1+E2; 海南、云南和重慶非靶向代謝物的數(shù)據(jù)集E1+E2+E3; 全部非靶向代謝物和重慶靶向代謝物的數(shù)據(jù)集 E1+E2+E3+CQ; 所有非靶向代謝物和靶向代謝物的數(shù)據(jù)集E1+E2+E3+CQ+YN。

1.2 基因組預(yù)測的統(tǒng)計模型

1.2.1 嶺回歸最佳線性無偏估計(rrBLUP) 該方法是基因組預(yù)測最常用的模型之一, 它通過在訓(xùn)練群體中估計標(biāo)記效應(yīng), 結(jié)合預(yù)測群體的基因型信息將標(biāo)記效應(yīng)累計, 最終獲得預(yù)測群體的個體估計育種值。該方法在構(gòu)建標(biāo)記效應(yīng)矩陣時, 賦予部分標(biāo)記較大權(quán)重, 從而放大有效標(biāo)記的貢獻, 弱化無效標(biāo)記的作用, 以使預(yù)測結(jié)果更符合實際[39]。本研究使用基于R 語言的“rrBLUP”軟件包實現(xiàn)模型構(gòu)建,此軟件包含3 個函數(shù): A.mat 函數(shù)主要用于構(gòu)建加性效應(yīng)矩陣; mixed.solve 函數(shù)主要用于混合模型的參數(shù)求解, 可以用于預(yù)測標(biāo)記效應(yīng)或育種值; kin.blup函數(shù)使用加性效應(yīng)關(guān)系矩陣預(yù)測基因型值, 進而預(yù)測表型。本研究首先使用A.mat 函數(shù)計算加性效應(yīng)矩陣, 再通過kin.blup 函數(shù)進行表型預(yù)測, 計算預(yù)測結(jié)果與真實值的皮爾遜相關(guān)系數(shù)(Pearson correlation efficiency,r), 即預(yù)測精度。

1.2.2 最小絕對收縮選擇算子(LASSO) 該模型通過構(gòu)造一個懲罰函數(shù)得到一個較為精煉的模型,使它收縮一些系數(shù), 同時設(shè)定一些系數(shù)為零, 以實現(xiàn)對變量的選擇和對模型復(fù)雜程度的降低。因其對數(shù)據(jù)的要求極低, 無論變量是連續(xù)還是離散的, 都能用LASSO 進行處理。本研究利用R語言中“l(fā)attic”軟件包進行LASSO 預(yù)測, 首先使用createFolds 函數(shù)將數(shù)據(jù)隨機分為5 份, 再基于訓(xùn)練集使用train 函數(shù)中的LASSO 參數(shù)進行模型構(gòu)建, 然后基于測試集的組學(xué)數(shù)據(jù)使用predict 函數(shù)進行表型預(yù)測。模型的預(yù)測精度, 為預(yù)測結(jié)果與真實值的皮爾遜相關(guān)系數(shù)。

1.2.3 基因組預(yù)測的交叉驗證 本研究使用5 折交叉驗證計算模型的預(yù)測值與真實值的相關(guān)性, 并將5 次重復(fù)交叉驗證結(jié)果的均值來評估模型的性能。該方法將數(shù)據(jù)集隨機切分為5 個互不相交且大小相同的子集, 使用其中4 個子集訓(xùn)練模型, 剩下的1 個子集當(dāng)作測試集測試模型。將上述步驟重復(fù)5 次, 每次挑選不同的子集作為測試集, 訓(xùn)練得到5個模型, 共得到5 個預(yù)測值與真實值的皮爾遜相關(guān)系數(shù), 對這5 次的皮爾遜相關(guān)系數(shù)取平均值得到一個交叉驗證的結(jié)果。該方法會使所有數(shù)據(jù)參與到訓(xùn)練和預(yù)測中, 在最大化利用數(shù)據(jù)訓(xùn)練模型的同時避免過擬合的現(xiàn)象。

2 結(jié)果與分析

2.1 基于基因組數(shù)據(jù)分析對農(nóng)藝性狀和品質(zhì)性狀的預(yù)測差異

基于基因組數(shù)據(jù)和rrBLUP 模型, 研究發(fā)現(xiàn), 品質(zhì)性狀的平均預(yù)測精度(r=0.628)明顯優(yōu)于預(yù)測農(nóng)藝性狀(r=0.504) (圖1-A)。在農(nóng)藝性狀中, 僅有5 個性狀的預(yù)測精度高于0.6, 分別是粒長(kernel length)、散粉期(pollen shed)、雄穗分枝數(shù)(tassel branch number)、穗位葉長(ear leaf length)、抽雄期(heading date), 其中粒長的預(yù)測效果最好, 達到0.689 (圖1-B)。在品質(zhì)性狀中, 60%以上的性狀預(yù)測精度大于0.6,其中籽粒含油量(OIL)的預(yù)測效果最好, 達到0.887,預(yù)測效果最佳的前5 個性狀分別是含油量(OIL)、油酸(C18_1)、亞油酸(C18_2)、花生烯酸(C20_1)、花生酸(C20_0), 預(yù)測精度均高于0.7 (圖1-C)。

圖1 基于基因組數(shù)據(jù)分析對農(nóng)藝性狀和品質(zhì)性狀的預(yù)測差異Fig.1 Prediction differences in agronomic traits and quality traits based on genomic data analysis

2.2 不同組學(xué)數(shù)據(jù)對性狀預(yù)測的差異

對比分析基于基因組、轉(zhuǎn)錄組、代謝組預(yù)測不同類型性狀的結(jié)果, 研究發(fā)現(xiàn), 對于農(nóng)藝性狀,基因組數(shù)據(jù)的預(yù)測精度都顯著高于轉(zhuǎn)錄組的預(yù)測精度(r=0.504 vsr=0.459,P=4E-3)。此外, 轉(zhuǎn)錄組數(shù)據(jù)的預(yù)測精度顯著高于代謝組的預(yù)測精度(r=0.459 vsr=0.333,P=5.39E-6)。品質(zhì)性狀也遵循該趨勢(圖2-A)。

圖2 不同組學(xué)數(shù)據(jù)對性狀預(yù)測的差異Fig.2 Differences in trait prediction between omics data

此外, 相比于其他組學(xué)數(shù)據(jù), 基因組數(shù)據(jù)對27個性狀的預(yù)測精度最高。在這27 個性狀中, 花生烯酸(C20_1)利用基因組數(shù)據(jù)的預(yù)測精度高達 0.832,明顯高于該性狀基于轉(zhuǎn)錄組數(shù)據(jù)和代謝組數(shù)據(jù)的預(yù)測結(jié)果(0.776 和0.638)。其中, 軟脂酸與棕櫚油酸的比例(C160.C161)呈現(xiàn)出基因組數(shù)據(jù)和轉(zhuǎn)錄組數(shù)據(jù)的預(yù)測差值最大(rdiff.=0.251), 雄穗主軸長(Tassel main axis length)呈現(xiàn)出基因組數(shù)據(jù)和代謝組數(shù)據(jù)的預(yù)測差值最大(rdiff.=0.419), 同時該性狀也表現(xiàn)出轉(zhuǎn)錄組數(shù)據(jù)和代謝組數(shù)據(jù)的最大預(yù)測差值(rdiff.=0.275)。我們發(fā)現(xiàn), 有21 個性狀使用基因組數(shù)據(jù)和轉(zhuǎn)錄組數(shù)據(jù)的預(yù)測精度沒有明顯差異(rdiff.<0.05),但均優(yōu)于代謝組數(shù)據(jù)的預(yù)測精度。其中, 籽粒含油量(OIL)表現(xiàn)出基因組數(shù)據(jù)和轉(zhuǎn)錄組數(shù)據(jù)的最高的預(yù)測精度(r=0.887 和r=0.862), 明顯高于代謝組對該性狀的預(yù)測精度(r=0.734)。油酸含量(C18_1P)、亞油酸含量(C18_2P)、油酸與亞油酸的比例(C181.C182)的基因組和代謝組數(shù)據(jù)預(yù)測精度沒有明顯差異, 但均顯著高于轉(zhuǎn)錄組數(shù)據(jù)的預(yù)測精度, 同時,我們發(fā)現(xiàn)代謝組預(yù)測精度在各組學(xué)中達到最高。吐絲期(silking stage)使用轉(zhuǎn)錄組數(shù)據(jù)的預(yù)測精度達到0.606, 明顯優(yōu)于基因組和代謝組數(shù)據(jù)的預(yù)測精度(0.555 和0.526)。此外, 粒長(kernel length)被基因組、轉(zhuǎn)錄組和代謝組數(shù)據(jù)預(yù)測精度都非常高(r>0.68)(圖2-B)。

通過分析不同染色體標(biāo)記的預(yù)測差異, 我們發(fā)現(xiàn), 農(nóng)藝性狀、品質(zhì)性狀都呈現(xiàn)出了相同的預(yù)測趨勢。1 號染色體(n=7976)、5 號染色體(n=6093)和6號染色體(n=4098)的標(biāo)記預(yù)測精度高, 9 號染色體(n=3703)和10 號染色體(n=3480)的標(biāo)記預(yù)測精度低,這可能說明9 號、10 號染色體預(yù)測精度低是因為染色體上標(biāo)記數(shù)目少(圖2-C)。為了進一步驗證這個猜想, 本研究以10 條染色體中最少的標(biāo)記數(shù)目為參考,對染色體的標(biāo)記進行隨機挑選, 并分別進行表型預(yù)測, 將結(jié)果進行比較, 發(fā)現(xiàn)仍為1 號染色體、5 號染色體、6 號染色體上的標(biāo)記預(yù)測效果最佳, 9 號染色體和10 號染色體上的標(biāo)記預(yù)測效果較差。將隨機預(yù)測結(jié)果與染色體上所有標(biāo)記預(yù)測結(jié)果相比, 發(fā)現(xiàn)所有標(biāo)記預(yù)測結(jié)果略微優(yōu)于挑選標(biāo)記預(yù)測結(jié)果, 同時仍呈現(xiàn)出1 號染色體、5 號染色體、6 號染色體上的挑選標(biāo)記預(yù)測效果最佳, 9 號染色體和10 號染色體上的挑選標(biāo)記預(yù)測效果最差, 這說明1 號染色體、5號染色體和6 號染色體上可能存在對表型變異貢獻最大的標(biāo)記。因此, 研究初步推測在訓(xùn)練基于基因組數(shù)據(jù)的基因組預(yù)測模型時, 加大1 號染色體、5 號染色體、6 號染色體上的標(biāo)記數(shù)目, 可以提升全基因組數(shù)據(jù)的預(yù)測效果。

基于基因高表達的玉米組織, 本研究將轉(zhuǎn)錄組數(shù)據(jù)分為根高表達基因集Root (n=6322)、葉片高表達基因集 Leaf (n=3101)和穗高表達基因集 Ear(n=2577)三大類, 以及所有基因集(n=18,548)。結(jié)果發(fā)現(xiàn), 使用3 種組織數(shù)據(jù)預(yù)測2 種性狀的結(jié)果表現(xiàn)趨勢略微不同, 對于農(nóng)藝性狀而言, 穗高表達基因集的預(yù)測結(jié)果最好(r=0.435), 葉片高表達基因集的預(yù)測結(jié)果其次(r=0.434), 根高表達基因集的預(yù)測效果最差(r=0.428); 對于品質(zhì)性狀而言, 葉片高表達基因集的預(yù)測結(jié)果最好(r=0.534), 穗高表達基因集的預(yù)測結(jié)果其次(r=0.523), 根高表達基因集的預(yù)測效果最差(r=0.518) (圖2-D)。研究進一步對比基于不同組織預(yù)測2 種類型性狀的結(jié)果, 發(fā)現(xiàn)3 種組織數(shù)據(jù)集中均呈現(xiàn)出品質(zhì)性狀的預(yù)測結(jié)果優(yōu)于農(nóng)藝性狀。

通過逐步整合不同環(huán)境的代謝物數(shù)據(jù), 得到5個數(shù)據(jù)集, 分別是 E1、E1+E2、E1+E2+E3、E1+E2+E3+CQ、E1+E2+E3+CQ+YN, 分別用以上數(shù)據(jù)集分別預(yù)測農(nóng)藝性狀與品質(zhì)性狀, 發(fā)現(xiàn)呈現(xiàn)不同的預(yù)測趨勢。對于農(nóng)藝性狀, 數(shù)據(jù)集增大伴隨著預(yù)測精度的下降, 對于品質(zhì)性狀, 數(shù)據(jù)集增大伴隨著預(yù)測精度的上升(圖2-E)。為了解釋更多環(huán)境的代謝物會降低農(nóng)藝性狀預(yù)測精度的現(xiàn)象, 我們分析了5 個環(huán)境的代謝物與農(nóng)藝和品質(zhì)性狀的相關(guān)性。統(tǒng)計結(jié)果發(fā)現(xiàn), 與農(nóng)藝性狀正相關(guān)和負(fù)相關(guān)的代謝物數(shù)目比較接近, 因此整合不同代謝數(shù)據(jù)可能會存在貢獻抵消的情況, 而導(dǎo)致預(yù)測精度下降; 對品質(zhì)性狀而言, 不同環(huán)境代謝物負(fù)相關(guān)數(shù)目遠小于正相關(guān), 因此整合數(shù)據(jù)預(yù)測有更大潛力提升預(yù)測精度。

2.3 不同模型對性狀預(yù)測的影響

總體來說, rrBLUP 相比于LASSO 模型對不同性狀的預(yù)測精度更高, 但對于某些品質(zhì)性狀, LASSO模型的預(yù)測效果更佳。具體來說, 基于基因組數(shù)據(jù),發(fā)現(xiàn)對于55 種性狀rrBLUP 均為最佳模型。基于轉(zhuǎn)錄組數(shù)據(jù), rrBLUP 對于20 個農(nóng)藝性狀是最優(yōu)預(yù)測模型, 對品質(zhì)性狀中的大多數(shù)性狀為最優(yōu)模型, 而軟脂酸與棕櫚油酸的比例(C160.C161)和生育酚比值(ratio)則是LASSO 為最佳預(yù)測模型。基于代謝組數(shù)據(jù), LASSO 對于4 個農(nóng)藝性狀和8 個品質(zhì)性狀是最優(yōu)模型。在農(nóng)藝性狀中LASSO 預(yù)測效果最佳的性狀是株高(Plant height) (r=0.331), 品質(zhì)性狀中預(yù)測效果最佳的性狀是亞麻酸含量(C18_3P) (r=0.622)(圖3-A)?；谏鲜鼋Y(jié)果, 我們初步推測rrBLUP 適合用于預(yù)測以農(nóng)藝性狀為代表的數(shù)量性狀, LASSO更適合用于預(yù)測以品質(zhì)性狀為代表的質(zhì)量性狀。這可能是由于LASSO 模型本質(zhì)是進行變量選擇, 該方法很難捕獲到微效基因的效應(yīng), 而對于由主效基因控制的質(zhì)量性狀, 該方法可以放大主效基因的作用,以實現(xiàn)更好的預(yù)測效果。

圖3 不同模型和數(shù)據(jù)組合對性狀預(yù)測的整合評估Fig.3 Integrated evaluation of models and omic-data combinations on trait prediction

為探究不同模型之間精度變化規(guī)律, 本研究設(shè)定預(yù)測精度為0～0.3 為低精度、0.3～0.6 為中精度、0.6～1.0 為高精度(圖3-B)。利用基因組數(shù)據(jù), 發(fā)現(xiàn)使用rrBLUP 模型時所有性狀的預(yù)測精度均高于0.3,其中27 個性狀預(yù)測精度較高。使用LASSO 模型時,41 個性狀預(yù)測精度中等或較高。其中, 4 個性狀被rrBLUP 和LASSO 預(yù)測均較高, 均為品質(zhì)性狀, 分別是油酸(C18_1)、亞油酸(C18_2)、油分濃度(OIL)、生育酚總量(total)。利用轉(zhuǎn)錄組數(shù)據(jù), 發(fā)現(xiàn)使用rrBLUP 模型時 18 個性狀預(yù)測精度較高。使用LASSO 模型, 有39 個性狀預(yù)測精度中等或較高。其中1 個性狀-gamma 生育酚(gamma)被rrBLUP 和LASSO 預(yù)測精度均較高。利用代謝組數(shù)據(jù), 發(fā)現(xiàn)使用rrBLUP 模型時40 個性狀預(yù)測精度中等或較高,其中, 10 個性狀預(yù)測精度較高。使用LASSO 模型時,37 個性狀預(yù)測精度中等或較高, 其中, 6 個性狀預(yù)測精度較高(圖3-B)。這些結(jié)果表明, 針對不同類型性狀, 不同組學(xué)數(shù)據(jù)在rrBLUP 和LASSO 模型下表現(xiàn)出不同的精度變化模式。

通過整合分析, 因此我們鑒定出不同性狀的最佳預(yù)測的數(shù)據(jù)和模型組合(圖3-C)。結(jié)果發(fā)現(xiàn), 45 個性狀的最優(yōu)預(yù)測組合為基因組數(shù)據(jù)和rrBLUP 模型(Geno_rrBLUP), 其中22 個性狀預(yù)測精度大于0.6,其中預(yù)測精度最高的性狀為籽粒含油量(OIL) (r=0.887)。3 個性狀的最優(yōu)預(yù)測組合為轉(zhuǎn)錄組數(shù)據(jù)和rrBLUP 模型(Tran_rrBLUP), 預(yù)測精度均大于0.6,其中預(yù)測精度最高的性狀為散粉期(Pollenshed)(r=0.651)。生育酚比值(ratio)最優(yōu)預(yù)測組合為轉(zhuǎn)錄組數(shù)據(jù)和 LASSO 模型(Tran_LASSO), 預(yù)測精度為0.459; 3 個性狀的最優(yōu)預(yù)測組合為代謝組數(shù)據(jù)和rrBLUP 模型(Meta_rrBLUP), 預(yù)測精度均大于0.6,精度最高的性狀為油酸含量(C18_1P) (r=0.667); 3 個性狀的最優(yōu)預(yù)測組合為代謝組數(shù)據(jù)和LASSO 模型(Meta_LASSO), 其中預(yù)測精度最高的性狀是軟脂酸與硬脂酸的比例(C160.C180) (r=0.522)。此外, 結(jié)果顯示, 僅在農(nóng)藝性狀中出現(xiàn)轉(zhuǎn)錄組數(shù)據(jù)為最優(yōu)預(yù)測組合的數(shù)據(jù)來源, 僅在品質(zhì)性狀中出現(xiàn)代謝組數(shù)據(jù)為最優(yōu)預(yù)測組合的數(shù)據(jù)來源。這說明, 轉(zhuǎn)錄組數(shù)據(jù)對農(nóng)藝性狀預(yù)測更有價值, 代謝組數(shù)據(jù)對品質(zhì)性狀預(yù)測作用更大。

2.4 材料異質(zhì)性對不同組學(xué)預(yù)測的影響

本研究的368 個自交系根據(jù)系譜來源可分為溫帶(temperate, TEM)材料和熱帶(tropic/sub-tropic,TST)材料?；诨蚪M、轉(zhuǎn)錄組和代謝組數(shù)據(jù), 利用rrBLUP 模型對溫?zé)釒Р牧线M行預(yù)測分析。研究發(fā)現(xiàn),訓(xùn)練集和測試集的材料類型相同時(即溫帶預(yù)測溫帶, TEM-TEM; 或熱帶預(yù)測熱帶, TST-TST; 精度0.211～0.582), 3 種組學(xué)數(shù)據(jù)預(yù)測效果均優(yōu)于訓(xùn)練集和測試集類型不同的情況(即溫帶預(yù)測熱點, TEMTST; 或熱帶預(yù)測溫帶, TST-TEM; 精度0.120～0.197)(圖4-A)。同時, 研究發(fā)現(xiàn), 使用TEM 預(yù)測TEM 的效果明顯優(yōu)于使用TST 預(yù)測TST 的效果, 而使用TST 預(yù)測TEM 的效果略優(yōu)于使用TEM 預(yù)測TST 的效果(圖4-B), 這可能是由于熱帶材料的變異程度比溫帶材料的變異程度豐富, 從而使得用TST 訓(xùn)練集很難覆蓋TST 測試集的遺傳變異, 但是包含了TEM中的大部分變異。除此之外, 研究發(fā)現(xiàn)使用歸屬于TST 的代謝組數(shù)據(jù)訓(xùn)練模型時, 預(yù)測TEM 與預(yù)測TST 的精度相差無幾, 僅有0.015, 這個結(jié)果說明在使用親緣關(guān)系較遠的材料進行預(yù)測時, 可以加入代謝組學(xué)數(shù)據(jù)進行模型的構(gòu)建, 可能會提升預(yù)測精度。

圖4 材料系譜差異對不同組學(xué)預(yù)測的影響Fig.4 Influence of material genealogy on prediction by omics data

對農(nóng)藝和品質(zhì)性狀, 材料異質(zhì)性對不同組學(xué)預(yù)測進行分析。對比發(fā)現(xiàn), 基于基因組數(shù)據(jù)(圖4-C)和代謝組數(shù)據(jù)(圖4-E)的預(yù)測結(jié)果, 在4 種不同訓(xùn)練集和測試集的組合中, 均表現(xiàn)出品質(zhì)性狀的預(yù)測效果優(yōu)于農(nóng)藝性狀的預(yù)測結(jié)果; 基于轉(zhuǎn)錄組數(shù)據(jù)的預(yù)測結(jié)果(圖4-D), 發(fā)現(xiàn)除了TEM-TEM 的組合之外, 其余3 種組合均表現(xiàn)出農(nóng)藝性狀的預(yù)測效果優(yōu)于品質(zhì)性狀的預(yù)測結(jié)果。此外, 研究發(fā)現(xiàn)基于代謝組數(shù)據(jù)的預(yù)測結(jié)果也表現(xiàn)出 TST-TST 的預(yù)測結(jié)果與TST-TEM 的結(jié)果相當(dāng), 而對于品質(zhì)性狀, TST-TEM的預(yù)測均值(r=0.255)略高于TST-TST (r=0.238)。

3 討論

本研究基于基因組、轉(zhuǎn)錄組和代謝組數(shù)據(jù), 使用rrBLUP 和LASSO 兩種統(tǒng)計模型, 對玉米的農(nóng)藝性狀和品質(zhì)性狀進行了基因組預(yù)測分析。結(jié)果發(fā)現(xiàn),基于基因組數(shù)據(jù), 品質(zhì)性狀的預(yù)測能力整體高于農(nóng)藝性狀。對比不同組學(xué)數(shù)據(jù), 2 種類型性狀均呈現(xiàn)基因組數(shù)據(jù)預(yù)測精度高于轉(zhuǎn)錄組數(shù)據(jù)預(yù)測精度, 轉(zhuǎn)錄組數(shù)據(jù)預(yù)測精度高于代謝組數(shù)據(jù)預(yù)測精度的預(yù)測趨勢。對比不同模型, rrBLUP 在大多數(shù)農(nóng)藝性狀中表現(xiàn)出了良好的預(yù)測效果, 而LASSO 在一些品質(zhì)性狀中預(yù)測效果更好。整合分析發(fā)現(xiàn), 轉(zhuǎn)錄組數(shù)據(jù)對農(nóng)藝性狀預(yù)測更有價值, 代謝組數(shù)據(jù)對品質(zhì)性狀預(yù)測作用更大。此外, 使用親緣關(guān)系較遠的材料進行預(yù)測時, 加入代謝組學(xué)數(shù)據(jù)進行模型構(gòu)建, 可提升預(yù)測精度。

由于數(shù)據(jù)和模型的局限性, 本研究還存在值得改進和深入分析的地方。對于不同組學(xué)數(shù)據(jù), 由于轉(zhuǎn)錄組和代謝組數(shù)據(jù)維度要遠小于基因組數(shù)據(jù), 同時, 本研究使用的轉(zhuǎn)錄組數(shù)據(jù)僅來自玉米籽粒, 可能忽略了某些組織特異性基因的表達差異。此外,Zheng 等[40]曾證明禾本科的組蛋白修飾在植物生長和脅迫適應(yīng)中發(fā)揮著重要作用, 這表明表觀遺傳修飾在植物生長中的重要性, 而在本研究中未涉及到表觀組學(xué)的預(yù)測分析。因此, 未來可通過不斷補充組學(xué)數(shù)據(jù)類型和擴充組學(xué)數(shù)據(jù)維度, 進一步提高預(yù)測精度。對于不同模型, 本研究僅選擇了2 種穩(wěn)健的統(tǒng)計模型rrBLUP 與LASSO 進行分析和比較, 發(fā)現(xiàn)rrBLUP 更適合用于預(yù)測以農(nóng)藝性狀為代表的數(shù)量性狀, LASSO 更適合用于預(yù)測以品質(zhì)性狀為代表的質(zhì)量性狀。同時, Campos 等[41]曾驗證相同的數(shù)據(jù)集預(yù)測相同的目標(biāo)性狀可能會出現(xiàn)多個同樣好的預(yù)測模型, 很難直接確定GS 的最優(yōu)統(tǒng)計方法。因此,即使rrBLUP 和LASSO 在GS 中表現(xiàn)出更優(yōu)的性能,未來仍可以嘗試更多機器學(xué)習(xí)或神經(jīng)網(wǎng)絡(luò)算法, 以更好地學(xué)習(xí)多組學(xué)數(shù)據(jù)內(nèi)部的復(fù)雜非線性互作關(guān)系[24]。此外, 本研究僅是基于單個組學(xué)的預(yù)測效果分析, 而已有報道利用遺傳標(biāo)記作為變量并結(jié)合轉(zhuǎn)錄組和代謝組數(shù)據(jù), 進行 3 層遺傳特征學(xué)習(xí)的MLLASSO (multilayered least absolute shrinkage and selection operator)模型, 產(chǎn)量的預(yù)測精度(R2)可由0.1588 提升至0.2451[33], 表明基于中心法則, 借助機器學(xué)習(xí)算法優(yōu)勢, 建立時序性的多模態(tài)多組學(xué)預(yù)測模型[33], 是合理利用多組學(xué)數(shù)據(jù)精確預(yù)測玉米復(fù)雜性狀的潛在解決方案。

4 結(jié)論

本研究基于rrBLUP 和LASSO 兩種統(tǒng)計模型,分別利用基因組、轉(zhuǎn)錄組、代謝組數(shù)據(jù)進行基因組預(yù)測分析, 系統(tǒng)測試了多個玉米性狀在不同組學(xué)數(shù)據(jù)下預(yù)測能力的差異, 為多維組學(xué)數(shù)據(jù)進行基因組預(yù)測分析提供了數(shù)據(jù)基礎(chǔ)和參考依據(jù)。