多位點(diǎn)關(guān)聯(lián)分析方法學(xué)的研究進(jìn)展

溫陽(yáng)俊,馮建英,張瑾,*

(1.南京農(nóng)業(yè)大學(xué)理學(xué)院,江蘇 南京 210095;2.南京農(nóng)業(yè)大學(xué)作物遺傳與種質(zhì)創(chuàng)新國(guó)家重點(diǎn)實(shí)驗(yàn)室,江蘇 南京 210095)

近年來,多位點(diǎn)模型(multi-locus model)在全基因組關(guān)聯(lián)分析(genome-wide association study,GWAS)中引起了廣泛關(guān)注[1-2],它在多重檢驗(yàn)、群體結(jié)構(gòu)和多基因背景控制方面都體現(xiàn)了一定優(yōu)勢(shì)。其中,以混合線性模型為框架的多位點(diǎn)GWAS方法,比如,多位點(diǎn)隨機(jī)SNP效應(yīng)混合線性模型(multi-locus random-SNP-effect mixed linear model,mrMLM)[3]、快速多位點(diǎn)隨機(jī)SNP效應(yīng)有效混合模型關(guān)聯(lián)(FAST multi-locus random-SNP-effect efficient mixed model association,FASTmrEMMA)[4-5]以及全基因組復(fù)合區(qū)間作圖(genome-wide composite interval mapping,GCIM)[5-7]能夠進(jìn)一步提高數(shù)量性狀核苷酸(quantitative trait nucleotide,QTN)的檢測(cè)功效,降低假陽(yáng)性率,特別適合小效應(yīng)和連鎖QTN的檢測(cè)。多位點(diǎn)方法在高粱[8]、水稻[9]、玉米[10-11]、棉花[12]、大麥[13]、大豆[14]和桃[15]等植物中已廣泛應(yīng)用,并挖掘出新的顯著QTN及候選基因。為此,筆者對(duì)多位點(diǎn)關(guān)聯(lián)分析方法學(xué)的研究進(jìn)展進(jìn)行綜述。

1 多位點(diǎn)GWAS方法學(xué)的研究進(jìn)展

1.1 GWAS混合線性模型方法學(xué)的建立

混合線性模型(mixed linear model,MLM)或線性混合模型(linear mixed model,LMM)是指同時(shí)含有固定效應(yīng)和隨機(jī)效應(yīng)的線性模型,其中隨機(jī)效應(yīng)和誤差具有方差,因此也稱之為方差分量模型[16]。

2005年以前,GWAS以基于Case-Control數(shù)據(jù)和家系數(shù)據(jù)的算法為主。雖然有人建議采用MLM方法,但相關(guān)成果未見報(bào)道。為此,Zhang等[17]利用品種資源群體品種系譜、分子標(biāo)記和數(shù)量性狀表型信息,通過品種系譜計(jì)算品種間數(shù)量性狀基因座(quantitative trait locus,QTL)的后裔同樣(identity-by-descent,IBD)矩陣以及多基因背景的加性親緣關(guān)系矩陣,建立了GWAS的混合線性模型方法。此方法把QTL效應(yīng)看成隨機(jī)效應(yīng),模型中包含了QTL效應(yīng)方差、多基因背景方差和誤差方差共3個(gè)方差分量。

在自然群體研究中,由于品種系譜信息缺失而導(dǎo)致結(jié)果不可靠。針對(duì)這一問題,Yu等[18]改進(jìn)了Zhang等[17]的方法,提出由分子標(biāo)記信息計(jì)算的品種間親緣關(guān)系K矩陣代替品種間IBD矩陣,并引入群體結(jié)構(gòu)Q矩陣。這稱為關(guān)聯(lián)分析的Q+K混合線性模型方法。它把QTN效應(yīng)看成固定效應(yīng),模型中僅包含多基因背景方差和誤差方差2個(gè)方差分量。

Kang等[19]建立了有效混合模型關(guān)聯(lián)(efficient mixed-model association,EMMA)方法,被視作混合線性模型GWAS算法的黃金標(biāo)準(zhǔn)。對(duì)方差分量估計(jì)來說屬于精確算法,即檢測(cè)每一個(gè)SNP都需要重新估計(jì)多基因背景方差和誤差方差的比值。它利用譜分解將似然函數(shù)和參數(shù)估計(jì)量都表示為特征值、特征向量、標(biāo)記信息和表型值的標(biāo)量運(yùn)算形式,極大優(yōu)化了似然函數(shù)的求解,速度比早期Yu等[18]提出的方法要快很多。盡管與現(xiàn)在的算法相比,速度不算快,但需要放在歷史背景下進(jìn)行比較。

1.2 GWAS混合線性模型方法學(xué)的發(fā)展

在MLM方法中,如果SNP標(biāo)記數(shù)太多和群體容量太大,則計(jì)算的運(yùn)行時(shí)間長(zhǎng)。為了降低運(yùn)行時(shí)間,提高QTN檢測(cè)功效,涌現(xiàn)了許多以單位點(diǎn)(single-locus)分析為主的快速M(fèi)LM方法。標(biāo)準(zhǔn)的單位點(diǎn)混合線性模型為:

y=Wα+Xβ+Zu+ε

(1)

Zhang等[22]在Q+K的MLM關(guān)聯(lián)分析基礎(chǔ)上,將品種聚類成幾組(壓縮方法),用組間親緣關(guān)系矩陣代替?zhèn)€體間親緣關(guān)系矩陣,這稱為壓縮MLM(compressed mixed linear model,CMLM)方法;同時(shí),在全基因組標(biāo)記掃描時(shí)固定多基因背景方差與誤差方差的比值,即P3D(population parameters previously determined)算法,它避免了多基因背景方差分量的重復(fù)估計(jì),節(jié)約了計(jì)算時(shí)間。P3D對(duì)方差分量估計(jì)來說屬于近似算法。后來Li等[23]從8種聚類算法(非加權(quán)組平均、非加權(quán)組質(zhì)心、最長(zhǎng)距離、可變類平均、加權(quán)的類間平均、加權(quán)的類間重心、最短距離和離差平方和)與3種類間親緣關(guān)系算法(平均數(shù)、最大值和中位數(shù))的24種組合中尋找最優(yōu)組合,以實(shí)施CMLM方法。這稱為增強(qiáng)CMLM(enriched CMLM,ECMLM),它能提高檢測(cè)功效10%左右。

Kang等[24]同期也提出用P3D算法提高EMMA的計(jì)算速度,這稱為EMMA eXpedited(EMMAX)算法。隨后,Zhou等[25]提出了全基因組有效混合模型關(guān)聯(lián)(genome-wide efficient mixed-model association,GEMMA)精確算法。GEMMA算法有3個(gè)特點(diǎn):1)對(duì)親緣關(guān)系矩陣只需1次譜分解。2)似然函數(shù)、待估參數(shù)的估計(jì)量及其涉及的一階、二階導(dǎo)數(shù)中的行列式、矩陣跡、向量-矩陣-向量乘積,都表示為標(biāo)量運(yùn)算形式。它顯著提高了EMMA的運(yùn)算速度,對(duì)群體容量較大的情形尤其突出。3)與EMMA相比,GEMMA不允許基因型數(shù)據(jù)有缺失,需要完整的或者imputed SNP。

Svishcheva等[26]整合了基于得分檢驗(yàn)快速關(guān)聯(lián)(fast association score test-based analysis,FASTA)[27]和基于混合模型與回歸的全基因組快速關(guān)聯(lián)(genome-wide rapid association using mixed model and regression,GRAMMAR)[28]2種方法的優(yōu)點(diǎn),發(fā)展了GRAMMAR-Gamma兩步法。第1步,先估計(jì)衡量樣本關(guān)系的群體參數(shù)以及GRAMMAR-Gamma因子,通過K矩陣變換得到新的表型值;第2步,提出一種含有GRAMMAR-Gamma因子的得分檢驗(yàn)統(tǒng)計(jì)量,來檢驗(yàn)每一個(gè)SNP標(biāo)記是否與性狀顯著關(guān)聯(lián),其運(yùn)算復(fù)雜度非常接近理論上的最小值。

Lippert等[29]提出了譜分解變換線性混合模型(factored spectrally transformed linear mixed model,FaST-LMM)。針對(duì)大規(guī)模數(shù)據(jù)集,它不是利用所有SNP構(gòu)建遺傳相似矩陣(genetic similarity matrix),而是均勻地隨機(jī)抽取部分SNP構(gòu)建真實(shí)親緣關(guān)系矩陣(realized relationship matrix,RRM),并且這部分SNP數(shù)量小于樣本個(gè)體數(shù)量,從而在保證沒有信息損失的前提下,計(jì)算仍然準(zhǔn)確,而運(yùn)行速度卻得到了顯著提高。同時(shí),它在檢測(cè)所有標(biāo)記過程中,只需1次譜分解;利用K矩陣的譜分解變換表型、標(biāo)記信息向量以及協(xié)變量,可使變換后的數(shù)據(jù)不相關(guān)。

Listgarten等[30]發(fā)展了FaST-LMM-Select方法,通過線性回歸獲得每個(gè)SNP的P值,當(dāng)遺傳控制因子(genomic control factor)達(dá)到最小值時(shí),確定P值作為閾值,選擇滿足此閾值的一小部分SNP及其附近(比如2 cM以內(nèi))的標(biāo)記,不參與親緣關(guān)系矩陣構(gòu)建,再利用FaST-LMM方法進(jìn)行檢測(cè)。這種方法能夠提高檢測(cè)功效,節(jié)約計(jì)算成本。

Wang等[31]提出了SUPER(settlement of MLM under progressively exclusive relationship)方法,包括如下基本步驟:1)先進(jìn)行單標(biāo)記的一般線性模型和混合線性模型分析。2)將基因組劃分為若干bin區(qū)間,每個(gè)bin區(qū)間選出最顯著的SNP作為代表。3)用限制性最大似然方法確定最優(yōu)的bin區(qū)間大小以及數(shù)目,選出潛在關(guān)聯(lián)QTN來控制背景效應(yīng)。4)檢測(cè)每個(gè)標(biāo)記時(shí),去除潛在關(guān)聯(lián)QTN中與待檢測(cè)標(biāo)記高連鎖不平衡(linkage disequilibrium,LD)的SNP,再用FaST-LMM方法檢測(cè)顯著關(guān)聯(lián)標(biāo)記。這樣通過剔除與待檢測(cè)標(biāo)記高LD的潛在關(guān)聯(lián)QTN,避免親緣關(guān)系矩陣對(duì)待檢測(cè)標(biāo)記的過度矯正,提高檢測(cè)功效。

Loh等[32]假設(shè)多基因效應(yīng)服從高斯混合分布,利用快速變分近似計(jì)算表型殘差,憑借表型預(yù)測(cè)貝葉斯模型與經(jīng)典關(guān)聯(lián)檢驗(yàn)方法有機(jī)結(jié)合的追溯(retrospective)得分統(tǒng)計(jì)量檢驗(yàn)殘差與檢測(cè)標(biāo)記間的相關(guān)性,這稱為BOLT-LMM方法。它具有快速計(jì)算的特點(diǎn),表現(xiàn)在:計(jì)算類似相關(guān)個(gè)體混合模型關(guān)聯(lián)得分檢驗(yàn)(mixed-model association score test on related individuals,MASTOR)[33]的擬似然得分檢驗(yàn)統(tǒng)計(jì)量;利用類似于GRAMMAR-Gamma方法進(jìn)行校正;利用預(yù)處理共軛梯度來估計(jì)方差分量,代替矩陣特征值分解。它也具有高功效特點(diǎn),表現(xiàn)在:利用快速變分近似計(jì)算后驗(yàn)概率;利用LD得分回歸校正檢驗(yàn)統(tǒng)計(jì)量。

Jiang等[34]提出了基于混合線性模型全基因組關(guān)聯(lián)(MLM-based GWA,fastGWA),它是基于MLM框架的GWAS工具,利用SNP衍生的主成分來控制群體結(jié)構(gòu),利用系譜信息或稀疏遺傳關(guān)系矩陣來控制親緣關(guān)系,對(duì)biobank-scale數(shù)據(jù)進(jìn)行GWAS。它利用稀疏矩陣的Cholesky分解,提出了基于格子搜索求解限制性最大似然函數(shù)(fastGWA-REML)算法,代替了矩陣行列式和逆運(yùn)算,采用近似GRAMMAR-Gamma得分檢驗(yàn)統(tǒng)計(jì)量。模擬研究表明:fastGWA具有可靠、穩(wěn)健和高效的特點(diǎn)。

上述的所有快速方法中,都是1次僅對(duì)1個(gè)SNP進(jìn)行檢驗(yàn),屬于單位點(diǎn)分析,且待檢測(cè)SNP效應(yīng)都是固定效應(yīng)。

1.3 多位點(diǎn)模型方法的發(fā)展

雖然單位點(diǎn)分析應(yīng)用廣泛,但是復(fù)雜性狀是由多個(gè)微效基因共同控制,所以單位點(diǎn)分析并不符合數(shù)量性狀的真實(shí)模型。因此,多位點(diǎn)模型方法得到了快速發(fā)展。標(biāo)準(zhǔn)的多位點(diǎn)遺傳模型為:

(2)

式中:y、W、α和ε同模型(1);p表示關(guān)聯(lián)標(biāo)記個(gè)數(shù);Xi和βi分別表示第i個(gè)n×1關(guān)聯(lián)標(biāo)記基因型向量和效應(yīng)。

多位點(diǎn)分析一般以懲罰壓縮方法或貝葉斯方法或這兩者相結(jié)合的方式為主,例如彈性網(wǎng)[35]、經(jīng)驗(yàn)貝葉斯[36]和經(jīng)驗(yàn)貝葉斯LASSO(empirical Bayesian least absolute shrinkage and selection operator,EBLASSO)[37]。

Cho等[35]提出的借助彈性網(wǎng)變量選擇懲罰壓縮方法構(gòu)建多位點(diǎn)模型,是一種多步策略。第1步,對(duì)每一個(gè)SNP進(jìn)行線性回歸,得到部分顯著的標(biāo)記;第2步,對(duì)這些標(biāo)記利用彈性網(wǎng)懲罰回歸,進(jìn)一步得到顯著相關(guān)標(biāo)記;第3步,利用bootstrap抽樣,對(duì)上述顯著的SNP進(jìn)行顯著性檢驗(yàn)。

Lü等[36]提出了主效QTN、上位性互作和環(huán)境互作檢測(cè)的經(jīng)驗(yàn)貝葉斯方法。它將上述所有效應(yīng)放入同一個(gè)模型中,利用Xu[38]提出的經(jīng)驗(yàn)貝葉斯方法估計(jì)所有效應(yīng),并對(duì)非零效應(yīng)進(jìn)行似然比檢驗(yàn),達(dá)到同時(shí)檢測(cè)主效QTN、上位性互作和環(huán)境互作的目的。若標(biāo)記數(shù)目大于樣本容量50倍,經(jīng)驗(yàn)貝葉斯方法不能有效估計(jì)模型中包含的效應(yīng)。為克服這一問題,Wen等[39]將具有3層先驗(yàn)分布、采用內(nèi)積和外積迭代方式的EBLASSO算法[37]推廣到部分NCII交配設(shè)計(jì)上位性關(guān)聯(lián)作圖中,以剖析雜種優(yōu)勢(shì)的遺傳基礎(chǔ),構(gòu)建的遺傳模型同時(shí)包含加性效應(yīng)、顯性效應(yīng)、加性×加性效應(yīng)、加性×顯性效應(yīng)、顯性×加性效應(yīng)和顯性×顯性效應(yīng)。

如果標(biāo)記數(shù)量是樣本容量的幾百甚至數(shù)千倍時(shí),上述幾種方法將會(huì)失效。因此,應(yīng)該考慮如何在多位點(diǎn)模型中有效減少待估效應(yīng)數(shù)量。

Zhou等[40]將LMM和稀疏回歸模型相結(jié)合,在貝葉斯框架下構(gòu)建了多位點(diǎn)混合模型。假定每個(gè)SNP效應(yīng)是服從2個(gè)正態(tài)分布的混合分布的隨機(jī)效應(yīng),以滿足稀疏性,即大部分SNP效應(yīng)為零;同時(shí)利用親緣關(guān)系矩陣控制多基因背景。利用馬爾科夫鏈蒙特卡洛(Markov chain Monte Carlo,MCMC)方法,從觀測(cè)數(shù)據(jù)中獲得參數(shù)后驗(yàn)分布的近似樣本,用來估計(jì)超參數(shù)。這稱為BSLMM(Bayesian sparse linear mixed model)方法。由于稀疏性,其運(yùn)行速度較快,能處理至少1萬個(gè)體(每個(gè)體50萬SNP)的數(shù)據(jù)集。Moser等[41]考慮了貝葉斯混合分布模型(Bayesian mixture model,BayesR),將SNP效應(yīng)看成是隨機(jī)效應(yīng),服從 4個(gè)不同數(shù)量級(jí)方差正態(tài)分布的混合分布,需要估計(jì)的方差分量數(shù)不多。經(jīng)過幾輪Gibbs抽樣后,大約有500個(gè)效應(yīng)最終保留在模型當(dāng)中。計(jì)算時(shí)間過長(zhǎng)是多位點(diǎn)貝葉斯方法的主要缺點(diǎn)。

Tamba等[42]提出多位點(diǎn)模型兩步ISIS EM-BLASSO(iterative modified-sure independence screening EM-Bayesian LASSO)算法。第1步,利用SIS-SCAD[43]算法篩選潛在關(guān)聯(lián)標(biāo)記。其做法是:檢測(cè)每個(gè)標(biāo)記與數(shù)量性狀表型值相關(guān)的顯著性,選擇概率P≤0.01的SNP進(jìn)行SCAD[44]壓縮估計(jì),以篩選潛在關(guān)聯(lián)的p1個(gè)SNP。在剩余的標(biāo)記中,再用SIS-SCAD算法篩選出另一部分潛在關(guān)聯(lián)的p2個(gè)SNP;第2步,將第1步篩選的p1+p2個(gè)潛在關(guān)聯(lián)SNP標(biāo)記放入同一模型,用EM-BLASSO算法[38]和似然比檢驗(yàn),獲得與數(shù)量性狀顯著關(guān)聯(lián)的SNP及其效應(yīng)與貢獻(xiàn)大小。此算法的最大優(yōu)點(diǎn)是計(jì)算速度快和檢測(cè)功效高。

Wang等[45]提出了稀疏貝葉斯學(xué)習(xí)(sparse Bayesian learning,SBL)方法。該方法采用坐標(biāo)下降算法,在所有其他參數(shù)當(dāng)前值的基礎(chǔ)上,一次更新1個(gè)參數(shù)來估計(jì)標(biāo)記效應(yīng)。相比LASSO算法,它利用L2類型懲罰函數(shù),可以處理超過100 000的樣本量數(shù)據(jù)。同時(shí),Wang等[46]還構(gòu)建了壓縮嶺回歸(deshrinking ridge regression,DDR)檢驗(yàn)統(tǒng)計(jì)量,進(jìn)一步提高多位點(diǎn)模型QTN的檢測(cè)功效。

1.4 多位點(diǎn)GWAS混合模型方法的發(fā)展

植物數(shù)量性狀一般情況是受少數(shù)主基因和大量微效基因控制,因此考慮品種間親緣關(guān)系矩陣的多位點(diǎn)混合模型方法更符合其生物學(xué)特點(diǎn)。

Wang 等[47]將壓縮算法自適應(yīng)LASSO(adaptive LASSO)與線性混合模型LMM相結(jié)合,利用LASSO得到初始顯著的q個(gè)SNP,對(duì)每個(gè)k(0—q),考慮多基因遺傳背景,構(gòu)建LMM,利用最大似然估計(jì)方差分量,獲得協(xié)方差矩陣并進(jìn)行特征值分解,得到矩陣變換,并變換原始模型,使變換后的表型彼此不相關(guān),同時(shí)計(jì)算貝葉斯信息準(zhǔn)則(Bayesian information criteria,BIC)、拓展BIC(extended BIC,EBIC)和赤池信息準(zhǔn)則(Akaike’s information criteria,AIC)。對(duì)k重復(fù)以上步驟,以BIC、EBIC或AIC最小為準(zhǔn)則,找到最優(yōu)子集模型。

Segura等[48]提出了一種多位點(diǎn)混合模型關(guān)聯(lián)分析方法(multi-locus mixed-model,MLMM),它利用所有遺傳標(biāo)記構(gòu)建了親緣關(guān)系矩陣,將向前、向后的逐步回歸思想應(yīng)用到線性混合模型中,在回歸的每一步之前先對(duì)多基因背景方差和誤差方差進(jìn)行估計(jì),然后利用廣義最小二乘來估計(jì)SNP效應(yīng),進(jìn)行F檢驗(yàn)并獲得P值。隨后,將最顯著的SNP作為協(xié)變量放入模型中進(jìn)行下一步逐步回歸,所有協(xié)變量的P值與方差分量一起重新估計(jì)。重復(fù)這一過程,直到符合終止條件。它利用Gram-Schmidt正交化過程,以及對(duì)協(xié)變量矩陣進(jìn)行QR分解來提高計(jì)算效率。

Liu等[49]提出了固定與隨機(jī)交替概率統(tǒng)一模型(fixed and random model circulating probability unification,FarmCPU),同時(shí)結(jié)合MLMM[48]和FaST-LMM-Select[30]的優(yōu)點(diǎn),交替使用固定效應(yīng)模型和隨機(jī)效應(yīng)模型,盡量避免群體結(jié)構(gòu)、親緣關(guān)系和候選QTN間的混雜,達(dá)到降低假陽(yáng)性率的目的。在固定效應(yīng)模型中,一次僅檢驗(yàn)1個(gè)標(biāo)記,可能關(guān)聯(lián)SNP作為協(xié)變量來控制假陽(yáng)性率;從顯著關(guān)聯(lián)SNP中選取可能關(guān)聯(lián)SNP來計(jì)算親緣關(guān)系矩陣,并用隨機(jī)效應(yīng)模型進(jìn)行預(yù)測(cè)優(yōu)化。在隨機(jī)效應(yīng)模型中,計(jì)算該親緣關(guān)系矩陣所解釋的表型變異是否達(dá)到極大似然值,以此來防止固定效應(yīng)模型的過擬合問題。固定和隨機(jī)效應(yīng)模型一直交替使用直到?jīng)]有新的可能關(guān)聯(lián)SNP進(jìn)入模型中。

上述3種方法以及1.2節(jié)的所有快速混合線性模型方法都將待檢測(cè)SNP效應(yīng)視為固定效應(yīng)。Goddard等[50]認(rèn)為將SNP視為隨機(jī)效應(yīng),構(gòu)建隨機(jī)SNP效應(yīng)模型,其效果優(yōu)于固定SNP效應(yīng)模型。例如,隨機(jī)模型可以把SNP效應(yīng)壓縮為零[3,51]。然而,Goddard等[50]并沒有給出估計(jì)隨機(jī)SNP效應(yīng)的快速算法。

Wang等[3]提出了基于兩步法的mrMLM算法。首先,對(duì)每一標(biāo)記進(jìn)行全基因組單標(biāo)記掃描,以一種較為寬松的顯著標(biāo)準(zhǔn)選擇潛在關(guān)聯(lián)SNP;其次,將所選擇的潛在關(guān)聯(lián)標(biāo)記放入多位點(diǎn)模型,通過經(jīng)驗(yàn)貝葉斯[38]估計(jì)和似然比檢驗(yàn)獲得顯著關(guān)聯(lián)的QTN。它結(jié)合了P3D、FaST-LMM的模型變換,檢測(cè)標(biāo)記的隨機(jī)效應(yīng),并用固定效應(yīng)和誤差方差聯(lián)合估計(jì)這些技術(shù)。計(jì)算機(jī)模擬和真實(shí)數(shù)據(jù)分析表明,與單位點(diǎn)模型及固定SNP效應(yīng)模型方法相比,mrMLM的檢測(cè)功效和參數(shù)估計(jì)值精度更高,假陽(yáng)性率和假陰性率得到了有效控制。Zhang等[52]利用GEMMA[25]算法思想、矩陣變換以及Miller矩陣等式[53]來化簡(jiǎn)矩陣行列式和逆運(yùn)算,提高了mrMLM[3]方法的運(yùn)行速度,這稱為FASTmrMLM算法。

Wen等[4]提出了FASTmrEMMA算法。它是一種快速高效多位點(diǎn)兩步GWAS方法。第1步,將SNP效應(yīng)看成是隨機(jī)效應(yīng),且采用3種加速技術(shù)來減少運(yùn)行時(shí)間。這3種加速技術(shù)分別是:1)對(duì)原始MLM進(jìn)行特殊的矩陣變換,將多基因背景和誤差變異變換成標(biāo)準(zhǔn)正態(tài)離差;2)固定多基因背景方差與誤差方差的比值,對(duì)全基因組每一個(gè)SNP進(jìn)行單標(biāo)記掃描;3)檢測(cè)標(biāo)記信息矩陣的非零特征值數(shù)降至1個(gè)。第2步,對(duì)全基因組掃描概率P≤0.005的SNP放入多位點(diǎn)模型,利用經(jīng)驗(yàn)貝葉斯算法[38]進(jìn)行參數(shù)估計(jì),利用似然比檢驗(yàn)選出與目標(biāo)性狀顯著關(guān)聯(lián)的QTN。Wen等[5]進(jìn)一步優(yōu)化了FASTmrEMMA算法,利用Woodbury矩陣恒等式和性質(zhì),替換了特征向量的計(jì)算,適用于關(guān)聯(lián)分析大群體,計(jì)算速度比原算法提高了至少60%,其他性能與原算法一致。該算法具有QTN檢測(cè)功效高、精度高、運(yùn)算速度快、假陽(yáng)性率和假陰性率低的特點(diǎn)。

目前,FASTmrEMMA算法的矩陣變換已有5個(gè)拓展。Zhang等[54]將其他染色體多基因背景與誤差轉(zhuǎn)換成標(biāo)準(zhǔn)正態(tài)離差,然后通過最小角回歸在目標(biāo)染色體上選擇與性狀潛在關(guān)聯(lián)的標(biāo)記,最后將所有潛在關(guān)聯(lián)標(biāo)記放入多位點(diǎn)模型以檢測(cè)顯著關(guān)聯(lián)QTN,這就是pLARmEB方法。與此類似的還有pKWmEB[55]、TSLRF[56]、FastRR[57]以及F2群體GCIM[7]算法。其中,pKWmEB算法結(jié)合了非參數(shù)檢驗(yàn)Kruskal-Wallis;TSLRF算法結(jié)合了隨機(jī)森林;FastRR算法結(jié)合了壓縮嶺回歸DRR[46]檢驗(yàn)統(tǒng)計(jì)量。

1.5 軟件平臺(tái)的重要發(fā)展

目前有很多GWAS軟件包可供利用(表1)。這里僅對(duì)幾種主流的MLM方法軟件平臺(tái)作簡(jiǎn)單介紹。

表1 全基因組關(guān)聯(lián)分析混合線性模型方法常用軟件平臺(tái)Table 1 The common software platform for genome-wide association study(GWAS)based on mixed linear model(MLM)

TASSEL(trait analysis by association,evolution and linkage)[58]是一款早期的關(guān)聯(lián)分析軟件(https://www.maizegenetics.net/tassel),能夠執(zhí)行廣義線性模型(general linear model,GLM)和MLM方法[18],有效控制群體結(jié)構(gòu)和多基因遺傳背景對(duì)關(guān)聯(lián)分析的影響。

GEMMA(genome-wide efficient mixed model association)[25]是一款實(shí)現(xiàn)標(biāo)準(zhǔn)MLM的關(guān)聯(lián)分析軟件(http://www.xzlab.org/software.html),適用于大規(guī)模數(shù)據(jù),計(jì)算效率高,可利用免費(fèi)開源數(shù)字庫(kù)。包括了 4種算法:1)單變量線性混合模型(GEMMA)[25],用于單性狀分析,解釋群體分層和樣本結(jié)構(gòu),估計(jì)遺傳率;2)多變量線性混合模型(mvLMM)[59],用于多性狀分析,估計(jì)復(fù)雜性狀之間遺傳相關(guān)性;3)貝葉斯稀疏線性混合模型(BSLMM)[40],用于貝葉斯框架下的多位點(diǎn)分析;4)利用個(gè)體水平數(shù)據(jù)或匯總統(tǒng)計(jì)量(summary statistics)對(duì)方差分量進(jìn)行HE、REML和MQS估計(jì)[60]。

GAPIT(genome association and prediction integrated tool)[61]自開發(fā)以來,目前已更新至v3.0版本(http://zzlab.net/GAPIT),可執(zhí)行多種MLM關(guān)聯(lián)分析方法,包括Q+K[18]、CMLM[22]、ECMLM[23]、SUPER[31]、MLMM[48]、FarmCPU[49]和BLINK[62]。同時(shí)還包括了2種預(yù)測(cè)方法:sBLUP和cBLUP[63]。

mrMLM v4.0[52]是章元明教授團(tuán)隊(duì)研發(fā)的R語(yǔ)言包(https://cran.r-project.org/web/packages/mrMLM/index.html)。它能夠執(zhí)行mrMLM[3]、FASTmrMLM[52]、FASTmrEMMA[4-5]、ISIS EM-BLASSO[42]、pLARmEB[54]和pKWmEB[55]6種多位點(diǎn)GWAS方法。

2 GWAS方法學(xué)研究的影響因素

2.1 假陽(yáng)性率

GWAS面臨著較高假陽(yáng)性率(inflated false positive rates)的問題,這是因?yàn)槿后w結(jié)構(gòu)(population structure)[64-65]、品種間的家系結(jié)構(gòu)(family structure)和未知關(guān)聯(lián)(cryptic relatedness)不準(zhǔn)確[20],從而給GWAS帶來了新的挑戰(zhàn)。

目前,基于群體結(jié)構(gòu)控制假陽(yáng)性率的主要方法有基因組控制(genomic control,GC)[66]、結(jié)構(gòu)關(guān)聯(lián)(structured association,SA)[67]、主成分分析(principal components analysis,PCA)[68]和多維尺度(multidimensional scaling,MDS)[69]。這些方法能夠較好控制群體分層,但是并沒有考慮到所有品種的完整系譜,不適用于復(fù)雜的結(jié)構(gòu),例如品種間如果表現(xiàn)出復(fù)雜的親緣關(guān)系,就不適合采用上述方法[48]。混合線性模型(MLM)方法能有效控制由大量微小效應(yīng)多基因遺傳背景和群體結(jié)構(gòu)導(dǎo)致的偏差以及QTN檢測(cè)的假陽(yáng)性率,因此廣泛應(yīng)用于GWAS。

最近,Klasen等[70]認(rèn)為群體結(jié)構(gòu)校正過于嚴(yán)格,會(huì)影響真實(shí)關(guān)聯(lián)的鑒定,從而提出了數(shù)量性狀聚類關(guān)聯(lián)檢驗(yàn)(quantitative trait cluster association test,QTCAT)。該檢驗(yàn)在考慮標(biāo)記間相關(guān)性的同時(shí),進(jìn)行多標(biāo)記關(guān)聯(lián)。因此,QTCAT不需要進(jìn)行群體結(jié)構(gòu)校正,卻比單標(biāo)記方法更能反映出復(fù)雜性狀的多基因特性。模擬數(shù)據(jù)分析表明QTCAT明顯優(yōu)于MLM方法[70]。

2.2 模型選擇

2.2.1 單位點(diǎn)遺傳模型的問題在GWAS中,單標(biāo)記快速檢測(cè)算法運(yùn)算速度快,但它降低了檢測(cè)功效。這是因?yàn)槿绻麡?biāo)記數(shù)量很多,利用單標(biāo)記掃描算法,一次僅檢驗(yàn)1個(gè)標(biāo)記與性狀間的關(guān)聯(lián),則需要進(jìn)行Bonferroni校正,而這種校正過于嚴(yán)格,會(huì)導(dǎo)致一些小效應(yīng)QTN達(dá)不到顯著標(biāo)準(zhǔn),即檢測(cè)不到[3-4]。在基于混合模型的GWAS中,考慮了多基因背景控制,為了提高運(yùn)行速度,每次對(duì)QTN進(jìn)行單標(biāo)記掃描,采用了固定多基因背景方差與誤差方差的比值這種P3D算法,這是方差分量的近似算法。這無疑會(huì)使檢測(cè)功效降低,不同群體降低程度有所不同。此外,當(dāng)樣本較大時(shí),單標(biāo)記掃描算法需要大規(guī)模求解矩陣逆運(yùn)算,時(shí)間復(fù)雜度較高;并且對(duì)性狀方差總的遺傳貢獻(xiàn)率估計(jì)困難。

2.2.2 多位點(diǎn)遺傳模型的提出多位點(diǎn)模型比單位點(diǎn)模型更能解釋植物復(fù)雜性狀的遺傳基礎(chǔ),它可以考慮相鄰位點(diǎn)之間的潛在關(guān)系,具有更強(qiáng)的判別能力。在數(shù)學(xué)和統(tǒng)計(jì)上它是一種變量子集選擇的過程,但它的難點(diǎn)是,可能的子集個(gè)數(shù)隨SNP數(shù)的增加呈指數(shù)增長(zhǎng)。若將海量標(biāo)記全部放入同一多位點(diǎn)模型,則模型中變量(SNP數(shù))遠(yuǎn)遠(yuǎn)大于樣本容量,常用的線性回歸分析無法解決這類問題[43]。雖然當(dāng)前的壓縮估計(jì)方法和Bayesian方法可以部分解決這類問題,但是以當(dāng)前的CPU等硬件條件,運(yùn)算量很大,運(yùn)算時(shí)間很長(zhǎng)。目前,懲罰壓縮估計(jì)和Bayesian估計(jì)是多位點(diǎn)遺傳模型的主流算法,有時(shí)也將這2種估計(jì)結(jié)合使用。但是,這2種方法很少有考慮到多基因背景或群體結(jié)構(gòu)。大多數(shù)植物數(shù)量性狀是由少數(shù)大效應(yīng)和多數(shù)小效應(yīng)基因所控制的,那么需要考慮多基因背景控制和群體結(jié)構(gòu)的影響。對(duì)于微效QTN,采用多位點(diǎn)模型進(jìn)行懲罰壓縮估計(jì),容易出現(xiàn)這些小效應(yīng)都被壓縮至零的現(xiàn)象[3]。不同模型的選擇,會(huì)影響GWAS的結(jié)果。

2.3 計(jì)算與運(yùn)行速度

數(shù)據(jù)量龐大帶來了計(jì)算負(fù)擔(dān)。除開發(fā)GPU等硬件設(shè)施之外,從算法本身進(jìn)行創(chuàng)新與提速才是最重要的。充分利用降維方法,結(jié)合大規(guī)模矩陣分解與變換,是提高速度的有效途徑之一。目前快速M(fèi)LM方法的加速技術(shù)具體體現(xiàn)在:

1)單標(biāo)記掃描事先固定多基因背景方差與誤差方差的比率。CMLM[22]利用P3D技術(shù)加速M(fèi)LM方法。EMMAX[24]利用P3D技術(shù)加速EMMA[19]方法。mrMLM[3]、FASTmrMLM[52]、FASTmrEMMA[4-5]都利用P3D技術(shù)加速。P3D技術(shù)對(duì)估計(jì)多基因背景方差和誤差方差來說屬于近似算法,它可以避免單標(biāo)記掃描過程中重復(fù)估計(jì)這2個(gè)方差分量,從而提高計(jì)算速度。這是因?yàn)榇蟛糠諷NP對(duì)性狀不顯著,因此在檢測(cè)每一個(gè)標(biāo)記時(shí),這個(gè)比率近似相等。

2)化簡(jiǎn)大規(guī)模矩陣行列式和逆運(yùn)算。MLM方法的似然函數(shù)或限制性似然函數(shù)對(duì)參數(shù)進(jìn)行估計(jì)時(shí)都會(huì)涉及大規(guī)模的矩陣行列式和逆運(yùn)算。比如,GEMMA[25]利用一次親緣關(guān)系矩陣的譜分解,將行列式、矩陣跡以及向量-矩陣-向量乘積,都用標(biāo)量乘積的形式表示出來。FASTmrMLM[52]采用GEMMA算法思想以及Miller矩陣等式[53],將逆運(yùn)算化簡(jiǎn)成向量外積。FASTmrEMMA[4-5]利用矩陣特殊變換后,將矩陣行列式化簡(jiǎn)成向量?jī)?nèi)積運(yùn)算以及一個(gè)向量-矩陣-向量乘法運(yùn)算,二次型化簡(jiǎn)為特征向量的加權(quán)內(nèi)積運(yùn)算,單個(gè)標(biāo)記信息矩陣的非零特征值只是一個(gè)正數(shù)(轉(zhuǎn)換為向量?jī)?nèi)積運(yùn)算)。利用Woodbury矩陣等式,將二次型進(jìn)一步化簡(jiǎn)為與特征向量無關(guān)的向量?jī)?nèi)積運(yùn)算。fastGWA[34]利用稀疏矩陣的Cholesky分解,提出基于格子搜索求解限制性最大似然函數(shù)(fastGWA-REML)算法,代替矩陣行列式和逆運(yùn)算。充分利用數(shù)理統(tǒng)計(jì)知識(shí)、矩陣不同的分解形式等進(jìn)行化簡(jiǎn),化簡(jiǎn)成標(biāo)量乘積、向量?jī)?nèi)積、向量-矩陣-向量乘法等運(yùn)算,來提高運(yùn)算速度。

3)采取有效的降維策略。在大數(shù)據(jù)時(shí)代,利用單一方法,比如單位點(diǎn)分析或多位點(diǎn)分析,有可能導(dǎo)致計(jì)算不可行或檢測(cè)不準(zhǔn)確。而采取多步降維方法是一種有效策略。比如,在進(jìn)行多位點(diǎn)分析之前,可以將標(biāo)記數(shù)量降到接近樣本數(shù)量。由于LD使標(biāo)記之間存在共線性,是高維數(shù)據(jù)的主要問題,可以利用LD進(jìn)行SNP修剪[71],刪除一些高度相關(guān)的標(biāo)記,選擇部分標(biāo)記進(jìn)行多位點(diǎn)分析?；趩伪缎偷腟NP tagging[72]與SNP binning[51]也是高維數(shù)據(jù)降維的有效方式。此外,信念獨(dú)立篩選(sure independence screening,SIS)[43]以及迭代信念獨(dú)立篩選(iterative SIS)[43]算法可以處理超高維數(shù)據(jù),它們主要利用標(biāo)記與性狀之間的邊緣線性相關(guān)性來達(dá)到快速降維的目的。

3 小結(jié)與展望

近年來,多位點(diǎn)關(guān)聯(lián)分析在剖析動(dòng)植物的重要性狀和人的復(fù)雜疾病中引起了學(xué)者們廣泛關(guān)注,特別是以混合線性模型為框架構(gòu)建的多位點(diǎn)GWAS方法,既能很好控制群體結(jié)構(gòu)與多基因背景,降低假陽(yáng)性率,又能提高小效應(yīng)的檢測(cè)功效。深入研究多位點(diǎn)GWAS方法不僅有助于深入了解復(fù)雜性狀的遺傳結(jié)構(gòu),還可以為分子設(shè)計(jì)育種和分子生物學(xué)研究提供新思路。

稀有等位基因(minor allele frequency,MAF<1%)檢測(cè)一直是GWAS研究方法學(xué)的挑戰(zhàn)[73]。在數(shù)據(jù)進(jìn)行關(guān)聯(lián)分析之前,通常是刪除MAF<5%的位點(diǎn),因?yàn)檫@些位點(diǎn)基本檢測(cè)不到,此外由于Bonferroni校正,會(huì)導(dǎo)致閾值增加,其他顯著位點(diǎn)統(tǒng)計(jì)功效降低[74]。但是沒有理由認(rèn)為稀有等位基因在生物學(xué)上不重要。事實(shí)上,由于純化選擇,許多有害等位基因都將以低頻率出現(xiàn)[75]。目前,已有新的統(tǒng)計(jì)模型及軟件用來分析人的稀有等位基因[76]。這為動(dòng)植物數(shù)量性狀稀有等位基因的挖掘提供了思路。

實(shí)際數(shù)據(jù)分析中,針對(duì)復(fù)雜性狀的遺傳結(jié)構(gòu)和群體,如何選擇最好的方法,并沒有一個(gè)統(tǒng)一標(biāo)準(zhǔn)。不同GWAS方法往往會(huì)得到相似但不一致的結(jié)果[1-2,70]。例如某些顯著QTN已經(jīng)由一種方法檢測(cè)到并通過功能驗(yàn)證具有生物學(xué)意義,但是其他方法卻檢測(cè)不到。這是因?yàn)槔昧司哂胁煌卣骱蛢?yōu)勢(shì)的統(tǒng)計(jì)模型,比如單位點(diǎn)模型或多位點(diǎn)模型。我們建議用幾種方法同時(shí)分析同一組數(shù)據(jù),多種方法同時(shí)檢測(cè)到的顯著QTN應(yīng)該是可靠的。同時(shí),利用基因注釋、表達(dá)、京都基因和基因組百科全書(Kyoto Encyclopedia of Genes and Genomes,KEGG)以及網(wǎng)絡(luò)分析等現(xiàn)代組學(xué)手段,是比較容易發(fā)現(xiàn)顯著QTN附近的可靠候選基因。

確定合適的顯著性閾值一直是GWAS研究的熱點(diǎn)[77]。在人的數(shù)據(jù)分析中,5×10-8是常用閾值。單位點(diǎn)分析中,Bonferroni校正是通用的方法,即0.05/p(p表示標(biāo)記數(shù))。但它假設(shè)了標(biāo)記之間互相獨(dú)立,并且對(duì)于作物遺傳數(shù)據(jù)來說過于嚴(yán)格,以至丟失許多重要的小效應(yīng)位點(diǎn)[1]。Xu[51]提出利用標(biāo)記有效數(shù)來代替標(biāo)記數(shù)p,修正了Bonferroni閾值。假發(fā)現(xiàn)率(false discovery rate,FDR)[78]也假設(shè)了獨(dú)立性,使遺傳區(qū)域的SNP由于LD具有相似的檢驗(yàn)統(tǒng)計(jì)量,不太適合GWAS。置換檢驗(yàn)被認(rèn)為是建立顯著性閾值的標(biāo)準(zhǔn)方法[79]。該方法在無效假設(shè)條件下,從檢驗(yàn)統(tǒng)計(jì)量的分布中直接抽樣。通過置換表型,同時(shí)保持基因型數(shù)據(jù)不變來計(jì)算每個(gè)標(biāo)記的檢驗(yàn)統(tǒng)計(jì)量。但該方法非常耗時(shí),計(jì)算量龐大,不適合大數(shù)據(jù)的GWAS。同時(shí),在MLM框架下,置換表型可能會(huì)破壞源于遺傳關(guān)系的協(xié)方差結(jié)構(gòu)[80]。多位點(diǎn)分析中,由于所有潛在關(guān)聯(lián)標(biāo)記的個(gè)數(shù)及效應(yīng)能夠在一個(gè)線性模型中同時(shí)確定并估計(jì)出來,所以無需進(jìn)行Bonferroni校正[1-2]?？梢栽O(shè)定LOD=3,對(duì)應(yīng)的p閾值為0.000 2,比通用的0.05嚴(yán)格,能很好控制假陽(yáng)性率,同時(shí)獲得較高的統(tǒng)計(jì)功效[3-5]。這為確定合適的顯著性閾值提供了方向。

隨著其他類型的組學(xué)數(shù)據(jù),如轉(zhuǎn)錄組學(xué)、蛋白質(zhì)組學(xué)、代謝組學(xué)以及表觀遺傳學(xué)數(shù)據(jù)的開發(fā)利用[81-82],GWAS方法可以拓展到組學(xué)關(guān)聯(lián)分析(omic-wide association study,OWAS)[83]。多組學(xué)變量在GWAS中被認(rèn)為是性狀的一種,可以在一定程度上彌補(bǔ)基因型和表型之間的未知關(guān)聯(lián),最終有利于選擇性育種。例如,組學(xué)變量(在不同層次)被映射到與農(nóng)藝性狀相同的基因組位置,可以為復(fù)雜的遺傳結(jié)構(gòu)和潛在的生物途徑提供多方面闡釋。GWAS也能進(jìn)行多性狀分析,它可以利用個(gè)體水平的表型與基因型數(shù)據(jù)或者整合單性狀GWAS分析的前期結(jié)果進(jìn)行元分析[84]。在多組學(xué)關(guān)聯(lián)分析研究中應(yīng)用GWAS方法同時(shí)分析多水平和多性狀,結(jié)合多組學(xué)的多位點(diǎn)關(guān)聯(lián)分析方法能夠?yàn)檫z傳學(xué)大數(shù)據(jù)分析提供新的途徑。