佟 良，周 影，趙衛(wèi)績(jī)，丁 慧，鄒大偉，李成鳳，馬春華

（1.綏化學(xué)院 信息工程學(xué)院，黑龍江 綏化152061；2.黑龍江大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，黑龍江 哈爾濱150080）

數(shù)量性狀位點(diǎn)（QTL）泛指影響數(shù)量性狀的基因。將包含有遺傳信息的染色體片段量化，使人們可以利用更多的統(tǒng)計(jì)手段來探求影響生物性狀的QTL的數(shù)目、位置以及效應(yīng)的大小，從而為遺傳育種、疾病防治等工作提供支持。

分子標(biāo)記被廣泛應(yīng)用到動(dòng)植物的遺傳研究中。Lander和 Botstein［1］，Zeng［2］，Kao［3］在這方面做出了突出貢獻(xiàn)。但是上面所有方法都依賴于高精度的基因型測(cè)量技術(shù)，這一點(diǎn)由于基因型記分軟件的某些缺陷以及生物化學(xué)的異?，F(xiàn)象，所得到的基因型數(shù)據(jù)往往都帶有一定測(cè)量誤差。實(shí)際研究發(fā)現(xiàn)，即使是很小的基因型測(cè)量誤差也會(huì)給遺傳疾病的研究帶來嚴(yán)重的負(fù)面影響，例如連鎖分析（linkage analysis）研究和遺傳距離的估計(jì)［4］以及連鎖不平衡參數(shù)的估計(jì)［5］等。Lebrec等人［6］討論了基因型誤差對(duì)復(fù)雜性狀連鎖定位的影響。Hou等人［7］對(duì)QTL位于標(biāo)記位點(diǎn)上的情況進(jìn)行了分析，Tong［8］等人就回交群體給出了分析，但對(duì)F2群體沒有做進(jìn)一步的研究。

本文基于F2群體在基因型數(shù)據(jù)含有誤差的情形下考慮了每個(gè)個(gè)體所有可能的基因型（基因圖譜）。根據(jù)所有可能的基因型，在模型框架下給出了模型參數(shù)的估計(jì)方法。文中假設(shè)所有標(biāo)記位點(diǎn)的基因型錯(cuò)誤率相同，用EM 算法［9］和加權(quán)的EM 算法［10］進(jìn)行參數(shù)模擬。研究表明，考慮F2群體基因型帶有誤差能夠減少基因型誤差給QTL定位帶來的影響。

1 理論與方法

1.1 記號(hào)和背景知識(shí)

考慮N個(gè)自交個(gè)體，緊密連鎖的M＋1個(gè)標(biāo)記位點(diǎn)構(gòu)成M個(gè)標(biāo)記區(qū)間。假設(shè)標(biāo)記位點(diǎn)均為兩個(gè)等位基因Mj和mj，令

其中：Xij和?Xij（i＝1，…，N，j＝1，…，M＋1）分別表示第i個(gè)個(gè)體，第j個(gè)標(biāo)記位點(diǎn)真實(shí)基因型和含有誤差的基因型。Yi表示第i個(gè)個(gè)體的表型值，X＊ij表示第i個(gè)個(gè)體，第j個(gè)標(biāo)記位點(diǎn)內(nèi)的潛在QTL基因型，令

令γj表示第j個(gè)標(biāo)記區(qū)間兩側(cè)標(biāo)記之間的重組率，γj1為第j個(gè)標(biāo)記與第j個(gè)標(biāo)記區(qū)間內(nèi)的QTL之間的重組率，其中γj已知，假設(shè)一個(gè)標(biāo)記區(qū)間內(nèi)至多存在一個(gè)QTL，令p（）為第i個(gè)個(gè)體第j個(gè)標(biāo)記區(qū)間基因型給定條件下的條件概率，它可由QTL位置以及QTL所在區(qū)間的雙側(cè)標(biāo)記基因型計(jì)算求得。對(duì)于F2群體來說，其條件概率列于表1。

表1 F2群體標(biāo)記基因型已知情況下QTL基因型的條件概率

假設(shè)每個(gè)標(biāo)記位點(diǎn)的等位基因Mj和mj以概率為θ等可能被誤編。這里θ＝p（＝k｜Xij≠k）表示基因型誤差率（k＝－1，0，1），φi是第i個(gè)個(gè)體的聯(lián)合誤差率。ki表示M＋1個(gè)標(biāo)記位點(diǎn)的錯(cuò)誤編碼個(gè)數(shù)，它是可以計(jì)算的。假定不同標(biāo)記位點(diǎn)標(biāo)記基因型是否有誤差是相互獨(dú)立的，可得到

式中：α為群體均值，βj為第j個(gè)標(biāo)記區(qū)間內(nèi)QTL的效應(yīng)值，εi為隨機(jī)誤差，這里X＊ij與εi是相互獨(dú)立的，且εi～N（0，σ2）。

1.2 基因型含有誤差時(shí)模型參數(shù)的EM方法

模型的參數(shù)向量 Ω＝（α，β1…βM，γ，θ，σ2），其中γ＝（γ11…γM1）。由于觀測(cè)到的基因型數(shù)據(jù)是含有誤差的，因此，與通常的區(qū)間定位不同的是標(biāo)記位點(diǎn)真實(shí)的基因型和QTL基因型均為潛在數(shù)據(jù)。這里采用EM算法進(jìn)行參數(shù)估計(jì)。針對(duì)可能帶有誤差的標(biāo)記基因型數(shù)據(jù)，借助前面的基因型圖譜考慮每個(gè)個(gè)體所有可能的基因型；針對(duì)QTL基因型將依據(jù)相應(yīng)的傳遞概率（見表1）對(duì)其補(bǔ)值。用標(biāo)記位點(diǎn)真實(shí)的基因型和QTL基因型補(bǔ)全觀測(cè)數(shù)據(jù)｛Yi），i＝1，…，N｝，進(jìn)而得到完全數(shù)據(jù)｛，Yi，Xi，），i＝1，…，N｝。因而第i個(gè)個(gè)體的完全似然函數(shù)為

對(duì)于E 步來說，考慮到?X，Y，Ω（k），計(jì)算lc（Ω）的條件概率，有

其中，Ω（k）為Ω當(dāng)前的估計(jì)值，?X＝（?X1…?XN）表示N個(gè)個(gè)體中含有誤差的標(biāo)記基因型，令

對(duì)于M 步來說，方程（1）的第1項(xiàng)只與α，β1，…，βM，σ2有關(guān)，第2項(xiàng)只與γ有關(guān)，最后1項(xiàng)只與θ有關(guān)，對(duì)于方程（1）的第1項(xiàng)，為了簡(jiǎn)單，考慮其矢量和矩陣形式。

令

E1可以通過下面公式推導(dǎo)求得

這里?P＝P（?Xc｜Xc，Ω（k）），，s，t＝1…M，k＝－1，1；l＝－1，1。

為了求得γ的迭代值，引入一個(gè)隨機(jī)變量

對(duì)式（2）進(jìn)行極大化，得到迭代公式

2 模擬與結(jié)果

模擬比較考慮基因型誤差情形和不考慮基因型誤差情形兩種參數(shù)估計(jì)方法的好壞。為了簡(jiǎn)單起見，只考慮兩個(gè)數(shù)量性狀位點(diǎn)。標(biāo)記密度對(duì)QTL效應(yīng)的影響不大，同時(shí)，隨著標(biāo)記間距的減小，位置估計(jì)的準(zhǔn)確性也有所提高，且趨向一致。所以只討論標(biāo)記間距為10cm的情況（標(biāo)記間距可轉(zhuǎn)化為重組率［12］），樣本容量選定為500、1 000，兩種方法的表型值＋εi，其中α是群體均值，εi～N（0，σ2）表示隨機(jī)誤差，，j＝1，2為QTL基因型的值，取值為－1，0，1。為了體現(xiàn)基因型誤差θ，

2.1 模擬設(shè)計(jì)

重組率γ和QTL效應(yīng)β1，β2之間的影響，在緊密連鎖情形下模擬了θ＝0.00、0.01、0.05、0.10和γ11＝0.02，γ21＝0.03，參數(shù)β1，β2選取使得遺傳力在0.05、0.1、0.2附近。對(duì)每個(gè)參數(shù)用本文提出的方法（記為PM）和不考慮標(biāo)記位點(diǎn)誤差的方法（記為QM）進(jìn)行估計(jì)，并將整個(gè)過程重復(fù)500次來計(jì)算參數(shù)估計(jì)的均值。為了評(píng)價(jià)參數(shù)估計(jì)的偏性，給出了每個(gè)參數(shù)的均方誤差（MSE）。

2.2 模擬結(jié)果

表2給出4種基因型誤差率（0，0.01，0.05，0.1）和遺傳力h2＝0.1情形下兩種方法的參數(shù)估計(jì)值和MSE。容易看出，θ＝0時(shí)，所估計(jì)的參數(shù)值是等價(jià)的，這是合理的，因?yàn)檫@時(shí)條件期望E1，E2等于它們本身。θ≠0時(shí)隨著基因型誤差率的增加，兩種方法的參數(shù)估計(jì)值逐漸偏離參數(shù)真值，但PM方法所估計(jì)的參數(shù)值更接近參數(shù)真值。這說明PM方法能夠降低基因型誤差率對(duì)參數(shù)估計(jì)的影響。同時(shí)看出，兩種方法的MSE隨著基因型誤差率的增加而增加，但是PM方法的MSE小于QM方法的MSE，這再次表明考慮基因型誤差能夠降低測(cè)量誤差帶來的影響。

表2 樣本量N＝500，h2＝0.1不同誤差率情形下兩種方法的模擬結(jié)果

PM方法能夠估計(jì)誤差率θ，表3給出遺傳力N＝500，h2＝0.01、0.05、0.1時(shí)，誤差率θ＝0的估計(jì)值，從表中可以看出，在遺傳力和效應(yīng)不同的條件下，誤差率θ＝0.01的估計(jì)值隨著θ＝0.05真值的增加，偏離程度逐漸增大，MSE隨著基因型誤差率的增加而增加。

表3 遺傳力和QTL效應(yīng)不同時(shí)誤差率的估計(jì)

3 檢驗(yàn)統(tǒng)計(jì)量

事實(shí)上，獲得參數(shù)向量Ω的估計(jì)之后，可進(jìn)一步討論區(qū)間中是否顯著存在QTL，用似然比統(tǒng)計(jì)量進(jìn)行研究。原假設(shè)為區(qū)間內(nèi)沒有QTL，即

假設(shè)區(qū)間內(nèi)至少有一個(gè)QTL，lc（Ω）表示對(duì)數(shù)似然函數(shù)，則相應(yīng)的似然比統(tǒng)計(jì)量為

其中，＾Ω是在H1下Ω的 MLE，＾Ω0是在H0下Ω的MLE。通過簡(jiǎn)單的推導(dǎo)很容易求得Ω在H0和H1下的MLE。

4 結(jié)論與討論

本文提出了基于F2群體基因型帶有誤差情形下區(qū)間定位的參數(shù)估計(jì)方法，模擬結(jié)果表明，新方法能較好地定位QTL位置，比不考慮基因型誤差的方法更優(yōu)。在相同條件下做樣本容量為500，1 000的參數(shù)估計(jì)時(shí)，發(fā)現(xiàn)用于QTL定位的群體規(guī)模對(duì)QTL定位效果有明顯的影響，隨著群體規(guī)模的擴(kuò)大，所有參數(shù)估計(jì)的準(zhǔn)確度均隨之提高。因此，在實(shí)際QTL定位中，讓群體達(dá)到一定的規(guī)模，有助于提高QTL定位的可靠性，這與Jeon等人［13］的結(jié)論是一致的。對(duì)不同遺傳力水平下QTL定位效果的分析表明，目標(biāo)數(shù)量性狀的遺傳力高低對(duì)QTL定位的準(zhǔn)確性也有明顯的影響。當(dāng)對(duì)高遺傳力性狀實(shí)施QTL定位時(shí)，QTL位置估計(jì)的準(zhǔn)確性也相對(duì)較高。本文提出的方法可應(yīng)用于誤差率不同時(shí)QTL的區(qū)間定位。當(dāng)誤差率不同時(shí)，只需適當(dāng)?shù)卣{(diào)整基因圖譜，同樣可以得到各參數(shù)的迭代公式。

當(dāng)然，本文提出的方法也存在不足之處。由于標(biāo)記基因型帶有誤差和QTL基因型的未知性，當(dāng)標(biāo)記位點(diǎn)數(shù)目比較多時(shí)，算法中運(yùn)算量會(huì)非常大。另外，EM算法也有自身的局限之處，它收斂速度較慢，并且收斂速度對(duì)初始值的選擇有較大的依賴性。鑒于QTL多區(qū)間定位在基因病的研究中起到至關(guān)重要的作用，在未來的工作中對(duì)這類問題將會(huì)作進(jìn)一步考慮，尋求更高效的算法，以便適應(yīng)標(biāo)記位點(diǎn)更多的情況。

［1］ LANDER E S，BOTSTEIN D B.Mapping mendelian factors underlying quantitative traits using RFLP linkage mals.Genetics，1989，121：185－199.

［2］ ZENG Z B.Precision mapping of quantitative trait loci.Genetics，1994，135：1457－1468.

［3］ KAO C H，ZENG Z B，TEASDALE R D.Multiple interval mapping for quantitative trait loci.Genetics，1999，152：1203－1216.

［4］ GOLDSTEIN D R，ZHAO H，SPEED T P.The effects of genotyping errors and interference on estimation of genetics distance.Hum Hered，1997，47：86－100.

［5］ AKEY J M，ZHANG K，XIONG M，et al.The effect that genotyping errors have on the robustness of common linkage disequilibrium measures.Am J Hun Genet，2001，68：1447－1456.

［6］ LEBREC J J，PUTTER H，HOUWING－DUISTERMAAT J J，et al.Influence of Genotyping Errors in Linkage Mapping for Complex Traits－An Analytic Study.BMC Genetics，2008，9：57.

［7］ HOU Y J，MA W J，ZHOU Y，et al.Parameter estimation in quantitative trait loci mapping when using data with genotyping errors.The Proceedings of 2010 International Conference on Probability and Statistics of the International Institute for General Systems Studies，2010，1：236－240.

［8］ TONG L，MA W J，ZHOU Y，et al.Simultaneous estimation of QTL effects and positions when using genotype data with errors.J.Genet.2015，94：27 –34.

［9］ DEMPSTER A P，LAIRD N M，RUBIN D B.Maximum likelihood from incomplete.data via the EM algorithm.Journal of the Royal Statistical.Society，Series B，1977，39（1）：1－38.

［10］IBRAHIM J G.Incomplete data in generalized linear models.J.Am.Statist.Association，1990，85（411）：765－769.

［11］MA W J，ZHOU Y，et al.A two－step method for estimating QTL effects and positions.Genet Res.，2011，93：115－124.

［12］周影，韓國(guó)牛，史寧中，等.約束下多子女家系數(shù)據(jù)重組率的最大似然估計(jì)［J］.中國(guó)科學(xué)：數(shù)學(xué)，2010，40（10）：971－984.

［13］JEON G J.The effects of population size and dominance of quatitative trait loci（QTL）on the detection of linkage between markers and QTL for livestock.AJAS，1995，8：651－655.