基于全基因組關(guān)聯(lián)研究的tag SNP 預(yù)測方法的研究與分析*

劉高偉

（中國石油大學(xué)（華東）計算機(jī)與通信工程學(xué)院 青島 266580）

1 引言

全基因組關(guān)聯(lián)研究［1］（Genome-Wide Association Study，GWAS）在全基因組層面上，通過大樣本、多中心、反復(fù)驗(yàn)證的方法來尋找大規(guī)模群體的各種脫氧核糖核酸（Deoxyribonucleic acid，DNA）遺傳標(biāo)記［2］，如單核苷酸多態(tài)性（Single Nucleotide Polymorphism，SNP）［3］或 基 因 拷 貝 數(shù) 變異（Copy Number Variations，CNV）［4］等與復(fù)雜疾病相關(guān)的遺傳因素，從而全面揭示與人類疾病發(fā)生、發(fā)展和治療相關(guān)的遺傳基因。全基因組關(guān)聯(lián)研究中的單核苷酸多態(tài)性（single-nucleotide polymorphism，SNP）是指基因組DNA 序列同一位置單個核苷酸變異所引起的多態(tài)性，其占人體基因組遺傳多態(tài)性的90%以上［5～7］。SNP 集分析方法在識別致病SNP 時有較高的功效，且越來越受到生物醫(yī)學(xué)界的關(guān)注［8～9］。

理論上，在尋找遺傳疾病相關(guān)基因的過程中，是要通過對所有SNP 位點(diǎn)進(jìn)行鑒定，然而，人體基因中一共有300 萬個SNP，對所有SNP 位點(diǎn)進(jìn)行基因分析，這樣將消耗大量成本［10～13］。研究表明，許多SNP 位點(diǎn)之間具有一定的關(guān)聯(lián)關(guān)系，其中部分SNP 位點(diǎn)信息就可以代表所有SNP 位點(diǎn)所攜帶的信息，我們將具有代表性的SNP 成為標(biāo)簽SNP（tag SNP）集［14～16］。通過對tag SNP 集進(jìn)行基因分型，tag SNP 集的質(zhì)量與基因分型的準(zhǔn)確率有密切的關(guān)系，因此tag SNP 集評估方法已經(jīng)成為GWAS 研究領(lǐng)域的熱點(diǎn)研究內(nèi)容［17］。

評估tag SNP 集質(zhì)量是SNP 集分析中的一個重要內(nèi)容，高質(zhì)量的tag SNP 集可以更好節(jié)約基因分型成本，并且提高基因分型的準(zhǔn)確率。2005 年，Halperin 等［18］提出了STAMPA 方法評估tag SNP 集質(zhì)量，該方法在預(yù)測過程中使用了容易獲得的基因型數(shù)據(jù)，而不是單倍體數(shù)據(jù)，并且該方法預(yù)測過程不受單倍體塊劃分限制，對每一個非tag SNP 位點(diǎn)進(jìn)行預(yù)測，以此達(dá)到評估tag SNP集質(zhì)量的目的，該方法在預(yù)測非tag SNP時只考慮了與所需預(yù)測的非tag SNP物理距離最鄰近的tag SNPs，而沒有考慮到SNP 之間的生物距離，這樣也就導(dǎo)致預(yù)測效率不高；Ilhan 等［19］使用支持向量機(jī)（Support Vector Machine，SVM）方法對所選擇的tag SNP 進(jìn)行預(yù)測，通過粒子群優(yōu)化方法對支持向量機(jī)中的參數(shù)進(jìn)行優(yōu)化，使用遺傳算法（Genetic Algorithm，GA）選擇tag SNP，通過每一次迭代對參數(shù)進(jìn)行優(yōu)化，以使得SVM 方法更加有效，該方法在預(yù)測tag SNP集時，通過每一次訓(xùn)練調(diào)整參數(shù)，在預(yù)測時需要所有tag SNPs 位點(diǎn)信息，并且需要迭代的次數(shù)要足夠大，這樣才可以有較好的評估tag SNP 集；Mauawad 等［20］基于SNP 之間的連鎖不平衡度，提出Linkage Coverage 方法評估tag SNP 集質(zhì)量，使用所有tag SNPs預(yù)測non-tag SNP，充分考慮了SNP 之間的關(guān)聯(lián)關(guān)系，然而該方法只計算了tag SNP 集與tag SNPs 集之間的連鎖不平衡，沒有考慮SNPs之間物理距離，這樣也就降低tag SNP集的評估質(zhì)量。

基于上文所提出的問題，本文提出一種基于連鎖不平衡度的Genetic-majority dominance（GMD）tag SNP 預(yù)測方法評估tag SNP 質(zhì)量，GMD 通過使用遺傳方法迭代選擇出tag SNP 集，然后以SNP 之間的連鎖不平衡為基礎(chǔ)，充分考慮了SNP之間的關(guān)聯(lián)關(guān)系以及實(shí)際生物距離，在使用tag SNP 集預(yù)測非tag SNP集時，并不需要使用所有tag SNP，只需要找出與所需預(yù)測的non-tag SNP 集關(guān)系最為鄰近的tag SNPs，所使用的也是易獲得的基因型數(shù)據(jù)，GMD 可以定量分析tag SNP 集質(zhì)量且也保證了評估效果。

2 基于連鎖不平衡度的GMD方法

2.1 連鎖不平衡度

連鎖不平衡是指不同位點(diǎn)上兩個等位基因出現(xiàn)在同一個染色體上的頻率要高于預(yù)期的隨機(jī)頻率的現(xiàn)象［21］。由于Human Leucocyte Antigen（HLA）不同基因位點(diǎn)的某些等位基因經(jīng)常連鎖在一起遺傳，而連鎖的基因并非完全隨機(jī)地組成單體型，有些基因總是較多地在一起出現(xiàn)，致使某些單體型在群體中呈現(xiàn)較高的頻率，從而引起連鎖不平衡。

設(shè)群體有n 個個體，每個個體有p個SNP 位點(diǎn)，即原始SNP 集中有p 個SNPs，記原始SNP 集為L=｛1，2，…，p｝。因?yàn)槿祟愂嵌扼w生物，所以共有2n 個 單 體 型，令zi=(zi1，zi2，…，zip) ，i=1，2，…，2n，Rij表示SNPi與SNPj的連鎖不平衡度，通過下式可得［22］：

2.2 遺傳算法-多數(shù)占優(yōu)方法（GMD）

遺傳算法-多數(shù)占優(yōu)方法（GMD）是將遺傳算法與多數(shù)占優(yōu)原則結(jié)合的一種tag SNP 集預(yù)測方法，該方法使用遺傳算法選擇出tag SNP集，然后以多數(shù)占優(yōu)原則為基礎(chǔ)評估tag SNP集的質(zhì)量。

遺傳算法是一種基于“適者生存”的高度并行、隨機(jī)和自適應(yīng)的優(yōu)化算法，通過復(fù)制、交叉、變異將問題解編碼表示的“染色體”群一代代不斷進(jìn)化，最終收斂到最適應(yīng)的群體，從而求得問題的最優(yōu)解或滿意解［23～24］。在本文中，使用遺傳算法從原始SNP集中選擇出tag SNP 集。使用tag SNP 集對非tag SNP 進(jìn)行預(yù)測時，對樣本基因數(shù)據(jù)進(jìn)行訓(xùn)練，得到數(shù)據(jù)之間的關(guān)系，然后根據(jù)多數(shù)占優(yōu)原則使用所選擇的tag SNPs預(yù)測非tag SNP集。

3 算法描述

3.1 符號定義

設(shè)群體有n 個個體，每個個體有p個SNP 位點(diǎn)，即原始SNP 集中有p 個SNPs，記原始SNP 集為L=｛1，2，…，p｝。設(shè)一個p 維向量gi=｛0，1，2｝p表示一個基因片段，每一個位點(diǎn)上取值為0，1，2。0，1分別表示次等位基因和主等位基因，2 表示雜合子基因，prei表示在i位置上得到的預(yù)測值。每一個基因片段是由兩個單體型組成，單體型中每一個位點(diǎn)上的取值為0 或1。設(shè)gij所表示的是第i 個基因片段上第j個位置上的基因型信息，使用hij1與hij2表示組成該基因的兩個單體型對應(yīng)位點(diǎn)上的單體型值，hijp∈0｛0，1｝，如果hij1=hij2=1，gij=1，如果hij1=hij2=0，gij=0，如果hij1≠hij2，則gij=2。

T 表示tag SNP 集，|T|表示tag SNP 集中元素的個數(shù)，即tag SNP的個數(shù)。

3.2 遺傳算法選擇tag SNP集

3.2.1 染色體編碼方案

使用二元向量表示tag SNP 選擇問題的染色體，使用Ci=｛0，1｝p表示一個染色體，0 表示non-tag SNP，即非tag SNP，1 表示tag SNP。所使用的樣本數(shù)量一共是n 個，所以有Ci=｛ci1，ci2，…，cip｝，cij的取值為0或1，i=1，2，…，n，j=1，2，…，p，設(shè)每一個染色體中有k個tag SNP，即有k個值為1。

3.2.2 選擇

計算每一個染色體的適應(yīng)值，即tag SNP 集預(yù)測non-tag SNP 集的預(yù)測精度值pv（具體預(yù)測方法在下節(jié)介紹），這樣每一個染色體就有一個對應(yīng)的pv，按照pv從小到大進(jìn)行排序，C'=｛C1'，C2'，…，Cp'｝。設(shè) 兩 個 參 數(shù)Wup（Wup＞1）和Wdown（Wdown＜1），d=（Wup-Wdown）/p-1，令Rank（Ci'）=Rank（C1'）+（i-1）*d。設(shè)閾值t0，如果Rank（Ci'）＞t0，則選擇Ci'，且將Rank（Ci'）=Rank（Ci'）-1；在0-1 之間隨機(jī)產(chǎn)生一個分?jǐn)?shù)f1，如果Rank（Ci'）＞f1，則將Ci'選擇，否則不選擇，一直進(jìn)行該操作，一直選擇出n個染色體。

3.2.3 交叉

選擇操作完成之后，進(jìn)行交叉操作，再本文中使用的是兩點(diǎn)交叉，設(shè)交叉率為Xrate，判斷染色體Ci'與Ci+1'之間是否需要進(jìn)行交叉，需要在0～1 之間隨機(jī)產(chǎn)生一個分?jǐn)?shù)f2，如果Xrate＜f2，則不進(jìn)行交叉操作，否則，將染色體Ci與Ci+1進(jìn)行交叉操作，交叉過程為：隨機(jī)產(chǎn)生兩個整數(shù)r1與r2，1 ≤r1≤r2≤p，將染色體Ci與Ci+1的區(qū)間［r1，r2］之間的片段進(jìn)行交叉操作。

3.2.4 突變

突變操作：對于每一個染色體的每一個SNP進(jìn)行突變操作，設(shè)突變率為Mrate，判斷SNP 是否需要突變，在0～1 之間隨機(jī)產(chǎn)生一個分?jǐn)?shù)f3，若f3＞Mrate，則不發(fā)生突變，否則將SNP 位點(diǎn)進(jìn)行突變操作，即該位點(diǎn)信息本來為0，則變?yōu)?，如果是1，則突變?yōu)?。

3.2.5 修復(fù)

事先設(shè)置的tag SNP 個數(shù)為k，在經(jīng)過選擇，交叉，變異操作之后，為了保證tag SNP個數(shù)為k，需要進(jìn)行修復(fù)操作，對于每一個染色體有t 個tag SNP，如果t＞k，則隨機(jī)選擇t-k 個tag SNP 將其值改為0，如果t＜k，則隨機(jī)選擇t-k 個non-tag SNP 的值改為1。

在修復(fù)操作完成之后，使用預(yù)測算法對每一個染色體中的non-tag SNP 進(jìn)行預(yù)測精度。若迭代次數(shù)未達(dá)到總次數(shù)，則繼續(xù)迭代，否則停止。

3.3 tag SNP預(yù)測算法

4 結(jié)果與討論

4.1 數(shù)據(jù)描述

在實(shí)驗(yàn)中，一共使用了四組數(shù)據(jù)集：包括了HTR2A 基因、ASAH1 基因、OLMF4 基因和PHGDH基因數(shù)據(jù)集。

HTR2A 是13 號染色體上與精神分裂癥、強(qiáng)迫癥等疾病有關(guān)的基因，該基因位于13 號染色體位置q14-21，一共有169個SNPs在基因HTR2A上。

ASAH1 是8 號染色體上與法伯脂肪肉芽腫病疾病有關(guān)的基因，該基因位于8 號染色體位置8p22，一共有177個SNPs。

OLMF4 是13 號染色體上與胃癌、胰腺癌等疾病有關(guān)的基因，該基因位于13 號染色體位置q14.3，一共有56個SNPs。

PHGDH 是1 號染色體上與癲癇疾病有關(guān)的基因，該基因位于1 號染色體位置1p12，一共有64 個SNPs在基因PHGDH。

實(shí)驗(yàn)中，對于每一個數(shù)據(jù)集使用HAPGEN2 軟件產(chǎn)生1000 組模擬數(shù)據(jù)，每一組數(shù)據(jù)中有500 個case個體，仿真數(shù)據(jù)的連鎖不平衡結(jié)構(gòu)是基于Hap-Map計劃產(chǎn)生的。

4.2 實(shí)現(xiàn)環(huán)境

在實(shí)現(xiàn)本方法時，所使用的編程語言是R 語言。所運(yùn)行的環(huán)境是在Linux 操作系統(tǒng)下運(yùn)行的，所使用的PC 的設(shè)施要求4GHz的內(nèi)存以及500G 的硬盤存儲空間。

4.3 預(yù)測精度比較

在實(shí)驗(yàn)中，設(shè)置遺傳算法的最大迭代次數(shù)為50，參數(shù)Wup=1.2，Wdown=0.8，交叉率Xrate=0.8，突變率Mrate=0.1，在對HTR2A 和ASAH1 數(shù)據(jù)集進(jìn)行實(shí)驗(yàn)時，tag SNP 的 個 數(shù) 從5 遞 增 到30，對OLMF4 和PHGDH 數(shù)據(jù)集進(jìn)行實(shí)驗(yàn)時，tag SNP 的個數(shù)從3 遞增到10，每一次遞增數(shù)量為1。

在本文中，首先使用遺傳方法選擇出每一次迭代所產(chǎn)生的tag SNP 集，然后通過使用不同tag SNP預(yù)測方法衡量每一次迭代出的tag SNP 集的質(zhì)量，得到對應(yīng)的適應(yīng)值，將本文中所提出的GMD 方法與STAMPA預(yù)測方法進(jìn)行比較，比較結(jié)果圖1所示。

圖1 GMD方法與STAMPA方法分別對HTR2A基因、ASAH1基因、OLMF4基因、PHGDH基因進(jìn)行預(yù)測比較圖

通過圖1 可以看出當(dāng)tag SNP 的個數(shù)比較少時，STAMPA 方法所預(yù)測的精度較小，隨著tag SNP個數(shù)不斷增加，預(yù)測精度也在逐漸提高，且可以看出本文提出的GMD 預(yù)測方法在HTR2A、ASAH1、OLMF4 和PHGDH 數(shù)據(jù)集中，預(yù)測結(jié)果都優(yōu)于STAMPA 方法。STAMPA 方法在預(yù)測non-tag SNP時，只考慮了物理距離最短的tag SNPs，而沒有考慮到所需要預(yù)測的non-tag SNP 與tag SNP 的生物意義，因?yàn)樵诨蚱沃写嬖诙鄠€單體型塊，使用STAMPA 所選擇的tag SNPs 與所需要預(yù)測的non-tag SNP不在同一個單體型塊，這樣tag SNPs與non-tag SNP 之間的關(guān)聯(lián)關(guān)系比較弱，使用不在同一個單體型塊中的tag SNP 預(yù)測non-tag SNP，這樣大大降低了預(yù)測精度，所以該方法不能很好評估tag SNP 集。而本文所提出的預(yù)測方法則是可以充分考慮到了tag SNP 與non-tag SNP 之間的生物距離，而且也保證了與所需要預(yù)測的non-tag SNP 關(guān)聯(lián)關(guān)系最為緊密的tag SNP，這樣在預(yù)測時，提高了預(yù)測精度，可見與STAMPA 方法比較，GMD 方法可以更好去評估所選擇出的tag SNP。

在本文中，也將GMD 方法與Linkage Coverage方法進(jìn)行了比較，因?yàn)長inkage Coverage 方法在評估tag SNP 質(zhì)量時，是通過計算tag SNP 集與non-tag SNP 集之間的連鎖不平衡度之和的覆蓋率評估tag SNP 集質(zhì)量，并沒有對具體SNP 位點(diǎn)進(jìn)行預(yù)測，所以不能與GMD 方法進(jìn)行直接比較，在實(shí)驗(yàn)中，通過將Linkage Coverage 方法與GMD 方法作為評估tag SNP 集方法與遺傳方法選擇tag SNP 集結(jié)合，遺傳方法迭代達(dá)到50 次之后，對所選擇的tag SNP 集使用STAMPA 方法預(yù)測，然后再對所預(yù)測的結(jié)果進(jìn)行比較，比較結(jié)果圖2所示。

圖2 GMD方法與Linkage Coverage方法分別對HTR2A基因、ASAH1基因、OLMF4基因、PHGDH基因進(jìn)行tag SNP 預(yù)測比較圖

通過圖2 可以看出，在HTR2A 和ASAH1 數(shù)據(jù)集中，當(dāng)tag SNP數(shù)量較少時，兩者之間的預(yù)測精度差還是比較少的，但是隨著tag SNP數(shù)量不斷增大，本文所提出的GMD 方法與Linkage Coverage 方法之間的差距越來越大，且一直優(yōu)于Linkage Coverage 方法。在OLMF4 和PHGDH 數(shù)據(jù)集，通過圖2（c）和圖2（d）可以看出本文所提出的方法預(yù)測tag SNP 集都是優(yōu)于Linkage Coverage 方法。因?yàn)镚MD方法不僅考慮了SNPs 之間的連鎖不平衡，而且也使用到了實(shí)際樣本中的基因型數(shù)據(jù)，充分考慮了基因型數(shù)據(jù)之間的關(guān)系，這樣也就使得GMD 的預(yù)測精度比較高。

5 結(jié)語

本文提出了一種新的tag SNP 預(yù)測方法——GMD 預(yù)測方法，該方法通過使用遺傳算法選擇出tag SNP，然后確定tag SNP 集，通過訓(xùn)練基因型數(shù)據(jù)，使用多數(shù)占優(yōu)的原則對non-tag SNP 進(jìn)行預(yù)測，這樣實(shí)現(xiàn)了對tag SNP 的定量分析，將tag SNP 集的質(zhì)量進(jìn)行量化分析，該方法在預(yù)測tag SNP 集時不僅考慮了SNPs 之間的連鎖不平衡度，而且考慮了樣本的基因型數(shù)據(jù)關(guān)系，所以GMD 方法可以比較好的評估tag SNP 集。通過將該GMD 與Linkage Coverage、STAMPA 方法在HTR2A、ASAH1、OLMF4和PHGDH 數(shù)據(jù)集上進(jìn)行比較，實(shí)驗(yàn)結(jié)果顯示GMD相比于其他兩種方法可以更好評估tag SNP 集質(zhì)量。

盡管GMD 方法可以較好評估tag SNP 集質(zhì)量，然而，在評估非tag SNP 集時只是使用了與非tag SNP 關(guān)聯(lián)關(guān)系最強(qiáng)的兩個tag SNPs，沒有考慮其他的tag SNPs，盡管其他tag SNPs影響較小，但是其中可能存在潛在的生物意義，所以在接下來的工作中，我們除了使用與所預(yù)測的非tag SNP 關(guān)聯(lián)關(guān)系最近的tag SNPs，還需要考慮到與需要預(yù)測的非tag SNP有潛在生物意義的tag SNPs。