趙迎利，朱 旭

(1.西安交通大學(xué)城市學(xué)院數(shù)學(xué)教研室,陜西 西安 710018; 2.西安交通大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院,陜西 西安 710049)

0 引 言

近年來，影像遺傳學(xué)(Imaging Genetics)成為一種新興的交叉學(xué)科，其主要是結(jié)合影像學(xué)和遺傳學(xué)方法檢測腦結(jié)構(gòu)與神經(jīng)精神疾病、認(rèn)知和情緒調(diào)節(jié)等行為相關(guān)腦功能的遺傳變異研究[1]。醫(yī)學(xué)影像的功能性磁共振成像(Functional Magnetic Resonance Imaging, FMRI)數(shù)據(jù)和遺傳學(xué)的單核苷酸多態(tài)性(Single Nucleotide Polymorphisms, SNP)數(shù)據(jù)提供了精神類疾病不同但互補(bǔ)的信息，聯(lián)合并挖掘這2種類型的數(shù)據(jù)對于精神類疾病的探究可以提供實質(zhì)性的幫助，然而此類數(shù)據(jù)往往具有樣本較少和維度較高的特點，因此，從大量的FMRI數(shù)據(jù)和SNP數(shù)據(jù)中聯(lián)合分析出與精神類疾病相關(guān)的生物標(biāo)記物成為眾多研究學(xué)者面對的主要難題。

典型相關(guān)分析[2](Canonical Correlation Analysis, CCA)和偏最小二乘[3](Partial Least Squares, PLS)被提出來分析2組數(shù)據(jù)間的關(guān)系，其主要是通過最大化2組數(shù)據(jù)集變量間的線性組合之間的相關(guān)性進(jìn)行關(guān)聯(lián)分析，但是由于SNP的小樣本高維度問題，即生物標(biāo)志物的數(shù)量遠(yuǎn)大于樣本數(shù)量，這種問題通常會導(dǎo)致多重共線性問題和過擬合的問題出現(xiàn)[4]。為了解決該問題，眾多研究學(xué)者致力于使用lasso、group lasso等以l1范數(shù)和l2,1范數(shù)為稀疏懲罰項的正則化方法進(jìn)行高維基因表達(dá)數(shù)據(jù)的稀疏，例如Parkhomenko等人[4-6]應(yīng)用lasso懲罰項提出稀疏典型相關(guān)分析(Sparse Canonical Co rrelation Analysis, SCCA)，其主要用來分析雙多變量間的關(guān)系以便從高維的基因表達(dá)數(shù)據(jù)中選擇出與疾病相關(guān)的潛在致病機(jī)制?？紤]到精神類疾病重要的先驗信息通常以感興趣腦區(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)的生物標(biāo)志物。

考慮到上述問題，本文針對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)分析，并保持?jǐ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)的易感基因和風(fēng)險腦區(qū)。

1 數(shù)學(xué)模型

1.1 典型相關(guān)分析

典型相關(guān)分析是一種經(jīng)典的多元統(tǒng)計分析方法，經(jīng)常被用于機(jī)器學(xué)習(xí)領(lǐng)域，其主要是為了2組多維變量X∈Rn×p和Y∈Rn×q分別尋找基向量u∈Rp和v∈Rq，使得隨機(jī)變量uTX和vTY之間的相關(guān)性達(dá)到最大。模型的表示具體如下：

s.t.uTXTXu=1vTYTYv=1

(1)

其中u和v分別表示X和Y的基向量，也被稱為X和Y的典型向量，Xu和Yv被稱為X和Y的典型變量。

1.2 稀疏典型相關(guān)分析

由于現(xiàn)有的基因表達(dá)數(shù)據(jù)通常具有樣本數(shù)目遠(yuǎn)小于特征數(shù)目的特點，直接使用CCA易出現(xiàn)過擬合現(xiàn)象。為了解決此問題，Parkhomenko等人提出了稀疏典型相關(guān)分析，其主要思想是使用lasso的懲罰項，使不重要的特征的權(quán)重很小或者為0，重要特征的權(quán)重較大，從而稀疏掉不重要的特征，以便從大量的特征中找到與精神類疾病相關(guān)的潛在致病機(jī)制，具體的模型表示如下：

s.t.‖u‖2≤1,‖v‖2≤1,‖u‖1≤c1,‖v‖1≤c2

(2)

其中，‖u‖1和‖v‖1是lasso懲罰，c1和c2控制典型向量u和v的稀疏程度。

1.3 組稀疏典型相關(guān)分析

神經(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)

1.4 組聯(lián)合稀疏典型相關(guān)分析模型的提出

影像學(xué)中，腦區(qū)之間存在著一定的聯(lián)系，比如杏仁核活性一旦升高，下丘腦、海馬、前扣帶回以及前額葉皮層等都會得到相應(yīng)的反饋。遺傳學(xué)上，具有相似功能的SNP位點通常位于同一基因，基因組中的基因之間也是存在功能聯(lián)系的，而組稀疏典型相關(guān)分析中的懲罰項均不會保留具有相關(guān)性的多組群組特征和組內(nèi)的多個特征，只會選擇出和精神類疾病相關(guān)的一種風(fēng)險腦區(qū)或者一條易感基因，這對于未知樣本的特征選擇預(yù)測問題的效果易產(chǎn)生影響?？紤]到此問題，本文引入exclusive group lasso的懲罰項l1,2范數(shù)并將其與l2,1范數(shù)組合共同懲罰典型相關(guān)分析，可以同時選擇出具有相關(guān)性的組間特征和組內(nèi)特征，構(gòu)造新的目標(biāo)函數(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，則目標(biāo)函數(shù)在0點處不可微，可以通過分別給ui,vj加上一個很小的正數(shù)η來改善。然后利用極值存在的必要條件得到：

(6)

(7)

即

(8)

(9)

最后求解上式得到u,v，即

(10)

(11)

算法1 模型算法的偽代碼

輸入: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ù)進(jìn)行分組，設(shè)置組信息；

3:While未達(dá)到終止條件do

4:根據(jù)公式(10)更新ut+1；

6:end

7:While未達(dá)到終止條件do

8:根據(jù)公式(11)更新vt+1;

10:end

11:更新迭代次數(shù)t=t+1;

12:end while

2 模擬數(shù)據(jù)的實驗分析

2.1 模擬數(shù)據(jù)的構(gòu)造

為了驗證GJSCCA的正確性，本文生成一組模擬數(shù)據(jù)集進(jìn)行測試[13]，設(shè)置X∈R80×100和Y∈R80×120，具體包含以下5個步驟：1)隨機(jī)產(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個坐標(biāo)；4)生成滿足yi～N(ξiv,Σy)和(Σy)pl=exp-|vp-vl|的數(shù)據(jù)集Y=(y1,y2,…,yq)，其中vρ和vl分別為v的第ρ和l個坐標(biāo)；5)根據(jù)產(chǎn)生的典型變量進(jìn)行分組，主要分為3種不同類型的組別，第1類組特征的權(quán)重部分為0，部分不為0，第2類組特征的權(quán)重均為0，第3類組特征的權(quán)重均不為0。具體分法如表1所示。

表1 模擬數(shù)據(jù)的分組

2.2 參數(shù)設(shè)置

從公式(10)和公式(11)中可以看到有4個可調(diào)參數(shù)，分別是λ1、λ2、α、β，采用基于CCA觀察的策略，根據(jù)文獻(xiàn)[14],在一個合適的范圍內(nèi)調(diào)整參數(shù)，可從[10-2,10-1,100,101,102]中產(chǎn)生最優(yōu)的參數(shù)，所有的參數(shù)都通過五折交叉驗證進(jìn)行調(diào)整。

(12)

其中，Xtest和Ytest表示測試集，u和v表示從訓(xùn)練集得到的典型向量或者特征權(quán)重，本文選擇滿足CV(λ,β,α)的參數(shù)作為最優(yōu)參數(shù)。

2.3 模擬實驗結(jié)果與分析

為了說明本文模型的有效性，模擬實驗分別從特征選擇的準(zhǔn)確性和數(shù)據(jù)集相關(guān)系數(shù)的準(zhǔn)確性2個方面與SCCA和GSCCA進(jìn)行對比評價。由于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ù)

3 精神分裂癥數(shù)據(jù)集上的實驗結(jié)果

精神分裂癥數(shù)據(jù)集來源于MCIC聯(lián)盟(Mind Clinical Imaging Consortium)數(shù)據(jù)庫中的208個樣本，其中包括92名精神分裂癥患者和116名正常人，共含有41236個FMRI體素和777365個SNP位點。由于實驗的目的是找出與精神分裂癥相關(guān)的易感基因和風(fēng)險腦區(qū)，因此本文使用p_value<0.5刪除精神分裂癥患者和正常人之間沒有關(guān)系的生物標(biāo)志物，最終得到36900個SNP位點和2741個FMRI體素，并且將位點位于同一基因的視為一組，體素位于同一腦區(qū)的視為一組。

3.1 易感基因的選擇

表3分別給出3種方法選中的風(fēng)險基因，經(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ì)的受體，被認(rèn)為是精神分裂癥的潛在候選基因[20]，DRD4的-120-bp重復(fù)多態(tài)性與精神分裂癥相關(guān)[21]，DTNBP1是被認(rèn)為在調(diào)節(jié)谷氨酸能系統(tǒng)中起關(guān)鍵作用的基因之一，它的變異導(dǎo)致對精神分裂癥的易感性[22]，GRIN2B會導(dǎo)致精神病患者的谷氨酸能功能障礙[23]，MTHFR C677T多態(tài)性是精神分裂癥的一個危險因素[24]，IL1B有助于精神分裂癥的遺傳風(fēng)險[25]。對比3種方法，本文方法可選擇出更多的與精神分裂癥相關(guān)的易感基因，證實本文方法的有效性。

表3 不同方法選中的基因

3.2 風(fēng)險腦區(qū)的選擇

表4給出3種方法選中的與精神分裂癥相關(guān)的風(fēng)險腦區(qū)，SCCA選中2個風(fēng)險腦區(qū)，分別是腦島和丘腦，研究表明精神分裂癥患者腦島的皮質(zhì)厚度減少[26]，丘腦的腦灰質(zhì)體積明顯減少[27]。GSCCA選中2個風(fēng)險腦區(qū)，分別是海馬和顳中回，其中海馬體的損傷與神經(jīng)發(fā)育有直接的聯(lián)系[28]，精神分裂癥患者的學(xué)習(xí)和記憶選擇性損害與內(nèi)側(cè)顳葉功能障礙一致[29]。GJSCCA選中7個風(fēng)險腦區(qū)，分別是海馬、丘腦、額中回、顳上回、左頂葉、杏仁核和小腦，其中丘腦是SCCA和GJSCCA共同選中的風(fēng)險腦區(qū)，海馬體是GSCCA和GJSCCA共同選中的風(fēng)險腦區(qū)，據(jù)研究，精神分裂癥患者在情緒處理的過程中小腦活動較弱[30]，左頂葉結(jié)構(gòu)連通性較降[31]，額葉的體積發(fā)生變化[32]，杏仁核體積減少[33]，左顳上回較小[34]。對比3種方法，本文方法可選擇出更多的與精神分裂癥相關(guān)的風(fēng)險腦區(qū)，再次證實本文方法的有效性。

表4 不同方法選中的腦區(qū)

為了更直觀地觀察風(fēng)險腦區(qū)在大腦皮層的具體位置，采用BrainNetViewer繪制3種方法的腦區(qū)全方位可視化圖如圖2～圖4。3種方法選中的與精神分裂癥相關(guān)的風(fēng)險腦區(qū)以不同的灰度被標(biāo)注在相應(yīng)位置。根據(jù)AAL116圖譜，腦區(qū)被分為116個腦區(qū)，其中前90個腦區(qū)均屬于大腦結(jié)構(gòu)，后26個腦區(qū)均屬于小腦結(jié)構(gòu)，從3種方法可視化圖的對比清晰看出，一方面，SCCA和GSCCA選中的風(fēng)險腦區(qū)主要集中在大腦，GJSCCCA選中的風(fēng)險腦區(qū)在大腦和小腦中均有體現(xiàn)，而小腦主要控制肢體等行為活動，精神分裂癥患者一般會出現(xiàn)遲發(fā)性運(yùn)動障礙。另一方面，SCCA選中的風(fēng)險腦區(qū)主要分布在左側(cè)大腦半球和右側(cè)大腦半球，GSCCA和GJSCCCA選中的風(fēng)險腦區(qū)均集中在左側(cè)大腦半球，但是GJSCCCA選中的風(fēng)險腦區(qū)更多。

圖2 SCCA選中的風(fēng)險腦區(qū)可視化圖

圖3 GSCCA選中的風(fēng)險腦區(qū)可視化圖

圖4 GJSCCA選中的風(fēng)險腦區(qū)可視化圖

4 結(jié)束語

本文針對具有群組結(jié)構(gòu)的影像遺傳學(xué)數(shù)據(jù)提出了一種在大量基因和腦影像數(shù)據(jù)之間找到潛在相關(guān)特征的方法，利用稀疏正則化項的稀疏和平滑作用，通過交替最小二乘法選擇出潛在特征，在仿真數(shù)據(jù)和精神分裂癥真實數(shù)據(jù)集上進(jìn)行實驗，仿真實驗表明此模型不但提高數(shù)據(jù)間相關(guān)系數(shù)的準(zhǔn)確性，而且保證重要的組間特征和組內(nèi)特征選取的相關(guān)性，真實數(shù)據(jù)表明本文模型可以選擇出更多的與精神分裂癥的易感基因和風(fēng)險腦區(qū)。