一種基于連鎖不平衡的tagSNPs選擇算法

王鈞峰，王新贈(zèng)

(山東科技大學(xué)數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院，山東青島266590)

一種基于連鎖不平衡的tagSNPs選擇算法

王鈞峰，王新贈(zèng)

(山東科技大學(xué)數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院，山東青島266590)

進(jìn)行全基因組關(guān)聯(lián)研究(genome－wide association studies，簡(jiǎn)記為GWAS)時(shí)，我們需要獲得一個(gè)足夠密集的單核苷酸多態(tài)性(single－nucleotide polymorphism，簡(jiǎn)記為SNP)標(biāo)記集來(lái)解釋常見(jiàn)疾病遺傳風(fēng)險(xiǎn)的一部分.候選基因中SNP的數(shù)量是有限的，但是直接分析所有現(xiàn)存的SNPs是無(wú)效的，因?yàn)樵谶@些位點(diǎn)上的基因型有很強(qiáng)的關(guān)聯(lián)性，會(huì)導(dǎo)致大量的冗余信息，并且會(huì)造成基因分型成本的增加，消耗大量的時(shí)間.所以我們?cè)谶M(jìn)行關(guān)聯(lián)檢驗(yàn)時(shí)，沒(méi)有必要對(duì)所有的SNPs進(jìn)行基因分型，只需要選擇出具有代表性，并且數(shù)量很少的SNPs進(jìn)行分型，并對(duì)這些SNPs進(jìn)行關(guān)聯(lián)檢驗(yàn).這里選擇出的SNPs稱為標(biāo)簽SNPs(記為tagSNPs)，它為SNPs的一個(gè)小的子集，在每個(gè)單體型區(qū)域中足以捕獲單體型的信息.選擇tagSNPs的方法有很多，本文我們提出了一種新的tagSNPs選擇方法，通過(guò)使用基于連鎖不平衡(linkage disequilibrium，簡(jiǎn)記為L(zhǎng)D)的兩兩r2準(zhǔn)則對(duì)單體型分組，分成不相交的組，并在每個(gè)組中選擇標(biāo)簽SNPs.與基于原始SNPs集的檢驗(yàn)方法比較，我們的方法產(chǎn)生了更少的tagSNPs，在最大化所選標(biāo)記提供信息含量的同時(shí)，降低了基因分型成本，提高了效率.

全基因關(guān)聯(lián)研究;SNP;標(biāo)簽SNP;連鎖不平衡

1 引言

在人類基因組中，由于SNP高的豐富性，低突變率，易于高通量基因分型［1］，所以SNP在疾病關(guān)聯(lián)研究中起到了很重要的作用.TagSNPs的選擇已經(jīng)成為一個(gè)非?；钴S的研究課題.如果能在SNP數(shù)據(jù)集中選擇出tagSNP集，就可以減少用于關(guān)聯(lián)檢驗(yàn)所使用的SNP的數(shù)量，這樣也就降低了基因分型的成本和計(jì)算的時(shí)間復(fù)雜度.目前已經(jīng)被提出的tagSNP選擇方法主要有兩種，分別是基于單體塊結(jié)構(gòu)的識(shí)別方法［2－5］和基于LD的識(shí)別方法［6－8］.很多算法也已經(jīng)被提出來(lái)檢測(cè)單體型塊和標(biāo)記的選擇.Patil［9］等定義一個(gè)單體型塊，單體型塊是一個(gè)區(qū)域，在其中所有觀察的a%的單體型是常見(jiàn)單體型，常見(jiàn)單體型是單體型頻率大于某一閾值的單體型，通過(guò)這個(gè)方法選擇所有可能的單體型塊，然后通過(guò)一個(gè)貪婪優(yōu)化算法選擇出連續(xù)沒(méi)有重疊的單體型塊和tagSNPs.Johnson［2］等是基于連鎖不平衡計(jì)算兩兩SNPs間的連鎖不平衡程度，如果連鎖不平衡程度大于某個(gè)閾值，那么其中一個(gè)就可以作為tagSNP.Zhang［3］等使用一個(gè)動(dòng)態(tài)編程方法來(lái)進(jìn)行單體型塊的劃分和tagSNP的選擇.Cardon［10］的策略是選擇一個(gè)具有代表性的SNPs集，不考慮剩余SNPs，目的是保留大多數(shù)原始集信息的同時(shí)，所選擇的SNPs沒(méi)有信息的重疊.Carlson［8］利用貪婪算法識(shí)別tagSNP，從所有超過(guò)某閾值的SNPs開(kāi)始，與最大數(shù)量SNPs之間的連鎖不平衡程度都大于某個(gè)閾值的SNP作為tagSNP.Zhang［5］等人，介紹基于LD的方法，這些方法搜索一個(gè)小的SNP集，并與其他不被選擇的SNPs有強(qiáng)的連鎖不平衡.Stram［11］等給出了一種統(tǒng)計(jì)方法，在其中多個(gè)tagSNPs可以被用來(lái)代表每個(gè)沒(méi)有被標(biāo)記的SNP.

在應(yīng)用我們提出的方法選擇tagSNPs之后，我們要驗(yàn)證所選擇的tagSNPs是否可以驗(yàn)證與疾病的關(guān)聯(lián)性，所以我們要進(jìn)行假設(shè)檢驗(yàn).研究者已經(jīng)提出了很多檢驗(yàn)方法.CAST［12］是對(duì)于每個(gè)個(gè)體把在一個(gè)區(qū)域內(nèi)(例如，一個(gè)基因的外顯子)所有罕見(jiàn)變體信息重疊成一個(gè)二分變量，通過(guò)判斷個(gè)體是否有任何罕見(jiàn)變體在這區(qū)域內(nèi)，然后運(yùn)用一個(gè)單變量檢驗(yàn)［13］.CMC［14］擴(kuò)展了CAST方法，通過(guò)在一個(gè)等位基因頻率的基礎(chǔ)上，把罕見(jiàn)變體重疊在分組的區(qū)域內(nèi)，重疊所分的組就如CAST方法一樣，并對(duì)所分的組用一個(gè)多元檢驗(yàn).加權(quán)和檢驗(yàn)(WST)［15］考慮一組病例對(duì)照，把一個(gè)SNP集重疊成一個(gè)罕見(jiàn)等位基因數(shù)量的單一的加權(quán)平均.Wu［16］提出了一種基于邏輯核機(jī)器的邏輯核檢驗(yàn)方法(sequence kernel association test，簡(jiǎn)記為SKAT)，SKAT假定標(biāo)記的回歸系數(shù)的一個(gè)分布，其方差取決于靈活的權(quán)重.SKAT執(zhí)行一個(gè)基于得分的方差分量檢驗(yàn)，它的計(jì)算只需要擬合空模型通過(guò)單獨(dú)在協(xié)變量上回歸表型和解析計(jì)算P值.SKAT能夠直接獲得一個(gè)P值而不需要排列求P值.SKAT的一個(gè)重要特性是它允許結(jié)合靈活的加權(quán)函數(shù)來(lái)提高分析功效.所以在進(jìn)行假設(shè)檢驗(yàn)時(shí)，我們使用高效靈活的SKAT.

2 方法

2.1 tagSNP的選擇

考慮包含P個(gè)雙等位基因SNP標(biāo)記a1，a2，…，αp的一個(gè)集.進(jìn)一步假設(shè)所有這些標(biāo)記次要等位基因頻率(MAF)超過(guò)一個(gè)特定的閾值(在這使用0.05).首先，計(jì)算兩兩LD測(cè)量r2［17］.如果兩個(gè)標(biāo)記ai和aj間的r2大于一個(gè)特定的閾值r0，那么就說(shuō)它們兩個(gè)有強(qiáng)的LD，表示為r2(ai，aj)≥r0(在這個(gè)研究中r0=0.8)，兩個(gè)都可以被考慮作為對(duì)方的tagSNPs，在其中ai可以用來(lái)作為aj的一個(gè)替代，反之亦然.

我們的目的是找到一個(gè)tagSNP集，對(duì)于基因分型，我們開(kāi)發(fā)了一種算法來(lái)識(shí)別tagSNPs子集，從超過(guò)一個(gè)給定MAF閾值的所有SNPs中選擇.從超過(guò)MAF閾值的所有SNPs開(kāi)始，對(duì)我們所選擇的基因區(qū)域進(jìn)行分組，分成幾個(gè)SNP子集，SNP集中任意兩個(gè)SNPs同屬于一個(gè)組當(dāng)且僅當(dāng)這兩個(gè)SNPs之間的r2大于等于給定的閾值，也就是說(shuō)同一個(gè)組中的SNPs至少與同組中的一個(gè)SNP連鎖不平衡，分組的過(guò)程是迭代的，每一次循環(huán)分析所有未被分組的SNPs，直到所有的SNPs被分組，這樣，就分成了幾個(gè)組，但是會(huì)出現(xiàn)這樣一種特殊情況，有的SNPs和任何SNPs都不連鎖不平衡，那么我們就把這樣的單個(gè)SNP單獨(dú)作為一個(gè)組.

這樣，組中一個(gè)SNP被指定為“tagSNPs”，每個(gè)組只有一個(gè)tagSNP將會(huì)需要被基因分型.下面，我們就在每個(gè)組中選擇一個(gè)tagSNP，再把每個(gè)組中選擇的tagSNP組成總的tagSNP集，然后進(jìn)行假設(shè)檢驗(yàn).選擇tagSNP的方法具體如下，在一個(gè)已經(jīng)分好的組中進(jìn)行選擇，我們首先計(jì)算出組中最大r2值的兩個(gè)SNPs，然后再?gòu)倪@兩個(gè)SNPs中選擇其中一個(gè)作為tagSNP，分別計(jì)算這兩個(gè)SNPs與其他除去二者本身的組內(nèi)剩余SNPs的r2值的和，哪個(gè)值大，我們就選擇哪一個(gè)來(lái)作為tagSNP，如果相等那么我們就任選其中一個(gè)來(lái)作為tagSNP，我們選擇的研究對(duì)象沒(méi)有出現(xiàn)這種情況.

2.2 r2的計(jì)算

給出m個(gè)個(gè)體，m/2個(gè)病例和m/2個(gè)對(duì)照，所有P個(gè)位點(diǎn)上的單體型Zij∈{0，1}，i=1，2，…，2m，j =1，2，…，P.計(jì)算生物學(xué)中描述SNP間相關(guān)關(guān)系的連鎖不平衡系數(shù)［18］r2:

2.3 SKAT

現(xiàn)有n個(gè)獨(dú)立個(gè)體，對(duì)于給定的含s個(gè)SNPs的SNP集，設(shè)Zi1，Zi2，…，Zis，是第i個(gè)個(gè)體在這s個(gè)SNPs上的基因型值，i=1，2，…，n.顯然，Zij∈{0，1，2}.第i個(gè)個(gè)體的定性性狀用yi表示，若個(gè)體i患病，則yi=1，否則yi=0.

下式(1)給出了個(gè)體定性性狀和基因型值間的半?yún)?shù)模型

本文使用SKAT檢驗(yàn)方法，分別對(duì)基于原始SNP集和基于標(biāo)簽SNP集進(jìn)行檢驗(yàn).基于原始SNP集的檢驗(yàn)用SKAT表示，基于標(biāo)簽SNP集的檢驗(yàn)用SKAT－tag表示.然后比較二者的P值和功效.

3 仿真研究

3.1 仿真數(shù)據(jù)

為了計(jì)算對(duì)兩種SNP集檢驗(yàn)的P值和功效，我們進(jìn)行了大量的仿真實(shí)驗(yàn)，本文的仿真數(shù)據(jù)均由HAPGEN2［22］產(chǎn)生，并且事先假定所有原因SNPs都會(huì)增加致病風(fēng)險(xiǎn).第13號(hào)染色體攜帶許多與乳腺癌有關(guān)的基因，我們選擇其中的MTRF1基因，它包含62個(gè)HapMap［23］SNPs.這62個(gè)SNPs中的10個(gè)SNPs已經(jīng)由Illumina HumanHap 500 array給出了確定的基因型.我們使用HapMap上的CEU樣本，用HAPGEN2基于CEU樣本的連鎖不平衡結(jié)構(gòu)產(chǎn)生仿真數(shù)據(jù).

我們使用HAPGEN2在不同的參數(shù)下產(chǎn)生MTRF1基因上62個(gè)SNPs的1000組仿真數(shù)據(jù)，每組包含500個(gè)病例和500個(gè)對(duì)照.我們從這100組中隨機(jī)選擇1組并從該組中隨機(jī)選取50個(gè)病例和50個(gè)對(duì)照的200條單體型數(shù)據(jù)上，在這200條單體型數(shù)據(jù)上運(yùn)用我們的方法選取tagSNPs，之后所有的仿真實(shí)驗(yàn)和假設(shè)檢驗(yàn)都以現(xiàn)在選取的SNPs作為tagSNPs.

3.2 P值計(jì)算

我們使用來(lái)自HapMap計(jì)劃中174個(gè)CEU種族個(gè)體的真實(shí)單體型數(shù)據(jù)，分別對(duì)tagSNPs集和原始SNPs集進(jìn)行關(guān)聯(lián)檢驗(yàn)，我們?cè)陲@著水平a=0.05下使用SKAT檢驗(yàn)方法，求其P值，見(jiàn)表1.

表1 P值

使用原始SNPs集，求得的P值為0.037，遠(yuǎn)小于0.05，所以使用原始集在統(tǒng)計(jì)學(xué)上顯著關(guān)聯(lián).使用我們方法選擇的tagSNPs集，所求的P值為0.048，也小于0.05，所以我們的方法選擇的SNPs集在統(tǒng)計(jì)學(xué)上微弱顯著關(guān)聯(lián)，因?yàn)槲覀兊姆椒ㄟx擇了少數(shù)的SNPs，所以丟失很多信息，我們的方法微弱顯著性關(guān)聯(lián)，說(shuō)明我們選擇的tagSNPs具有代表性，也說(shuō)明我們方法選擇的tagSNPs可以用來(lái)進(jìn)行疾病關(guān)聯(lián)檢驗(yàn).雖然我們的方法不如使用原始集關(guān)聯(lián)顯著，但是我們方法計(jì)算速度明顯高于基于原始集的方法.

3.3 第I類錯(cuò)誤率估計(jì)

利用HAPGEN2產(chǎn)生空模型下的1000組仿真數(shù)據(jù)來(lái)估計(jì)第I類錯(cuò)誤率，每組方針數(shù)據(jù)包含500個(gè)病例和500個(gè)對(duì)照.如表2得到的第I類錯(cuò)誤率，顯著水平為a=0.05.SKAT方法和SKAT－tag方法的第I類錯(cuò)誤率分別為0.049和0.042.說(shuō)明SKAT檢驗(yàn)方法都能很好地控制第I類錯(cuò)誤率.

表2 第I類錯(cuò)誤率

3.4 功效估計(jì)

我們將基因MTRF1上由Illumina HumanHap 500 array給出的10個(gè)已確定基因型的SNPs每個(gè)輪流作為致病SNP，并有HAPGEN2進(jìn)行仿真，10個(gè)已確定基因型的SNPs每個(gè)輪流一次就仿真1000組，總共為10000組，我們假定雜合子致病風(fēng)險(xiǎn)為1.25，純合子致病風(fēng)險(xiǎn)為1.5.表3給出了MTRF1基因上已確定基因型的SNPs.

表3 MTRF1基因上已確定的10個(gè)SNPs

我們將表3給出的10個(gè)SNPs中每個(gè)SNPs輪流作為致病SNP，每一個(gè)SNPs作為致病SNPs，分別使用SKAT和SKAT－tag求一次功效，顯著水平a=0.05，然后比較兩種方法，10個(gè)SNPs各輪流作為致病SNPs的功效，如圖1所示.

圖1 10個(gè)致病SNP分別在使用SKAT和SKAT－tag方法時(shí)的功效

4 結(jié)果分析和討論

在本文中我們提出了一種tagSNPs的選擇方法，降低了基因分型的成本和計(jì)算的時(shí)間復(fù)雜度.與原始SNPs集相比，我們不需要對(duì)所有的SNPs進(jìn)行基因型的測(cè)定，也不需要對(duì)所有SNPs進(jìn)行關(guān)聯(lián)檢驗(yàn)，而只需要對(duì)我們選擇的tagSNPs進(jìn)行基因分型和關(guān)聯(lián)檢驗(yàn).通過(guò)我們的方法對(duì)MTRF1基因上62個(gè)SNPs進(jìn)行tagSNPs的選擇，最終我們選擇了rs666930，rs586650，rs550174，rs483180，rs616111，rs6668589，rs523395，rs2246410，rs512854這11個(gè)tagSNPs，數(shù)量大約為原始SNPs集的1/6，基因分型的成本也就降低了5/6，P值的計(jì)算時(shí)間也會(huì)減少.因?yàn)槲覀兪紫仍谝恍〔糠謹(jǐn)?shù)據(jù)中選擇tagSNPs，所以總的來(lái)說(shuō)，使用tagSNPs進(jìn)行關(guān)聯(lián)檢驗(yàn)的時(shí)間復(fù)雜度比使用原始SNPs集進(jìn)行關(guān)聯(lián)檢驗(yàn)的時(shí)間復(fù)雜度要小.

通過(guò)使用來(lái)自HapMap的174個(gè)CEU種族個(gè)體數(shù)據(jù)，我們對(duì)兩種SNPs集使用SKAT檢驗(yàn)方法求得P值，如表1所示，基于tagSNPs集檢驗(yàn)的P值小于基于原始SNPs集檢驗(yàn)的P值.但是二者的P值都小于顯著性水平，也就兩種SNPs集對(duì)疾病都顯著性關(guān)聯(lián)，所以我們的方法求得的tagSNPs可以用來(lái)進(jìn)行疾病關(guān)聯(lián)檢驗(yàn).從表2可以看出，SKAT對(duì)第I類錯(cuò)誤率都是可控的.圖1表明基于原始SNPs集檢驗(yàn)的功效與基于tagSNPs集檢驗(yàn)的功效在大部分情況下相差不大.但在某些情況下，基于tagSNPs集檢驗(yàn)的功效小于基于原始SNPs集檢驗(yàn)的功效.還有使用tagSNPs集得到P值小于使用原始SNPs集得到的P值，很大程度上是因?yàn)殛P(guān)聯(lián)檢驗(yàn)時(shí)沒(méi)有包含所有SNPs，所以造成了信息的大量丟失.某些情況我們方法功效大于原始SNPs集方法，可能因?yàn)镾NPs數(shù)量的減少，造成了自由度的降低.總的來(lái)說(shuō)，我們的方法計(jì)算復(fù)雜度低，基因分型成本低，是可以用來(lái)進(jìn)行疾病關(guān)聯(lián)檢驗(yàn).

［1］Collins F.S，et al.Variations on a theme:cataloging human DNA sequence variation［J］.Science，1997(278):1580－1581.

［2］G.C.Johnson，L.Esposito，B.J.Barratt，et al.Haplotype tagging for the identification of common disease genes［J］.Nature Genetics，2001(2):233－237.

［3］K.Zhang，M.Deng，T.Chen，M.S.Waterman，F(xiàn).Sun.A dynamic programming algorithm for haplotype block partitioning［J］.Proceedings of the National Academy of Sciences of the United States of America，2002(11):7335－7339.

［4］E.C.Anderson，J.Novembre.Finding haplotype block boundaries by using the minimum－description－length principle［J］.American Journal of Human Genetics，2003(2):336－354.

［5］K.Zhang，P.Calabrese，M.Nordborg，F(xiàn).Sun.Haplotype block structure and its applications to association studies:power and study designs［J］.American Journal of Human Genetics，2002(6):1386－1394.

［6］Weale M.E.，et al.Selection and evaluation of tagging SNPs in the neuronal－sodiumchannel gene SCN1A:implications for linkagedisequilibrium gene mapping［J］.Am.J.Hum.Genet，2003(73):551－565.

［7］Ke X.，Cardon L.R.Efficient selective screening of haplotype tag SNPs［J］.Bioinformatics，2003(19):287－288.

［8］Carlson C.S.，et al.Selecting a maximally informative set of single－nucleotide polymorphisms for association analyses using linkage disequilibrium［J］.Am.J.Hum.Genet，2004(74):106－120.

［9］N.Patil，A.J.Berno，D.A.Hinds，et al.Blocks of limited haplotype diversity revealed by high－resolution scanning of human chromosome 21［J］.Science，2001(5547):1719－1723.

［10］Cardon L.R.，Abecasis，G.R.Using haplotype blocks to map human complex trait loci［J］.Trends Genet.，2003(19):135－140.

［11］Stram D.O.，et al.Choosing haplotype－tagging SNPs based on unphased genotype data using preliminary sample of unrelated subjects with an example from the multiethic cohort study［J］.Hum.Hered.，2003(55):27－36.

［12］Morgenthaler S.，Thilly W.G.A strategy to discover genes that carry multi－allelic or mono－allelic risk for common diseases:a cohort allelic sums test(CAST)［J］.Mutat.Res，2007(615):28－56.

［13］Morris A.P.，Zeggini E.An evaluation of statistical approaches to rare variant analysis in genetic association studies［J］.Genet.Epidemiol，2010(34):188－193.

［14］Li B.，Leal S.M.Methods for detecting associations with rare variants for common diseases:application to analysis of sequence data［J］.Am.J.Hum.Genet，2008(83):311－321.

［15］Madsen B.E.，Browning S.R.A groupwise association test for rare mutations using a weighted sum statistic［J］.PLoS Genet，2009(5): e1000384.

［16］Wu M C，Kraft P，Epstein M P，et al.Powerful SNP－set analysis for case－control genome－wide association studies［J］.The Ameri-can Journal of Human Genetics，2010，86(6):929－942.

［17］Devlin B，Risch N.A comparison of linkage disequilibrium measures for fine－scale mapping［J］.Genomics，1995(29):311－322.

［18］Hill W G，Robertson A.Linkage disequilibrium in finite populations［J］.Theoretical and Applied Genetics，1968，38(6):226－231.

［19］Liu D，Ghosh D，Lin X.Estimation and testing for the effect of a genetic pathway on a disease outcome using logistic kernel machine regression via logistic mixed models［J］.BMC bioinformatics，2008，9(1):1－11.

［20］Zhang D，Lin X.Hypothesis testing in semiparametric additive mixed models［J］.Biostatistics，2003，4(1):57－74.

［21］Lin X，Cai T，Wu M C，et al.Kernel machine SNP‐set analysis for censored survival outcomes in genome‐wide association studies［J］.Genetic epidemiology，2011，35(7):620－631.

［22］Su Z，Marchini J，Donnelly P.HAPGEN2:simulation of multiple disease SNPs［J］.Bioinformatics，2011，27(16):2304－2305.

［23］The International HapMap Consortium.The International HapMap Project［J］.Nature，2003(426):789－796.

The Method of Selecting tagSNPs Based on Linkage Disequilibrium

WANG Jun－feng，WANG Xin－zeng
(School of Mathematics and Systems Science，Shandong University of Science and Technology，Qingdao，266590，China)

In genome－wide association studies，we need to have a sufficiently dense single nucleotide polymorphisms set to explain part of the genetic risk for common diseases.Within candidate genes，the number of common polymorphisms is finite，but direct assay of all existing common polymorphism is inefficient，because genotypes at many of these sites are strongly correlated，can lead to a large amount of redundant information，and will result in an increase in the cost of genotyping，consume large amounts of time.So when we test the association of markers with disease，typing all available SNP markers is inefficient and not necessary.We only need to select a representative，small number of SNPs for genotyping，and test the association between these SNPs and disease.The SNPs selected here called tagSNPs，it is a small subset of the SNPs，and enough to capture the haplotype information in every haplotype region.The selection of tagSNPs has become a very active research topic and many strategies have been proposed.In this paper，we put forward a new kind of tagSNPs selection method，by using measure based on pairwise LD to group the haplotype，divided into disjoint groups，and selected the tagSNPs in each group.Compared with testing method based on original SNPs sets，our method has produced less tagSNPs，while simultaneously maximizing the information content by selected markers，reducing the cost of genotyping，and improving the efficiency.

genome－wide association studies;SNP;tagSNPs;linkage disequilibrium

Q811.4

A

1672－2590(2016)03－0049－06

2016－04－03

國(guó)家自然科學(xué)基金資助項(xiàng)目(61572522)

王鈞峰(1990－)，男，河北滄州人，山東科技大學(xué)數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院碩士研究生.