基于融合模型的遺傳位點分析?

張繼榮 寇 磊

（西安郵電大學通信與信息工程學院 西安 710121）

1 引言

遺傳疾病致病位點或致病基因的提早發(fā)現(xiàn)，對于疾病的個性化、精準化治療［1］有著非常重要的意義。目前全基因組關(guān)聯(lián)分析，主要是借助于SNP遺傳標記，進行總體的關(guān)聯(lián)性分析，在全基因組范圍內(nèi)選擇遺傳性狀進行基因分類，比較異常性狀和對照組之間每個遺傳頻率的差異，總體統(tǒng)計分析每個異常與目標性狀之間的關(guān)聯(lián)性的大小，選出最具有相關(guān)性的遺傳變異進行驗證，并根據(jù)驗證結(jié)果最終確認相關(guān)性是否存在。但是，在分析過程中進行了多重假設(shè)檢驗，導致假陽性的關(guān)聯(lián)，降低定位的準確性，對進一步鑒別基因甚至是后續(xù)的個性化治療造成相當巨大的偏差［2］。因此，使用精準的方法提高關(guān)聯(lián)性檢驗的準確性［3］，降低甚至消除因多重檢驗造成的假陽性對全基因組關(guān)聯(lián)分析有著至關(guān)重要的作用。

2 模型的建立和求解

2.1 模型假設(shè)

假設(shè)一:含有致病位點的基因?qū)儆谥虏』蚍懂牐?/p>

假設(shè)二:性狀不受環(huán)境影響；

假設(shè)三:1000個樣本來自同一個群體；

假設(shè)四:提供的樣本為無關(guān)的個體樣本。

2.2 堿基編碼方式

DNA是由分別帶有A，T，C，G四種堿基的脫氧核苷酸鏈接組成的雙螺旋長鏈分子。所以堿基編碼形式有n2（種），但是因為SNP關(guān)聯(lián)分析其中有是會相互重復的，只需要保留一種［4］。所以建立數(shù)學模型:

使用直接替換的方法把原始數(shù)據(jù)中每個位點的堿基（A，T，C，G）編碼方式轉(zhuǎn)換數(shù)值編碼方式。具體對照碼表如表1所示。

表1 堿基編碼轉(zhuǎn)化數(shù)值編碼

2.3 SNP關(guān)聯(lián)分析

利用SNP位點集搜索算法——二階段蟻群尋找SNP位點集與遺傳疾病之間的聯(lián)系。此種方法使用位點集與類標間的卡方統(tǒng)計量作為評價函數(shù)，通過蟻群算法在可行解的空間中搜索。二階段蟻群算法的優(yōu)點是對邊緣效應有高的魯棒性，而且采用兩輪蟻群算法，對這兩輪的迭代過程分別設(shè)置參數(shù)［3］。

算法流程如圖1所示。

二階段蟻群算法需要指定所用螞蟻數(shù)目iAntCount，蒸發(fā)率 ρ ，初始信息素 τ0，尋找SNP子集中位點的數(shù)目n，第一輪搜索位點數(shù)目為N1，第二輪搜索位點數(shù)目為N2以及每輪迭代次數(shù)iter_max，其中 N1＞N2＞n。在一輪迭代中，每只人工螞蟻會按順序選擇個位點集，這組位點便是一個可行解。假設(shè)共有L個位點。某只螞蟻在決策時是否選擇位點k的概率公式為

其中 τ（ki）為第i輪搜索第k個位點上的信息素［7］。但是需要注意的是，本算法中ηkβ=1，即螞蟻在每一次選擇新解時只會受信息素的影響，而忽略評價函數(shù)［8］。這么做是因為如果 ηkβ≠1，那么在每只螞蟻選擇可行解時，都需要重新計算所有位點加入后的卡方統(tǒng)計量，運算復雜度過高，影響算法效率，占用大量內(nèi)存。在每一輪迭代結(jié)束后，每只螞蟻都從所有K個位點中選出了一個SNP位點集，每組中含有n個位點。設(shè)螞蟻m選擇了SNP位點集Sm，則Sm中的所有位點k上的信息素會被更新，更新公式如下:

其中Δτk(i)為位點集Sm與疾病之間的0.1?χ2。若某位點k未在任意Sm中出現(xiàn)，則Δτk(i)=0。

圖1 算法流程圖

二階段蟻群算法采用了兩個階段［2］的搜索策略。因為兩輪搜索后的計算結(jié)果中可能包含有相同的SNP位點，所以二階段蟻群算法對結(jié)果采用了最小化假陽率（Minimizing false positive）來進行處理。用EIall表示最小化假陽率前的結(jié)果，EIall初始化為兩輪搜索后的結(jié)果，用EIm表示最小化假陽率后的結(jié)果，EIm初始化為空。對EIall中的每一個元素Ii執(zhí)行以下操作:若EIm中沒有元素與Ii含有相同的SNP位點，則將移動 Ii到 EIm中:若EIm中有元素Ji與Ii含有至少一個相同的SNP位點，則比較Ji與Ii，取 ρ-Value值較大的元素保留在其中，將ρValue值較小的元素丟棄。

對于二階段蟻群的偽代碼為

二階段蟻群算法:

輸入

1 for（第一輪搜索子集中的SNP位點數(shù)目，第二輪搜索子集中的SNP位點數(shù)目）

2 if setsize==第一輪搜索子集中的位點數(shù)目

3 iItCount=第一輪蟻群算法迭代次數(shù)；

4 else iItCount=第二輪蟻群算法迭代次數(shù)；

5 設(shè)置每個位點的初始信息水平

6 i=0，開始進行迭代

7 根據(jù)式（1）為每只螞蟻選擇一個含有setsize序列個位點的SNP子集，

8 計算每只螞蟻選出的SNP子集的卡方檢測統(tǒng)計量，

9 更新SNP子集中所有的位點上的信息素，

10 記錄卡方統(tǒng)計量最高的SNP子集，清除螞蟻

11 i=i+1

12 if i=iItCount停止迭代

后期處理:窮舉搜索 χ2值最大的前iTopModel個子集和前iTopLoci個最大信息素子集，選出所有含有iEpiModl個位點而且閾值小于設(shè)定值的子集。

可選處理:最大限度地減少誤報

返回結(jié)果:以p-Values報告所有檢測到的上位相互作用。

注:其中iItCount為人工螞蟻數(shù)；Largesetsie，Smallsetsize為兩輪搜索子集中SNP位點數(shù)；iTopModel為SNP候選子集數(shù)目；iTopLoci為SNP候選位點數(shù)目；iEpiModl為SNP子集中位點數(shù)目。

2.4 Logistic回歸分析

問題以位點集合為單位分別與疾病性狀做關(guān)聯(lián)分析，既每個單核苷酸多態(tài)性集合（Single Nucleotide Polymorphisms set，SNPs）與性狀均需建立一次模型，以二階段蟻群算法構(gòu)建的SNP關(guān)聯(lián)分析模型為基礎(chǔ)，分析數(shù)據(jù)屬性可以得到以下結(jié)論:

1）反應變量為二分類的分類變量或是某事件的發(fā)生率。

2）自變量與Logit（π）之間為線性關(guān)系。

3）殘差合計為0，且服從二項分布。

4）各觀測值間相互獨立。

Logistic回歸模型可以很好地滿足對此類數(shù)據(jù)的建模需求。所以本文選用Logistic回歸方法，并分別檢測每個SNPs的顯著性。

本文中探討遺傳疾病的危險因素，將選擇兩組人群，一組是患病組，一組是非患病組，兩組人群肯定有不同的體征和生活方式等。這里的因變量就是是否患病，即“是”或“否”（數(shù)據(jù)表示為0或1），為兩分類變量，自變量就是基因包含的SNPs位點。自變量可以是連續(xù)的，也可以是分類的。通過Logistic回歸分析，就可以了解哪些因素是患病的危險因素。

邏輯回歸中的數(shù)據(jù)結(jié)構(gòu)如表2所示，其中 y表示的是否患病，Xp表示基因中的位點信息。

表2 Logistic回歸模型的數(shù)據(jù)結(jié)構(gòu)

參照線性回歸方程，建立下面形式的回歸模型:

對分類變量直接擬合，實質(zhì)上擬合的是發(fā)生概率，顯然，該模型可以描述當各自變量變化時，因變量的發(fā)生概率會怎樣變化，可以滿足分析的基本要求。

令y=1為患病，y=0為未患病，且將患病的幾率記為P，它與自變量x1，x2…，xp之間的Logistic回歸模型為

由此可知未患病的概率為

經(jīng)數(shù)學變化可以得到:

定義:

為Logistic變換，即:

2.5 模型求解

使用二階段蟻群算法SNP關(guān)聯(lián)分析模型搜索計算得到的致病位點，如表3所示。

表3 與疾病最有可能的致病位點

由二階段蟻群算法SNP關(guān)聯(lián)分析模型輸出表3，通過查找的 ρ-Value≤0.01［13］致病位點有10個，依據(jù)此10個致病位點查找盡300個基因，找出包含這些致病位點的基因。這10個基因分別是gene_102、gene_98、gene_218、gene_191、gene_260、gene_149、gene_24、gene_272，gene_84、gene_125。

計算有SNP位點構(gòu)成的基因矩陣的向量間的相似度，選取與致病位點相似度最大的位點，并構(gòu)成新的位點集，以該位點集來表示包含致病位點的基因，找出的位點集總共42個。通過二元邏輯回歸模型求出基因與疾病的關(guān)聯(lián)性。二元邏輯回歸函數(shù)為

其中x為各基因為包含的位點集合，用作自變量，yθ(x)表示患病的概率。

表4 疾病與基因的關(guān)聯(lián)性

3 模型檢測

3.1 隨機森林做模型檢測

為驗證二階段蟻群算法檢測到的致病位點的有效性，將這些致病位點在數(shù)據(jù)樣本集（1000×9445）中提取出來作為新的驗證樣本集（1000×10）。由于是關(guān)于致病位點與患病與否之間的關(guān)系，因此驗證致病位點有效性可以視為分類問題，本文選用隨機森林算法［8］進行驗證。

隨機森林［9］是一種由多個決策樹組合而成的分類器，它利用bootstrap重抽樣方法從原始樣本中抽取多個樣本，對每個bootstrap樣本進行決策樹建模，然后將這些決策樹組合在一起，通過投票或輸出的平均值得出最終分類或預測的結(jié)果。

將驗證樣本集（1000×10）以7:3的比例隨機劃分為樣本集和測試集，利用隨機森林算法進行識別分類，通過分析實驗結(jié)果，平均識別正確率達到82.15%，驗證了該方法的有效性。實驗結(jié)果如表5所示。

表5 隨機森林檢測結(jié)果

3.2 SPSS軟件做模型檢測

對回歸模型識別率下降排序，基因的排序為gene_98、gene_218、gene_102、gene_191、gene_260、gene_149、gene_24、gene_272，選取前六個基因為與疾病可能相關(guān)的基因集合。

利用SPSS軟件對邏輯回歸模型進行檢測［11］，將基因102、98、91、281、260、149的位點新集合輸入SPSS軟件中，輸出一系列的模型參數(shù)，采用模型系數(shù)混合檢驗作為檢測參數(shù)。

表6 模型系數(shù)的綜合檢測

模型系數(shù)的混合檢驗（Omnibus Tests ofModel Coefficients）是針對步驟、模塊和模型開展模型系數(shù)的綜合檢驗（如表6所示）。表6中給出卡方值及其相應的自由度、P值及Sig.值。模型卡方臨界值的計算步驟如下:取顯著性水平 Sig.為 0.05［13］，自由度數(shù)目為表中的df，通過查表法，可以查出卡方臨界值。當表中的卡方值小于卡方臨界值，Sig.值大于0.05時，模型通過檢驗。

表7 卡方臨界值

通過計算以上各基因的卡方臨界值得到表7，各基因的模型系數(shù)綜合檢驗的卡方值均小于各自基因的卡方臨界值，且各基因的顯著性水平大于0.05，因此各基因的模型通過檢驗，驗證了二元邏輯回歸模型尋出與疾病最相關(guān)基因的正確性。

4 模型評價分析

為了驗證模型的有效性，在本文模型的基礎(chǔ)上和文獻［7］中的傳統(tǒng)蟻群模型與文獻［12］中的Logistic回歸模型進行對比。采用庫都為國際標準pubmed庫。模型評價采用通用準確率P評價分析［12］

r代表檢索結(jié)果中識別正確的數(shù)目，u代表檢索結(jié)果中識別錯誤的數(shù)目。

為了驗證模型的有效性和算法的合理性，在上述庫中用上述模型在不同數(shù)量樣本下進行識別，其準確率如圖2所示。

圖2 五種樣本數(shù)目下的準確率

由上述圖像對比可知，相對于傳統(tǒng)的遺傳性疾病和性狀的遺傳位點分析模型，在大樣本情況下，尤其是樣本數(shù)目超過600個的時候，融合模型的識別準確率有顯著優(yōu)勢。

5 結(jié)語

本文提出了一種基于融合模型的遺傳性疾病和性狀的遺傳位點分析模型，這個模型針對基因數(shù)據(jù)的關(guān)聯(lián)性特點融合不同的模型進行數(shù)據(jù)分析并綜合運用各種方法得到有效的相關(guān)性分析。本文使用邏輯判斷和基于條件判斷的數(shù)據(jù)轉(zhuǎn)換算法進行編碼方式轉(zhuǎn)變，在面對遺傳疾病與位點的相關(guān)性，給出二階段蟻群算法模型，得到了理想的關(guān)聯(lián)性結(jié)果。在面對疾病、基因和位點多層次相關(guān)問題，通過分析之前的結(jié)果，使用邏輯回歸模型中的Logistic回歸模型，通過隨機森林算法和卡方檢測測試驗證，大大加強了融合模型的可靠性。在文中所運用的二階段蟻群算法是對基礎(chǔ)蟻群算法的優(yōu)化，在蟻群算法的信息素更新階段加入了順序更新的規(guī)則，分成兩個階段對螞蟻走過路徑上的信息素進行了更新，當然同時也提出了兩種信息素跟新階段的劃分策略，并且把每段路徑上的信息素的值限定在一個整數(shù)區(qū)間里，算法的收索能力明顯增強。