基于Cox半?yún)?shù)模型檢測(cè)生存性狀的印記數(shù)量性狀位點(diǎn)

周曉晶，方銘

（1.黑龍江八一農(nóng)墾大學(xué)理學(xué)院數(shù)學(xué)系，大慶163319；2.黑龍江八一農(nóng)墾大學(xué)生命科學(xué)技術(shù)學(xué)院）

印記基因是指子代中來(lái)源于親本的等位基因僅一方親本的基因表達(dá)，而來(lái)自另一親本的基因表達(dá)總是受到抑制，導(dǎo)致該基因呈現(xiàn)單一來(lái)源表達(dá)的遺傳現(xiàn)象。自從1991年印記現(xiàn)象被發(fā)現(xiàn)后[1]，研究者在植物、動(dòng)物及人類身上觀測(cè)到了許多類似于印記的現(xiàn)象，隨著基因印記研究的進(jìn)展，陸續(xù)有學(xué)者基于區(qū)間定位方法提出了檢測(cè)印記數(shù)量性狀位點(diǎn)的遺傳模型和統(tǒng)計(jì)方法[2-6]及貝葉斯定位方法[7-8]。

印記基因不但控制生長(zhǎng)發(fā)育性狀，也有可能控制生存性狀。這類帶有“右尾”偏態(tài)分布的性狀常見(jiàn)于植物的開(kāi)花期、動(dòng)物的存活時(shí)間和死亡時(shí)間。目前，大多數(shù)定位數(shù)量性狀位點(diǎn)的統(tǒng)計(jì)方法是在正態(tài)分布的假設(shè)下，顯然不適于分析生存性狀。因此，研究者將生存分析中非參數(shù)、參數(shù)及半?yún)?shù)的模型合并到檢測(cè)生存性狀位點(diǎn)的區(qū)間定位中[9-11]?；诩铀偈Ｐ?，Cheng和Tzeng提出了參數(shù)和半?yún)?shù)的區(qū)間定位方法[12]。而這兩種方法都涉及到誤差分布的選擇。采用Cox比例風(fēng)險(xiǎn)模型和偏似然函數(shù)方法的一個(gè)明顯優(yōu)勢(shì)就是避免了選擇經(jīng)常令研究者“討厭”的基線風(fēng)險(xiǎn)函數(shù)，其復(fù)雜的估計(jì)通常都比較耗時(shí)。Fang基于Cox比例風(fēng)險(xiǎn)模型，提出了在偏似然函數(shù)的框架內(nèi)表征數(shù)量性狀位點(diǎn)的方法[13]，但他采用的最小二乘算法估計(jì)QTL參數(shù)，即用QTL基因型指示變量的先驗(yàn)概率的數(shù)學(xué)期望去估計(jì)指示變量，而這個(gè)過(guò)程中會(huì)產(chǎn)生“誤差膨脹”。此外，如上提到的方法主要討論以孟德?tīng)柗绞竭z傳的QTL的相關(guān)問(wèn)題[14]，如QTL在基因組的位置，遺傳效應(yīng)、QTL之間的互作等?；蚪M印記現(xiàn)象，也被稱為父母起源效應(yīng)，是由每個(gè)父母基因組對(duì)后代表型的不等遺傳貢獻(xiàn)導(dǎo)致的。目前，還沒(méi)有文獻(xiàn)研究應(yīng)用半?yún)?shù)方法檢測(cè)生存性狀的印記數(shù)量性狀位點(diǎn)。

首先借助于Cox半?yún)?shù)比例風(fēng)險(xiǎn)模型構(gòu)造了iQTL對(duì)生存時(shí)間的效應(yīng)，此模型的好處就是不需要對(duì)性狀分布做任何假設(shè)。其次，為檢測(cè)iQTL，在偏似然函數(shù)的框架內(nèi)，研究參數(shù)的期望最大化算法，并實(shí)現(xiàn)了整個(gè)基因組范圍內(nèi)的統(tǒng)計(jì)顯著性的檢驗(yàn)。最后，利用一個(gè)公開(kāi)發(fā)表的老鼠模型數(shù)據(jù)集證實(shí)了此方法的統(tǒng)計(jì)行為。

1 方法

1.1 遺傳模型

假設(shè)在給定的群體中每個(gè)位點(diǎn)有四個(gè)可區(qū)分的基因型，記為QMQP，QMqP，qMQP和qMqP，這里下腳標(biāo)代表父母起源，觀測(cè)到帶有已知連鎖圖譜的n個(gè)個(gè)體的生存時(shí)間（T：t1，t2，L，tn），并標(biāo)記基因型。假設(shè)單個(gè)QTL位于相鄰的兩個(gè)標(biāo)記Mk和Mk+1之間，控制目標(biāo)性狀?；贑ox比例風(fēng)險(xiǎn)模型構(gòu)造QTL對(duì)生存時(shí)間的效應(yīng)，其中風(fēng)險(xiǎn)函數(shù)為

這里，λ0（tj）為基線風(fēng)險(xiǎn)函數(shù)，即當(dāng)所有相互獨(dú)立的變量值為零時(shí)個(gè)體的風(fēng)險(xiǎn)函數(shù)[15]。注意，如果采用參數(shù)方法時(shí)，那么就需要估計(jì)基線風(fēng)險(xiǎn)函數(shù)，但半?yún)?shù)方法不需要。Xj=（1 zjwjsj）為按照生存時(shí)間大小排列的第j個(gè)個(gè)體所對(duì)應(yīng)的協(xié)變量值構(gòu)成的向量，β=（μa d i）T為待估計(jì)系數(shù)構(gòu)成的向量。zj，wj和sj分別為對(duì)應(yīng)于加性效應(yīng)、顯性效應(yīng)和印記效應(yīng)的基因型指示變量，這些指示變量由Mantey等人于2005年給出[16]，詳細(xì)如下：

1.2 極大似然估計(jì)

記R（tj）為觀測(cè)到的個(gè)體的風(fēng)險(xiǎn)集，其生存時(shí)間不小于tj，從而根據(jù)貝葉斯理論，觀測(cè)到的個(gè)體的生存時(shí)間的條件概率密度為

依賴于QTL的四個(gè)基因型，記為f（tj｜l），l=1，2，3，4分別對(duì)應(yīng)著QTL基因型QM QP，QM qP，qM QP和qMqP。令pj｜l，（l=1，2，3，4）為相應(yīng)的條件概率，從而構(gòu)造混合模型為

假設(shè)性狀值相互獨(dú)立，那么QTL參數(shù)的似然函數(shù)為n個(gè)個(gè)體上相互獨(dú)立的混合模型的連乘積：

此公式稱為偏似然函數(shù)[17]，進(jìn)一步，偏似然對(duì)數(shù)函數(shù)為：

對(duì)參數(shù)求導(dǎo)得：

接下來(lái)，采用期望最大化算法[18]求解β的極大似然估計(jì)，迭代步驟如下：

（1）選取 β=（μa d i）T的初始值

（2）將 β（0）帶入方程（7）。

（4）用 β（1）代替初始值 β（0），回到（2）繼續(xù)迭代。

（5）繼續(xù)迭代運(yùn)算，直到結(jié)果符合設(shè)定的收斂標(biāo)準(zhǔn)為止。得到的參數(shù)值作為極大似然估計(jì)，記為

1.3 顯著性檢驗(yàn)

一般的，顯著性檢驗(yàn)步驟如下：（i）在整個(gè)基因組范圍內(nèi)檢測(cè)QTL；（ii）確定檢測(cè)到的QTL的遺傳方式；（iii）鑒別印記QTL的印記類型。

具體檢驗(yàn)過(guò)程如下：

首先，構(gòu)造影響生存性狀的QTL的存在性的假設(shè)

其次，當(dāng)H0被拒絕時(shí)，檢測(cè)到的QTL的遺傳模式的假設(shè)

如果接受H0，那么檢測(cè)到的QTL以孟德?tīng)柗绞竭z傳；否則是印記的。

第三，檢測(cè)到的印記QTL的印記類型取決于遺傳效應(yīng)a或d是否等于印記效應(yīng)i，印記類型和相應(yīng)的零假設(shè)參見(jiàn)參考文獻(xiàn)[7]和[19]。

確定相鄰兩個(gè)標(biāo)記區(qū)間內(nèi)的顯著效應(yīng)處的QTL是否存在，可以通過(guò)計(jì)算零假設(shè)和備擇假設(shè)構(gòu)成的對(duì)數(shù)似然比檢驗(yàn)統(tǒng)計(jì)量，即，這里和＾1分別代表零假設(shè)H0和備擇假設(shè)H1下的未知參數(shù)的極大似然估計(jì)。為實(shí)現(xiàn)顯著性檢驗(yàn)，采用Churchill和Doerge提出的置換檢驗(yàn)方法[20]進(jìn)行。

2 實(shí)際例子分析

我們應(yīng)用提出的方法分析由高氧性急性肺損傷引發(fā)的生存時(shí)間的表觀遺傳屬性。本研究數(shù)據(jù)集來(lái)自一個(gè)公開(kāi)發(fā)表的小鼠模型系統(tǒng)，來(lái)源于一對(duì)互惠的近交系。其中B系（C57BL/6J）小鼠對(duì)高氧性急性肺損傷敏感，而S系小鼠對(duì)高氧性急性肺損傷具有較高的抗性。雌性B系小鼠與雄性S系小鼠交配（BS）和雌性S系小鼠與雄性B系小鼠交配（SB）產(chǎn)生F1代小鼠群體[21-22]。F1代小鼠間平均生存時(shí)間的差異提供了存在親代效應(yīng)的強(qiáng)有力的證據(jù)。F1代小鼠再采用相互交配的方式（BS×BS，BS×SB，SB×BS和SB×SB）的到4個(gè)可能的F2群體，每個(gè)群體大約產(chǎn)生200個(gè)個(gè)體。這些雜交個(gè)體可以用于估計(jì)潛在印記效應(yīng)的遺傳研究。共有840個(gè)F2代個(gè)體的生存時(shí)間（小時(shí)）表型值被測(cè)定，并且記錄了分布于小鼠全基因組（包括X染色體在內(nèi)）上的97個(gè)多肽微衛(wèi)星標(biāo)記。為了消除由于雌親、雄親及性別等系統(tǒng)環(huán)境因素對(duì)小鼠存活時(shí)間的影響，對(duì)觀測(cè)數(shù)據(jù)取對(duì)數(shù)進(jìn)行校正。在全基因組所有基因位點(diǎn)估計(jì)4種基因型的條件概率，對(duì)于有BS×BS或SB×SB的雜交方式僅得到一種雜合子，并且雜合子基因型與其親代基因型相同，我們可以使用這種方法利用F2雜交群體估計(jì)三種可能的QTL基因型的條件概率，在給定相鄰標(biāo)記間的QTL基因型包括一種雜合子和兩個(gè)純合子，而且把另一種交互雜合子設(shè)置為0；對(duì)于來(lái)至BS×SB和SB×BS雜交的F2代個(gè)體，可能攜帶交互雜合子中的一個(gè)。假設(shè)雌雄重組率為1.25∶1.0，目的是估計(jì)那些由于缺乏可區(qū)分的交互雜合子基因型的測(cè)量標(biāo)記基因型的條件概率。

檢測(cè)到的iQTL的參數(shù)估計(jì)值及印記模式參見(jiàn)表。圖描述了整個(gè)基因組范圍內(nèi)的對(duì)數(shù)似然比（LR）檢驗(yàn)統(tǒng)計(jì)量曲線。在顯著性水平為5%下經(jīng)驗(yàn)臨界值為16.7?？梢?jiàn)，在LR圖上有5個(gè)峰值超過(guò)了臨界值，證實(shí)了控制小鼠急性肺損傷的生存時(shí)間的iQTL的存在性，它們分別位于1號(hào)、4號(hào)、9號(hào)、15號(hào)和17號(hào)染色體上，17號(hào)染色體的QTL以孟德?tīng)柗绞竭z傳。

3 討論

大多數(shù)現(xiàn)存的表征數(shù)量性狀位點(diǎn)的方法是在表型值服從正態(tài)分布及完全觀測(cè)到的假設(shè)條件下的，然而，這些方法不適于分析生存數(shù)據(jù)，因?yàn)檫@類數(shù)據(jù)通常不服從正態(tài)分布并且常帶有刪失數(shù)據(jù)。研究者將生存分析中的參數(shù)的、半?yún)?shù)的和加速失效模型合并到生存性狀位點(diǎn)的區(qū)間定位的框架內(nèi)。而這些方法主要關(guān)注以孟德?tīng)柗绞竭z傳的QTL的檢測(cè)中。在參數(shù)比例風(fēng)險(xiǎn)模型中，需要估計(jì)復(fù)雜的基線風(fēng)險(xiǎn)函數(shù)，并且涉及到如何選擇最優(yōu)的生存函數(shù)去擬合生存曲線。在參數(shù)的加速失效模型中，需要選擇對(duì)應(yīng)于生存分布的誤差分布，也有可能當(dāng)選擇錯(cuò)誤的誤差分布時(shí)，會(huì)導(dǎo)致QTL的效應(yīng)和位置的估計(jì)有較大的偏差。參數(shù)模型通常都會(huì)涉及到分布模型誤指定問(wèn)題。

表1 控制小鼠急性肺損傷的生存時(shí)間的iQTL的參數(shù)估計(jì)Table1 Parameter Estimates of survival time for acute lung injury survival time ofmice

圖1 Cox半?yún)?shù)區(qū)間定位的LR檢驗(yàn)統(tǒng)計(jì)量曲線Fig.1 The profile of LR test statistics obtained with the intervalmapping basedon Cox semi-parametricmodel for HALIsurvival time inmice

借助于Cox比例風(fēng)險(xiǎn)模型，給出了半?yún)?shù)的區(qū)間定位方法去表征控制生存時(shí)間的數(shù)量性狀位點(diǎn)，此方法不需要估計(jì)基線風(fēng)險(xiǎn)函數(shù)。然后采用期望最大化算法得到印記參數(shù)的極大似然估計(jì)值，這些參數(shù)包括描述iQTL的遺傳效應(yīng)和印記類型。此方法成功應(yīng)用于分析一個(gè)公開(kāi)發(fā)表的數(shù)據(jù)集，圖和表分別描述了采用本研究所給出的統(tǒng)計(jì)方法檢測(cè)到的顯著iQTL所在基因組的位置及印記性質(zhì)的剖析。

[1]DeChiara，T.M.，Robertson，E.J.，&Efstratiadis，A.Parental imprinting of themouse insulin-like growth factor II gene[J].Cell，1991，64（4），849-859.

[2]Knott，S.A.，Marklund，L.，Haley，C.S.，et al.Multiple marker mapping of quantitative trait loci in a cross between outbred wild boar and large white pigs[J].Genetics1998，149（2），1069-1080.

[3]de Koning，D.J.，Rattink，A.P.，Harlizius，B.，et al.Genomewide scan for body composition in pigs reveals important role of imprinting[J].Proc Natl Acad Sci U SA，2000，97（14），7947-7950.

[4]de Koning，D.J.，Bovenhuis，H.，&van Arendonk，J.A.On the detection of imprinted quantitative trait loci in experimental crosses of outbred species[J].Genetics，2002，161（2），931-938.

[5]Tuiskula-Haavisto，M.，de Koning，D.J.，Honkatukia，M.，et al.Quantitative trait loci with parent-of-origin effects in chicken[J].Genet Res，2004，84（1），57-66.

[6]Cui，Y.，Lu，Q.，Cheverud，J.M.，et al.Model for mapping imprinted quantitative trait loci in an inbred F2 desig[J].Genomics，2006，87（4），543-551.

[7]Yang，R.，Wang，X.，Wu，Z.，et al.Bayesian model selection for characterizing genomic imprinting effects and patterns[J].Bioinformatics，2010，26（2），235-241.

[8]Hayashi，T.，&Awata，T.A Bayesian method for simultaneously detecting Mendelian and imprinted quantitative trait loci in experimental crosses of outbred species[J].Genetics，2008，178（1），527-538.

[9]Symons，R.C.，Daly，M.J.，F(xiàn)ridlyand，J.，et al.Multiple genetic locimodify susceptibility to plasmacytoma-related morbidity in E（mu）-v-abl transgenic mice[J].Proc Natl Acad SciUSA，2002，99（17），11299-11304.

[10]Broman，K.W.Mapping quantitative trait loci in the case of a spike in the phenotype distribution[J].Genetics，2003，163（3），1169-1175.

[11]Diao，G.，Lin，D.Y.，&Zou，F(xiàn).Mapping quantitative trait loci with censored observations[J].Genetics，2004，168（3），1689-1698.

[12]Cheng，J.Y.，&Tzeng，S.Parametric and semiparametric methods for mapping quantitative trait loci[J].Computational Statistics and Data Analysis，2009，53，1843-1849.

[13]Fang，Y.A note on QTL detecting for censored traits[J].Genet Sel Evol，2006，38（2），221-229.

[14]Lander，E.S.，&Botstein，D.Mapping mendelian factors underlying quantitative traits using RFLP linkage maps[J].Genetics，1989，121（1），185-199.

[15]Kalbeisch，J.D.，&Prentice，R.L.The Statistical Analysis of Failure Time Data（Second ed.） [M].New York:Wiley，2002.

[16]Mantey，C.，Brockmann，G.A.，Kalm，E.，et al.Mapping and exclusion mapping of genomic imprinting effects in mouse F2 families[J].JHered，2005，96（4），329-338.

[17]Cox，D.R.Regression Models and Life-Tables[J].JR Statis Soc，1972，34（2），187-220.

[18]Dempster，A.P.，Laird，N.M.，&Rubin，D.B.Maximum likelihood from incomplet data via the EM algorithm[J].JR Statist Soc B，1977，39（1），1-38.

[19]Cheverud，J.M.，Hager，R.，Roseman，C.，et al.Genomic imprinting effects on adult body composition in mice[J].Proc Natl Acad SciUSA，2008，105（11），4253-4258.

[20]Churchill，G.A.，&Doerge，R.W.Empirical threshold values for quantitative trait mapping[J].Genetics，1994，138（3），963-971.

[21]Prows，D.R.，Hafertepen，A.P.，Gibbons，W.J.，et al.A genetic mouse model to investigate hyperoxic acute lung injury survival[J].Physiol Genomics，2007，30（3），262-270.

[22]Prows，D.R.，Hafertepen，A.P.，Winterberg，A.V.，et al.Genetic analysis of hyperoxic acute lung injury survival in reciprocal intercross mice[J].Physiol Genomics，2007，30（3），271-281.