使用肺癌GWAS數據進行遺傳風險預測的方法和策略研究

南京醫科大學公共衛生學院生物統計學系（211166）

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

使用肺癌GWAS數據進行遺傳風險預測的方法和策略研究

南京醫科大學公共衛生學院生物統計學系（211166）

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

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

肺癌 遺傳風險得分 支持向量機 隨機森林 最優預測子集 單核苷酸多態性

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

資料與方法

1．研究對象與質量控制

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

2．研究策略與統計分析

（1）研究策略

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

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

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

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

最后，在測試集中評價相關預測方法和上述策略的預測效果。

（2）統計分析

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

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

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

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

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

（3）遺傳風險預測評價指標

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

（4）研究參數設置

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

結 果

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

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

表1 訓練集自身驗證預測準確度

圖1 測試集中預測準確度變化圖

討 論

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

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

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

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

［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］陳峰，柏建嶺，趙楊，等．全基因組關聯研究中的統計分析方法．中華流行病學雜志，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．

（責任編輯：郭海強）

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

國家自然科學基金（81473070，81373102）

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