趙迎利,朱 旭
(1.西安交通大學(xué)城市學(xué)院數(shù)學(xué)教研室,陜西 西安 710018; 2.西安交通大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院,陜西 西安 710049)
近年來,影像遺傳學(xué)(Imaging Genetics)成為一種新興的交叉學(xué)科,其主要是結(jié)合影像學(xué)和遺傳學(xué)方法檢測腦結(jié)構(gòu)與神經(jīng)精神疾病、認知和情緒調(diào)節(jié)等行為相關(guān)腦功能的遺傳變異研究[1]。醫(yī)學(xué)影像的功能性磁共振成像(Functional Magnetic Resonance Imaging, FMRI)數(shù)據(jù)和遺傳學(xué)的單核苷酸多態(tài)性(Single Nucleotide Polymorphisms, SNP)數(shù)據(jù)提供了精神類疾病不同但互補的信息,聯(lián)合并挖掘這2種類型的數(shù)據(jù)對于精神類疾病的探究可以提供實質(zhì)性的幫助,然而此類數(shù)據(jù)往往具有樣本較少和維度較高的特點,因此,從大量的FMRI數(shù)據(jù)和SNP數(shù)據(jù)中聯(lián)合分析出與精神類疾病相關(guān)的生物標記物成為眾多研究學(xué)者面對的主要難題。
典型相關(guān)分析[2](Canonical Correlation Analysis, CCA)和偏最小二乘[3](Partial Least Squares, PLS)被提出來分析2組數(shù)據(jù)間的關(guān)系,其主要是通過最大化2組數(shù)據(jù)集變量間的線性組合之間的相關(guān)性進行關(guān)聯(lián)分析,但是由于SNP的小樣本高維度問題,即生物標志物的數(shù)量遠大于樣本數(shù)量,這種問題通常會導(dǎo)致多重共線性問題和過擬合的問題出現(xiàn)[4]。為了解決該問題,眾多研究學(xué)者致力于使用lasso、group lasso等以l1范數(shù)和l2,1范數(shù)為稀疏懲罰項的正則化方法進行高維基因表達數(shù)據(jù)的稀疏,例如Parkhomenko等人[4-6]應(yīng)用lasso懲罰項提出稀疏典型相關(guān)分析(Sparse Canonical Co rrelation Analysis, SCCA),其主要用來分析雙多變量間的關(guān)系以便從高維的基因表達數(shù)據(jù)中選擇出與疾病相關(guān)的潛在致病機制??紤]到精神類疾病重要的先驗信息通常以感興趣腦區(qū)(ROI)和富集區(qū)基因塊(LD)的形式存在,若直接使用SCCA會忽視本體潛在的生物群組特征[7-8],Lin等人[9]聯(lián)合應(yīng)用lasso懲罰項和group lasso懲罰項共同懲罰典型相關(guān)分析用來分析挖掘出全腦影像數(shù)據(jù)和全基因數(shù)據(jù)重要的組間變量信息和組內(nèi)變量信息。
雖然上述方法均可以識別出2組高維數(shù)據(jù)間的潛在特征或者潛在群組特征,但是仍然存在不足:一方面,為了求解方便,上述方法均假設(shè)特征之間正交,另一方面,基于lasso和group lasso的稀疏正則化項l1范數(shù)和l2,1范數(shù)往往傾向于只選擇出其中的單個特征或者一組組特征,但是基因之間和腦影像數(shù)據(jù)之間往往存在某種關(guān)聯(lián)[10],并且特征之間并不存在正交的特性,因此并不能完整地選擇出具有相關(guān)性的多個關(guān)聯(lián)特征或者多組組關(guān)聯(lián)特征,容易造成過稀疏的問題出現(xiàn),不能識別出更多的與精神類疾病相關(guān)的生物標志物。
考慮到上述問題,本文針對FMRI和SNP數(shù)據(jù)的特點,提出一種基于l1,2范數(shù)[11]和l2,1范數(shù)聯(lián)合懲罰典型相關(guān)分析的特征選擇的方法。該方法通過數(shù)據(jù)驅(qū)動獲得群組結(jié)構(gòu)并將其作為先驗信息,將可以實現(xiàn)組內(nèi)稀疏和保持組間特征平滑的正則化項l1,2范數(shù)和可以實現(xiàn)組間稀疏和保持組內(nèi)特征平滑的正則化項l2,1范數(shù)聯(lián)合共同懲罰典型相關(guān)分析,并保持數(shù)據(jù)間原有的不正交的特性,以互相制約具有相關(guān)性的群組特征和組內(nèi)特征的稀疏性能,有效地解決具有相關(guān)性的群組或者多個關(guān)聯(lián)特征之間的特征選擇問題。仿真實驗結(jié)果表明該方法不但可以估計出2組數(shù)據(jù)間的相關(guān)系數(shù),而且可以選擇出具有相關(guān)性的組間特征和組內(nèi)特征,在精神分裂癥數(shù)據(jù)集上的實驗表明該方法可以選擇出更多的與精神分裂癥相關(guān)的易感基因和風險腦區(qū)。
典型相關(guān)分析是一種經(jīng)典的多元統(tǒng)計分析方法,經(jīng)常被用于機器學(xué)習(xí)領(lǐng)域,其主要是為了2組多維變量X∈Rn×p和Y∈Rn×q分別尋找基向量u∈Rp和v∈Rq,使得隨機變量uTX和vTY之間的相關(guān)性達到最大。模型的表示具體如下:
s.t.uTXTXu=1vTYTYv=1
(1)
其中u和v分別表示X和Y的基向量,也被稱為X和Y的典型向量,Xu和Yv被稱為X和Y的典型變量。
由于現(xiàn)有的基因表達數(shù)據(jù)通常具有樣本數(shù)目遠小于特征數(shù)目的特點,直接使用CCA易出現(xiàn)過擬合現(xiàn)象。為了解決此問題,Parkhomenko等人提出了稀疏典型相關(guān)分析,其主要思想是使用lasso的懲罰項,使不重要的特征的權(quán)重很小或者為0,重要特征的權(quán)重較大,從而稀疏掉不重要的特征,以便從大量的特征中找到與精神類疾病相關(guān)的潛在致病機制,具體的模型表示如下:
s.t.‖u‖2≤1,‖v‖2≤1,‖u‖1≤c1,‖v‖1≤c2
(2)
其中,‖u‖1和‖v‖1是lasso懲罰,c1和c2控制典型向量u和v的稀疏程度。
神經(jīng)影像遺傳學(xué)數(shù)據(jù)中重要的特征通常存在于大腦中的感興趣區(qū)域和基因中的連鎖不平衡塊之間,SCCA忽略基因組和腦區(qū)的群組效應(yīng),為了解決該問題,Lin等人提出GSCCA,使用group lasso的懲罰項計算出群組特征的權(quán)重是否接近于0,從而稀疏掉不重要的群組,再次考慮到群組內(nèi)特征選取的重要性,繼續(xù)添加lasso的懲罰項計算群組內(nèi)特征的權(quán)重,從而可同時挖掘組間重要的特征和組內(nèi)重要的特征,此模型的表示具體如下:
(3)
影像學(xué)中,腦區(qū)之間存在著一定的聯(lián)系,比如杏仁核活性一旦升高,下丘腦、海馬、前扣帶回以及前額葉皮層等都會得到相應(yīng)的反饋。遺傳學(xué)上,具有相似功能的SNP位點通常位于同一基因,基因組中的基因之間也是存在功能聯(lián)系的,而組稀疏典型相關(guān)分析中的懲罰項均不會保留具有相關(guān)性的多組群組特征和組內(nèi)的多個特征,只會選擇出和精神類疾病相關(guān)的一種風險腦區(qū)或者一條易感基因,這對于未知樣本的特征選擇預(yù)測問題的效果易產(chǎn)生影響。考慮到此問題,本文引入exclusive group lasso的懲罰項l1,2范數(shù)并將其與l2,1范數(shù)組合共同懲罰典型相關(guān)分析,可以同時選擇出具有相關(guān)性的組間特征和組內(nèi)特征,構(gòu)造新的目標函數(shù)如下:
s.t.uTXTXu=1vTXTXv=1
(4)
1.4.1 模型的求解以及算法設(shè)計
稀疏典型相關(guān)分析和組稀疏典型相關(guān)分析模型均應(yīng)用l1稀疏約束,此類正則優(yōu)化問題都會假設(shè)特征之間正交,即XTX=I,YTY=I[12]。然而,影像遺傳學(xué)的數(shù)據(jù)往往不具有正交的特性,因此此假設(shè)將會限制識別有關(guān)聯(lián)信息的特征??紤]到此問題,本文采用拉格朗日乘子法,使用交替最小二乘法求解此優(yōu)化問題,使用拉格朗日數(shù)乘法的具體過程如下。
首先構(gòu)造拉格朗日方程
(5)
由于本文模型中使用的懲罰項中含有l(wèi)1范數(shù),如果|ui|=0,|vj|=0,則目標函數(shù)在0點處不可微,可以通過分別給ui,vj加上一個很小的正數(shù)η來改善。然后利用極值存在的必要條件得到:
(6)
(7)
即
(8)
(9)
最后求解上式得到u,v,即
(10)
(11)
輸入:X=[x1,x2,...,xn]T,Y=[y1,y2,...,yn]T
輸出:典型向量uandv
1:初始化:t=0,ut∈Rp×1,vt∈Rq×1,D1,F1;
2:根據(jù)數(shù)據(jù)特點將數(shù)據(jù)進行分組,設(shè)置組信息;
3:While未達到終止條件do
4:根據(jù)公式(10)更新ut+1;
6:end
7:While未達到終止條件do
8:根據(jù)公式(11)更新vt+1;
10:end
11:更新迭代次數(shù)t=t+1;
12:end while
為了驗證GJSCCA的正確性,本文生成一組模擬數(shù)據(jù)集進行測試[13],設(shè)置X∈R80×100和Y∈R80×120,具體包含以下5個步驟:1)隨機產(chǎn)生稀疏典型向量u=(u1,u2,…,up)和v=(v1,v2,…,vq);2)生成潛在變量d~N(0,In×n),d=(d1,d2,…,dn)T保證特征變量和組特征變量之間具有相關(guān)性;3)生成滿足xi~N(ξiu,Σx)和(Σx)pl=exp-|up-ul|的數(shù)據(jù)集X=(x1,x2,…,xp),其中uρ和ul分別為u的第ρ個和第l個坐標;4)生成滿足yi~N(ξiv,Σy)和(Σy)pl=exp-|vp-vl|的數(shù)據(jù)集Y=(y1,y2,…,yq),其中vρ和vl分別為v的第ρ和l個坐標;5)根據(jù)產(chǎn)生的典型變量進行分組,主要分為3種不同類型的組別,第1類組特征的權(quán)重部分為0,部分不為0,第2類組特征的權(quán)重均為0,第3類組特征的權(quán)重均不為0。具體分法如表1所示。
表1 模擬數(shù)據(jù)的分組
從公式(10)和公式(11)中可以看到有4個可調(diào)參數(shù),分別是λ1、λ2、α、β,采用基于CCA觀察的策略,根據(jù)文獻[14],在一個合適的范圍內(nèi)調(diào)整參數(shù),可從[10-2,10-1,100,101,102]中產(chǎn)生最優(yōu)的參數(shù),所有的參數(shù)都通過五折交叉驗證進行調(diào)整。
(12)
其中,Xtest和Ytest表示測試集,u和v表示從訓(xùn)練集得到的典型向量或者特征權(quán)重,本文選擇滿足CV(λ,β,α)的參數(shù)作為最優(yōu)參數(shù)。
為了說明本文模型的有效性,模擬實驗分別從特征選擇的準確性和數(shù)據(jù)集相關(guān)系數(shù)的準確性2個方面與SCCA和GSCCA進行對比評價。由于5折交叉驗證的5次實驗中訓(xùn)練集與測試集的結(jié)果存在一定的差異,因此本文最終的結(jié)果由訓(xùn)練集與測試集所得的相關(guān)系數(shù)之差最小的一組來決定,得到的結(jié)果具體如圖1所示。
圖1 不同方法下典型向量的估計值
圖1分別給出典型向量u和v的真實值和估計值,橫軸代表特征的索引,縱軸代表特征的權(quán)重,權(quán)重越大表明該特征越重要,第1行表示典型向量u和v的真實值,第2行至第4行分別表示SCCA、GSCCA和GJSCCA模型計算出的典型向量的估計值。
從第1行看出,重要的特征被依次分布在10
表2給出3種方法下5折交叉驗證后訓(xùn)練集與測試集的相關(guān)系數(shù)的均值及其平均估計誤差。模擬數(shù)據(jù)集的真實相關(guān)系數(shù)為0.79,表中可看出,相比于SCCA和GSCCA,訓(xùn)練集和測試集上的GJSCCA計算出的平均估計誤差最小,這表明本文方法計算出的相關(guān)系數(shù)更接近真實的相關(guān)系數(shù),優(yōu)于SCCA和GSCCA。
表2 五折交叉驗證下不同方法估計出的典型相關(guān)系數(shù)
精神分裂癥數(shù)據(jù)集來源于MCIC聯(lián)盟(Mind Clinical Imaging Consortium)數(shù)據(jù)庫中的208個樣本,其中包括92名精神分裂癥患者和116名正常人,共含有41236個FMRI體素和777365個SNP位點。由于實驗的目的是找出與精神分裂癥相關(guān)的易感基因和風險腦區(qū),因此本文使用p_value<0.5刪除精神分裂癥患者和正常人之間沒有關(guān)系的生物標志物,最終得到36900個SNP位點和2741個FMRI體素,并且將位點位于同一基因的視為一組,體素位于同一腦區(qū)的視為一組。
表3分別給出3種方法選中的風險基因,經(jīng)過生物醫(yī)學(xué)信息系統(tǒng)pubmed的檢索,SCCA成功選中一條與精神分裂癥相關(guān)的易感基因PLA2G6,據(jù)報道男性吸煙中的尼古丁有害物質(zhì)可引起該基因多態(tài)性發(fā)生病變[15]。GSCCA成功選中4條與精神分裂癥相關(guān)的易感基因,其中APOE與精神分裂癥和情緒障礙有關(guān)[16],DRD2中的rs6276位點與精神分裂癥遺傳易感性患者的缺陷綜合征之間存在潛在聯(lián)系[17],PLXNA2配體信號素在精神分裂癥患者的小腦中被顯示為上調(diào)[18],TPH1基因內(nèi)的序列變異與精神分裂癥有關(guān)[19]。GJSCCA選中7條與精神分裂癥相關(guān)的易感基因,其中APOE是GSCCA和GJSCCA共同選中的易感基因,DAO作為精神分裂癥相關(guān)神經(jīng)遞質(zhì)的受體,被認為是精神分裂癥的潛在候選基因[20],DRD4的-120-bp重復(fù)多態(tài)性與精神分裂癥相關(guān)[21],DTNBP1是被認為在調(diào)節(jié)谷氨酸能系統(tǒng)中起關(guān)鍵作用的基因之一,它的變異導(dǎo)致對精神分裂癥的易感性[22],GRIN2B會導(dǎo)致精神病患者的谷氨酸能功能障礙[23],MTHFR C677T多態(tài)性是精神分裂癥的一個危險因素[24],IL1B有助于精神分裂癥的遺傳風險[25]。對比3種方法,本文方法可選擇出更多的與精神分裂癥相關(guān)的易感基因,證實本文方法的有效性。
表3 不同方法選中的基因
表4給出3種方法選中的與精神分裂癥相關(guān)的風險腦區(qū),SCCA選中2個風險腦區(qū),分別是腦島和丘腦,研究表明精神分裂癥患者腦島的皮質(zhì)厚度減少[26],丘腦的腦灰質(zhì)體積明顯減少[27]。GSCCA選中2個風險腦區(qū),分別是海馬和顳中回,其中海馬體的損傷與神經(jīng)發(fā)育有直接的聯(lián)系[28],精神分裂癥患者的學(xué)習(xí)和記憶選擇性損害與內(nèi)側(cè)顳葉功能障礙一致[29]。GJSCCA選中7個風險腦區(qū),分別是海馬、丘腦、額中回、顳上回、左頂葉、杏仁核和小腦,其中丘腦是SCCA和GJSCCA共同選中的風險腦區(qū),海馬體是GSCCA和GJSCCA共同選中的風險腦區(qū),據(jù)研究,精神分裂癥患者在情緒處理的過程中小腦活動較弱[30],左頂葉結(jié)構(gòu)連通性較降[31],額葉的體積發(fā)生變化[32],杏仁核體積減少[33],左顳上回較小[34]。對比3種方法,本文方法可選擇出更多的與精神分裂癥相關(guān)的風險腦區(qū),再次證實本文方法的有效性。
表4 不同方法選中的腦區(qū)
為了更直觀地觀察風險腦區(qū)在大腦皮層的具體位置,采用BrainNetViewer繪制3種方法的腦區(qū)全方位可視化圖如圖2~圖4。3種方法選中的與精神分裂癥相關(guān)的風險腦區(qū)以不同的灰度被標注在相應(yīng)位置。根據(jù)AAL116圖譜,腦區(qū)被分為116個腦區(qū),其中前90個腦區(qū)均屬于大腦結(jié)構(gòu),后26個腦區(qū)均屬于小腦結(jié)構(gòu),從3種方法可視化圖的對比清晰看出,一方面,SCCA和GSCCA選中的風險腦區(qū)主要集中在大腦,GJSCCCA選中的風險腦區(qū)在大腦和小腦中均有體現(xiàn),而小腦主要控制肢體等行為活動,精神分裂癥患者一般會出現(xiàn)遲發(fā)性運動障礙。另一方面,SCCA選中的風險腦區(qū)主要分布在左側(cè)大腦半球和右側(cè)大腦半球,GSCCA和GJSCCCA選中的風險腦區(qū)均集中在左側(cè)大腦半球,但是GJSCCCA選中的風險腦區(qū)更多。
圖2 SCCA選中的風險腦區(qū)可視化圖
圖3 GSCCA選中的風險腦區(qū)可視化圖
圖4 GJSCCA選中的風險腦區(qū)可視化圖
本文針對具有群組結(jié)構(gòu)的影像遺傳學(xué)數(shù)據(jù)提出了一種在大量基因和腦影像數(shù)據(jù)之間找到潛在相關(guān)特征的方法,利用稀疏正則化項的稀疏和平滑作用,通過交替最小二乘法選擇出潛在特征,在仿真數(shù)據(jù)和精神分裂癥真實數(shù)據(jù)集上進行實驗,仿真實驗表明此模型不但提高數(shù)據(jù)間相關(guān)系數(shù)的準確性,而且保證重要的組間特征和組內(nèi)特征選取的相關(guān)性,真實數(shù)據(jù)表明本文模型可以選擇出更多的與精神分裂癥的易感基因和風險腦區(qū)。