自然群體多性狀表型缺失值預(yù)測方法的比較

王媛,溫陽俊,王艷萍,劉漢欽,馬若洵,吳清太,張瑾

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

缺失數(shù)據(jù)是影響數(shù)據(jù)分析準確性的重要原因之一,完整表型和基因型數(shù)據(jù)的分析對防治疾病、遺傳育種等研究具有重要意義。目前,在全基因組關(guān)聯(lián)研究(genome-wide association study,GWAS)中基因型數(shù)據(jù)的缺失已經(jīng)有了相對成熟的處理方法[1-5],而表型數(shù)據(jù)的缺失直接減少樣本量,削弱GWAS的統(tǒng)計推斷能力。表型的缺失值主要有3種處理方法:直接刪除法、重測量法和插補法。直接刪除法將數(shù)據(jù)中含有缺失值的樣本點全部刪除,從而得到一個不含缺失值的“完整”數(shù)據(jù)集。由于被刪除的樣本中可能包含重要信息,直接刪除會造成原始信息的損失,減少樣本量,降低后續(xù)統(tǒng)計推斷的功效。重測量法,即進行重復(fù)試驗,采集相同條件下的數(shù)據(jù),對原有缺失值進行插補。雖然重測量法能夠插補缺失的數(shù)據(jù),但耗費大量的人力和物力。數(shù)據(jù)插補法利用現(xiàn)有觀測數(shù)據(jù)信息,借助統(tǒng)計學(xué)方法預(yù)測缺失值以產(chǎn)生完整數(shù)據(jù)[6],插補表型缺失避免了原始信息的損失、增大了樣本量,具有省時高效的優(yōu)點。

精確的缺失表型預(yù)測可以提高后續(xù)統(tǒng)計分析的準確性,常用的方法有:均值插補法(mean)[7],利用觀測數(shù)據(jù)的均值插補缺失值,簡單易行,但是扭曲了樣本的分布,降低了數(shù)據(jù)的變異程度;回歸插補法[8],利用完整的觀測數(shù)據(jù)建立回歸模型,依據(jù)此模型預(yù)測表型缺失值;多重插補法[9],對缺失數(shù)據(jù)進行多次插補得到多組完整數(shù)據(jù)集,通過統(tǒng)計方法和綜合分析推斷出一組最佳完整數(shù)據(jù)集;EM算法[10],根據(jù)觀測數(shù)據(jù)的邊緣分布,利用EM算法進行極大似然估計,對含有缺失數(shù)據(jù)的不完整數(shù)據(jù)進行預(yù)測,但數(shù)據(jù)量大時收斂速度慢、計算繁瑣,效果不佳。

近年來,機器學(xué)習(xí)算法在缺失數(shù)據(jù)預(yù)測插補中也得到了應(yīng)用[11]。K近鄰(K-nearest neighbor,KNN)插補法[12],將最鄰近的k個樣本觀測值的加權(quán)平均值作為缺失數(shù)據(jù)的估計值,其插補效果易受異常值的影響,產(chǎn)生誤差較大;支持向量機[13]將完整數(shù)據(jù)集作為訓(xùn)練集,獲得預(yù)測模型后插補缺失數(shù)據(jù),大樣本情況下的插補精確度高;隨機森林[14]利用bootstrap抽取k個樣本集,組成k個決策樹,利用全部決策樹構(gòu)建的隨機森林插補缺失值;神經(jīng)網(wǎng)絡(luò)[15]將誤差反向傳播來訓(xùn)練多層網(wǎng)絡(luò),通過優(yōu)化網(wǎng)絡(luò)輸出不斷減小誤差,在隨機缺失數(shù)據(jù)插補中應(yīng)用廣泛。機器學(xué)習(xí)算法在處理數(shù)據(jù)缺失問題時速度快、泛化能力強,但由于訓(xùn)練數(shù)據(jù)對參數(shù)影響較大,插補能力還有待提高。此外,基于混合線性模型的缺失值插補方法[16-19]進一步提高了表型缺失值插補能力。

與單性狀插補方法相比,多性狀聯(lián)合插補分析可以利用性狀之間的遺傳結(jié)構(gòu)[20],并對多個性狀缺失值同時進行插補,省時高效。多重插補法(multiple imputation by chained equations,MICE)[9],利用觀測數(shù)據(jù)進行多變量回歸,以回歸的估計值插補缺失值,由于每個插補值都要進行多重插補計算,插補時間較長;PHENIX(phenotype imputation expediated)插補法[21]建立在貝葉斯多表型混合模型上,在親緣關(guān)系已知的情況下利用變分貝葉斯方法進行擬合;softImpute插補法[22]是一種利用核范數(shù)(或奇異值)正則化來插補矩陣缺失值的交替迭代算法,利用核范數(shù)定義損失函數(shù)并利用奇異值分解求解核范數(shù)使損失函數(shù)最小,精度高、速度快;多表型混合模型(multi-phenotype mixed models,MPMM)[16]是一種基于多個相關(guān)性狀的混合模型方法,利用其條件期望或最優(yōu)線性無偏預(yù)測對缺失的數(shù)據(jù)進行插補,具有精確度高的特點。多性狀基因組模型(multiple-trait genomic model,MTGM)[23]在單性狀模型的基礎(chǔ)上進行拓展,可應(yīng)用于估計育種值的預(yù)測問題。多性狀聯(lián)合插補方法是在利用變量之間關(guān)系的前提條件下對缺失值進行預(yù)測,結(jié)果往往比較準確,通常會明顯優(yōu)于單變量插補方法[24]。多表型插補技術(shù)已經(jīng)成功應(yīng)用于人類[20]、酵母、鼠、雞、小麥和大豆[25]等多個物種的缺失數(shù)據(jù)集上,有效提高了缺失數(shù)據(jù)的預(yù)測精度,進而提升全基因組關(guān)聯(lián)分析的計算效率和準確性。

表型缺失數(shù)據(jù)插補的準確性受到缺失率、樣本量等因素的影響,上述方法尚未針對這些因素對插補效果進行評判。本研究利用mean、KNN、決策樹(decision tree,DT)[26]、MICE、PHENIX和softImpute插補方法對多表型模擬缺失數(shù)據(jù)進行預(yù)測,以均方誤差(mean squared error,MSE)、皮爾遜相關(guān)系數(shù)(Pearson correlation coefficient,r)和計算時間為指標(biāo),比較在不同表型缺失率(phenotypic missing rate)、性狀數(shù)(D)、樣本量(n)和性狀相關(guān)性(cor)下的插補效果。使用6種方法對擬南芥實際數(shù)據(jù)的長日照花期(days to flowering under long day)、短日照花期(days to flowering under short day)、春化長日照花期(days to flowering under long day with vernalization)和春化短日照花期(days to flowering under short day with vernalization)的表型缺失值進行多性狀聯(lián)合插補,并對插補后的完整數(shù)據(jù)進行全基因組關(guān)聯(lián)分析,通過基因功能驗證插補數(shù)據(jù)的可靠性,查找出與目標(biāo)性狀關(guān)聯(lián)的已驗證基因。

1 材料與方法

1.1 模擬試驗材料

基于線性模型,生成等位基因頻率為[0.1,0.5]的10 000個單核苷酸多態(tài)性(single nucleotide polymorphism,SNP),隨機設(shè)置10個數(shù)量性狀核苷酸(quantitative trait nucleotide,QTN),其總遺傳率是60%,正態(tài)誤差方差設(shè)置為1.0。在此基礎(chǔ)上,設(shè)計了模擬試驗(表1)。

表1 4個模擬研究設(shè)計

利用表型數(shù)據(jù)插補方法對上述缺失數(shù)據(jù)進行預(yù)測,使用MSE、r和計算時間等指標(biāo)對插補數(shù)據(jù)的準確性進行評判。

1.2 真實試驗材料

真實數(shù)據(jù)集[27]來自擬南芥自然群體的199個個體,216 130個SNP(http://www.arabidopsis.usc.edu/),考慮長日照花期、短日照花期、春化長日照花期和春化短日照花期4個花期相關(guān)性狀(https://www.arabidopsis.org/portals/genAnnotation/index.jsp),其表型缺失率分別為16.0%、18.6%、15.6%和 20.1%。利用mean、KNN、DT、MICE、PHENIX和softImpute插補方法對缺失數(shù)據(jù)進行預(yù)測,并對插補后的完整數(shù)據(jù)集進行全基因組關(guān)聯(lián)分析得到顯著SNP,利用TAIR基因庫(https://www.arabidopsis.org)驗證位于顯著SNP附近20 kb的關(guān)聯(lián)基因,挖掘出與花期相關(guān)的已驗證基因。

1.3 插補方法

1)mean插補法[7]:用當(dāng)前變量的所有觀測值平均值代替其所有缺失值的方法。2)KNN插補法[12]:首先選出與缺失值距離最近的k個樣本觀測數(shù)據(jù),根據(jù)其距離加權(quán)計算出缺失數(shù)據(jù)的估計值,用該值來插補對應(yīng)的缺失值。KNN插補法基于R軟件程序包DMwR2實現(xiàn),k默認設(shè)為10。3)DT[26]插補法:利用觀測數(shù)據(jù)集構(gòu)建基于回歸樹的DT模型,然后將含有缺失值的數(shù)據(jù)代入DT模型進行插補。DT插補法基于R軟件程序包rpart實現(xiàn)。4)MICE插補法[9]:是一種重復(fù)計算缺失值的處理方法。它首先多次計算缺失值,并產(chǎn)生m個完整數(shù)據(jù)集;然后對這m個插補數(shù)據(jù)集進行統(tǒng)計分析,得到m個分析結(jié)果;最后綜合m個分析結(jié)果做出統(tǒng)計推斷,最終得到缺失值的最優(yōu)估計值。綜合起來,其主要步驟可總結(jié)為:插補、分析、綜合。MICE插補法基于R軟件程序包mice實現(xiàn)。5)PHENIX插補法[21]:充分利用表型之間的相關(guān)性和親緣關(guān)系矩陣,基于貝葉斯多表型混合線性模型,利用變分貝葉斯方法對模型進行擬合,根據(jù)觀測數(shù)據(jù)的后驗多變量正態(tài)分布得到缺失數(shù)據(jù)的近似估計值。PHENIX插補法基于R軟件程序包phenix實現(xiàn)。6)softImpute 插補法[22]:將矩陣完備化問題轉(zhuǎn)化為核范數(shù)最小二乘問題,利用軟閾值奇異值估算矩陣的缺失值。該算法的計算原理類似于EM算法,每一次迭代,首先通過對估計矩陣的核范數(shù)(或奇異值)的懲罰來估計低秩表型矩陣的近似矩陣,并利用估計值對缺失矩陣進行插補,然后利用軟閾值奇異值分解法,向目標(biāo)函數(shù)的最小值逼近,得到近似最優(yōu)解的集合,不斷迭代以對完整矩陣進行不斷優(yōu)化。softImpute插補法基于R軟件程序包softImpute實現(xiàn)。

1.4 評估指標(biāo)

采用均方誤差MSE和皮爾遜相關(guān)系數(shù)r來衡量基于不同表型缺失率、性狀數(shù)、樣本量和性狀相關(guān)性的插補準確性。MSE計算公式為:

(1)

(2)

2 結(jié)果與分析

2.1 模擬數(shù)據(jù)的缺失數(shù)據(jù)插補

利用mean、KNN、DT、MICE、PHENIX和softImpute插補方法對模擬缺失數(shù)據(jù)進行預(yù)測。針對不同表型缺失率、性狀數(shù)、樣本量和性狀相關(guān)性,比較各個方法插補效果的影響。

模擬研究1(表1)考察缺失率分別為5%、10%和20%時缺失數(shù)據(jù)的插補效果。結(jié)果表明,隨著表型缺失率的增大,6種方法的MSE(圖1-A)不斷增大,r不斷減小(圖1-B),這意味著缺失值插補的準確性隨著表型缺失率的增大而降低,這與事實相符。其中,PHENIX插補法準確性最優(yōu),隨著缺失率的增大,該方法在插補準確性上的優(yōu)勢變得越來越明顯,說明PHENIX插補法更適用于高缺失率的情形。KNN和DT插補法的準確性比PHENIX略微遜色,其次是softImpute和mean插補法,而MICE準確性受表型缺失率的影響最明顯,尤其在高缺失率的情況下表現(xiàn)不佳。

圖1 不同表型缺失率下mean、KNN、決策樹(DT)、MICE、PHENIX和softImpute插補值與觀測值的均方誤差(A)和皮爾遜相關(guān)系數(shù)(B)(重復(fù)100次)

模擬研究2(表1)考察性狀數(shù)(D)為2、3、5和10時缺失數(shù)據(jù)的插補效果。mean插補法的效果不受性狀數(shù)的影響,其他5種方法隨性狀數(shù)的增大,MSE(圖2-A)不斷減小,r不斷增大(圖2-B),說明上述多表型缺失數(shù)據(jù)插補方法均能夠利用多性狀之間的相關(guān)性和潛在信息,性狀數(shù)不斷增加,插補的準確性不斷提高。其中,PHENIX的插補效果最優(yōu),其次是DT和KNN,softImpute的插補效果隨性狀數(shù)的增加逐漸優(yōu)化,這4種插補方法的效果均優(yōu)于mean和MICE插補法。

圖2 不同性狀數(shù)下mean、KNN、DT、MICE、PHENIX和softImpute插補值與觀測值的均方誤差(A)和皮爾遜相關(guān)系數(shù)(B)(重復(fù)100次)

模擬研究3(表1)考察樣本量(n)為500、1 000和2 000時缺失數(shù)據(jù)的插補效果。在相同表型缺失率(10%)下,樣本量的增大對于插補精確度沒有顯著影響,但是MSE和r的箱線圖方差變小(圖3),說明樣本量越大,結(jié)果越穩(wěn)定。趨勢與上述模擬試驗相似,依然是PHENIX插補法的效果最好,其次是KNN和DT插補法,這3種插補方法的效果均優(yōu)于softImpute、mean和MICE插補法。

圖3 不同樣本量下mean、KNN、DT、MICE、PHENIX和softImpute插補值與觀測值的均方誤差(A)和皮爾遜相關(guān)系數(shù)(B)(重復(fù)100次)

模擬試驗4(表1)考察性狀相關(guān)性(cor)為0.1、0.3、0.5和0.8時缺失數(shù)據(jù)的插補效果。結(jié)果表明mean插補法的效果不受性狀之間相關(guān)性強弱的影響,而其他5種方法隨著性狀之間相關(guān)性的增加,MSE(圖4-A)不斷減小,r不斷增大(圖4-B)。這意味著性狀之間的相關(guān)性越強,KNN、DT、MICE、PHENIX和softImpute方法越能夠利用相似的遺傳結(jié)構(gòu),插補效果越好。在弱相關(guān)的情況下,softImpute方法和KNN、DT、PHENIX插補法的效果相似,具有較高的準確性,隨著性狀相關(guān)性的增大,softImpute插補準確性與它們差距變大,上述方法的準確性均優(yōu)于mean和MICE插補法。

圖4 不同性狀相關(guān)性下mean、KNN、DT、MICE、PHENIX和softImpute插補值與觀測值的均方誤差(A)和皮爾遜相關(guān)系數(shù)(B)(重復(fù)100次)

從模擬試驗的計算時間(表2)可見:mean和softImpute的計算速度最快,均小于2 s,隨后是KNN和DT插補法,計算時間為3～40 s,而MICE和PHENIX插補時間相對比較長。由表2也可以看出:數(shù)據(jù)的表型缺失率和性狀相關(guān)性對所有方法的運算時間影響均不顯著,而隨性狀數(shù)和樣本量的增大,運算時間呈現(xiàn)顯著增加的趨勢。MICE插補法對于性狀數(shù)較為敏感,隨性狀數(shù)增加,MICE的計算時間呈現(xiàn)出指數(shù)級的增長;而PHENIX插補法對于樣本量較為敏感,計算時間與樣本量呈現(xiàn)出指數(shù)級增長的關(guān)系。

表2 4種模擬研究100次模擬數(shù)據(jù)集的mean、KNN、DT、MICE、PHENIX和softImpute的計算時間

2.2 真實數(shù)據(jù)的缺失數(shù)據(jù)插補

真實擬南芥數(shù)據(jù)包括199個個體、216 130個SNP,4個與花期相關(guān)的性狀(LD、SD、LDV和SDV),其表型缺失率分別為16.0%、18.6%、15.6%和20.1%,4個性狀之間具有較強的相關(guān)性,相關(guān)系數(shù)為0.65～0.85。

采用mean、KNN、DT、MICE、PHENIX和softImpute插補方法對真實數(shù)據(jù)進行聯(lián)合插補,并對插補后的數(shù)據(jù)進行全基因組關(guān)聯(lián)分析。利用TAIR基因庫驗證位于顯著SNP位點(Bonferroni校正)附近20 kb的關(guān)聯(lián)基因,挖掘出與花期相關(guān)的已驗證基因。由圖5可以看出:PHENIX、DT和KNN檢測的已驗證基因數(shù)最多,分別為22、18和18,并且共同檢測到超過17個已驗證基因,說明3種方法能夠利用擬南芥真實數(shù)據(jù)的復(fù)雜遺傳結(jié)構(gòu),得到更接近真實情況的插補值。mean檢測的已驗證基因數(shù)為14,并且與其他方法共同檢測的基因較多,softImpute和MICE檢測的已驗證基因數(shù)較少,分別為7和2。

圖5 利用mean、KNN、DT、MICE、PHENIX和softImpute插補數(shù)據(jù)和觀測數(shù)據(jù)檢測以及共同檢測的已驗證基因數(shù)

6種插補方法中,mean、KNN、DT和softImpute插補法的計算時間在同一個數(shù)量級上,均小于0.1 s,PHENIX插補法計算時間略長,約0.4 s,而MICE插補法長于其他5種方法,為3.6 s,這與模擬數(shù)據(jù)的計算時間一致。由分析結(jié)果、檢測準確性和運算時間等指標(biāo)可以看出,PHENIX、DT和KNN多性狀聯(lián)合分析插補法可以充分利用性狀之間的相關(guān)性,插補效果優(yōu)于其余方法,并適宜在真實數(shù)據(jù)分析中進行表型缺失的插補。

2.3 數(shù)量性狀表型插補前后關(guān)聯(lián)分析結(jié)果的比較

表型數(shù)據(jù)插補能夠利用現(xiàn)有觀測數(shù)據(jù)信息,預(yù)測缺失值以產(chǎn)生完整數(shù)據(jù),避免直接刪除數(shù)據(jù)帶來的信息損失,增大樣本量,提高后續(xù)研究的準確性。因此,本研究對擬南芥實際數(shù)據(jù)中表型數(shù)據(jù)插補前、后的結(jié)果進行比較,利用關(guān)聯(lián)分析檢測到的顯著SNP進行已驗證基因篩選,原始觀測數(shù)據(jù)篩選的已驗證基因數(shù)為15,少于PHENIX、KNN和DT方法,比mean、softImpute和MICE方法檢測到的已驗證基因數(shù)多(圖5)。這進一步說明了PHENIX、DT和KNN插補方法得到的缺失數(shù)據(jù)估計值更加準確、更接近真實值,因此,其關(guān)聯(lián)分析結(jié)果功效更高、檢測的已驗證基因數(shù)更多,而mean、softImpute和MICE方法的結(jié)果不能令人滿意。

3 結(jié)論與討論

數(shù)據(jù)缺失是進行全基因組關(guān)聯(lián)分析時普遍存在且難以避免的問題,不但增大了統(tǒng)計分析難度,而且還會造成分析結(jié)果的偏倚,使結(jié)果無法準確反映總體的真實特征,甚至得出錯誤的分析結(jié)果。對于不同類型的缺失數(shù)據(jù),選擇合適的預(yù)測方法能夠為后續(xù)分析結(jié)果的準確性和可靠性奠定基礎(chǔ)。因此,比較不同插補方法處理不同缺失樣本的效果具有重要意義。

本研究針對不同研究目的生成4組模擬數(shù)據(jù)集,采用mean、KNN、DT、MICE、PHENIX和softImpute插補方法,分析多性狀的表型缺失率、性狀數(shù)、樣本量和性狀相關(guān)性對表型插補準確性的影響。模擬研究結(jié)果顯示,在不同的試驗中PHENIX、KNN和DT的插補效果較好,具有較小的MSE和較高的r;其次是softImpute插補方法;而mean插補法僅僅利用了表型的平均值,降低了表型的變異程度,表現(xiàn)不佳;MICE插補法的MSE最大、r最小,其結(jié)果不能令人滿意。同時,利用6種插補方法對實際缺失數(shù)據(jù)進行多性狀聯(lián)合插補,并通過全基因組關(guān)聯(lián)分析方法驗證了插補方法的可靠性。結(jié)果顯示多性狀聯(lián)合插補方法PHENIX、KNN和DT可以利用性狀之間的遺傳結(jié)構(gòu),檢測到的已驗證基因數(shù)最多,并且它們共同檢測到的基因較多。

從計算時間的角度來看,mean、KNN、DT和softImpute插補法在模擬分析和實際數(shù)據(jù)分析時均比較快,MICE和PHENIX插補法的計算時間相對較長,MICE插補法的計算時間受性狀數(shù)的影響比較大,而PHENIX插補法受樣本量的影響比較大。

不同的缺失數(shù)據(jù)處理方法有不同的特點,進行插補方法選擇時,應(yīng)盡量全面了解數(shù)據(jù)的背景和特征,盡可能利用各種輔助信息來進行缺失值的預(yù)測插補,如性狀之間的相關(guān)性等,有利于提高缺失數(shù)據(jù)的插補準確性。