基于協(xié)方差估計的多因變量回歸模型的多性狀QTL定位研究

張 慧，葉景山，申佳瑜，劉慧銘，尹 寧，李立婷，溫永仙

(1 福建農(nóng)林大學(xué)a計算機與信息學(xué)院，b統(tǒng)計及應(yīng)用研究所，福建 福州，350002； 2 漳州農(nóng)業(yè)發(fā)展集團有限公司，福建 漳州，363000；3 廈門華廈學(xué)院，福建 廈門 361021)

數(shù)量性狀基因座(quantitative trait loci,QTL)與連續(xù)變化的數(shù)量性狀表型有密切關(guān)系，常用DNA分子標(biāo)記技術(shù)對數(shù)量性狀基因遺傳位置進行標(biāo)記，QTL定位研究是遺傳學(xué)領(lǐng)域的一個重點。

早期的QTL定位方法是利用分子標(biāo)記與QTL之間的連鎖關(guān)系，定位出QTL在染色體上的位置，并估算出相應(yīng)QTL效應(yīng)值。但初期的單個性狀QTL定位存在一些問題。隨后，研究者提出了多性狀聯(lián)合定位分析方法。Jiang和Zeng[1]提出了一種多性狀的復(fù)合區(qū)間定位方法(composite interval mapping,CIM)，利用所考慮性狀的相關(guān)結(jié)構(gòu)進行定位，可以提高QTL檢測的準(zhǔn)確性。還有研究結(jié)果表明，同時檢測多個性狀比單獨檢測1個性狀更有效[2-4]。

多性狀QTL定位的實質(zhì)是在多因變量回歸模型的基礎(chǔ)上進行變量選擇。近年來，很多學(xué)者嘗試對多性狀QTL定位進行研究。Jansen 和Stam[5]提出了參數(shù)多變量回歸模型，研究多個性狀與分子標(biāo)記之間的關(guān)系，并通過極大似然比檢驗來找出與性狀相關(guān)的QTL位點。但這種方法計算量較大，為此，Lange和Whittaker[6]提出廣義估計方程，該方程不需要對具體分布進行假設(shè)，大大縮短了計算時間。肖靜等[7]和Xiao等[8]提出了多性狀主基因聯(lián)合分離分析方法(multivariate segregation analysis,MSA)，通過對比單個性狀和多個性狀聯(lián)合分析的模擬結(jié)果發(fā)現(xiàn)，多個性狀聯(lián)合分析效果較好，統(tǒng)計功效和效應(yīng)估計值的準(zhǔn)確度也較高。Banerjee等[9]在多性狀分析中引入貝葉斯模型，并結(jié)合馬爾科夫鏈蒙特卡洛(markov chain monte carlo,MCMC)算法進行模擬，建立相關(guān)表型和不相關(guān)表型兩個模型。Xu等[10]利用貝葉斯模型分析多個性狀與分子標(biāo)記之間的關(guān)系，通過壓縮系數(shù)方式來估計所有標(biāo)記區(qū)間內(nèi)的遺傳效應(yīng)。

關(guān)于多性狀聯(lián)合基因關(guān)聯(lián)分析，O′Reilly等[11]用MultiPhen方法，以可解釋的方式同時快速模擬了多種表型，提高了功效[11]。Bolormaa等[12]和Zhu等[13]用meta分析方法對多個性狀進行了基因關(guān)聯(lián)分析。Cheng等[14]用混合先驗貝葉斯回歸方法進行多性狀回歸分析發(fā)現(xiàn)，其效果優(yōu)于單性狀基因關(guān)聯(lián)分析。Tong等[15]結(jié)合期望最大化算法(expectation maximization,EM)，提出多性狀特征多區(qū)間下估計參數(shù)的方法(multiple trait multiple-interval mapping-new,MTMIM-NEW)。Montesinos- López等[16]用基于奇異值分解(singular value decomposition,SVD)的四階段分析方法進行多性狀基因關(guān)聯(lián)分析發(fā)現(xiàn)，其在參數(shù)估計和預(yù)測精度方面與使用貝葉斯多性狀多環(huán)境模型(bayesian multiple-trait and multiple-environment model,BMTME)獲得的結(jié)果類似。Yang等[17]提出了一個具有多性狀的全關(guān)聯(lián)的整合函數(shù)線性模型，利用懲罰函數(shù)解決了單核苷酸多態(tài)性(single nucleotide polymorphism,SNP)的高維性和多性狀相關(guān)性問題。Lin等[18]提出一種基于混合線性模型的多性狀聯(lián)合基因關(guān)聯(lián)分析方法，模擬結(jié)果表明，多性狀全基因組關(guān)聯(lián)研究(genome-wide association studies,GWAS)在檢測多效性位點的影響方面較單個性狀效果更好。Tran等[19]在繪制多數(shù)量性狀位點的統(tǒng)計方法中考慮到X染色體，擴展了一種多QTL模型選擇的懲罰似然方法。

此外，還有一些降維方法被應(yīng)用于解決多性狀基因關(guān)聯(lián)分析問題，包括主成分分析[20-21]、典型相關(guān)分析[22]、偏最小二乘法[23]和貝葉斯Lasso方法[24]。本研究采用Rothman等[25]提出的基于協(xié)方差估計的多因變量回歸(multivariate regression with covariance estimation,MRCE)模型，通過計算機模擬產(chǎn)生基因型數(shù)據(jù)和性狀表型數(shù)據(jù)，利用MRCE模型進行參數(shù)估計，探究基因位點解釋的方差比、表型相關(guān)系數(shù)、遺傳率對模擬效果的影響，并將此模型應(yīng)用于水稻群體標(biāo)記數(shù)據(jù)中，完成基因定位，估計其參數(shù)，以期為多性狀QTL定位研究提供參考。

1 研究模型的構(gòu)建

1.1 模型建立

假設(shè)一個遺傳群體包含n個個體，若不考慮群體結(jié)構(gòu)等因素，對第i個個體，在遺傳關(guān)聯(lián)分析中假設(shè)有p遺傳標(biāo)記為xi1,xi2,…,xip(i=1,2,…,n)，若有q個數(shù)量性狀，線性遺傳模型可以表示為：

1,2,…，q;k=1,2,…，p)。

式中：yij表示第i個個體第j個性狀表型值。xik表示第i個個體在第k個基因標(biāo)記位點的指示變量值，若A和a表示1對等位基因，當(dāng)基因型是AA時，xij取1；當(dāng)基因型是Aa時，xik取0；當(dāng)基因型是aa時，xik取-1。β0j代表第j個數(shù)量性狀的均值，βkj代表第k個基因標(biāo)記位點對第j個數(shù)量性狀所表現(xiàn)的遺傳效應(yīng)值。εij為隨機誤差，一般εij之間不是相互獨立的，假定它們服從均值均為0，協(xié)方差矩陣為∑的多元正態(tài)分布。當(dāng)q=1時，模型為經(jīng)典的單因變量回歸模型。將線性遺傳模型寫成矩陣形式，分別用X、Y、B、ε表示：

則有：

Y=XB+ε。

Rothman[25]提出了MRCE，B的稀疏估計量，該方法在負對數(shù)似然函數(shù)上加入了兩個懲罰項，求解B的稀疏估計值，具體形式為：

式中：Ω=∑-1=[ωj′j]，∑-1是協(xié)方差矩陣∑的逆矩陣，ωj′j是逆矩陣中的元素。λ1≥0,λ2≥0,二者均是調(diào)整參數(shù)，用k折交叉驗證來選擇參數(shù)λ1和λ2。

1.2 模型的假設(shè)檢驗

首先，原假設(shè)效應(yīng)系數(shù)都為0，通過基于Pillai-Bartlett跡、Hotelling-Lawley跡和Wilks’s Lambda的近似F分布檢驗進行模型檢驗[26-30]。其次，用Lβ=0的方法對基因標(biāo)記位點的遺傳效應(yīng)βij(i=1,2,…,p;j=1,2,…,q)進行檢驗[31]，其中L是c×p+1階的矩陣，用來識別檢驗假設(shè)的遺傳效應(yīng)。如對β1j的假設(shè)檢驗可以寫：

對于假設(shè)可采用F檢驗進行，F(xiàn)檢驗的形式為：

2 研究模型的模擬與實際應(yīng)用

2.1 模擬試驗設(shè)計

2.1.1 SNPs生成 參照黃楊岳等[32]的SNPs數(shù)據(jù)仿真方法，生成純合SNP模擬數(shù)據(jù)SNPs，包含500個個體和200個基因位點，其基因型為AA、aa。

2.1.2 數(shù)量性狀表型值生成的多元仿真框架 (1)給定截距b的合適值。

(2)按照Porter和O’Reilly[33]的方法，給定v值，v是遺傳變異所解釋的表型方差遺傳效應(yīng)向量。例如，當(dāng)v=(0,1,0.5)時，對應(yīng)于SNP，表示解釋了性狀1的 0.1%表型方差，性狀2的0.5%表型方差。

性狀表型相關(guān)矩陣為R：

式中：ρij(i,j=1,2,…,q)表示第i個性狀與第j個性狀的表型相關(guān)系數(shù)。

則協(xié)方差∑為：

數(shù)量性狀表型值y的計算公式為：

y=b+f(v)x+ε。

式中：x代表基因型指示變量，基因型為AA時取1，基因型為aa時取-1。

2.1.3QTL檢驗功效的計算 對于染色體上的某個基因位點需要對其進行參數(shù)估計及統(tǒng)計檢驗，若對于給定的顯著性水平，該位點的遺傳效應(yīng)值達到顯著，說明在該位點檢測到QTL。若假設(shè)情況的計算機模擬共重復(fù)m次，染色體基因位點能檢測到m0次，則該位點的QTL檢測功效為m/m0。效應(yīng)值是f(v)，估計值是用MRCE方法估計的f(v)。

(1)模擬1。給定3組v值，v1=(0.5,0.5)，v2=(0.5,0.1)，v3=(0.5,0.0)，計算相應(yīng)的f(v)：f(v1)=(0.1253,0.1253)，f(v2)=(0.1253,0.0559)，f(v3)=(0.1253,0.0000)，進而得到相應(yīng)性狀表型值，利用MRCE模型進行參數(shù)估計，根據(jù)功效，比較v值對QTL定位模擬效果的影響。

(2)模擬2。設(shè)v=(0.5,0.5)，相關(guān)系數(shù)從-0.9每次增加0.1直到0.9，分別產(chǎn)生相應(yīng)的性狀表型值，利用MRCE方法進行參數(shù)估計，探究相關(guān)系數(shù)對模擬效果的影響；同時設(shè)v=(0.5,0.1)、遺傳率為(0.05,0.05)時，研究相關(guān)系數(shù)對QTL定位模擬效果的影響。

(3)模擬3。設(shè)v=(0.5,0.5)，給定不同遺傳率組合(0.05，0.05)，(0.05，0.10)，(0.05，0.15)，(0.10，0.05)，(0.10，0.10)，(0.10，0.15)，(0.15，0.05)，(0.15，0.10)，(0.15，0.15)，對相應(yīng)性狀表型值進行功效模擬，分析遺傳率對QTL定位模擬效果的影響。

2.2 實例數(shù)據(jù)收集

2.2.1實例1數(shù)據(jù) 實例1數(shù)據(jù)選自qtlnetwork軟件，是一個水稻DH群體，包含12條水稻染色體中的3條染色體，共54個標(biāo)記，每條染色體上標(biāo)記數(shù)量不等，99個個體，2個環(huán)境(1998年和1999年)。由于水稻DH群體數(shù)據(jù)中存在缺失數(shù)據(jù)，需通過相鄰平均值方法進行填補，再將1998年與1999年數(shù)據(jù)進行整合，最終得到的數(shù)據(jù)集為54個標(biāo)記和198個樣本量，提取ph6、ph8作為性狀表型值，兩性狀的相關(guān)系數(shù)為0.9464。

2.2.2實例2數(shù)據(jù) 實例2數(shù)據(jù)是包含12條染色體的水稻永久F2群體試驗數(shù)據(jù)[37-40]，該群體由來自珍汕97×明恢63，含有210個株系的重組自交系(RIL)群體隨機交配生成，共產(chǎn)生278個雜種株系，其遺傳圖譜共有1619個標(biāo)記序號(Bin1～Bin1619)，包含單株產(chǎn)量、分蘗數(shù)、穗粒數(shù)、粒質(zhì)量4個性狀，本研究僅對穗粒數(shù)和粒質(zhì)量進行聯(lián)合分析，剔除缺失數(shù)據(jù)后獲得2組完整數(shù)據(jù)，其中1998年有246個，1999年有276個，為了簡化考慮，本研究僅考慮其加性效應(yīng)。

3 結(jié)果與分析

3.1 MRCE模型對QTL定位效果的模擬試驗結(jié)果

3.1.1 模擬試驗1 通過MRCE模型對不同v值情況下的QTL進行定位發(fā)現(xiàn)，任意給定一個固定的相關(guān)系數(shù)和遺傳率時，不同v值對QTL定位的影響規(guī)律大致相同，所以本研究選擇其中3個相關(guān)系數(shù)(0.1，0.5，0.9)且遺傳率為(0.05，0.05)時進行分析，結(jié)果見表1。從表1可以看出，當(dāng)相關(guān)系數(shù)分別為0.1，0.5，0.9，遺傳率為(0.05，0.05)時，v值越大，功效越大；v為0時，功效為0或接近0。所以，如果遺傳變異所解釋的方差比大小合適，則利用MRCE模型進行QTL定位是可行的。

表1 不同v值情況下QTL定位的模擬結(jié)果Table 1 Simulation of QTL mapping with different v values

3.1.2 模擬試驗2 圖1表明，當(dāng)v相同時，兩端功效略高于中間部分，說明相關(guān)系數(shù)絕對值越大，其功效越高。

圖1 v=(0.5,0.5)時不同相關(guān)系數(shù)對MRCE模型用于QTL定位模擬效果的影響Fig.1 Simulation of QTL mapping based on MRCE for different correlation coefficients and v=(0.5,0.5)

表2是當(dāng)v=(0.5,0.1)、遺傳率為(0.05,0.05)時相關(guān)系數(shù)對模擬效果的影響。從表2可以看出，相關(guān)系數(shù)絕對值越大，QTL1和QTL2估計值越接近效應(yīng)值，功效也越高?？梢奙RCE模型可用于QTL定位。

表2 v=(0.5,0.1)時不同相關(guān)系數(shù)情況下QTL定位的模擬結(jié)果Table 2 Simulation of QTL mapping with different correlation coefficients and v=(0.5,0.1)

表2(續(xù)) ConutinuedTable 2

3.1.3 模擬試驗3 分析v=(0.5,0.5)時遺傳率對模擬結(jié)果的影響,結(jié)果見表3。

表3 v=(0.5,0.5)時不同遺傳率下QTL定位的模擬結(jié)果Table 3 Simulation of QTL mapping for different heritability and v=(0.5,0.5)

從表3可以看出，遺傳率越高，其效應(yīng)估計值越接近真值，功效也越好，在其他不同遺傳率假設(shè)下也有上述相似結(jié)果。綜上可知，利用MRCE模型進行QTL定位分析是可行的，同時遺傳變異所占方差比越大，相關(guān)系數(shù)絕對值越大，遺傳率越大，則模擬效果越好。

3.2 MRCE模型對QTL定位效果的應(yīng)用實例

3.2.1 應(yīng)用實例1 表4和表5分別為用qtlnetwork軟件和MRCE模型得出的水稻DH群體數(shù)據(jù)QTL定位結(jié)果，定位到的QTL均通過了顯著性檢驗。

表4 基于qtlnetwork的水稻DH群體數(shù)據(jù) QTL定位結(jié)果Table 4 QTL mapping for DH population of rice by qtlnetwork

表5 基于MRCE模型的水稻DH群體數(shù)據(jù)QTL定位結(jié)果Table 5 QTL mapping for DH population of rice based on MRCE

從表5可以看出，通過MRCE模型發(fā)現(xiàn)，8個標(biāo)記MK6、MK15、MK16、MK18、MK31、MK32、MK52、MK54與ph6性狀有關(guān)，6個標(biāo)記MK15、MK16、MK18、MK31、MK52、MK54與ph8性狀有關(guān)。

由表5還可以看出，與qtlnetwork軟件定位的QTL結(jié)果對比，基于MRCE模型選出的標(biāo)記中，有6個標(biāo)記與真實結(jié)果一致，尤其是MK6這個標(biāo)記僅與ph6有關(guān)；此外，還多定位到了3個標(biāo)記，分別為MK18、MK32、MK52；且這些標(biāo)記與qtlnetwork軟件定位的QTL結(jié)果相鄰。MK18與qtlnetwork軟件定位的MK15-MK16相鄰，MK32與MK30-MK31相鄰，MK52與MK53-MK54相鄰，多定位到的標(biāo)記可能與鄰近QTL效應(yīng)的影響以及在qtlnetwork軟件定位過程中的閾值設(shè)定有關(guān)，由此可知,基于MRCE模型的QTL定位與用qtlnetwork軟件的定位結(jié)果基本相符，進一步說明MRCE模型應(yīng)用于QTL定位是可行的。

3.2.2 應(yīng)用實例2 由于MRCE模型不能用于樣本量(n)小于遺傳標(biāo)記個數(shù)(p)的情況，所以本研究計算了遺傳標(biāo)記與性狀表型值的邊際相關(guān)系數(shù)，且邊際相關(guān)系數(shù)越高，則該遺傳標(biāo)記與對應(yīng)性狀表型值的相關(guān)性越高，最終選取邊際相關(guān)系數(shù)絕對值較大的200個標(biāo)記數(shù)據(jù)進行初步降維，QTL定位結(jié)果見表6。表6表明，利用MRCE模型檢測到1998年穗粒數(shù)在第3、第6和第7條染色體上各有1個QTL；粒質(zhì)量在第1和第5條染色體上各有3個QTL，第3和第7條染色體上各有2個QTL。利用MRCE模型檢測到1999年穗粒數(shù)在第3條染色體有1個QTL，第7條染色體有2個QTL，粒質(zhì)量第1和第3條染色體各有1個QTL，第5和第7條染色體各有2個QTL。對比1998和1999年的定位結(jié)果可知，穗粒數(shù)都定位到Bin436，粒質(zhì)量都定位到Bin65、Bin439、Bin699、Bin769和Bin1008，是因為1998年穗粒數(shù)與粒質(zhì)量2個性狀之間的相關(guān)系數(shù)(0.15)大于1999年(0.05)，故1998年定位出更多QTL。

綜合實例1和實例2的結(jié)果，說明MRCE模型不僅可以用于模擬QTL定位，而且在實際定位中也同樣適用，結(jié)果較好。

表6 基于MRCE模型的水稻永久F2群體穗粒數(shù)和粒質(zhì)量的QTL定位結(jié)果Table 6 QTL mapping of grains per panicle and grain weight for immortalized F2 population of rice by MRCE

4 討 論

采用不同的QTL定位方法和不同數(shù)據(jù)檢測到的QTL數(shù)目和位置可能有差異，若能定位到更多的QTL，在一定程度上可以彌補用其他方法未找到的備選QTL，但是否是真實的QTL，需用生物檢測方法進行驗證。Yu等[39]用超高密度SNP圖譜，檢測出穗粒數(shù)性狀在第1、第3和第7條染色體上各有1個QTL，粒質(zhì)量在第1和第3條染色體上各有2個QTL，第5和第9條染色體上各有1個QTL。對比本研究結(jié)果可以發(fā)現(xiàn)，利用MRCE模型檢測出的穗粒數(shù)、粒質(zhì)量QTL更多；其中1998年的數(shù)據(jù)中多檢測出穗粒數(shù)第6條染色體上的1個QTL，粒質(zhì)量多檢測出第7條染色體上的2個QTL，且與Yu等[39]檢測到的QTL位置大致相近；但本研究利用MRCE模型檢測時丟失了穗粒數(shù)第1條染色體上的1個QTL和粒質(zhì)量性狀第9條染色體上的1個QTL、第1條染色體Bin172位置附近的1個QTL。原因可能是只考慮了加性效應(yīng)而沒有考慮顯性效應(yīng)，或丟失QTL的LOD值都比較小，剛好超過給定的閾值[39]。

本研究僅驗證了MRCE模型定位QTL的可行性和優(yōu)勢，即表型性狀聯(lián)合定位時相關(guān)系數(shù)越大，效果越好。且MRCE模型不適用樣本量小于維度的情況，對此情況可利用降維手段，先將高維數(shù)據(jù)降為低維。對于水稻永久F2群體數(shù)據(jù)分析中定位到的QTL少的問題，可以增加顯性效應(yīng)進一步研究。