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

熊思燦,胡桂開,阮周生

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

1.引言

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

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

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

2.模型建立

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

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

令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]中的研究中,被用作比較模型.正如作者所述,當不存在親緣效應時,該模型的表現是與含親源效應的模型表現是相似的.因此,在只為探測主基因效應時,模型(2.2)不失為一個理想模型.接下來,本文重點從組LASSO懲罰函數法的角度來求解模型(2.2).

3.組LASSO方法

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

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

其相應的條件概率密度函數為

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

對隨機效應項α1,··· ,αL進行積分,可得觀測數據的似然函數為

注意到

從而,似然函數可改寫為

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

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

盡管最小化問題(3.3)對所有參數是一個非凸(non-convex)最優化問題,難于求解.不過,當固定μ,σ2e以及σ2a時,該問題為一個易于求解的凸最優化問題.注意到我們的主要目的是進行數量性狀位點定位,其效應值的估計則變為次級目標.因此,我們將對輪廓對數似然函數(profile log-likelihood function),并采用迭代加權最小二乘法進行求解.

4.迭代加權最小二乘法

為了進行數量性狀位點定位,我們對似然函數(3.1)取對數,并去掉其中的常數項,得到如下的以負輪廓對數似然函數表示的目標函數：

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

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

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

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

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

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

綜合上述分析,當給定調節參數λ后,迭代加權最小二乘法的更新過程如下：

算法1(IRLS 算法)

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

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

步3 令k=k+1,計算

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

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

5.迭代加權最小二乘法的收斂性

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

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

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

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

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

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

注意到,如下不等式

對任意的k ≥0成立,故迭代加權最小二乘法產生的序列是一個使得函數值單調不

再由如下事實,

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

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

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

這將導出式(5.2).

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

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

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

6.調節參數的選取

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

給定調節參數λ,記迭代加權最小二乘法收斂所得β的最優估計為

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

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

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

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

步3 當i ≤N時,返回步2.否則,轉向步4.

步4 使得由式(6.1)計算所得的BIC值最小的對應λ,不妨記為λs,即為最優的調節參數.與之對應的活動集={i1,··· ,iS},1≤i1<···

7.模擬計算

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

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

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

表7.1 混淆矩陣

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

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

表7.2 不同模擬設置下的模型平均表現

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

8.結論

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