高級醫(yī)學統(tǒng)計學：交互作用(二)

交互作用(二)

DetectingInteractionviaMachineLearningMethods研究生《醫(yī)學統(tǒng)計學》課程2美國人類基因組研究所/GWAStudies截止2013年3月，共有1,552篇公開發(fā)表的GWAS文獻。幾乎所有的研究僅僅局限于單位點的研究，而沒有深入探討基因—基因交互作用。主要內(nèi)容窮舉式檢索BOOST多因子降維隨機式檢索隨機森林啟發(fā)式檢索螞蟻優(yōu)化算法3GWAS中交互作用分析方法分類窮舉式檢測考慮變量間所有組合。若變量總數(shù)為p，檢測k-1階交互作用，檢驗次數(shù)為。計算負擔極重，多重比較校正力度極大。即使僅有100個變量：一階交互作用有4,950種組合；二階交互作用有161,700種組合；三階交互作用有3,921,225種組合。隨機式檢測從所有變量中隨機挑選部分變量，考察這些變量的重要性。重復上述步驟，盡可能覆蓋所有變量。變量數(shù)很多時，效率和效能難以平衡：需要重復極多次，才能覆蓋所有變量。此時，計算負擔仍然很重；若重復次數(shù)有限，無法覆蓋所有變量。此時，容易遺漏交互作用。啟發(fā)式檢測通過指定的統(tǒng)計學或生物學標準，過濾掉絕大部分噪音變量，保留少數(shù)信號變量。在剩余子集中進行分析。不恰當?shù)臉藴嗜菀走z漏交互作用：統(tǒng)計學：過濾掉無主效應的位點，僅針對有主效應的位點分析；生物學：僅保留具有蛋白編碼、轉(zhuǎn)錄調(diào)控、等生物學功能的位點。GWAS中交互作用分析方法匯總窮舉邏輯運算掃描檢測法(BOOST)、GBOOST、BiForce多因子降維(MDR)、MDRGPU、PLINK、FastEpistasis基于樹的上位效應關聯(lián)地圖(TEAM)、SHEsisEPI、EpiBlaster、iLOCi、EpiBlaster、等隨機隨機森林(randomforest)、randomjungle、epiForest貝葉斯上位效應關聯(lián)地圖法(BEAM)、PBEAM、epiMODE、BEAM2、BEAM3MegaSNPHunter、TRM、等啟發(fā)螞蟻算法(AntEpiSeeker)極大熵條件概率模型(MECPM)INTERSNP、SNPRuler、Gene-MDRSNPHarvester、SIXPAC、等檢索交互作用的方式一窮舉式檢索exhaustivesearch遍歷所有變量組合的情況，對交互作用進行檢測。設位點總數(shù)為p，檢測(k-1)階交互作用時，假設檢驗次數(shù)達到次；遍歷所有組合的交互作用，計算負擔極大。實際應用時需要極高的計算效率(硬件+算法)。GWAS中常采用保守的Bonferroni法校正多重比較，以控制一類錯誤。6檢索交互作用的方式二隨機式檢索stochasticsearch隨機地從總變量中挑選變量子集，考察這一部分變量的重要性。重復以上步驟，以保證能夠覆蓋到所有的變量。對同一變量計算其平均重要性獲得交互作用的估計值。在變量總數(shù)很低時，在可接受的重復次數(shù)內(nèi)，能夠覆蓋所有變量。若變量總數(shù)極大時，較低的重復次數(shù)，有降低把握度之嫌。7檢索交互作用的方式三啟發(fā)式檢索heuristicsearch首先通過某種判斷標準過濾掉絕大部分的噪音變量，比如無邊際效應的變量。僅對剩余變量探索交互作用。啟發(fā)式檢索可以大大降低計算負擔和統(tǒng)計學校正力度。然而，被過濾的變量可能會與其他變量存在交互作用。該檢索方式容易遺漏無邊際效應但存在交互作用的變量。89窮舉式檢索的方法10隨機式檢索的方法11啟發(fā)式檢索的方法GWAS中交互作用分析的挑戰(zhàn)(1)計算速度以50萬位點的數(shù)據(jù)為例。僅分析一階交互作用，需檢驗1.25×1011次。高性能工作站每次分析0.001秒，需日夜計算近4年。高階交互作用分析，負擔更為驚人！目前絕大多數(shù)方法或相應軟件仍然沒有勝任GWAS一階交互作用分析。多核并行計算。GPU取代CPU，利用顯卡計算，成為軟件開發(fā)的熱點。統(tǒng)計方法基因-環(huán)境-疾病，三者間關系復雜。參數(shù)模型難以準確刻畫三者關系。機器學習多采用非參數(shù)模型。其用于假設檢驗、避免過擬合、特征選擇的permutation、bootstrap或交叉驗證手段，反而加重計算負擔。多重比較采用Bonferroni法校正一類錯誤。位點間存在連鎖不平衡，過于保守。采用permutation控制總一類錯誤(FWER)或錯誤發(fā)現(xiàn)率(FDR)獲得可靠的檢驗水準，該手段加重計算負擔。GWAS中交互作用分析的挑戰(zhàn)(2)降維分析先篩選有交互作用的基因，再定位有交互作用的位點，可大幅降低計算負擔和校正力度。多以位點-位點交互作用分析為主。鮮有以基因為單位，進行基因-基因交互作用分析的方法。暫未發(fā)現(xiàn)基因水平上，交互作用分析的軟件。以基因為單位，如何評價基因-基因交互作用的效應？控制混雜不控制混雜因素可能導致假陽性或降低把握度。暫未發(fā)現(xiàn)控制混雜因素的GWAS交互作用分析方法。方法改進大多數(shù)方法僅支持一階交互作用分析，高階交互作用分析方法待完善。不少方法并非專門檢測交互作用，而是檢測聯(lián)合作用。GWAS中交互作用分析方法的評價(1)統(tǒng)計性能（文獻綜述和課題組模擬試驗的結(jié)論）多個研究者先后采用模擬試驗，對不同交互作用方法的一類錯誤和把握度進行了評價，涉及：MDR、隨機森林、BEAM、SNPHarvester、TEAM、SNPRuler、MECPM、SIXPAC、等數(shù)十個方法。BOOST表現(xiàn)突出！BOOST檢測的交互作用而非聯(lián)合作用。BOOST能夠一類錯誤控制。BOOST把握度相對最高。ZhangY(NatGenet,2007)Motsinger-ReifAA(GenetEpidemiol,2008)SuchestonL(BMCGenomics,2010)WangY(Bioinformatics,2011)ChenL(BMCGenomics,2011)ShagJ(BMCBioinformatics,2011)張汝陽(博士論文,2013)GWAS中交互作用分析方法的評價(2)計算速度（5000個樣本，1萬個位點的計算時間）BOOST表現(xiàn)出色。2011年，推出了GBOOST(GPU計算)。WanX(Bioinformatics,2010a,2010b)WangY(Bioinformatics,2011)535倍202倍65倍一、BOOSTBOOST的算法介紹BOOST的軟件應用16500,000個位點50萬次檢驗/10分鐘程序需要跑4.76年！171.1BOOST的算法介紹BOOST基于對數(shù)線性模型(log-linearmodel)推導出衡量交互作用的統(tǒng)計量采取先篩選后檢驗(screeningbefortesting)的策略檢測交互作用。篩選期：利用相對熵(kullback-leiblerdivergence)快速篩選潛在的交互作用。檢驗期：利用極大似然法，根據(jù)似然比卡方檢驗考察交互作用。181.1BOOST的算法介紹病例、對照人群中兩個位點不同基因型組合下頻數(shù)nijk的如表所示，總樣本量為n。191.1BOOST的算法介紹設表中單元格中的頻數(shù)服從Poisson分布，每個格子期望頻率是則每個格子的期望頻數(shù)是可以構(gòu)建似然函數(shù)其對數(shù)似然函數(shù)為：201.1BOOST的算法介紹含交互作用項的飽和模型(saturatedmodel,MS)極大似然值為不含交互作用項的均質(zhì)關聯(lián)模型(homogeneousassociationmodel,MH)極大似然值為：211.1BOOST的算法介紹上述兩個模型的似然值之差為22有無交互項對數(shù)線性模型似然值之差πijk與pijk的相對熵πijk

=nijk/nKirkwood疊加近似法(KSA)估計pijk23KSA法估計24二維表格一維表格計算示例1.1BOOST的算法介紹篩選期利用KSA法快速定義交互作用假設檢驗的統(tǒng)計量該統(tǒng)計量近似服從自由度為4的卡方分布。251.1BOOST的算法介紹利用KSA方法存在兩點優(yōu)勢通過列聯(lián)表中的觀察頻數(shù)可直接估計

，進而獲得交互作用統(tǒng)計量，無需迭代，計算速度快。與極大似然法相比，KSA法高估交互作用統(tǒng)計量，如此保證不會遺漏潛在的交互作用?？稍O置較低的閾值初篩交互作用。BOOST軟件中默認的閾值為30，對應交互作用的P為4.89×10-6。261.1BOOST的算法介紹271.1BOOST的算法介紹檢驗期針對初篩所得交互作用，利用廣義迭代法獲得

的極大似然估計值。通過似然比卡方檢驗獲得交互作用統(tǒng)計量，進一步檢驗交互作用。281.1BOOST的算法介紹優(yōu)點BOOST的最大特色是在篩選期利用不需要迭代的KSA法快速獲得衡量交互作用的一個統(tǒng)計量；在檢驗期進一步進行精確的假設檢驗。軟件采用C語言編寫，執(zhí)行效率高。缺點BOOST僅能檢測一階交互作用無法調(diào)整協(xié)變量，如年齡、性別、人群分層現(xiàn)象等。291.2BOOST的軟件應用下載BOOST軟件：

http://bioinformatics.ust.hk/BOOST.html數(shù)據(jù)文件名simdata.txt第一列是否病例986病例、1014對照其余列基因分析1000個SNP301.2BOOST的軟件應用參數(shù)文件filenamelist.txt指定數(shù)據(jù)文件的名稱simdata.txt可以同時讀入多個數(shù)據(jù)文件CMD窗口運行命令311.2BOOST的軟件應用結(jié)果輸出MarginalAssoc.txt存儲單位點邊際效應的結(jié)果InteractionRecords.txt存儲交互作用的結(jié)果653號SNP與929號SNP存在交互作用卡方統(tǒng)計量33.49，自由度為432二、多因子降維多因子降維的算法介紹多因子降維的軟件應用332.1多因子降維的算法介紹多因子降維(Multifactordimensionalityreduction,MDR)最早由RitchieMD于2001年提出，并于2003年推出相應的軟件。MDR是一種非參數(shù)(nonparametric)，無遺傳模型假設(geneticmodel-free)的數(shù)據(jù)挖掘技術之一。非參數(shù)是指MDR不涉及任何參數(shù)估計。無模型假設是指MDR不指定任何遺傳模型，即變量的線性或非線性形式。MDR廣泛應用在病例-對照研究中，用于檢測多變量間的交互作用。34352.1多因子降維的算法介紹將樣本分為10等份，用于10倍交叉驗證(10-foldcrossvalidation)。9份作為訓練集(trainingset)建立模型1份作為測試集(testingset)用于預測。從訓練集N個位點中選擇2個位點。基于兩個位點取值的組合構(gòu)建一個新變量。如此，便將變量從多維降低為一維。計算新變量每個取值下病例-對照比值。362.1多因子降維的算法介紹定義高、低危險組在平衡的(balanced)病例-對照研究中(病例、對照樣本量之比為1)，將病例-對照比值≥1者定義為高危險組(high-riskgroup)，將比值<1者定義為低危險組(low-riskgroup)。在非平衡的病例-對照研究中，可將該比值設定為總樣本中病例-對照樣本量之比。忽略無法計算比值的單元格。372.1多因子降維的算法介紹計算預測錯誤率高、低危險組分別被預測為病例、對照。利用剩余1份測試集計算當前模型的預測錯誤率(predictionerror)。依次利用其他9份訓練集重復步驟(1)至(5)，記錄10次交叉驗證中該模型預測錯誤率。計算10個預測錯誤率的均值作為該模型的最終預測錯誤率。382.1多因子降維的算法介紹選擇最佳模型重復步驟(1)至(5)獲得所有2位點模型的預測錯誤率。選擇預測錯誤率最低者作為最佳模型。據(jù)此，在總?cè)巳褐薪⒆罴涯Ｐ?92.1多因子降維的算法介紹對最佳模型的假設檢驗可借助Permutation完成。依次利用10個訓練集重復步驟(1)至(5)，記錄10次交叉驗證中該模型被定義為最佳模型的次數(shù)為m。隨機排列樣本是否患病的狀態(tài)，生成1000個Permutation數(shù)據(jù)集。對每個Permutation數(shù)據(jù)集，重復上述步驟再次考察10次交叉驗證中該模型被定義為最佳模型的次數(shù)m*。模型假設檢驗的P值為1000個Permutation數(shù)據(jù)集中m*≥m的比例。402.1多因子降維的算法介紹MDR遍歷所有2位點組合后，可進一步遍歷所有3位點組合、4位點組合，直至p位點組合，確定每種組合下的最佳模型。通過交叉驗證聯(lián)合Permutation的方式從p-1個最佳組合、環(huán)境變量中，確定最終模型。變量納入最終模型的原則：能夠使得預測錯誤率變低。412.1多因子降維的算法介紹MDR的優(yōu)點：MDR將多位點整合為一個新變量，同時檢測多個位點。MDR是一種非參數(shù)算法。若采用logistic回歸模型同時檢測多個交互作用，交互作用項的個數(shù)隨變量個數(shù)呈指數(shù)增長。HosmerandLemeshow(2000)建議logistic回歸模型中的變量個數(shù)不宜超過min(n1,n0)/10-1。其中n1，n0分別是病例、對照的樣本量。MDR不事先假設任何遺傳模型。實際上位點與疾病間存在復雜、未知的關系。MDR通過僅對最佳模型或最終模型進行假設檢驗，從而盡量最小化假陽性率(falsepositiverate)。422.1多因子降維的算法介紹MDR的弊端：若位點數(shù)極多，如GWAS的海量數(shù)據(jù)。MDR窮舉任意階的交互作用則是巨大的挑戰(zhàn)?？刹捎貌⑿杏嬎惴绞?，加速檢測速度。MDR所得最佳模型可能難以對應專業(yè)領域的解釋。若樣本量較小，高階交互模型的預測錯誤率將很大。采用最近鄰法(nearestneighbormethod)將空單元格歸類為大多數(shù)單元格所屬的類別。比如多數(shù)單元格都歸類為高危險組，則空單元格也歸類為高危險組。432.1多因子降維的算法介紹VelezDR(2007)對MDR進行改進，以適應不平衡病例-對照研究。兩種常用的校正不平衡的方式是再加權(quán)(reweighting)和再抽樣(resampling)。若當前觀察樣本中病例少、對照多。過抽樣有放回地在病例中重復抽樣從而不斷增加病例的樣本量；欠抽樣則在對照中隨機地抽取與病例等樣本量的對照。兩者均通過特定的方式使得病例對照樣本量相等。442.1多因子降維的算法介紹VelezDR的研究顯示當處理不平衡病例-對照數(shù)據(jù)時，閾值應該調(diào)整為原始樣本中病例-對照樣本量的比值，但是效果不佳。若調(diào)整閾值結(jié)合再抽樣技術，則顯著改善效果，且欠抽樣的效果要優(yōu)于過抽樣。而欠抽樣技術又不及利用測試集判斷訓練集模型時直接利用平衡精度函數(shù)(balancedaccuracyfunction)。平衡精度函數(shù)定義為(靈敏度+特異度)/2。452.1多因子降維的算法介紹平衡精度函數(shù)(靈敏度+特異度)/2462.1多因子降維的算法介紹MDR受位點主效應的影響，并非交互作用特異性檢測方法！472.2多因子降維的軟件應用安裝MDR包install.pakcages("MDR")library(MDR)載入數(shù)據(jù)集data(mdr1)482.2多因子降維的軟件應用利用平衡精度函數(shù)，通過交叉驗證獲得K個變量的最佳模型fit<-mdr.cv( mdr1, #分析數(shù)據(jù)集

K=2, #組成交互作用的位點數(shù)

cv=10, #10倍交叉驗證

ratio=1, #劃分高低危險組的閾值

equal=“HR“, #比值為1時歸為高危險組

genotype=c(0,1,2)) #基因型性編碼492.2多因子降維的軟件應用繪制最佳模型的高低危險組的圖plot(fit,data=mdr1)502.2多因子降維的軟件應用保存最佳模型的預測精度accuracy<-fit$'finalmodelaccuracy'保存最佳模型的變量組合loc<-fit$'finalmodel'512.2多因子降維的軟件應用通過permutation獲得最佳模型的經(jīng)驗P值perm<-permute.mdr( accuracy, #最佳模型的預測精度loci=loc, #最佳模型變量的位置N.permute=100,#Permutation次數(shù)method="CV“, #交叉驗證控制過擬合data=mdr1, #分析數(shù)據(jù)集cv=10, #10倍交叉驗證次數(shù)K=2, #至多2個位點形成交互作用equal="HR“, #病例對照之比為1時歸為高危險組genotype=c(0,1,2), #基因型的編碼LRT=FALSE) #不對交互項采用logistic回歸似然比卡方檢驗522.2多因子降維的軟件應用顯示經(jīng)驗P值perm$'PermutationP-value‘假設檢驗P值為0.01。532.2多因子降維的軟件應用GreeneCS等(2010)實現(xiàn)了MDR算法GPU的并行計算，推出MDRGPU。20萬點窮舉一階交互作用分析，耗時100小時/projects/mdr/files/mdrgpu/54三、隨機森林隨機森林的算法介紹隨機森林的軟件應用5556年齡基因年齡基因基因病例對照高年齡：突變基因保護作用低年齡：突變基因危險因素3.1隨機森林的分類樹3.1隨機森林的算法介紹隨機森林(randomforest,RF)最早由BreimanL于2001年提出。隨機森林算法的基本思想是隨機抽取部分樣本和部分變量作為訓練集建立多個分類樹，多個樹形成森林。將剩余數(shù)據(jù)作為測試集，逐個考察回歸樹中每個變量的重要性評分。變量在所有分類樹中重要性評分的平均值作為該變量的重要性指標。57583.1隨機森林的算法介紹以病例-對照設計下GWAS為例，詳細介紹隨機森林算法及變量篩選過程。假設有L個SNP。5960124533.1隨機森林的算法介紹(1)從總樣本中以bootstrap方式有放回地(withreplacement)或無放回地(withoutreplacement)隨機抽取部分樣本的部分變量作為測試集用于建立分類樹，將剩余樣本作為測試集，稱之為袋外數(shù)據(jù)(out-of-bagsample,OOBsample)。一般袋外數(shù)據(jù)的比例約為占總樣本量的1/3。設總變量數(shù)為L，分類樹上的變量數(shù)mtry設為sqrt(L)效果較好。613.1隨機森林的算法介紹(2)重復以上過程nTree次，得到多個分類樹，形成森林。對于GWAS而言，nTree一般設置至少10萬，盡量保證森林中的樹能夠覆蓋到所有變量。623.1隨機森林的算法介紹設根節(jié)點中病例、對照的頻率分別為f1，f0?；岵患兌?theginiimpurity)可定義為根據(jù)某變量生成左右兩個子節(jié)點，兩個子節(jié)點樣本量比例分別是

fL、

fR

。左右兩個子節(jié)點與根節(jié)點間不純度的差異可定義為633.1隨機森林的算法介紹子節(jié)點的上級被定義為其父節(jié)點(fathernode)。進一步，針對每個子節(jié)點重復上述過程進行劃分，直至子節(jié)點中純度不再增加。(如：不純度為0，即僅包含病例或者對照)。最后一級子節(jié)點被定義為終結(jié)點(terminalnode)。決策樹從根節(jié)點開始叉生為子節(jié)點，子節(jié)點再次叉生直至終節(jié)點。該過程類似于樹生長的過程。643.1隨機森林的算法介紹記

為變量

是否為節(jié)點

的叉生變量。每顆決策樹T上節(jié)點數(shù)為Nnode。對于變量

而言在單顆樹上的基尼重要性(giniimportance)是在整個森林中的基尼重要性是653.1隨機森林的算法介紹另一個衡量變量重要性的方式。利用袋外數(shù)據(jù)獲得當前決策樹的預測錯誤率(out-of-bagerrorrate,OOB-ER)，以衡量決策樹的預測精度。對決策樹上的待考察變量隨機排列取值后permutation，再次構(gòu)建決策樹，并計算袋外數(shù)據(jù)的預測錯誤率。兩個袋外數(shù)據(jù)預測錯誤率之差反映了該變量在決策樹中的重要性。森林中所有決策樹的重要性的均值被定義為該變量的原始重要性(rawimportance)。BreimanL研究指出基尼重要性評分和原始重要性評分效果近似，但前者計算負擔較輕。663.1隨機森林的算法介紹(3)通過隨機森林獲得所有變量或者特征(feature)的重要性評分(variableimportancemeasure,VIM)，按大小順序排序。重要性評分越大，表示變量越重要性，越有可能和結(jié)果存在關聯(lián)。一般利用隨機森林初篩重要變量。后期僅對重要變量進一步深入分析。673.1隨機森林的算法介紹(4)通過滑動窗口序貫向前特征選擇法(slidingwindowsequentialforwardfeatureselection,SWSFS)從所有變量中篩選最有可能和結(jié)果存在關聯(lián)的子集。按重要性評分從大到小依次選擇1個、2個、……、i個、……、L個變量，建立隨機森林，并利用袋外數(shù)據(jù)獲得當前森林的平均預測錯誤率。選擇袋外數(shù)據(jù)平均預測錯誤率最低的隨機森林所含變量。683.1隨機森林的算法介紹缺點分類樹中一個位點只有成為父節(jié)點，才能與其子節(jié)點形成交互作用。因此，只有位點存在主效應時才能被隨機森林發(fā)現(xiàn)與其他位點存在交互作用。YoshidaM等(2011)針對此缺陷，拓展算法，提出了專門用于篩選交互作用的方法：隨機交互森林(SNPInterForest)。其思想是將多個位點整合為一個新變量，對此新變量進行檢驗(聯(lián)合作用)。無論哪種方法均不是交互作用特異性的檢測方法！693.2隨機森林的軟件應用針對GWAS的專用工具RandomJunglehttp://www.imbs-luebeck.de/imbs/de/node/227R軟件Install.packages(“randomForest”)library(randomForest)703.2隨機森林的軟件應用讀入數(shù)據(jù)data(mdr1)指定預測變量和結(jié)局變量predictor<-mdr1[,2:26]pheno<-as.factor(mdr1[,1])713.2隨機森林的軟件應用隨機森林rf<-randomForest( x=predictor, #預測變量

y=pheno, #結(jié)果變量

mtry=5, #單顆樹上的變量

ntree=5000, #樹的數(shù)目

replace=TRUE, #有放回地bootstrap nPerm=1, #permutation次數(shù)

importance=TRUE)#輸出變量重要性評分723.2隨機森林的軟件應用獲得重要性評分vim<-importance(rf)vim<-vim[order(vim[,1],decreasing=T),]繪圖查看變量重要性評分varImpPlot(rf)73四、螞蟻優(yōu)化算法螞蟻優(yōu)化算法的算法介紹螞蟻優(yōu)化算法的軟件應用744.1螞蟻優(yōu)化算法的算法介紹螞蟻算法(antcolonyalgorithm)最早由DorigioM和GambardellaLM提出。螞蟻優(yōu)化算法(antcolonyoptimizationalgorithm,ACO)是優(yōu)化的螞蟻算法，用于病例-對照設計下的全基因組關聯(lián)研究中檢索交互作用。754.1螞蟻優(yōu)化算法的算法介紹自然界中螞蟻通過化學物質(zhì)-信息素(pheromones)進行交流溝通。螞蟻在尋在食物過程中，所經(jīng)之處會留下信息素。蟻群在前進的過程中選擇通過信息素較多的路徑，進而又留下更多的信息素，形成一個正反饋(positivefeedback)。螞蟻優(yōu)化算法則模擬蟻群尋找食物的機制，在海量數(shù)據(jù)中尋找交互作用。764.1螞蟻優(yōu)化算法的算法介紹(1)螞蟻優(yōu)化算法的機制。螞蟻優(yōu)化算法初始化時，設置多個人工螞蟻形成蟻群。螞蟻交流的信息素則是概率密度函數(shù)(probabilitydensityfunction,PDF)。變量集的交互作用越強，則會被賦予更大的概率(更多的信息素)。通過反復迭代的方式，更新變量的概率，實現(xiàn)信息素的更新。第i次迭代，選擇第k個變量的概率可已定義：774.1螞蟻優(yōu)化算法的算法介紹其中，

是第i次迭代時，第k個變量的信息素。

是變量k的先驗信息(priorinformation)，一般設為1，且信息素

均取相同值，意指初始化時同等對待所有變量。

用于指定螞蟻留下信息素的權(quán)重。根據(jù)上述概率公式，每個人工螞蟻選擇一個變量集。該變量集的聯(lián)合效應通過卡方檢驗進行估計，將此聯(lián)合效應的卡方統(tǒng)計量作為信息素

。784.1螞蟻優(yōu)化算法的算法介紹在第i+1次迭代時，變量集中的每個變量的信息素通過以下公式進行更新

是信息素的蒸發(fā)率

是第i+1次變量對應卡方統(tǒng)計量與第i次卡方統(tǒng)計量的差值。794.1螞蟻優(yōu)化算法的算法介紹(2)AntSeeker算法。AntSeeker算法由兩階段螞蟻優(yōu)化算法構(gòu)成。第一階段，AntSeeker算法根據(jù)信息素篩選潛在的變量集，其中所含變量數(shù)目高于構(gòu)成交互作用變量的數(shù)目。第二階段，AntSeeker算法在潛在變量集以及高信息素變量中窮舉變量的交互作用。804.1螞蟻優(yōu)化算法的算法介紹第一階段具體算法：設置兩種不同大小的變量集：大變量集、小變量集。依次對兩種類型變量集進行迭代操作。迭代次數(shù)分別是iItCountLarge，iItCountSmall。設定多個人工螞蟻，利用上述概率公式同時迭代信息素。每次迭代時，記錄蟻群中的對應聯(lián)合效應卡方統(tǒng)計量最大的變量集。迭代結(jié)束后可以分別獲得iItCountLarge，iItCountSmall個大變量集、小變量集，以及信息素較大的變量(大于指定閾值)。814.1螞蟻優(yōu)化算法的算法介紹第二階段具體算法：針對iItCountLarge+iItCountSmall個變量集，分別在變量集內(nèi)窮舉計算交互作用。另外，在信息素較大的變量內(nèi)部同樣窮舉計算交互作用。報告小于界值的交互作用。824.1螞蟻優(yōu)化算法的算法介紹(3)最小化假陽性率的方式。針對上步所篩選的變量集及高信息素變量，保留無重疊的變量集，或有重復單聯(lián)合卡方效應P值較低者。834.1螞蟻優(yōu)化算法的算法介紹缺點檢測聯(lián)合作用，并非交互作用特異性檢測方法以統(tǒng)計量大小來設置選擇概率，這是降低計算負擔的方式，但是容易使得檢索范圍限制在局部范圍。844.2螞蟻優(yōu)化算法的軟件應用軟件下載/~romdhane/AntEpiSeeker/index.html數(shù)據(jù)文件data0.txt85參數(shù)文件parameters.txtiAntCount 1000 //numberofantsiItCountLarge150//numberofiterationsforthelargehaplotypesiItCountSmall300 //numberofiterationsforthesmallhaplotypesalpha 1//weightgiventopheromonedepositedbyantsiTopModel1000//numberoftoprankinghaplotypesinthefirststageiTopLoci 200//numberoflociwithtoprankingpheromoneinthefirststagerou 0.05 //evaporationrateinAntColonyOptimizaionphe 100 //initialpheromonelevelforeachlocusLargehapsize6 //sizeofthelargehaplotypesSmallhapsize3 //sizeofthesmallhaplotypesiEpiModel 2 //numberofSNPsinanepistaticinteractionPvalue0.01//pvaluethreshold(afterBonferronicorrection)INPFILE data0.txt//inputfilenameforcase-controlgenotypedataOUTFILE result.txt//outputfilenamefordetectedepistaticinteractions864.2螞蟻優(yōu)化算法的軟件應用結(jié)果輸出文件result.txt結(jié)果實際上是聯(lián)合效應，并非真正意義的交互作用87中國人群肺癌研究交互作用分析實例全基因組基因-基因交互作用研究ChuM#,ZhangR#,ZhaoY#,WuC,GuoH,ZhouB,LuJ,ShiY,DaiJ,JinG,MaH,DongJ,WeiY,WangC,GongJ,SunC,ZhuM,QiuY,WuT,HuZ,LinD,ShenH*,ChenF*.Agenome-widegene-geneinteractionanalysisidentifiesanepistaticgenepairforlungcancersusceptibilityinHanChinese.Carcinogenesis.2014Mar;35(3):572-7.全基因組基因-環(huán)境交互作用研究ZhangR#,ChuM#,ZhaoY#,WuC,GuoH,ShiY,DaiJ,WeiY,JinG,MaH,DongJ,YiH,BaiJ,GongJ,SunC,ZhuM,WuT,HuZ,LinD,ShenH*,ChenF*.Agenome-widegene–environmentinteractionanalysisfortobaccosmokeandlungcancersusceptibility.Carcinogenesis.2014Jul;35(7):1528-35.實例1：基因-基因交互作用研究1.1.1研究設計（三階段病例-對照設計）一致結(jié)果合計:6377病例和7015對照樣本來自南京和北京1473病例1962對照858病例1115對照南京子研究

北京子研究2331病例3077對照AffymetrixGenome-WideHumanSNPArray6.01534病例1489對照一期驗證樣本來自北京SequenomMassARRAY二期驗證2512病例2449對照樣本來自沈陽和廣州TaqManassay1.1.2樣本收集樣本收集病例組對照組健康個體性別、年齡和居住地區(qū)與病例頻數(shù)匹配新發(fā)病例病理組織學或細胞學確診的肺癌新發(fā)病例1.1.3流行病學調(diào)查對調(diào)查對象使用統(tǒng)一的調(diào)查表進行面訪問卷調(diào)查。內(nèi)容包括一般情況、個人健康狀況、吸煙等生活方式。吸煙定義一生中吸煙平均每天<1支且延續(xù)不到1年者為不吸煙者，否則為吸煙者；在納入研究前戒煙超過1年的吸煙者為過去吸煙者。累積吸煙劑量“包-年”按(每日吸煙量/20×吸煙年數(shù))進行計算。目前仍在吸煙者和過去吸煙者根據(jù)對照人群累積吸煙量中位數(shù)進一步劃分為輕度吸煙者和重度吸煙者。1.2.1質(zhì)量控制（1）-刪除樣本2383病例和3160對照13例樣本基因分型成功率<95%2375病例和3155對照2375病例和3148對照2345病例和3089對照2334病例和3078對照2331病例和3077對照7例樣本染色體性別與調(diào)查性別不一致89例樣本存在重復、污染或親緣關系22例樣本雜合率異常4例樣本人群分層離群值1.2.1質(zhì)量控制（2）-刪除位點906703個位點905119個位點872522個位點591625個位點591370個位點1584個位點不在常染色體區(qū)域32597個位點子研究或總樣本分型率<95%280897個位點MAF<0.05或PHWE

<1e-05255個位點分型圖不合格1.3統(tǒng)計分析流程圖北京子研究PLRM≤1.00×10-4沈陽及廣州地區(qū)4961獨立樣本：1對交互作用位點總?cè)巳篜LRM≤2.86×10-13合計13392例樣本：rs2562796-rs16832404，PLRM

=1.03×10-13

北京地區(qū)3023獨立樣本：4對交互作用位點logistic回歸(LRM)檢驗493177對交互作用驗證要求：效應一致、分型清晰、無LDBOOST軟件篩選1.75×1011交互作用位點經(jīng)質(zhì)控的南京、北京兩個子研究：5408個病例、591370個位點

PBOOST≤1.00×10-6南京子研究PLRM≤1.00×10-4第2階段驗證GWAS篩選第1階段驗證1.4.1統(tǒng)計分析：基本信息描述變量GWAS篩選期驗證I驗證II南京子研究北京子研究北京地區(qū)沈陽、廣州地區(qū)病例對照病例對照病例對照病例對照(n=1473)(n=1962)(n=858)(n=1115)(n=1534)(n=1489)(n=2512)(n=2449)年齡60.1±10.359.4±9.760.0±10.262.5±9.258.1±9.753.9±13.259.6±12.458.5±12.5性別男10571214654872102999215391357女4167482042435054979731092吸煙狀況吸煙741636511447726553819726戒煙168838614318619421226不吸564124326152562291712721497吸煙量41.4±26.931.0±20.344.8±29.632.4±19.939.6±25.325.9±18.435.8±17.125.5±16.4組織類型鱗狀細胞癌421401587910腺癌8964087801095小細胞癌12949100370其他27067137表1.1位點rs16832404、rs2562796在篩選、驗證階段中交互作用分析的結(jié)果1.4.2統(tǒng)計分析：本研究識別的交互作用研究rs16832404A/Gars2562796G/Ta交互作用

OR(95%CI)bPb

OR(95%CI)bPb

OR(95%CI)bPbGWAS0.39(0.32-0.46)9.56×10-250.60(0.54-0.68)2.38×10-172.58(2.24-2.97)1.37×10-39驗證I0.96(0.80-1.16)6.87×10-10.95(0.82-1.09)4.50×10-11.17(0.99-1.38)6.37×10-2驗證II0.88(0.75-1.02)8.75×10-20.97(0.87-1.08)6.22×10-11.21(1.06-1.38)4.61×10-3綜合分析0.81(0.73-0.88)3.67×10-60.90(0.85-0.96)2.24×10-31.33(1.23-1.43)1.03×10-13研究rs16832404rs2562796

OR(95%CI)aPa

OR(95%CI)aPaGWAS1.06(0.97-1.17)2.12×10-11.03(0.95-1.12)4.94×10-1驗證I1.08(0.96-1.21)1.98×10-11.00(0.90-1.12)9.66×10-1驗證II1.01(0.92-1.10)8.67×10-11.05(0.97-1.14)2.13×10-1綜合分析1.05(0.99-1.11)9.53×10-21.05(1.00-1.10)5.62×10-2表1.2位點rs16832404、rs2562796在篩選、驗證階段中交互作用分析的結(jié)果表1.3位點rs16832404、rs2562796在篩選、驗證階段中單位點分析的結(jié)果a主要/次要等位基因；b調(diào)整年齡、性別、吸煙“年-包”數(shù)、人群分層特征向量。a調(diào)整年齡、性別、吸煙“年-包”數(shù)、人群分層特征向量。1.4.3統(tǒng)計分析：交互作用模式森林圖展示5.65×10-1GWASGGGTTT驗證IGGGTTT驗證IIGGGTTT綜合分析GGGTTT攜帶位點rs2562796不同基因型的人群0.36(0.28,0.44)1.04(0.90,1.21)2.55(2.10,3.09)0.94(0.76,1.17)1.14(0.95,1.37)1.31(0.99,1.73)0.91(0.76,1.08)1.06(0.92,1.22)1.21(0.97,1.52)0.80(0.72,0.88)1.08(1.00,1.18)1.41(1.26,1.57)位點rs16832404的效應1.21×10-194.36×10-216.04×10-11.58×10-16.06×10-22.61×10-14.33×10-18.70×10-21.71×10-56.43×10-21.29×10-9P值OR(95%CI)1.4.4統(tǒng)計分析：區(qū)域交互作用信號圖選擇rs2562796、rs16832404位點上下游30kb區(qū)域內(nèi)170個高質(zhì)量填補位點。將交互作用P值進行-log10對數(shù)變換，負對數(shù)值超過12者被截斷。P值由大到小依次用“黃、綠、紅”三色標記。本研究識別的交互作用位點用藍色方框標記。提示：所識別的交互作用位點附件區(qū)域存在強烈的交互作用信號！rs16832404rs25627961.4.5統(tǒng)計分析：分層分析分層變量rs2562796rs16832404交互作用OR(95%CI)a

PaOR(95%CI)a

PaOR(95%CI)a

Pa年齡≤60歲0.91(0.83-0.99)3.31×10-20.87(0.76-0.98)2.23×10-21.26(1.14-1.39)5.91×10-6>60歲0.90(0.81-0.99)2.44×10-20.74(0.64-0.85)1.58×10-51.42(1.27-1.59)1.57×10-9性別男0.86(0.80-0.94)4.08×10-40.78(0.69-0.88)3.33×10-51.40(1.27-1.54)6.37×10-12女0.98(0.88-1.09)6.98×10-10.84(0.72-0.98)2.64×10-21.22(1.08-1.38)1.18×10-3吸煙狀況吸煙0.82(0.74-0.91)1.37×10-40.69(0.60-0.80)8.65×10-71.58(1.40-1.79)2.90×10-13戒煙0.76(0.62-0.94)1.23×10-20.77(0.56-1.04)8.86×10-21.46(1.13-1.89)3.61×10-3不吸0.98(0.89-1.07)5.90×10-10.88(0.77-0.99)4.05×10-21.18(1.06-1.30)1.68×10-3吸煙量≤25年-包0.95(0.88-1.03)1.96×10-10.83(0.75-0.93)7.66×10-41.26(1.15-1.37)2.95×10-7>25年-包0.79(0.70-0.90)1.88×10-40.73(0.61-0.87)2.99×10-41.57(1.35-1.82)1.84×10-9病理組織類型鱗狀細胞癌0.96(0.87-1.05)3.70×10-10.85(0.74-0.98)2.10×10-21.31(1.18-1.46)1.17×10-6腺癌0.85(0.79-0.92)6.21×10-50.77(0.69-0.86)3.92×10-61.41(1.29-1.54)1.16×10-13小細胞癌0.92(0.79-1.06)2.48×10-10.78(0.62-0.97)2.43×10-21.25(1.04-1.50)1.58×10-2其他1.22(0.96-1.55)1.06×10-11.19(0.85-1.66)3.09×10-10.83(0.63-1.11)2.13×10-1表1.4位點rs16832404、rs2562796在不同亞組人群中交互作用分析的結(jié)果1.4.6統(tǒng)計分析：交互作用改善預測效果模型1：年齡、性別、吸煙“年-包”數(shù)模型2：年齡、性別、吸煙“年-包”數(shù)、位點rs2664241與rs2562796評價指標：category-freenetreclassificationimprovement(NRI)新模型將病例預測概率提高的比例新模型將對照預測概率降低的比例本研究中增加交互作用位點模型預測效果改善了10.52%。對象病例的預測概率合計降低提高病例3,1283,0956,223對照3,8283,0666,8941.5.1本研究的特點：病例對照匹配：新發(fā)確診病例與健康對照特征匹配，控制混雜因素。質(zhì)量控制嚴格：僅保留符合標準的樣本和基因位點，避免虛假關聯(lián)。研究尺度寬泛：首次從全基因組水平探索基因-基因一階交互作用。檢索范圍全面：采用窮舉式分析策略，覆蓋所有可能的基因組合。先初篩后檢驗：BOOST快速初篩，多因素logistic回歸模型檢驗。校正多重比較：保守的Bonferroni法校正一類錯誤：2.86×10-13。獨立人群驗證：多個獨立人群驗證結(jié)果，再次控制假陽性。內(nèi)部驗證：南京子研究+北京子研究；外部驗證：第1階段驗證；外部驗證：第2階段驗證。樣本量相對大：總樣本量達13392：病例6377+對照7015。研究結(jié)論可信：列聯(lián)表中的頻數(shù)均較大，結(jié)果穩(wěn)??；分層分析亞組人群中交互作用仍有統(tǒng)計學意義；所識別位點附近區(qū)域存在強烈的交互作用信號。1.5.2本研究的不足高階交互作用：本研究僅考慮一階交互作用。高階交互作用值得關注。BOOST僅能分析一階交互作用。我們改進BOOST，從理論上實現(xiàn)了高階交互作用的快速初篩方法。相乘交互作用：本研究關注的是常見的相乘交互作用。相加交互作用值得進一步探索。多重比較校正：由于位點間并非完全獨立，Bonferroni法校正往往過于嚴格，易遺漏潛在的交互作用。適當放松校正力度，可能獲得更多交互作用信號。復雜關聯(lián)模式：一般假定存在線性關聯(lián)。如，森林圖中OR值線性變化。非線性的關聯(lián)值得后續(xù)研究。實例2：基因-環(huán)境交互作用研究2.1研究設計（二階段病例-對照設計）一致結(jié)果合計:3865病例和4566對照樣本來自南京和北京1473病例1962對照858病例1115對照南京子研究

北京子研究2331病例3077對照AffymetrixGenome-WideHumanSNPArray6.01534病例1489對照一期驗證樣本來自北京SequenomMassARRAY基因-環(huán)境交互作用研究所采用的數(shù)據(jù)同基因-基因交互作用研究。很遺憾！沈陽和廣州標本不足以再次DNA分型。因此，無第2階段驗證。2.2統(tǒng)計分析流程圖北京子研究PLRM≤1.00×10-2總?cè)巳篜LRM≤5.00×10-5合計8431：例樣本：成功驗證2個與吸煙存在交互作用的位點rs1316298-吸煙，PLRM

=6.73×10-6rs4589502-吸煙，PLRM

=3.84×10-6

北京地區(qū)3023獨立樣本：6個位點與吸煙存在交互作用logistic回歸(LRM)檢驗591370對交互作用驗證要求：效應一致、分型清晰、無LD經(jīng)質(zhì)控的南京、北京兩個子研究：5408個病例、591370個位點南京子研究PLRM≤1.00×10-2獨立人群驗證GWAS篩選2.3.1統(tǒng)計分析：交互作用P值曼哈頓圖將591,370個基因-吸煙交互作用P值，進行l(wèi)og10對數(shù)變換，繪制曼哈頓圖。紅線表示P=5.0×10-5。染色體2、6、14、15、18等區(qū)域存在交互作用信號。2.3.2統(tǒng)計分析：交互作用-log10P值QQ圖繪制591,370個交互作用-log10P值實際分位數(shù)與理論分位數(shù)的QQ圖。實際數(shù)與理論數(shù)較為符合，貼近45度對角線?；蚺蛎浺蜃?genomic-controlinflation,λ)為1.021。結(jié)果提示：人群分層控制較好。陽性結(jié)果并非人群分層混雜因素導致。2.3.3統(tǒng)計分析：基因-環(huán)境交互作用分析研究基因吸煙交互作用OR(95%CI)aPaOR(95%CI)aPaOR(95%CI)aPars1316298-吸煙b南京子研究1.35(1.12,1.65)2.17×10-034.03(3.12,5.22)2.46×10-260.65(0.50,0.85)1.28×10-03北京子研究1.20(0.92,1.55)1.78×10-014.04(2.95,5.52)2.73×10-180.61(0.44,0.84)3.09×10-03GWAS總樣本1.25(1.08,1.45)3.53×10-033.81(3.15,4.61)3.42×10-430.66(0.55,0.81)4.15×10-05驗證階段1.05(0.88,1.27)5.74×10-014.34(3.45,5.47)1.53×10-350.64(0.50,0.83)8.87×10-04綜合分析1.11(1.00,1.23)5.67×10-023.99(3.48,4.58)6.71×10-860.71(0.62,0.83)6.73×10-06rs4589502-吸煙

c南京子研究0.68(0.53,0.88)3.39×10-032.80(2.22,3.55)6.24×10-181.69(1.20,2.36)2.33×10-03北京子研究0.75(0.54,1.05)9.10×10-022.69(2.04,3.54)1.69×10-121.76(1.16,2.68)8.37×10-03GWAS總樣本0.70(0.58,0.85)3.95×10-042.66(2.24,3.16)2.40×10-291.72(1.34,2.22)2.61×10-05驗證階段0.78(0.62,0.98)2.93×10-023.26(2.63,4.03)1.17×10-271.39(1.01,1.91)4.40×10-02綜合分析0.74(0.64,0.84)1.14×10-053.01(2.66,3.42)4.07×10-661.55(1.29,1.87)3.84×10-06表2.1肺癌基因-吸煙交互作用分析的所識別的2個位點a調(diào)整年齡、性別、人群分層特征向量。b位于14q22.1；c位于15q22.32。2.3.4統(tǒng)計分析：交互作用模式森林圖展示CT+TTI2=94.3%I2=79.8%TTCCTC+CC攜帶不同基因型的人群4.25(3.37,5.37)吸煙的OR(95%CI)4.21(3.60,4.93)3.14(2.74,3.59)2.56(2.15,3.04)2.56×10-34P3.91×10-711.38×10-621.49×10-26rs1316298rs45895023.00以下人群吸煙的危害更大，應控煙、禁煙預防肺癌發(fā)生攜帶位點rs1316298TT基因型攜帶位點rs4589502CT/TT基因型2.3.5統(tǒng)計分析：基因-吸煙聯(lián)合作用分析研究基因突變a吸煙rs1316298-吸煙rs4589502-吸煙OR(95%CI)bP

bOR(95%CI)bP

bGWAS否否參考組參考組否是3.95(3.25,4.80)2.95×10-432.68(2.26,3.18)2.09×10-29是否1.30(1.08,1.57)6.25×10-030.68(0.55,0.84)5.05×10-04是是2.96(2.42,3.62)4.21×10-263.20(2.58,3.96)1.69×10-26Pinteractionb0.58(0.45,0.74)9.48×10-061.76(1.33,2.33)8.94×10-05驗證階段否否參考組參考組否是4.26(3.36,5.39)2.88×10-333.28(2.65,4.07)1.21×10-27是否1.02(0.82,1.27)8.83×10-010.76(0.59,0.97)2.83×10-02是是2.79(2.19,3.56)1.70×10-163.41(2.58,4.50)6.62×10-18Pinteractionb0.65(0.47,0.88)5.66×10-031.39(0.98,1.96)6.60×10-02綜合分析否否參考組參考組否是4.09(3.55,4.71)1.31×10-843.03(2.67,3.44)4.48×10-66是否1.13(0.99,1,29)7.36×10-020.72(0.62,0.83)1.78×10-05是是2.99(2.58,3.46)1.87×10-483.39(2.89,3.99)1.13×10-49Pinteractionb0.65(0.54,0.78)2.23×10-061.56(1.27,1.92)2.64×10-05表2.2肺癌