南京醫(yī)科大學(xué)公共衛(wèi)生學(xué)院生物統(tǒng)計學(xué)系（211166）

段巍巍 趙 楊 張麗偉 胡志斌 陳 峰△

使用肺癌GWAS數(shù)據(jù)進行遺傳風(fēng)險預(yù)測的方法和策略研究

南京醫(yī)科大學(xué)公共衛(wèi)生學(xué)院生物統(tǒng)計學(xué)系（211166）

段巍巍 趙 楊 張麗偉 胡志斌 陳 峰△

目的探討基于肺癌全基因組關(guān)聯(lián)研究數(shù)據(jù)的遺傳風(fēng)險預(yù)測方法和策略。方法將肺癌GWAS數(shù)據(jù)中的南京子樣本和北京子樣本分別作為訓(xùn)練集和測試集，分別使用預(yù)測全集和最優(yōu)預(yù)測子集兩種策略，比較三種預(yù)測方法在不同連鎖不平衡結(jié)構(gòu)（LD）和初篩檢驗水準(zhǔn)（α）下的預(yù)測準(zhǔn)確度。結(jié)果wGRS在高LD結(jié)構(gòu)下，隨著-log（α）增大，預(yù)測準(zhǔn)確度呈現(xiàn)上升趨勢；RF和SVM對LD結(jié)構(gòu)不如wGRS敏感，但三種方法在低LD結(jié)構(gòu)（r2＜0．2）下預(yù)測準(zhǔn)確度優(yōu)于高LD結(jié)構(gòu)；wGRS方法下最優(yōu)預(yù)測子集效果略優(yōu)于預(yù)測全集效果，SVM下子集效果與全集近似，但略遜于全集，RF下子集效果則不如全集，且差距較大。結(jié)論基于LD結(jié)構(gòu)修剪SNP位點和選擇適當(dāng)?shù)某鹾Y水準(zhǔn)可以提高遺傳風(fēng)險預(yù)測準(zhǔn)確度，此時wGRS方法預(yù)測效果優(yōu)于SVM和RF。

肺癌 遺傳風(fēng)險得分 支持向量機 隨機森林 最優(yōu)預(yù)測子集 單核苷酸多態(tài)性

近些年來，全基因組關(guān)聯(lián)研究（genome-w ide association study，GWAS）蓬勃發(fā)展。截止2013年，全球的研究者們累計發(fā)現(xiàn)了與600多種表型（疾?。┫嚓P(guān)的15000多個單核苷酸多態(tài)性（single nucleotide polymorphism，SNP）位點［1］，而利用這些已發(fā)現(xiàn)的位點進行臨床個體化醫(yī)療實踐、疾病預(yù)防等成為后GWAS（post-GWAS）時代的主要目標(biāo)之一。而建立一個準(zhǔn)確的遺傳風(fēng)險預(yù)測模型則成為GWAS研究成果轉(zhuǎn)化最為關(guān)鍵的一步。早期的遺傳風(fēng)險預(yù)測研究表明，利用GWAS獲得的關(guān)聯(lián)位點進行預(yù)測并不成功［2－4］，其中的一個主要原因在于GWAS研究為了控制假陽性關(guān)聯(lián)位點，多階段驗證過于嚴(yán)格，具體到某一類疾病（表型）而言，關(guān)聯(lián)位點數(shù)量較少［5－7］。本研究基于肺癌病例-對照GWAS研究數(shù)據(jù)，選擇靈活的關(guān)聯(lián)位點初篩水準(zhǔn)，采用新的策略，即以能否增加預(yù)測效果作為進一步納入SNP位點進行預(yù)測的標(biāo)準(zhǔn)，并采用多種預(yù)測方法進行肺癌的遺傳風(fēng)險預(yù)測效果分析和比較。

資料與方法

1．研究對象與質(zhì)量控制

數(shù)據(jù)來自于一項非小細(xì)胞肺癌的病例-對照GWAS研究［7］，包含南京樣本（1473個病例和1962個對照）和北京樣本（858個病例和1115個對照）兩個部分，詳細(xì)的樣本人群信息參見文獻［7］。所有研究個體的DNA樣本來自于全血，使用Affymetrix Genome-W ide Human SNP Array 6．0芯片進行基因分型，經(jīng)過嚴(yán)格的質(zhì)量控制后［7－8］，用于分析的共570373個SNP位點。

2．研究策略與統(tǒng)計分析

（1）研究策略

將研究數(shù)據(jù)集分為訓(xùn)練集和測試集。首先在訓(xùn)練集中進行位點初篩，即每一個位點都與表型變量做logistic回歸。隨后僅對P值小于初篩水準(zhǔn)α的位點進行分析。對初篩后的位點考慮兩種研究策略：①全部納入初篩時有統(tǒng)計學(xué)意義（P＜α）的位點進行預(yù)測分析，簡稱“預(yù)測全集”，即Performance［phenotype，score（SNP）］；其中Performance表示預(yù)測準(zhǔn)確度，phenotype表示表型變量，連續(xù)表型時使用決定系數(shù)（R2），二分類表型時使用受試者工作曲線（receiver operating characteristic curve，ROC）下面積（area under curve，AUC）；score表示個體在候選位點上的得分值。②僅納入①中位點的最優(yōu)預(yù)測子集，選擇方法如下：

a．對于所有初篩后的m個位點，選擇其中具有最佳預(yù)測準(zhǔn)確度的位點。對于連續(xù)表型，有標(biāo)準(zhǔn)化后基因型矩陣G和表型向量y，令相關(guān)系數(shù)向量corr＝Gy則向量中最大值對應(yīng)的位點即為m個位點中的最佳預(yù)測位點，此處為每個位點回歸系數(shù)組成的向量；同理，對于二分類表型，則選擇m個位點中AUC最大的位點。

b．從剩下的m-1個位點中挑選一個位點納入預(yù)測模型，使得新模型的預(yù)測能力提升最大。選擇第i個位點：表示現(xiàn)有模型的遺傳得分表示第i個位點的遺傳得分。

c．重復(fù)過程b，直到納入某一個位點后模型的預(yù)測能力不再提高時，則停止搜索。此時納入模型的SNP位點即組成了最優(yōu)預(yù)測子集。

最后，在測試集中評價相關(guān)預(yù)測方法和上述策略的預(yù)測效果。

（2）統(tǒng)計分析

用南京樣本數(shù)據(jù)作為訓(xùn)練集，北京樣本數(shù)據(jù)作為測試集。位點初篩時使用logistic回歸，采用相加生物學(xué)模式，并校正了如下因素：年齡、性別、吸煙（包／年）以及人群分層的前兩個主成分。本文使用的遺傳預(yù)測方法有：

①效應(yīng)加權(quán)遺傳風(fēng)險得分（effect-weighted genetic risk score，wGRS），即此時Gi表示第i個位點的次等位基因計數(shù)（取值為0、1、2）向量，wi為第i個位點的權(quán)重，wGRS的方法是現(xiàn)行比較通用的方法［9－12］。

②支持向量機（support vectormachine，SVM），在SVM的參數(shù)設(shè)置中，考慮到位點初篩后的數(shù)據(jù)結(jié)構(gòu)（較低的維度和可能的非線性關(guān)系）以及常用核函數(shù)的相關(guān)性質(zhì)，本文擬選用高斯核函數(shù)（radial basis function，RBF），即與此同時，將其懲罰參數(shù)Cost設(shè)置為1，而核參數(shù)Gamma設(shè)置為位點數(shù)的倒數(shù)。SVM的相關(guān)原理和遺傳風(fēng)險預(yù)測可以參考其他研究［13－15］。

③隨機森林（random forest，RF），主要參數(shù)中森林中決策樹的數(shù)量（ntree）、內(nèi)部節(jié)點隨機選擇的特征數(shù)（m try）及終末節(jié)點最小樣本數(shù)（nodesize）分別設(shè)置為500、位點數(shù)開方和1，RF的原理和分類應(yīng)用可以參見早期的研究［16－18］。在應(yīng)用相關(guān)預(yù)測方法之前我們對缺失數(shù)據(jù)進行了均值填補。

數(shù)據(jù)處理和分析使用plink 1．7軟件和R v2．15語言包，其中SVM算法使用R語言e1071包，RF算法使用random Forest包。

（3）遺傳風(fēng)險預(yù)測評價指標(biāo)

本研究選擇AUC作為預(yù)測效果評價指標(biāo)，其理論取值范圍為0．5～1，取值越大，表明預(yù)測能力越強。AUC及其方差的估計使用Mann-Whitney檢驗的U統(tǒng)計量［19］。

（4）研究參數(shù)設(shè)置

考慮到GWAS研究中多重比較的檢驗水準(zhǔn)（α）較為嚴(yán)格，我們在預(yù)測研究中設(shè)定了靈活的檢驗水準(zhǔn)，即-log（α）設(shè)定為4～9；為了考慮連鎖不平衡（linkage disequilibrium，LD）結(jié)構(gòu)對預(yù)測結(jié)果的影響，我們對候選的SNP位點進行了修剪，分別設(shè)定LD參數(shù)為r2＜1、r2＜0．5和r2＜0．2，分別表示不修剪任何SNP位點、修剪掉r2大于0．5的位點、r2大于0．2的位點。

結(jié) 果

在兩種研究策略下分別比較wGRS、SVM和RF三種分類方法的預(yù)測能力。表1為wGRS方法下訓(xùn)練集自身驗證（self-validation）的結(jié)果。結(jié)果表明，最優(yōu)預(yù)測子集的預(yù)測效果始終是優(yōu)于全集，并具有較高的AUC值。LD結(jié)構(gòu)是遺傳數(shù)據(jù)分析中不可忽視的影響因素，隨著r2閾值的降低，納入各個預(yù)測集合的位點數(shù)降低。對于預(yù)測全集而言，r2＜0．5與r2＜0．2情況下預(yù)測能力優(yōu)于r2＜1，而它們之間的差距隨著初篩閾值逐漸嚴(yán)格（此時組成各自預(yù)測全集的位點集逐漸相同）而變小；在相同的α值下，高LD閾值下構(gòu)建的預(yù)測全集包含了低LD下的預(yù)測全集位點，而類似情況下，預(yù)測子集則不存在這種包含關(guān)系。隨著-log（α）的逐漸增大，AUC基本呈現(xiàn)一個緩慢下降的趨勢。

測試集驗證的結(jié)果見圖1。SVM隨著-log（α）增大，預(yù)測全集與最優(yōu)預(yù)測子集保持著相同的變化趨勢，RF在-log（α）＜7時AUC值保持穩(wěn)定，隨后下降明顯，而SVM則在5＜-log（α）＜7時穩(wěn)定，過高或過低的α值都會引起AUC值的下降；考慮到SVM和RF在共線性問題上相較wGRS更加穩(wěn)健、SVM對于噪聲位點的耐受性［20］以及RF可以納入位點間可能的交互作用等因素，在較低的-log（α）下，其結(jié)果優(yōu)于wGRS，但隨著集合內(nèi)位點減少，其效果亦出現(xiàn)下降；RF在-log（α）＝4處呈現(xiàn)領(lǐng)先優(yōu)勢，可能原因在于交互作用處理上的優(yōu)勢。由于wGRS對高LD的敏感性，其在高LD結(jié)構(gòu)下，隨著-log（α）增大，保持著上升趨勢（需要指出的是，本研究數(shù)據(jù)的特征是：-log（α）越大，全集和子集內(nèi)位點數(shù)減少，集合內(nèi)高相關(guān)的位點也逐漸減少）；而在低LD結(jié)構(gòu)下，隨著-log（α）增大，集合內(nèi)位點數(shù)目成為主要影響因素，因此呈現(xiàn)出與高LD不同的變化趨勢。r2＜0．2情況下預(yù)測結(jié)果基本上優(yōu)于前兩者；wGRS方法下預(yù)測子集效果略優(yōu)于預(yù)測全集效果，SVM下子集效果與全集近似，但略遜于全集，RF下子集效果則不如全集，且差距較大，可能的原因在于構(gòu)建最優(yōu)預(yù)測子集時位點采用線性結(jié)合方式，因此不適用于SVM和RF。在r2＜0．2和-log（α）＝6處獲得參數(shù)組合下最大AUC，估計值為0．7149，95%CI為0．6912～0．7387。

表1 訓(xùn)練集自身驗證預(yù)測準(zhǔn)確度

圖1 測試集中預(yù)測準(zhǔn)確度變化圖

討 論

本研究采用預(yù)測全集和最優(yōu)預(yù)測子集兩種研究策略，在不同LD結(jié)構(gòu)參數(shù)和初篩水準(zhǔn)α下比較wGRS、SVM和RF這三種常用預(yù)測方法的預(yù)測準(zhǔn)確度?？傮w而言，最優(yōu)子集下的wGRS方法在r2＜0．2和-log（α）＝6的參數(shù)下獲得了較為滿意預(yù)測準(zhǔn)確度（AUC＝0．715）。而對比同類研究，Li等人［12］利用其他GWAS研究發(fā)現(xiàn)的4個肺癌關(guān)聯(lián)位點進行預(yù)測，其AUC值僅為0．555。Wei等人［13］在1型糖尿病的預(yù)測研究中也發(fā)現(xiàn)，設(shè)定適當(dāng)寬松的α值有助于提高預(yù)測準(zhǔn)確度。如果假設(shè)SNP位點的效應(yīng)服從于一個均數(shù)為0的正態(tài)分布［21］，那么大部分的關(guān)聯(lián)研究都僅發(fā)現(xiàn)了分布尾巴處＂數(shù)量少、效應(yīng)大＂的位點。由此可見，忽略中間大部分低效應(yīng)位點必然影響預(yù)測模型的準(zhǔn)確度，而寬松的α值可以在一定程度上納入這些潛在的關(guān)聯(lián)位點，當(dāng)然本研究和Wei等人的研究也都展示了這種“寬松”是有限度的，即受制于噪聲位點的大量增多，預(yù)測準(zhǔn)確度會逐漸降低。

高LD結(jié)構(gòu)對線性結(jié)合的wGRS方法影響較大，因此參照其他研究對SNP位點進行修剪（LD-pruning）［22－24］。RF和SVM對LD結(jié)構(gòu)不如wGRS敏感，但三種方法在低LD結(jié)構(gòu)下表現(xiàn)均優(yōu)于高LD。對于兩種預(yù)測位點集合，考慮到訓(xùn)練集與測試集可能存在數(shù)據(jù)結(jié)構(gòu)上的差異，最優(yōu)預(yù)測子集在測試數(shù)據(jù)中并不是完全優(yōu)于預(yù)測全集，而構(gòu)建最優(yōu)預(yù)測子集的方法采用了wGRS的得分相加方式，因此最優(yōu)子集并不適用于RF和SVM方法。綜合而言，在低LD結(jié)構(gòu)的情況下，wGRS方法優(yōu)勢顯著，此時可以選擇最優(yōu)預(yù)測子集；而在高LD和寬松的初篩閾值下，由于過多的假陽性位點和共線性結(jié)構(gòu)，所有方法均表現(xiàn)不佳，但RF相對穩(wěn)定。

在同類的遺傳風(fēng)險預(yù)測研究中，有相當(dāng)一部分的結(jié)果［25－27］不甚理想，而這些預(yù)測研究的共同特點在于使用的風(fēng)險位點均是經(jīng)過多階段嚴(yán)格驗證。由此可見，受限于研究設(shè)計和相關(guān)技術(shù)，利用苛刻條件下發(fā)現(xiàn)的少量常見變異、強效應(yīng)關(guān)聯(lián)位點很難獲得令人振奮的預(yù)測能力。但是我們不應(yīng)該以消極的態(tài)度來看待關(guān)聯(lián)研究的轉(zhuǎn)化應(yīng)用。值得注意的是，我們認(rèn)識到遺傳度（heritability）低的疾病會限制遺傳預(yù)測方法效果的提升，以本研究的肺癌疾病為例，國外的一項研究表明［28］，瑞典人群肺癌的遺傳度僅在8%左右。對于本研究中使用的預(yù)測指標(biāo)AUC，很多學(xué)者提出異議［4，29］，他們發(fā)現(xiàn)數(shù)個較強的預(yù)測指標(biāo)線性組合后的AUC值相較某單一指標(biāo)提高幅度有限，于是對AUC作為預(yù)測評價指標(biāo)產(chǎn)生了質(zhì)疑，并提出了再分類優(yōu)值（net reclassification improvement，NRI）、綜合判別優(yōu)值（integrated discrimination improvement，IDI）等指標(biāo)［29］，但考慮到AUC指標(biāo)的綜合評價能力、廣泛的通用性以及研究目的，本文依然選擇AUC作為研究的預(yù)測評價指標(biāo)。這里討論的方法均為成熟或經(jīng)典的預(yù)測方法，最近有研究者［30－31］在探討使用混合效應(yīng)模型進行預(yù)測，亦有研究者整合多組學(xué)信息進行預(yù)測研究［32］，但這些方法還有待大量的實例研究去驗證。

（致謝：感謝南京醫(yī)科大學(xué)公共衛(wèi)生學(xué)院分子流行病實驗室和北京子研究的老師、同學(xué)給予的支持）

［1］Welter D，MacArthur J，Morales J，et al．The NHGRIGWAS Catalog，a curated resource of SNP-trait associations．Nucleic Acids Res，2014，42（Database issue）：D1001-1006．

［2］van der Net JB，Janssens AC，Sijbrands EJ，et al．Value of genetic profiling for the prediction of coronary heart disease．Am Heart J，2009，158（1）：105-110．

［3］M ihaescu R，Meigs J，Sijbrands E，etal．Genetic risk profiling for prediction of type2 diabetes．PLoSCurr，2011，3：RRN1208．

［4］Cook NR．Use and misuse of the receiver operating characteristic curve in risk prediction．Circulation，2007，115（7）：928-935．

［5］McCarthy MI，Abecasis GR，Cardon LR，etal．Genome-wide association studies for complex traits：consensus，uncertainty and challenges．Nat Rev Genet，2008，9（5）：356-369．

［6］Manolio TA．Genomewide association studies and assessment of the risk of disease．N Engl JMed，2010，363（2）：166-176．

［7］Hu Z，Wu C，Shi Y，etal．A genome-w ide association study identifies two new lung cancer susceptibility loci at13q12．12 and 22q12．2 in Han Chinese．Nature genetics，2011，43（8）：792-796．

［8］陳峰，柏建嶺，趙楊，等．全基因組關(guān)聯(lián)研究中的統(tǒng)計分析方法．中華流行病學(xué)雜志，2011，32（4）：400-404．

［9］Speliotes EK，W iller CJ，Berndt SI，et al．Association analyses of 249，796 individuals reveal18 new lociassociated with body mass index．Nat Genet，2010，42（11）：937-948．

［10］Ripatti S，Tikkanen E，Orho-Melander M，et al．A multilocus genetic risk score for coronary heart disease：case-control and prospective cohort analyses．Lancet，2010，376（9750）：1393-1400．

［11］Gui L，Wu F，Han X，et al．A multilocus genetic risk score predicts coronary heart disease risk in a Chinese Han population．Atherosclerosis，2014，237（2）：480-485．

［12］Li H，Yang L，Zhao X，et al．Prediction of lung cancer risk in a Chinese population using amultifactorial geneticmodel．BMCMed Genet，2012，13：118．

［13］Wei Z，Wang K，Qu HQ，et al．From disease association to risk assessment：an optim istic view from genome-w ide association studies on type 1 diabetes．PLoSGenet，2009，5（10）：e1000678．

［14］Becker N，Toedt G，Lichter P，et al．Elastic SCAD as a novel penalization method for SVM classification tasks in high-dimensional data．BMC Bioinformatics，2011，12：138．

［15］Cortes C，Vapnik V．Support-vector networks．Machine learning，1995，20（3）：273-297．

［16］Breiman L．Random forests．Machine learning，2001，45（1）：5-32．

［17］Zhang H，Yu CY，Singer B．Cell and tumor classification using gene expression data：construction of forests．Proc Natl Acad Sci USA，2003，100（7）：4168-4172．

［18］Yoon D，Kim YJ，Park T．Phenotype prediction from genome-w ide association studies：application to smoking behaviors．BMC Syst Biol，2012，6（2）：S11．

［19］Mason SJ，Graham NE．Areas beneath the relative operating characteristics（ROC）and relative operating levels（ROL）curves：Statistical significance and interpretation．Quarterly Journal of the Royal Meteorological Society，2002，128（584）：2145-2166．

［20］Xing EP，Jordan MI，Karp RM．Feature selection for high-dimensional genom ic m icroarray data．in ICML，2001：Citeseer．

［21］Speed D，Hemani G，Johnson MR，et al．Improved heritability estimation from genome-w ide SNPs．Am JHum Genet，2012，91（6）：1011-1021．

［22］Purcell SM，W ray NR，Stone JL，et al．Common polygenic variation contributes to risk of schizophrenia and bipolar disorder．Nature，2009，460（7256）：748-752．

［23］Consortium SWGO．Biological insights from 108 schizophrenia-associated genetic loci．Nature，2014，511（7510）：421-427．

［24］Chatterjee N，Wheeler B，Sampson J，et al．Projecting the performance of risk prediction based on polygenic analyses of genome-w ide association studies．Nat Genet，2013，45（4）：400-5，405e1-3．

［25］Zheng SL，Sun J，Wiklund F，etal．Genetic variantsand family history predict prostate cancer sim ilar to prostate-specific antigen．Clin Cancer Res，2009，15（3）：1105-1111．

［26］Mealiffe ME，Stokowski RP，Rhees BK，et al．Assessment of clinical validity of a breast cancer risk model combining genetic and clinical information．JNatl Cancer Inst，2010，102（21）：1618-1627．

［27］Wacholder S，Hartge P，Prentice R，et al．Performance of common genetic variants in breast-cancer risk models．N Engl JMed，2010，362（11）：986-993．

［28］Czene K，Lichtenstein P，Hemminki K．Environmental and heritable causes of cancer among 9．6million individuals in the Swedish family -cancer database．International Journal of Cancer，2002，99（2）：260-266．

［29］Pencina MJ，D’Agostino RS，D’Agostino RJ，et al．Evaluating the added predictive ability of a new marker：from area under the ROC curve to reclassification and beyond．StatMed，2008，27（2）：157-72；discussion 207-212．

［30］Golan D，Rosset S．Effective genetic-risk prediction using mixed models．Am JHum Genet，2014，95（4）：383-393．

［31］Speed D，Balding DJ．MultiBLUP：improved SNP-based prediction for complex traits．Genome Res，2014，24（9）：1550-1557．

［32］Wheeler HE，Aquino-M ichaels K，Gamazon ER，et al．Poly-om ic prediction of complex traits：OmicKriging．Genet Epidemiol，2014，38（5）：402-415．

（責(zé)任編輯：郭海強）

Strategies of Genetic Risk Prediction with Lung Cancer GWASData

Duan Weiwei，Zhao Yang，Zhang Liwei，et al．（DepartmentofBiostatistics，NanjingMedicalUniversity（211166），Nanjing）

ObjectiveTo investigate the performance of three genetic risk prediction methods，weighted genetic risk score（wGRS），supportvector machine（SVM）and random forest（RF），applied to high dimensional data of lung cancerwith two strategies．MethodsThis study served Nanjing and Beijing samples of GWAS data as training set and testing set respectively．Wemade use of the two strategies of Full predictive subset（FS）and Best predictive subset（BS）and compared the prediction accuracy within the threemethodsmentioned above with the combination of Linkage Disequilibrium（LD）and hypothesis testing levels（α）．ResultsUnder a high LD structure，the prediction accuracy of wGRSwas on the rise with the increasing-log（α）．RF and SVM were not sensitive to LD structures as wGRS，but the predictive accuracy of each method applied with a low LD structure（r2＜0．2）wasmainly better than itself with a high LD structure．Moreover，the performance of BS was slightly better than，approximately equal to or tiny less than and worse than FSwhen themethodswere respectively wGRS，SVM and RF．ConclusionThe prediction accuracy could be improved with the condition of LD-pruning and adopting a properα-value，meanwhile，wGRSwas better than SVM and RF in that condition．

Lung cancer；Genetic risk score；Support vectormachine；Random forest；Best predictive subset；Single nucleotide polymorphism

國家自然科學(xué)基金（81473070，81373102）

△通信作者：陳峰，Email：fengchen＠njmu．edu．cn