基于FASTmrEMMA、最小角回歸和隨機(jī)森林的全基因組選擇新算法

孫嘉利,吳清太,溫陽俊,張瑾

(南京農(nóng)業(yè)大學(xué)理學(xué)院,江蘇 南京 210095)

許多動(dòng)植物的重要性狀和人類的復(fù)雜疾病都是數(shù)量性狀,檢測(cè)、發(fā)掘與數(shù)量性狀關(guān)聯(lián)的基因并進(jìn)行全基因組選擇,對(duì)動(dòng)植物的性狀分析及人類疾病的預(yù)防和治療都具有重要意義[1-5]。全基因組選擇(genomic selection,GS)又稱全基因組預(yù)測(cè)(genomic prediction,GP),它是Meuwissen等[6]提出的一種利用全基因組高密度分子標(biāo)記數(shù)據(jù)及原始訓(xùn)練群體表型數(shù)據(jù)建立預(yù)測(cè)模型,估計(jì)每個(gè)分子標(biāo)記遺傳效應(yīng),從而預(yù)測(cè)個(gè)體全基因組育種值(genomic breeding value,GBV)的方法。由于GS無需鑒定與目標(biāo)性狀顯著關(guān)聯(lián)的位點(diǎn),即使微效標(biāo)記也都能被捕獲,只需個(gè)體遺傳信息即可計(jì)算其育種值,這大大縮短了育種周期,目前已廣泛應(yīng)用于動(dòng)植物的全基因組選擇[7-10]。

近年來,研究人員陸續(xù)提出了許多GS方法,例如,應(yīng)用廣泛的最佳線性無偏預(yù)測(cè)(best linear unbiased prediction,BLUP),它是基于系譜信息來定義個(gè)體間親緣關(guān)系矩陣,校正環(huán)境和非隨機(jī)交配造成的偏差,從而提供個(gè)體育種值的無偏估計(jì)值。GBLUP(genomic BLUP)[11]是利用全基因組遺傳標(biāo)記計(jì)算遺傳關(guān)系矩陣來預(yù)測(cè)表型缺失值。RR-BLUP(ridge regression BLUP)[6]是利用壓縮估計(jì)獲得標(biāo)記效應(yīng)來進(jìn)行全基因組預(yù)測(cè),但可能會(huì)被過度壓縮。BayesA和BayesB[6]可以有效彌補(bǔ)RR-BLUP的不足,它們假設(shè)每個(gè)標(biāo)記效應(yīng)服從先驗(yàn)分布,其中BayesB的先驗(yàn)分布為混合分布。此外,還有BayesC[12]、Bayesian LASSO[13]等方法。上述GS大多基于壓縮估計(jì)方法,其計(jì)算速度會(huì)隨著數(shù)據(jù)維度的增加而受到限制,不適用于現(xiàn)代測(cè)序技術(shù)產(chǎn)生的海量遺傳標(biāo)記數(shù)據(jù),且尚未考慮群體結(jié)構(gòu)和多基因背景,從而影響全基因組預(yù)測(cè)的準(zhǔn)確度。

植物交配設(shè)計(jì)相對(duì)復(fù)雜,考慮群體結(jié)構(gòu)和多基因背景可以有效提高全基因組預(yù)測(cè)的準(zhǔn)確度。TSLRF[14]是一種結(jié)合了最小角回歸(least angle regression,LARS)[15]和隨機(jī)森林(random forest,RF)[16]方法的兩階段關(guān)聯(lián)分析算法。它首先應(yīng)用FASTmrEMMA方法[17]進(jìn)行群體結(jié)構(gòu)和多基因背景控制,再通過LARS方法選擇與目標(biāo)性狀潛在關(guān)聯(lián)的單核苷酸多態(tài)(single nucleotide polymorphism,SNP),并用此SNP子集建立RF模型預(yù)測(cè)個(gè)體的GBV,即全基因組預(yù)測(cè)值。本研究以擬南芥自然群體的模擬數(shù)據(jù)和真實(shí)數(shù)據(jù)為研究對(duì)象,利用TSLRF方法針對(duì)植物數(shù)量性狀進(jìn)行全基因組預(yù)測(cè),其結(jié)果與TSRF[18]、RF和GBLUP進(jìn)行比較,驗(yàn)證TSLRF方法對(duì)植物數(shù)量性狀全基因組預(yù)測(cè)的準(zhǔn)確性和可行性,為后續(xù)的分子育種和優(yōu)異親本組合預(yù)測(cè)提供理論依據(jù)和現(xiàn)實(shí)參考。

1 材料與方法

1.1 模擬試驗(yàn)材料

以擬南芥數(shù)據(jù)[19]為模擬試驗(yàn)材料。按照文獻(xiàn)[20]的方法選取其中的10 000個(gè)SNP標(biāo)記,它們分別來自于5條染色體片段,其位置分別介于染色體1至5的11 226 256～12 038 776 bp、5 045 828～6 412 875 bp、1 916 588～3 196 442 bp、2 232 796～3 143 893 bp和19 999 868～21 039 406 bp[20]。選取6個(gè)QTN(quantitative trait nucleotide),分別設(shè)置在染色體1的11 298 364、11 655 607 bp標(biāo)記上(遺傳率設(shè)置為10%、5%)和染色體2的5 066 968、5 134 228、5 464 675、6 137 189 bp標(biāo)記上(遺傳率設(shè)置為5%、15%、5%、5%),其等位基因頻率均為0.30;總體均值和殘差均設(shè)置為10.0,樣本容量為199,隨機(jī)生成1 000次重復(fù)表型數(shù)據(jù)。該模擬數(shù)據(jù)集已應(yīng)用于多項(xiàng)研究[17,20-21]。

為考察不同表型缺失率的預(yù)測(cè)效果,針對(duì)上述模擬數(shù)據(jù)集的199個(gè)個(gè)體,分別以5%、10%、15%和20%的缺失率隨機(jī)刪除個(gè)體表型數(shù)據(jù),構(gòu)造缺失數(shù)據(jù)集,用來比較TSLRF、TSRF、RF和GBLUP這4種方法在不同缺失率下的全基因組預(yù)測(cè)準(zhǔn)確度(MAE和MSE)、預(yù)測(cè)模型擬合度和計(jì)算時(shí)間。

1.2 真實(shí)試驗(yàn)材料

真實(shí)數(shù)據(jù)集[19]來自擬南芥自然群體的199個(gè)個(gè)體,216 130個(gè)SNP(http://www.arabidopsis.usc.edu/),考慮長(zhǎng)日照花期(LD)、春化長(zhǎng)日照花期(LDV)和短日照花期(SD)3個(gè)花期相關(guān)性狀(https://www.arabidopsis.org/portals/genAnnotation/index.jsp),其表型缺失率分別為16.0%、15.6%和18.6%。利用TSLRF、TSRF、RF和GBLUP預(yù)測(cè)表型值,并利用預(yù)測(cè)值和已有觀測(cè)值進(jìn)行全基因組關(guān)聯(lián)分析,挖掘與性狀關(guān)聯(lián)的基因。

1.3 TSLRF方法

TSLRF先將模型中的群體結(jié)構(gòu)Q矩陣設(shè)為固定效應(yīng)來校正群體結(jié)構(gòu)效應(yīng),并借助FASTmrEMMA[17]將個(gè)體間多基因背景和誤差變異轉(zhuǎn)換成標(biāo)準(zhǔn)正態(tài)離差來校正多基因背景,便于后續(xù)進(jìn)行變量選擇;然后用LARS選擇出最可能與目標(biāo)性狀關(guān)聯(lián)的SNP子集;再用此SNP子集建立隨機(jī)森林模型預(yù)測(cè)個(gè)體的GBV。其過程如圖1所示。

圖1 基于最小角回歸和隨機(jī)森林的兩階段算法 (TSLRF)的流程圖Fig.1 The flow chart of two-stage algorithm based on least angle regression and random forest(TSLRF) LARS:最小角回歸Least angle regression;SNP:單核苷酸多態(tài)Single nucleotide polymorphism;RF:隨機(jī)森林Random forest;GBV:全基因組育種值Genomic breeding value.

1.3.1 遺傳模型遺傳模型如下:

y=1μ+Wα+Zγ+u+ε

(1)

式中:y=(y1,…,yn)T為自然群體的表型向量,n為樣本容量;1表示分量為1的n×1向量;μ為總體均值;α表示固定的群體結(jié)構(gòu)效應(yīng);γ表示p×1的標(biāo)記隨機(jī)效應(yīng)向量,p是標(biāo)記數(shù)量;W和Z分別是α和γ對(duì)應(yīng)的設(shè)計(jì)矩陣;u表示n×1的多基因隨機(jī)效應(yīng)向量;ε為n×1的誤差向量。

(2)

(3)

yc=1μ+Wcα+Zcγ+εc

(4)

式中:yc=Cy;Wc=CW;Zc=CZ;εc=Cu+Cε～MVNn(0,σ2In)。

Y=Xβ+εc

(5)

式中:Y為校正后的表型值;X為標(biāo)準(zhǔn)化后的基因型;β為包含固定效應(yīng)和隨機(jī)效應(yīng)的向量;εc為服從MVNn(0,σ2In)的誤差向量。

1.3.3 LARS算法采用LARS算法選擇出與表型Y具有絕對(duì)相關(guān)最大的變量,如xj1,則當(dāng)前路徑為s1*xj1,其中s1是xj1為與當(dāng)前殘差之間相關(guān)系數(shù)的符號(hào)。按照這一路徑前進(jìn),直到出現(xiàn)另一變量與當(dāng)前殘差的相關(guān)系數(shù)更大,記為xj2,將其加入回歸模型中并重復(fù)上述過程。設(shè)置LARS算法迭代到第k步,選擇出k個(gè)與目標(biāo)性狀潛在關(guān)聯(lián)的SNP來構(gòu)建RF模型。

1.3.4 RF算法針對(duì)LARS選擇出的變量集{xi1,xi2,…,xik}(i=1,2,…,n)以及原始表型觀測(cè)值,RF先使用Bagging算法生成500個(gè)自采樣樣本數(shù)據(jù)集,每次未被抽到的樣本組成500個(gè)袋外(out of bag,OOB)樣本數(shù)據(jù)集;然后通過節(jié)點(diǎn)隨機(jī)分裂技術(shù),從每棵樹的每個(gè)節(jié)點(diǎn)處的k個(gè)變量中隨機(jī)抽取mtry個(gè)變量作為備選分支變量(回歸樹一般取mtry=k/3),按照節(jié)點(diǎn)分裂原則,選擇最優(yōu)分支進(jìn)行分裂。生成500棵CART樹構(gòu)建RF模型,并預(yù)測(cè)所有個(gè)體的全基因組育種值。其中,mtry選取能使RF模型獲得最佳模型擬合度的參數(shù)值。本文中,RF將所有CART樹的平均值輸出作為預(yù)測(cè)值。

TSLRF可通過R軟件實(shí)現(xiàn),其中LARS和RF可分別通過程序包lars和randomForest實(shí)現(xiàn)。

1.4 評(píng)估指標(biāo)

采用平均絕對(duì)誤差(MAE)和均方誤差(MSE)來評(píng)估預(yù)測(cè)準(zhǔn)確度,用皮爾森相關(guān)系數(shù)(r)評(píng)估預(yù)測(cè)模型的擬合度。MAE計(jì)算公式為:

(6)

(7)

與MAE相比,MSE取誤差值的平方,增大了誤差的作用,即對(duì)誤差的估計(jì)更加敏感,是預(yù)測(cè)效果評(píng)估中具有代表性的指標(biāo)。MSE越大表示準(zhǔn)確度越低,越小表示準(zhǔn)確度越高。

皮爾森相關(guān)系數(shù)(r)計(jì)算公式為:

(8)

2 結(jié)果與分析

2.1 模擬數(shù)據(jù)的全基因組選擇

為了驗(yàn)證TSLRF的全基因組預(yù)測(cè)能力,利用TSLRF、TSRF、RF和GBLUP對(duì)不同表型缺失率(5%、10%、15%和20%)的模擬數(shù)據(jù)集進(jìn)行分析。這4種方法都是基于R軟件實(shí)現(xiàn)的,其中,3種隨機(jī)森林方法(TSLRF、TSRF和RF)通過程序包randomForest實(shí)現(xiàn);GBLUP通過程序包GAPIT實(shí)現(xiàn)。

計(jì)算1 000次重復(fù)的全基因組預(yù)測(cè)值與觀測(cè)值之間的MAE和MSE。結(jié)果表明,隨著表型缺失率的增大,4種方法的MAE(圖2-A)和MSE(圖2-B)增大。這意味著隨著表型缺失率的增大,全基因組預(yù)測(cè)準(zhǔn)確度不斷降低,并呈現(xiàn)出越來越顯著的差異,這與事實(shí)相符。在4種缺失率下,TSLRF的MAE顯著小于其他方法(方差分析檢驗(yàn),P=1.28×10-6),其次為TSRF和RF;GBLUP的MAE最大,在缺失率為5%和10%的情況下,大約是TSLRF的2倍。這說明TSLRF預(yù)測(cè)值與實(shí)際觀測(cè)值間的MAE最小,即預(yù)測(cè)準(zhǔn)確度最高。4種方法MSE的總體趨勢(shì)與MAE的情況相似。在4種缺失率下,RF的r最大,位于0.88～0.95;TSLRF在5%、10%、15%表型缺失率下僅次于RF,位于0.90～0.94,在20%缺失率時(shí),TSLRF(0.89)要稍優(yōu)于RF(0.88);而TSRF和GBLUP的全基因組預(yù)測(cè)能力要低于RF和TSLRF,TSRF的r位于0.85～0.91,GBLUP的r最小,位于0.84～0.91(圖2-C)?？傮w來說,RF的預(yù)測(cè)值與觀測(cè)值之間的皮爾森相關(guān)系數(shù)最大,TSLRF僅次于RF,而TSRF和GBLUP的模型擬合度不能令人滿意。

圖2 不同表型缺失率下TSLRF、TSRF、RF和GBLUP全基因組選擇的平均絕對(duì)誤差(A)、 均方誤差(B)和皮爾森相關(guān)系數(shù)(C)(重復(fù)1 000次)Fig.2 The mean absolute error(A),mean squared error(B)and Pearson correlation coefficient(C)of TSLRF,TSRF, RF and GBLUP in the genome-wide selection under various phenotypic missing rates(1 000 repeats) TSLRF:基于最小角回歸和隨機(jī)森林的兩階段算法Two-stage algorithm based on least angle regression and random forest;TSRF:基于隨機(jī)森林的兩階段變量選擇Two-stage stepwise variable selection based on random forest;RF:隨機(jī)森林Random forest;GBLUP:全基因組最佳線性無偏預(yù)測(cè)Genomic best linear unbiased prediction.

由表1可見:隨著表型缺失率的增大,4種方法的計(jì)算時(shí)間逐漸降低。在相同表型缺失率下,TSLRF和RF的計(jì)算時(shí)間均少于100 min,其中RF的計(jì)算時(shí)間最短,TSLRF僅次于RF,這是由于TSLRF進(jìn)行全基因組的群體結(jié)構(gòu)和多基因背景控制需要花費(fèi)一定的時(shí)間;TSRF的計(jì)算時(shí)間約為200 min;GBLUP的計(jì)算時(shí)間最長(zhǎng),約為400 min。

表1 不同表型缺失率下1 000次模擬數(shù)據(jù)集的TSLRF、TSRF、RF和GBLUP計(jì)算時(shí)間

2.2 真實(shí)數(shù)據(jù)的全基因組選擇

真實(shí)擬南芥數(shù)據(jù)包括199個(gè)個(gè)體,216 130個(gè)SNP,3個(gè)與花期相關(guān)的LD、LDV和SD性狀,其表型缺失率分別為16.0%、15.6%和18.6%。

采用TSLRF、TSRF、RF和GBLUP方法預(yù)測(cè)表型值,利用表型預(yù)測(cè)值和觀測(cè)值進(jìn)行關(guān)聯(lián)分析,得到顯著SNP并挖掘基因,考察這4種方法的全基因組預(yù)測(cè)能力。分別將TSLRF、TSRF和RF方法中重要性排名前20及GBLUP方法中效應(yīng)估計(jì)值排名前20的SNP視為QTN。利用TAIR基因庫(kù)(https://www.arabidopsis.org)查找位于QTN附近20 kb的關(guān)聯(lián)基因。從表2可見:TSLRF、GBLUP、TSRF、RF預(yù)測(cè)到被證實(shí)的基因數(shù)分別為40、30、24和21個(gè)。TSLRF基因預(yù)測(cè)能力最強(qiáng),GBLUP次之。

表2 4種方法預(yù)測(cè)的擬南芥LD、LDV和SD性狀前20個(gè)顯著QTN附近的已知基因數(shù)

TSLRF預(yù)測(cè)出與同一性狀顯著關(guān)聯(lián)的多個(gè)基因簇(表3),如, 位于染色體1的2 005 921 bp標(biāo)記附近的基因AT1G06515和AT1G06520被預(yù)測(cè)到與LD相關(guān);位于染色體5的26 004 094 bp標(biāo)記附近的基因AT5G65050、AT5G65060、AT5G65070和AT5G65080被預(yù)

表3 TSLRF預(yù)測(cè)的擬南芥LD、LDV和SD性狀前20個(gè)顯著QTN附近的已知基因Table 3 The known genes in the vicinity of the top 20 significant QTN for each of the LD,LDV and SD traits in A.thaliana using the TSLRF method

測(cè)到與LDV相關(guān)。TSLRF預(yù)測(cè)到的一些顯著關(guān)聯(lián)的基因,同時(shí)被其他方法預(yù)測(cè)到,例如,與LDV顯著相關(guān)的位于染色體2的9 965 680、9 964 693、9 966 813和9 967 491 bp標(biāo)記附近的基因AT2G23380可被TSLRF和GBLUP預(yù)測(cè)到;與SD顯著相關(guān)的位于染色體4的13 996 527、13 997 089和13 998 150 bp標(biāo)記附近的基因AT4G28190和AT4G28280可被TSLRF和TSRF同時(shí)預(yù)測(cè)到。由此可見,與其他3種方法相比,TSLRF具有更有效、更準(zhǔn)確的基因預(yù)測(cè)能力。

4種方法對(duì)真實(shí)數(shù)據(jù)分析的計(jì)算時(shí)間如表4所示。結(jié)果表明,TSLRF和RF的計(jì)算時(shí)間最短,均小于2 min,這與模擬研究的結(jié)果一致;GBLUP的計(jì)算時(shí)間約為25 min;TSRF的計(jì)算時(shí)間最長(zhǎng),至少需要60 min。因此,從運(yùn)行效率來看,利用TSLRF進(jìn)行全基因組預(yù)測(cè)是可行的。

表4 利用4種全基因組預(yù)測(cè)擬南芥LD、LDV和SD性狀的計(jì)算時(shí)間Table 4 Computing times for the LD,LDV and SD traits in A.thaliana using 4 methods min

3 結(jié)論與討論

針對(duì)動(dòng)植物重要性狀和人類復(fù)雜疾病進(jìn)行全基因組預(yù)測(cè),對(duì)作物改良和疾病防治具有重要意義。本文利用TSLRF方法進(jìn)行全基因組預(yù)測(cè),其先校正群體結(jié)構(gòu)和多基因背景,再結(jié)合變量選擇方法LARS和機(jī)器學(xué)習(xí)方法RF進(jìn)行預(yù)測(cè),既能顯著提高模型準(zhǔn)確度與擬合度,又能顯著提升全基因組的預(yù)測(cè)能力與效率。

模擬研究和真實(shí)數(shù)據(jù)分析均表明,表型缺失率是影響預(yù)測(cè)準(zhǔn)確度和后續(xù)關(guān)聯(lián)分析檢測(cè)能力的重要原因之一。隨表型缺失率增大,計(jì)算時(shí)間減少,預(yù)測(cè)值與觀測(cè)值之間的MAE和MSE增大,全基因組預(yù)測(cè)準(zhǔn)確度顯著降低,模型擬合度也降低。真實(shí)數(shù)據(jù)分析中,根據(jù)隨機(jī)森林相關(guān)研究[22-23],本文以經(jīng)驗(yàn)排名選擇顯著SNP個(gè)數(shù),分別選取前20、30、50的SNP。排名20以后(21～50)的SNP重要性得分明顯下降,與非關(guān)聯(lián)SNP之間差距不顯著;且被TAIR基因庫(kù)證實(shí)的基因數(shù)也未顯著增加。因此,本文選取排名前20的SNP作為顯著關(guān)聯(lián)標(biāo)記。

從全基因組預(yù)測(cè)運(yùn)行速度來看,TSLRF和RF比其他2種方法(TSRF和GBLUP)快。GBLUP在計(jì)算遺傳關(guān)系矩陣上需要花費(fèi)一定的時(shí)間,因此它的總計(jì)算時(shí)間最長(zhǎng);TSRF由于重復(fù)采樣、多次計(jì)算重要性得分,計(jì)算時(shí)間也較長(zhǎng),大約是TSLRF的60倍。綜上所述,結(jié)合全基因組預(yù)測(cè)的準(zhǔn)確度、擬合度以及計(jì)算時(shí)間等指標(biāo),TSLRF是高效、可行的全基因組選擇方法,這為機(jī)器學(xué)習(xí)方法應(yīng)用于海量數(shù)據(jù)的全基因組預(yù)測(cè)和全基因組關(guān)聯(lián)分析提供了新的途徑。