第二軍醫(yī)大學(xué)衛(wèi)生統(tǒng)計(jì)學(xué)教研室 (200433)　

宋嘉麒　金志超△　賀　佳△

基因-疾病關(guān)聯(lián)研究常利用單核苷酸多態(tài)性(single nucleotide polymorphisms,SNPs)篩查與疾病相關(guān)的基因突變位點(diǎn),然而由于樣本量和統(tǒng)計(jì)效能的問(wèn)題,常常導(dǎo)致假陰性,利用meta分析對(duì)原始研究進(jìn)行定量合并,能夠有效地提高統(tǒng)計(jì)效能。

在進(jìn)行基因-疾病關(guān)聯(lián)研究meta分析時(shí),常常遇到遺傳模型的假定問(wèn)題。遺傳模型的假定問(wèn)題既存在于原始文獻(xiàn)中,也存在于meta分析定量合并時(shí)。如果原始文獻(xiàn)中研究者主觀地進(jìn)行了遺傳模型的假定,如假定為顯性或者隱形遺傳模型,在其結(jié)果中有可能只給出特定遺傳模型下的效應(yīng)量(OR值),從而使meta分析作者無(wú)法同時(shí)獲得三種基因型的頻數(shù),最終導(dǎo)致其在定量合并時(shí)也不得不采用原始文獻(xiàn)中的遺傳模型,或者將該研究剔除。Minelli,Salanti等人對(duì)遺傳模型的假定問(wèn)題進(jìn)行了較為深入的研究[1-2]。Salanti等人的研究已經(jīng)能夠較好的處理上述兩種遺傳模型假定的情況[2],利用貝葉斯方法分析遺傳模型假定而存在的基因型合并問(wèn)題。但是在實(shí)際研究當(dāng)中,由于Salanti等人提出的貝葉斯方法較為復(fù)雜,且對(duì)于常用的回顧性似然未提供代碼(回顧性似然更適合采用病例對(duì)照研究的基因-疾病關(guān)聯(lián)研究[3]),使得此方法在實(shí)際應(yīng)用中使用的并不廣泛。同時(shí),Salanti等人在其研究中并未對(duì)遺傳模型的選擇提供結(jié)論。

本研究將在Minelli和Salanti等人研究的基礎(chǔ)上,利用回顧性似然和貝葉斯方法處理基因-疾病關(guān)聯(lián)研究meta分析中存在的原始文獻(xiàn)基因型合并的問(wèn)題,構(gòu)建一個(gè)遺傳模型選擇指示概率,以數(shù)據(jù)驅(qū)動(dòng)估計(jì)各種典型遺傳模型存在的相應(yīng)概率,最后提供相應(yīng)的JAGS代碼和R代碼用以實(shí)現(xiàn)上述過(guò)程,為其他類(lèi)似研究提供方法學(xué)支持。

方　　法

1.基本原理

常見(jiàn)的進(jìn)行基因關(guān)聯(lián)研究的原始文獻(xiàn)會(huì)提供如表1所示的病例組和對(duì)照組的三種基因型的頻率,假定三種基因型分別時(shí)AA、Aa和aa。

表1　常見(jiàn)基因關(guān)聯(lián)研究的數(shù)據(jù)表現(xiàn)形式

表1中cni(n=1,2,3)表示納入的某原始研究病例組三種基因型頻率,病例組樣本量為ci,tni(n=1,2,3)為對(duì)照組三種基因型頻率,對(duì)照組樣本量為ti。病例組和對(duì)照組基因型頻率服從多項(xiàng)分布,可用(1)式和(2)式表示:

(1)

(2)

(3)

式中系數(shù)λ表示遺傳模型,λ=0表示隱性遺傳模型,λ=0.5表示共顯性遺傳模型,λ=1表示顯性遺傳模型。式(3)是在病例組和對(duì)照組所有基因型都可以從納入文獻(xiàn)中獲取的情況下使用,但在納入的原始文獻(xiàn)中,作者在分析時(shí)有可能對(duì)遺傳模型進(jìn)行假定,得到的是一個(gè)對(duì)基因型頻數(shù)進(jìn)行合并了的四格表,如表2所示。

表2　假定遺傳模型下合并基因型的基因關(guān)聯(lián)研究數(shù)據(jù)表現(xiàn)形式

(4)

在估計(jì)φi時(shí),假設(shè)感興趣的等位基因A的頻率為κi,根據(jù)哈代-溫伯格平衡(Hardy-Weinberg equilibrium),φi=2(1-κi)/(2-κi)。相比利用極大似然估計(jì)等方法,貝葉斯方法在估計(jì)合并效應(yīng)量θ時(shí)更方便,能夠方便地借助外部信息,如哈代-溫伯格平衡。

2.實(shí)例

實(shí)例數(shù)據(jù)來(lái)自Huang等發(fā)表的關(guān)于ACE I/D基因多態(tài)性與腎移植預(yù)后關(guān)系的系統(tǒng)綜述[4],本實(shí)例選取了腎移植所致的急慢性排斥反應(yīng),共納入11項(xiàng)研究,其中9項(xiàng)提供了完全的基因型頻數(shù),2個(gè)研究存在合并的基因型頻數(shù)(II和ID基因型存在合并)[5-6],其基本信息如表3所示。

表3　 ACE I/D基因多態(tài)性對(duì)腎移植所致的急慢性排斥反應(yīng)納入文獻(xiàn)基本信息

3.軟件實(shí)現(xiàn)

本研究采用貝葉斯常用軟件JAGS,并利用R軟件中的R2jags包調(diào)用執(zhí)行。由于文章篇幅所限,本例中只提供關(guān)鍵的JAGS代碼,R代碼可通過(guò)郵件向通訊作者索取。JAGS完整代碼如下:

model{

#Complete studies

for(i in 1:9){

#Multinomial Likelihoods.matrix ′ca′ and ′co′ is 9*3

ca[i,1:3]～dmulti(p.case[i,1:3],tcase[i])

co[i,1:3]～dmulti(p.cont[i,1:3],tcont[i])

# 公式(3)

sum1[i]<-p.cont[i,1]+exp(lambda*theta[i])*p.cont[i,2]+exp(theta[i])*p.cont[i,3]

p.case[i,1]<-p.cont[i,1]/sum1[i]

p.case[i,2]<-exp(lambda*theta[i])*p.cont[i,2]/sum1[i]

p.case[i,3]<-exp(theta[i])*p.cont[i,3]/sum1[i]

for(k in 1:3){

p.cont[i,k]～dunif(0,1)

}

}

#Merged studies

for (i in 1:2){

#Multinomial Likelihoods.matrix ′ca.merg′ and ′co.merg′ is 9*2

ca.merg[i,1:2]～dmulti(p.casemerg[i,1:2],ca.n.merge[i])

co.merg[i,1:2]～dmulti(p.contmerg[i,1:2],co.n.merge[i])

#公式(4)

p.contr[i,1]<-p.contmerg[i,1]

p.contr[i,2]<-p.contmerg[i,2]*f[i]

p.contr[i,3]<-p.contmerg[i,2]*(1-f[i])

sum2[i]<-p.contmerg[i,1]+exp(lambda*theta[i+9])*f[i]*p.contmerg[i,2]+exp(theta[i+9])*(1-f[i])*p.contmerg[i,2]

p.case1[i,1]<-p.contr[i,1]/sum2[i]

p.case1[i,2]<-exp(lambda*theta[i+9])*p.contr[i,2]/sum2[i]

p.case1[i,3]<-exp(theta[i+9])*p.contr[i,3]/sum2[i]

for(k in 1:2){

p.contmerg[i,k]～dunif(0,1)

p.casemerg[i,k]～dunif(0,1)

}

#H-W平衡

f[i]<-2*(1-pa[i])/(2-pa[i])

pAA[i]<-(1-pa[i])*(1-pa[i])

n1[i]～dbin(pAA[i],co.n.merge[i])

pa[i]～dbeta(1,1)

}

for(i in 1:11){

theta[i]～dnorm(mean,prec)

}

# Priors information for effects and heterogeneity

mean～dnorm(0,0.0001)

tau ～ dnorm(0,1)I(0,) #半正態(tài)分布

prec=1/(tau*tau)

#Probability of genetic model,0:recessive model,0.5 co-dominant,1:dominant

d<-c(0,0.5,1)

p<-c(1/3,1/3,1/3)

K～dcat(p[])

lambda<-d[K]

#Results of interest

OR2<-exp(mean*lambda)

OR3<-exp(mean)

#probmodel[1]:recessive,probmodel[2]:co-dominant,probmodel[3]:dominant

for(g in 1:3){

probmodel[g]<-equals(K,g)

}

}

上述代碼的關(guān)鍵部分為涉及合并基因型的回顧性似然構(gòu)建和編碼,對(duì)于τ還可以有其先驗(yàn)分布,如dgamma(0.001,0.001)和dunif(0,10)等。

結(jié)　　果

上述實(shí)例中,相對(duì)于基因型ACE D/D,攜帶有基因型I/I和基因型I/D的腎移植患者發(fā)生排斥反應(yīng)風(fēng)險(xiǎn)無(wú)統(tǒng)計(jì)學(xué)意義,其OR值和95%置信區(qū)間分別為0.68 (0.41,1.02)和0.50(0.19,1.04),用以表示遺傳模型的系數(shù)λ及其95%置信區(qū)間為0.53 (0.50,1.00),其后驗(yàn)概率密度分布如圖1所示。三種遺傳模型的概率分別為隱性遺傳模型為0%,顯性遺傳模型6.2%,共顯性遺傳模型為93.8%。研究間變異τ2=0.98。

圖1　表示遺傳模型的系數(shù)λ的后驗(yàn)概率密度

另外,我們又選取了Ragland等發(fā)表的關(guān)于UCHL1 S/Y基因多態(tài)性與帕金森氏病發(fā)病風(fēng)險(xiǎn)關(guān)系的系統(tǒng)綜述[7],共納入18項(xiàng)研究,全部提供了完全的基因型頻數(shù)。首先用無(wú)遺傳模型約束的貝葉斯Meta分析方法[8-9]估計(jì)OR值,然后無(wú)放回地隨機(jī)抽取2個(gè)研究,合并S/Y和Y/Y基因型頻數(shù),并用本研究方法估計(jì)OR值,并重復(fù)9次。結(jié)果如表4所示,相對(duì)于基因型S/S,基因型Y/Y的OR值為0.876～0.912,均值為0.895,與納入全部研究的估計(jì)結(jié)果0.894相近,基因型S/Y的OR值為0.649～0.734,均值為0.704,與納入全部研究的估計(jì)結(jié)果0.704相近。

表4　隨機(jī)合并S/Y和Y/Y基因頻數(shù)及全部研究OR值估計(jì)結(jié)果

討　　論

本研究采用了貝葉斯方法處理了基因-疾病關(guān)聯(lián)研究meta分析中原始文獻(xiàn)中存在遺傳模型假定而帶來(lái)基因型合并問(wèn)題。在meta分析過(guò)程中不假定基因

型,利用完整的研究和哈代-溫伯格平衡估計(jì)合并基因型中每個(gè)基因型的頻率,由于納入了所有原始文獻(xiàn),未剔除存在合并基因型的研究,能夠提高統(tǒng)計(jì)檢驗(yàn)效能。利用Ragland等人實(shí)例提供的完全數(shù)據(jù)表明本文方法估計(jì)OR值的穩(wěn)定性、準(zhǔn)確性以及精確性均較高。

同時(shí),利用本研究的方法,在無(wú)充分生物學(xué)機(jī)制的基礎(chǔ)上,能夠利用數(shù)據(jù)驅(qū)動(dòng)估計(jì)幾種經(jīng)典遺傳模型的概率。利用meta分析的方法估計(jì)遺傳模型對(duì)于構(gòu)建慢性傳染性疾病的發(fā)病或預(yù)后模型具有重要意義[2]。在宏觀環(huán)境危險(xiǎn)因素的基礎(chǔ)上,利用合理的遺傳模型加入基因?qū)用娴挠绊懸蛩?能夠使預(yù)測(cè)模型更為精確,能夠促進(jìn)精準(zhǔn)醫(yī)學(xué)的開(kāi)展。

本研究存在的主要缺陷是哈代-溫伯格平衡在少量的研究中并不一定成立,因此需要考慮在哈代-溫伯格平衡不成立條件下的似然函數(shù)形式[10]。另外,本研究未對(duì)影響合并結(jié)果的其他因素如種族、檢測(cè)方法、性別等考慮在內(nèi),可以進(jìn)一步考慮在貝葉斯框架下使用meta回歸,使得結(jié)果更為合理。