王鈞峰，王新贈

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

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

王鈞峰，王新贈

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

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

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

1 引言

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

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

2 方法

2.1 tagSNP的選擇

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

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

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

2.2 r2的計算

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

2.3 SKAT

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

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

本文使用SKAT檢驗方法，分別對基于原始SNP集和基于標簽SNP集進行檢驗.基于原始SNP集的檢驗用SKAT表示，基于標簽SNP集的檢驗用SKAT－tag表示.然后比較二者的P值和功效.

3 仿真研究

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

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

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

3.2 P值計算

我們使用來自HapMap計劃中174個CEU種族個體的真實單體型數(shù)據(jù)，分別對tagSNPs集和原始SNPs集進行關聯(lián)檢驗，我們在顯著水平a=0.05下使用SKAT檢驗方法，求其P值，見表1.

表1 P值

使用原始SNPs集，求得的P值為0.037，遠小于0.05，所以使用原始集在統(tǒng)計學上顯著關聯(lián).使用我們方法選擇的tagSNPs集，所求的P值為0.048，也小于0.05，所以我們的方法選擇的SNPs集在統(tǒng)計學上微弱顯著關聯(lián)，因為我們的方法選擇了少數(shù)的SNPs，所以丟失很多信息，我們的方法微弱顯著性關聯(lián)，說明我們選擇的tagSNPs具有代表性，也說明我們方法選擇的tagSNPs可以用來進行疾病關聯(lián)檢驗.雖然我們的方法不如使用原始集關聯(lián)顯著，但是我們方法計算速度明顯高于基于原始集的方法.

3.3 第I類錯誤率估計

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

表2 第I類錯誤率

3.4 功效估計

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

表3 MTRF1基因上已確定的10個SNPs

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

圖1 10個致病SNP分別在使用SKAT和SKAT－tag方法時的功效

4 結果分析和討論

在本文中我們提出了一種tagSNPs的選擇方法，降低了基因分型的成本和計算的時間復雜度.與原始SNPs集相比，我們不需要對所有的SNPs進行基因型的測定，也不需要對所有SNPs進行關聯(lián)檢驗，而只需要對我們選擇的tagSNPs進行基因分型和關聯(lián)檢驗.通過我們的方法對MTRF1基因上62個SNPs進行tagSNPs的選擇，最終我們選擇了rs666930，rs586650，rs550174，rs483180，rs616111，rs6668589，rs523395，rs2246410，rs512854這11個tagSNPs，數(shù)量大約為原始SNPs集的1/6，基因分型的成本也就降低了5/6，P值的計算時間也會減少.因為我們首先在一小部分數(shù)據(jù)中選擇tagSNPs，所以總的來說，使用tagSNPs進行關聯(lián)檢驗的時間復雜度比使用原始SNPs集進行關聯(lián)檢驗的時間復雜度要小.

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

國家自然科學基金資助項目(61572522)

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