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

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

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

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

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

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

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

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

1.1 模型建立

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

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

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

則有：

Y=XB+ε。

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

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

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

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

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

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

2.1 模擬試驗(yàn)設(shè)計(jì)

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

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

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

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

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

則協(xié)方差∑為：

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

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

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

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

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

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

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

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

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

3 結(jié)果與分析

3.1 MRCE模型對(duì)QTL定位效果的模擬試驗(yàn)結(jié)果

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

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

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

圖1 v=(0.5,0.5)時(shí)不同相關(guān)系數(shù)對(duì)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)時(shí)相關(guān)系數(shù)對(duì)模擬效果的影響。從表2可以看出，相關(guān)系數(shù)絕對(duì)值越大，QTL1和QTL2估計(jì)值越接近效應(yīng)值，功效也越高?？梢奙RCE模型可用于QTL定位。

表2 v=(0.5,0.1)時(shí)不同相關(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 模擬試驗(yàn)3 分析v=(0.5,0.5)時(shí)遺傳率對(duì)模擬結(jié)果的影響,結(jié)果見表3。

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

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

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

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

表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個(gè)標(biāo)記MK6、MK15、MK16、MK18、MK31、MK32、MK52、MK54與ph6性狀有關(guān)，6個(gè)標(biāo)記MK15、MK16、MK18、MK31、MK52、MK54與ph8性狀有關(guān)。

由表5還可以看出，與qtlnetwork軟件定位的QTL結(jié)果對(duì)比，基于MRCE模型選出的標(biāo)記中，有6個(gè)標(biāo)記與真實(shí)結(jié)果一致，尤其是MK6這個(gè)標(biāo)記僅與ph6有關(guān)；此外，還多定位到了3個(gè)標(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é)果基本相符，進(jìn)一步說明MRCE模型應(yīng)用于QTL定位是可行的。

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

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

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