基于組LASSO的重組近親雜交合作小鼠基因位點(diǎn)定位研究

熊思燦,胡桂開,阮周生

(東華理工大學(xué)理學(xué)院,江西 南昌330013)

1.引言

雜交合作(Collaborative Cross,CC)小鼠計劃于2004年正式提出[1],它由具有不同性狀的八個基本(founder)小鼠品系近親雜交約20代形成.理論上,CC小鼠的基因型應(yīng)該是純合的(homozygous).實(shí)際中,當(dāng)CC小鼠基因型中的純合子比例達(dá)到98%及其以上時,便停止繁殖[2].此外,理論上,CC小鼠的每個位點(diǎn)具有等可能,即1/8的概率繼承8個祖先中的任何一個.實(shí)際上,每個祖先的貢獻(xiàn)可能不同[3].

重組近親雜交(recombinant inbred intercrosses,RIX)實(shí)驗(yàn)是ZOU等于2005年為了克服重組近親(recombinant inbred,RI)品系的不足之處,如樣本量小等,而提出的.[4]基于RIX實(shí)驗(yàn)所得品系的數(shù)量性狀位點(diǎn)定位(quantitative trait loci,QTL)功效會得到提高.當(dāng)RIX實(shí)驗(yàn)的雜交小鼠選擇為CC品系時,將會得到CC-RIX品系.CC-RIX品系具有很多CC品系以及RIX品系不具有的優(yōu)點(diǎn),它是一個可重復(fù)(reproducible)產(chǎn)生的雜合子(heterozygous)多親群體,其基因變異與人類是相似的.[5]基于CC-RIX品系,已有部分研究,如GONG和ZOU[6]考慮了系數(shù)隨時間變化而變化的時變非參數(shù)數(shù)量性狀位點(diǎn)定位方法,Giusti-Rodríguez等[7]考慮了抑制精神病藥物的副作用的遺傳基礎(chǔ)問題,Graham等[5]考慮了細(xì)胞免疫表型的遺傳位點(diǎn)定位問題[5],LIU等[8]考慮了帶親源效應(yīng)(parent-of-origin,PoO)和主基因效應(yīng)的聯(lián)合建模問題,并采用了分塊Gibbs算法對模型進(jìn)行求解.分塊Gibbs算法屬于貝葉斯(Bayesian)算法的一種,具有如整塊更新,可減少收斂時間等優(yōu)點(diǎn).但當(dāng)位點(diǎn)個數(shù)較大時,其計算量依然很大.

懲罰函數(shù)法作為高維數(shù)據(jù)處理的又一方法,往往具有運(yùn)算速度快,收斂時間短等特點(diǎn).如LASSO (least absolute shrinkage and selection operator)[9],SCAD (smoothly clipped absolute deviation)[10],adaptive LASSO[11],以及組LASSO(Group LASSO)[12]等.考慮到CCRIX的每個位點(diǎn)的基因型可能來自于八個祖先,若其中的某個或者某些祖先對應(yīng)的效應(yīng)非零,則表明該位點(diǎn)為數(shù)量性狀位點(diǎn).因此,基于CC-RIX品系的數(shù)量性狀位點(diǎn)定位問題,本質(zhì)上是一個組變量選擇(group variable selection)問題.本文在文[8]模型的基礎(chǔ)上,忽略親源效應(yīng),考慮僅有主基因效應(yīng)的數(shù)量性狀位點(diǎn)定位模型,以及組LASSO懲罰函數(shù)求解算法.注意到本問題中的非光滑性(non-smoothness)以及不可分割性(non-separable),以及設(shè)計矩陣不滿足部分正交性(partial orthogonal property),導(dǎo)致一些流行算法,如坐標(biāo)下降(coordinate descent)算法不能直接使用.因此,本文采用了迭代加權(quán)最小二乘法(Iteratively Re-weighted Least-Squares,IRLS)進(jìn)行模型求解,為基于CC-RIX品系的復(fù)雜性狀位點(diǎn)定位提供參考.

2.模型建立

假設(shè)共有L個CC父系,將他們按照RIX實(shí)驗(yàn)的繁殖辦法,可得n(≤L(L ?1)/2)個CCRIX樣本.對每個樣本,假定其總的位點(diǎn)個數(shù)為p,且記其相應(yīng)的表型(phenotype,即因變量)為yi,(i=1,··· ,n).在LIU等[8]所建立的混合線性模型(linear mixed model,LMM)基礎(chǔ)上,忽略親源效應(yīng)后,可得如下僅含主基因效應(yīng)的定位模型：

其中,μ,βj= (βj1,··· ,βj8)′以及αk分別表示總均值,主基因效應(yīng)(即八個祖先的等位效應(yīng),founder allelic effect),以及隨機(jī)多基因效應(yīng)(random polygenic effect).xij=(xij1,··· ,xij8)T,當(dāng)?shù)趇個CC-RIX樣本,即CC-RIXi(i=1,2,··· ,n)在第j(j=1,2,··· ,p)個位點(diǎn)繼承了第k(k=1,2,··· ,8)個祖先的0個,1個或者2個基因時,相應(yīng)的xijk取0,1或者2.當(dāng)?shù)趇個CC-RIX樣本,即CC-RIXi(i= 1,2,··· ,n)由0個,1個或者2個父系CCk(k= 1,2,··· ,L)雜交形成時,相應(yīng)的aik取0,1或者2.顯然有,因?yàn)槊總€CC-RIX樣本有且僅有2個父系.為方便描述,記矩陣A= (aik)n×L.按照文[8]中的假定,對隨機(jī)效應(yīng)項αk(k=1,2,··· ,L),我們依然假定αk ～N(0,σ2a).其中,σ2a為多基因效應(yīng)方差.同時,對隨機(jī)誤差項ei(i= 1,2,··· ,n),我們依然假定其獨(dú)立同分布于正態(tài)分布,即ei ～N(0,σ2e).其中,σ2e為隨機(jī)誤差方差.

令y= (y1,··· ,yn)T,μ=μ(1,··· ,1)T,x= (x1,··· ,xp),xj= (x1j,··· ,xnj)T,β=α=(α1,··· ,αL)T,以及e=(e1,··· ,en)T,則模型(2.1)可改寫為

模型(2.2)在文[8]中的研究中,被用作比較模型.正如作者所述,當(dāng)不存在親緣效應(yīng)時,該模型的表現(xiàn)是與含親源效應(yīng)的模型表現(xiàn)是相似的.因此,在只為探測主基因效應(yīng)時,模型(2.2)不失為一個理想模型.接下來,本文重點(diǎn)從組LASSO懲罰函數(shù)法的角度來求解模型(2.2).

3.組LASSO方法

由隨機(jī)多基因效應(yīng)α ～NL(0,σ2aIL),可得其概率密度函數(shù)為

因此,給定α的條件下,y的條件分布為

其相應(yīng)的條件概率密度函數(shù)為

從而,偽數(shù)據(jù)(pseudo-data){(y,x,A,α)}的全似然函數(shù)為

對隨機(jī)效應(yīng)項α1,··· ,αL進(jìn)行積分,可得觀測數(shù)據(jù)的似然函數(shù)為

注意到

從而,似然函數(shù)可改寫為

一般而言,因位點(diǎn)個數(shù)p遠(yuǎn)大于樣本量n,即p ?n,經(jīng)典的極大似然估計法,以及限制極大似然估計法(restricted maximum likelihood approach)[13]將不再適用.注意到,盡管位點(diǎn)個數(shù)較大,但是真正起作用的卻較少.因此,本文假定主基因效應(yīng),即β=(βT1,··· ,βTp)T滿足稀疏性條件,也即大量的βj(j= 1,··· ,p)是為零的.考慮到βj= (βj1,··· ,βj8)T是一個8×1維向量,如果其中的部分分量不為零,則其對應(yīng)的位點(diǎn)對表型存在顯著性影響.因此,對CC-RIX 品系的數(shù)量性狀位點(diǎn)定位問題本質(zhì)上是組變量選擇問題.為此,我們采用負(fù)對數(shù)似然函數(shù),外加組LASSO懲罰函數(shù)的方法來獲取主基因效應(yīng)的稀疏解.求解的目標(biāo)函數(shù)為：

這里,λ為非負(fù)的正則化參數(shù),∥βj∥= (βTj βj)1/2(j= 1,··· ,p)為βj的l2- 模.此處,我們采用l2-模表示主基因效應(yīng)的分量要么全部為零,要么不全為零,以此實(shí)現(xiàn)組變量選擇[12].當(dāng)某個位點(diǎn)的主基因效應(yīng)的所有分量全部為零時,相應(yīng)的位點(diǎn)不是數(shù)量性狀位點(diǎn).否則,相應(yīng)的位點(diǎn)為數(shù)量性狀位點(diǎn).參數(shù)的最優(yōu)解可表示為如下最小化問題：

盡管最小化問題(3.3)對所有參數(shù)是一個非凸(non-convex)最優(yōu)化問題,難于求解.不過,當(dāng)固定μ,σ2e以及σ2a時,該問題為一個易于求解的凸最優(yōu)化問題.注意到我們的主要目的是進(jìn)行數(shù)量性狀位點(diǎn)定位,其效應(yīng)值的估計則變?yōu)榇渭壞繕?biāo).因此,我們將對輪廓對數(shù)似然函數(shù)(profile log-likelihood function),并采用迭代加權(quán)最小二乘法進(jìn)行求解.

4.迭代加權(quán)最小二乘法

為了進(jìn)行數(shù)量性狀位點(diǎn)定位,我們對似然函數(shù)(3.1)取對數(shù),并去掉其中的常數(shù)項,得到如下的以負(fù)輪廓對數(shù)似然函數(shù)表示的目標(biāo)函數(shù)：

這里,且為了簡化計算,我們將γ以及σ2e當(dāng)成兩個固定常數(shù),這可以事先指定.這種等價替代(proxy)思想已被FAN等[14]研究過.他們指出,在一定的條件下,求解此等價替代問題依然能獲得正確的模型選擇結(jié)果,但可能會導(dǎo)致額外的估計偏差.本文,我們首先最小化(4.1)進(jìn)行數(shù)量性狀位點(diǎn)定位,然后采用限制極大似然估計法得到其他參數(shù)的估計值.為此,記γ以及σ2e的估計值為以及,矩陣D的相應(yīng)估計值為.不失一般性,假定μ≡0.否則,我們可以用=(1n,x)和=(μ,βT)T分別替換x和β.

這里,?>0為光滑化參數(shù).之所以進(jìn)行光滑化處理,其主要目的是可將等式(4.2)中的第二部分再次轉(zhuǎn)換成β的二次型形式,從而便于構(gòu)造迭代加權(quán)最小二乘法.

本文所研究的問題中,基因型矩陣x由0,1和2組成,容易導(dǎo)致矩陣xTx奇異.因此,標(biāo)準(zhǔn)的最優(yōu)化方法,如牛頓?拉斐遜算法(Newton-Raphson algorithm)將不再適合.為此,本文提出一種迭代加權(quán)最小二乘法來求解問題(4.2)的最小值.

迭代加權(quán)最小二乘法是一種用于求解特定優(yōu)化問題,如壓縮感知(sparse recovery)[15],穩(wěn)健回歸(robust regression)等的常用方法.該方法是一個逐漸迭代的過程,其每一步更新都會求解一個加權(quán)最優(yōu)化問題.假定在第k步迭代中,β的迭代值是β(k),其第k+1步的迭代值由下式給出：

最小化(β)的必要條件是(β)對β的偏導(dǎo)數(shù)為零,即

從而可得β(k+1)的顯示表達(dá),即

綜合上述分析,當(dāng)給定調(diào)節(jié)參數(shù)λ后,迭代加權(quán)最小二乘法的更新過程如下：

算法1(IRLS 算法)

步1 初始化參數(shù)：?= 10?6,γ= 1,β(0)= (0,··· ,0)T,以及k= 0.計算

步2 按前述公式計算Λ(k);

步3 令k=k+1,計算

步4 重復(fù)步2和步3直到∥β(k)?β(k+1)∥<δ(=10?6)時,停止迭代.

上述迭代加權(quán)最小二乘法,除了第三步每次都需要基于新的Λ(k)值,計算一個8p×8p矩陣的逆矩陣之外,其余步驟的計算都能快速高效地完成.第三步中矩陣逆矩陣的計算可以使用喬里斯基分解(Cholesky decomposition)[16]來進(jìn)行快速計算.因此,當(dāng)?shù)訖?quán)最小二乘法收斂時,記β的最優(yōu)值為β(?).而迭代加權(quán)最小二乘法的收斂性,將在接下來的一節(jié)中進(jìn)行討論.

5.迭代加權(quán)最小二乘法的收斂性

定理5.1假設(shè){β(k)}∞k=1是由迭代加權(quán)最小二乘法產(chǎn)生的β的估計序列,則該序列的極限存在,且該極限值使得目標(biāo)函數(shù)F?λ(β)達(dá)最小.

證眾所周知,最小化F?λ(β)的迭代加權(quán)最小二乘法,等價于最小化如下輔助(auxiliary)函數(shù)[15]：

這里,ω=(ω1,··· ,ωp)T.在第k步迭代中,ω的值ω(k)由下式給出

其中,β(k)為第k次迭代中β的估計值.易得,ω(k)的第j個分量

進(jìn)而,再由式(4.3)可得第k+1步β的估計值為

即最小化H(β,ω(k))和最小化等價,從而

注意到,如下不等式

對任意的k ≥0成立,故迭代加權(quán)最小二乘法產(chǎn)生的序列是一個使得函數(shù)值單調(diào)不

再由如下事實(shí),

事實(shí)上,β(km+1)滿足必要條件(4.5),即

從而,多次使用式(5.3)可得

這里,λmin(Λ(km))是對角正定矩陣Λ(km)的最小特征值.將上述不等式兩邊分別關(guān)于m ≥1求和,再注意到從而有

這將導(dǎo)出式(5.2).

對式(5.3)的兩邊分別關(guān)于m取極限,可得

這里,是一個8p ×8p對角矩陣,且是一個8×8對角矩陣.注意到,由式(4.1)所定義的目標(biāo)函數(shù)Fλ?(β)是一個光滑的嚴(yán)格凸函數(shù),其唯一的最小值點(diǎn)滿足如下必要條件

對比(5.4)和(5.5)可知,迭代加權(quán)最小二乘法所得的序列極限β(?)是最小化問題(4.1)的最優(yōu)解.證畢！

6.調(diào)節(jié)參數(shù)的選取

調(diào)節(jié)參數(shù)λ的選取對算法的表現(xiàn)十分重要,其常用的選擇方法有交叉驗(yàn)證(cross-validation),AIC (Akaike information criterion)[17]和BIC (Bayesian information criterion)[18]等等.本文,我們采用BIC準(zhǔn)則來選取調(diào)節(jié)參數(shù)λ.

給定調(diào)節(jié)參數(shù)λ,記迭代加權(quán)最小二乘法收斂所得β的最優(yōu)估計為

眾所周知,回歸模型可能存在將非零系數(shù)估計為零,或者將零系數(shù)估計為非零的可能.為此,我們需要選定一個合適的容許誤差,記為tols.若對某個1≤k ≤8,有(λ)(1≤j ≤p)≥tols,其相應(yīng)的位點(diǎn)j將被納入活動集(active set)?(0)(λ)中.否則,該位點(diǎn)將被納入非活動集(non-active set)((λ))C中,并將其對應(yīng)的效應(yīng)值壓縮為零.

為了選取最優(yōu)調(diào)節(jié)參數(shù)λ,我們需要針對其不同的可能取值,采用迭代加權(quán)最小二乘法.實(shí)際運(yùn)算中,我們將log(λmax)至log(λmin)等分得到從大到小排列的N個取值,并記為λ=(λ1,··· ,λN),然后按照下述步驟確定最優(yōu)調(diào)節(jié)參數(shù).

步1 令i=0,tols=0.001.記此時迭代加權(quán)最小二乘法所得的β的最優(yōu)估計為(λ).

步2 令i=i+1,λ=λi.取(λ)為β的初始估計,通過迭代加權(quán)最小二乘法可得β的第i步估計(λ).分別記此時的活動集和非活動集為(λ)和((λ))C,可按下式計算得BIC的值.

步3 當(dāng)i ≤N時,返回步2.否則,轉(zhuǎn)向步4.

步4 使得由式(6.1)計算所得的BIC值最小的對應(yīng)λ,不妨記為λs,即為最優(yōu)的調(diào)節(jié)參數(shù).與之對應(yīng)的活動集={i1,··· ,iS},1≤i1<···

7.模擬計算

為了評價本文所提方法的優(yōu)劣,我們進(jìn)行了模擬計算.這里,我們假定CC-RI父系共有L= 100個.按照ZOU等[4]所提的循環(huán)(loop)實(shí)驗(yàn)設(shè)計,將所有的L個CC-RI排成環(huán)形,然后每個CC-RI均與接下來的J= 3的CC-RI品系雜交,從而可產(chǎn)生n= 300個CC-RIX樣本.進(jìn)一步,假定有p=101和p=301個位點(diǎn)等間距分布在染色體上,其中有2個等間隔分布的數(shù)量性狀位點(diǎn).數(shù)量性狀位點(diǎn)所對應(yīng)的主基因效應(yīng)方差σ2j為1或者5,而非數(shù)量性狀位點(diǎn)所對應(yīng)的主基因效應(yīng)方差為0.主基因效應(yīng)βj由8維正態(tài)分布N(0,σ2jI8)產(chǎn)生.多基因效應(yīng)方差σ2a設(shè)置為0.1,而殘差方差σ2e設(shè)置為0.1或者1.從而,交叉組合可得8種不同的參數(shù)組合設(shè)置.

對每一種參數(shù)組合,隨機(jī)產(chǎn)生100個模擬數(shù)據(jù)集.其中,每個CC-RI品系的單個位點(diǎn)基因均有1/8的概率繼承8個祖先中的任何一個,按照RIX循環(huán)實(shí)驗(yàn)設(shè)計,即可得到基因型矩陣xj,以及親緣信息矩陣A.

對每一種參數(shù)組合下的每一個數(shù)據(jù)集,均采用迭代加權(quán)最小二乘法進(jìn)行數(shù)量性狀位點(diǎn)定位.為了衡量定位效果,我們將其處理成一個二分類的預(yù)測問題.如果某個位點(diǎn)為模擬設(shè)置的QTL位點(diǎn),則該位點(diǎn)標(biāo)記為P(positive,陽性),否則標(biāo)記為N(negative,陰性).如果一個QTL位點(diǎn)預(yù)測為模擬設(shè)置的QTL位點(diǎn),則稱為真陽性(true positive,TP),否則稱為假陽性(false positive,FP).如果一個預(yù)測非QTL位點(diǎn),為模型設(shè)置的非QTL位點(diǎn),稱之為真陰性(true negative,TN),否則稱之為假陰性(false negative,FN).從而可得表(7.1)所示的混淆矩陣.

表7.1 混淆矩陣

基于混淆矩陣,定義真陽性率(True Positive Rate,TPR),假陽性率(False Positive Rate,FPR),假發(fā)現(xiàn)率(False Discovery Rate,FDR)如下.

同時,對每一種參數(shù)組合下的每一個數(shù)據(jù)集,若預(yù)測的QTL恰好與模擬設(shè)置的QTL完全相同,稱之為“真實(shí)”(truth).綜合100次模擬,即100個數(shù)據(jù)集所得的總的“真實(shí)”頻率稱之為真實(shí)率(Truth Rate,TR).顯然,TPR以及TP越大越好,而FPR以及FDR越小越好.表(7.2)展示了不同模擬設(shè)置下,100次模擬所得的TPR,FPR以及FDR的平均值,以及真實(shí)率(TR)情況.

表7.2 不同模擬設(shè)置下的模型平均表現(xiàn)

從表(7.2)中,不難發(fā)現(xiàn),無論哪種模擬設(shè)置下,本文所提模型和方法的表現(xiàn)都非常不錯.例如,當(dāng)p= 101時,真陽性率都超過0.95,而對應(yīng)的假陽性率和假發(fā)現(xiàn)率卻十分接近于0.此外,真實(shí)率TR也十分接近于1.這充分表明本文所提模型和方法在進(jìn)行數(shù)量性狀位點(diǎn)定位中的有效性.若固定其他參數(shù)不變,當(dāng)非零數(shù)量性狀位點(diǎn)對應(yīng)的主基因效應(yīng)方差增大時,模型的表現(xiàn)會越來越好.因?yàn)?此時信號強(qiáng)度增強(qiáng),模型更容易識別真實(shí)的QTL位點(diǎn).若固定其他參數(shù),當(dāng)模型的殘差方差增大時,模型的表現(xiàn)會有所下滑.如情形1和情形2,當(dāng)σ2e由0.1增大到1時,相應(yīng)的TPR從0.995減少到0.95.因?yàn)榇藭r噪聲比例有所加大,從而模型表現(xiàn)整體變差.這一現(xiàn)象與我們的預(yù)期是一致的.整體而言,本文所提模型和方法在數(shù)量性狀位點(diǎn)識別方面具有較高的真陽性率,以及較低的假陽性率,是一種比較好的定位方法.

8.結(jié)論

本文在文[8]的模型的基礎(chǔ)上,忽略親源效應(yīng),考慮了一個僅含主基因效應(yīng)的數(shù)量性狀位點(diǎn)定位問題.考慮到研究背景,即CC-RIX品系的數(shù)據(jù)特點(diǎn),我們采用了組LASSO方法對模型進(jìn)行轉(zhuǎn)換,并設(shè)計了迭代加權(quán)最小二乘法求解模型,克服了設(shè)計矩陣容易奇異等計算難題.從模擬所得結(jié)果來看,本文所提模型和方法,在CC-RIX品系的基因位點(diǎn)定位中具有較好的表現(xiàn),能準(zhǔn)確識別出真實(shí)的數(shù)量性狀位點(diǎn),并且具有較高的真陽性率以及較低的假陽性率.同時,相比貝葉斯方法而言,本文所提模型和方法還具有計算量小,計算速度快等特點(diǎn).當(dāng)位點(diǎn)個數(shù)較多時,能體現(xiàn)出其計算上的顯著優(yōu)勢.