鄭州大學(xué)公共衛(wèi)生學(xué)院計(jì)算機(jī)與衛(wèi)生統(tǒng)計(jì)學(xué)教研室(450001)

羅子瀟 楊永利 賈曉燦 王玉平 施學(xué)忠△

全基因組關(guān)聯(lián)分析(genome-wide association study，GWAS)在識(shí)別疾病的常見(jiàn)變異方面取得了巨大進(jìn)展，目前已經(jīng)報(bào)道上萬(wàn)個(gè)單核苷酸多態(tài)性(single nucleotide polymorphism，SNP)位點(diǎn)與數(shù)百種復(fù)雜疾病存在關(guān)聯(lián)，這為表型變異的遺傳基礎(chǔ)提供了前所未有的視角[1-2]。但GWAS是基于個(gè)體水平基因型和表型數(shù)據(jù)的分析，因此需要更有效的方法基于匯總統(tǒng)計(jì)數(shù)據(jù)來(lái)識(shí)別復(fù)雜疾病中的罕見(jiàn)變異[3-4]。MGAS(multivariate gene-based association test by extended Simes procedure)和metaCCA(summary statistics-based multivariate meta-analysis of genome-wide association studies using canonical correlation analysis)是目前已知分析多元表型和基因型相關(guān)關(guān)系的有效方法，有不涉及候選基因等優(yōu)點(diǎn)，但同時(shí)在確定樣本量和識(shí)別罕見(jiàn)變異等方面也存在局限性[5-7]。對(duì)此，Wei Pan等[8]提出了動(dòng)力分?jǐn)?shù)測(cè)試(sum of powered score tests，SPU)，即利用不同的參數(shù)值來(lái)構(gòu)造SNP數(shù)據(jù)驅(qū)動(dòng)和變化的權(quán)重，從而適應(yīng)多個(gè)SNP之間未知的關(guān)聯(lián)強(qiáng)弱和關(guān)聯(lián)方向，該方法可以在基因和通路兩水平上進(jìn)行表型與基因型相關(guān)分析。本文將重點(diǎn)介紹不同關(guān)聯(lián)分析下SPU測(cè)試的方法、原理和實(shí)現(xiàn)，尤其是在基因和通路水平的多元表型與基因型相關(guān)分析，并探討其應(yīng)用前景。

不同關(guān)聯(lián)性分析下的SPU測(cè)試

1.多SNP-單性狀關(guān)聯(lián)性分析

GWAS是逐個(gè)檢測(cè)每個(gè)SNP，然后進(jìn)行多次檢測(cè)調(diào)整，選出符合要求的SNP。然而，由于罕見(jiàn)變異內(nèi)包含弱相關(guān)信號(hào)以及極小的等位基因頻率(minor allele frequency，MAF)導(dǎo)致GWAS效能較低[8]。Wei Pan等[9]提出了一類適用于匯總統(tǒng)計(jì)數(shù)據(jù)自適應(yīng)權(quán)值及其相應(yīng)的加權(quán)檢驗(yàn)法(adaptive sum of powered score tests，aSPUs)。該測(cè)試考慮了SNP之間的關(guān)聯(lián)性強(qiáng)弱及方向，適用于多SNP單性狀關(guān)聯(lián)性研究，其對(duì)多個(gè)GWAS的大規(guī)模meta分析以及單個(gè)GWAS或?qū)蝹€(gè)SNP的GWAS匯總統(tǒng)計(jì)具有實(shí)際可用性?；诼窂降腉WAS方法是利用基因功能方面的先驗(yàn)生物學(xué)知識(shí)來(lái)促進(jìn)對(duì)GWAS數(shù)據(jù)集的分析，在此基礎(chǔ)上aSPUs測(cè)試也擴(kuò)展到通路分析(aSPUspath)[10]。

2.單SNP-多性狀關(guān)聯(lián)性分析

aSPU測(cè)試的主要思想是在廣義估計(jì)模型(generalized estimation equation，GEE)的框架下構(gòu)建不同權(quán)重的測(cè)試，從這類加權(quán)測(cè)試中選擇最強(qiáng)大的一種，使其能夠在多種情況下保持高效能。aSPU測(cè)試旨在分析在個(gè)體水平或匯總統(tǒng)計(jì)下的單性狀多SNP相關(guān)性。因此，Wei Pan等[11]又將aSPU擴(kuò)展到適應(yīng)于匯總統(tǒng)計(jì)數(shù)據(jù)的多性狀-單SNP關(guān)聯(lián)性研究(the SPU and aSPU tests for multiple traits-single SNP association with GWAS summary statistics，MTaSPUs)。

3.多SNP-多性狀關(guān)聯(lián)性分析

基因與多種性狀的遺傳關(guān)聯(lián)研究已經(jīng)變得越來(lái)越重要，不僅因?yàn)樗袧摿μ岣呓y(tǒng)計(jì)效能，也因?yàn)槠淇紤]了復(fù)雜疾病不同表型之間的相關(guān)性。2016年II-Youp Kwak等[12]提出了基于基因的自適應(yīng)的測(cè)試(MTaSPUsSet)，用GWAS匯總統(tǒng)計(jì)對(duì)多個(gè)性狀進(jìn)行關(guān)聯(lián)分析并將其擴(kuò)展到通路水平(MTaSPUsSetPath)。與傳統(tǒng)多基因關(guān)聯(lián)分析相比，該方法在SNP和性狀層面上都是自適應(yīng)的，并考慮了SNP之間和不同性狀之間可能存在的各種關(guān)聯(lián)模式，使該測(cè)試在大部分情況下保持高功率。該方法也適用于從單個(gè)GWAS或多個(gè)GWAS的meta分析中獲得的Z統(tǒng)計(jì)量或P值的匯總數(shù)據(jù)[12-13]。

基因和通路水平SPU模型的建立

1.aSPU模型的建立

(1)構(gòu)建GEE模型：

假設(shè)對(duì)于每個(gè)目標(biāo)有i=1，…，n，有k個(gè)性狀Yi=(yi1，yi2，…，yik)′，xi=0，1或2為相關(guān)SNP的基因型得分，zi=(zi1，zi2，…，ziq)是協(xié)變量q的行向量。SNP和協(xié)變量通過(guò)邊際廣義線性模型(GLM)建模：

g(μi)=ηi=Ziφ+Xiβ=Hiθ

通過(guò)求解GEE，得到了兩個(gè)相容的漸近正態(tài)估計(jì)：

(2)建立零假設(shè)H0：β=(β1，…，βk)′=0；H1：β≠0；

為了構(gòu)建具有協(xié)變量Zi的基于分?jǐn)?shù)的測(cè)試，在假設(shè)特征具有獨(dú)立工作相關(guān)結(jié)構(gòu)的零假設(shè)下的得分向量為：

(3)構(gòu)造SPU測(cè)試：

由于SPU測(cè)試的冪函數(shù)曲線很難描述，所以我們用SPU測(cè)試的P值來(lái)估計(jì)它的冪函數(shù)。自適應(yīng)地選擇一個(gè)SPU測(cè)試，形成aSPU測(cè)試：

2.多SNP-單性狀關(guān)聯(lián)性分析的模型建立

(1)基因水平的關(guān)聯(lián)性分析——aSPUs的模型建立

首先假設(shè)個(gè)體的基因型和表型數(shù)據(jù)是可用的，采用如下廣義線性模型：

U的協(xié)方差矩陣可以估計(jì)為：

在H0成立的條件下，其均值為：

通過(guò)模擬SPU和aSPU測(cè)試的個(gè)體基因型和表型數(shù)據(jù)，僅利用匯總數(shù)據(jù)Z值就可以定義相應(yīng)的測(cè)試：

(2)通路水平的關(guān)聯(lián)性分析——aSPUsPath的模型建立

用Z分?jǐn)?shù)代替P值構(gòu)建模型的想法也可以擴(kuò)展到適應(yīng)性通路測(cè)試，定義基于基因和通路水平的SPU檢驗(yàn)為：

PathSPU(γ，γG：S)=∑g∈SSPU(γ：g)Γg

其中兩個(gè)整數(shù)γ>0和γG>0分別用于SNP和基因水平上的自適應(yīng)加權(quán)。例如，當(dāng)僅有較少的基因(或SNP)與性狀相關(guān)時(shí)，要想效率更高則需要更大的γG(或γ)。但由于(γ，γG)的最佳值是未知的，為了自適應(yīng)的選擇(γ，γG)，則提出：

3.多SNP-多性狀關(guān)聯(lián)性分析的模型建立

(1)基因水平的關(guān)聯(lián)性分析——MTaSPUsSet的模型建立

假設(shè)有d個(gè)SNP(例如在基因檢測(cè)中)其加性基因型分?jǐn)?shù)g=(g1，…，gd)′通過(guò)應(yīng)用廣義線性模型我們首先考慮表型Yh：

對(duì)于給定的數(shù)據(jù)集{(Yih，gi，ci)：i=1，…，n}有n個(gè)對(duì)象，βh的得分向量Uh=(Uh1，…，Uhd)′如下：

由于(γ1，γ2)的最優(yōu)值是未知的，從而提出了一個(gè)自適應(yīng)選擇(γ1，γ2)的方法：

(2)通路水平的關(guān)聯(lián)性分析——MTaSPUsSetPath的模型建立

對(duì)僅有GWAS匯總統(tǒng)計(jì)的案例進(jìn)行基于路徑的多性狀關(guān)聯(lián)測(cè)試，給出一個(gè)含有|S|基因的途徑S，在SNP水平上，對(duì)于基因g的第dg個(gè)SNP其Z值為：

Z(ig)=(Z(ig)1，Z(ig)2，…，Z(ig)dg)

將基于基因和路徑的測(cè)試定義為一個(gè)性狀，則多個(gè)性狀為：

SPUsPath(γ1，γ2；Z(i)，S)=(∑g∈SSPU(γ1；Z(ig)γ2)/|S|)1/γ2

MTSPUsSetPath(γ1，γ2，γ3；Z，S)=

其中γ1≥1，γ2≥1，γ3≥1分別對(duì)SNP、基因和性狀進(jìn)行加權(quán)以自適應(yīng)的選擇(γ1，γ2，γ3)，從而提出：

SPU模型的軟件實(shí)現(xiàn)

通過(guò)R軟件的aSPU軟件包實(shí)現(xiàn)(https：//cran.r-project.org/web/packages/aSPU/)，本研究中以MTaSPUsSet和MTaSPUsSetPath為例，介紹模型的實(shí)現(xiàn)。

1.MTaSPUsSet的模型實(shí)現(xiàn)

將原始GWAS匯總統(tǒng)計(jì)結(jié)果整理后得到包含在LCORL基因上的單個(gè)SNP對(duì)應(yīng)不同性狀的P值(Ps)或Z值(Zs)的數(shù)據(jù)集。MTaSPUsSet的軟件實(shí)現(xiàn)過(guò)程如下：

(1)利用Plink計(jì)算SNP之間的相關(guān)性(corSNP)，以歐洲后裔人群為例，輸人文件為SNP_id，代碼為：

plink2—file hapmap3—extract SNP_id—keepCEU_ hapmap—r2 inter-chr with-freqs—ld-window-r20-make-bed-out uppro

輸出結(jié)果為uppro，即為corSNP。

(2)下載aSPU安裝包，利用estcov函數(shù)計(jì)算表型之間的相關(guān)性(corPhe)，代碼為：

corPhe=estcov(Ps，Ps=True)

(3)使用以下命令執(zhí)行MTaSPUsSet：

library(aSPU)

(outFP<-MTaSPUsSet(PsF，corSNP=corSNPF，corPhe=corPheF，pow=c(1，2，4，8)，pow2=c(1，2，4，8)，n.perm=100，Ps=TRUE))

即可得出女性在不同SNP權(quán)重和不同性狀權(quán)重下用MTaSPUsSet計(jì)算的LCORL基因與身高、體重、BMI、腰圍、臀圍和腰臀圍這六個(gè)性狀相關(guān)性的P值，結(jié)果如表1。MTaSPUsSet=0.821782，即在不同權(quán)重取值下均適用的該基因-多性狀關(guān)聯(lián)測(cè)試的P值為0.821782。

表1 MTaSPUsSet的模型實(shí)現(xiàn)結(jié)果

(4)在LCORL基因上繪制SNP圖譜

plotG(someGs$LCORL[[1]]，main=“LCORL(P-values)”，zlim=c(0，18))

結(jié)果如圖1所示，呈現(xiàn)了基因LCORL上的SNP分別與身高、體重、BMI、腰圍、臀圍和腰臀圍相關(guān)關(guān)系的P值，與身高相關(guān)的SNP更多。

圖1 基因LCORL相關(guān)單核苷酸與六種性狀的對(duì)數(shù)轉(zhuǎn)換P值

2.MTaSPUsSetPath的模型實(shí)現(xiàn)

MTaSPUsPath的軟件實(shí)現(xiàn)過(guò)程如下：

(1)生成待測(cè)SNP的相關(guān)矩陣(corSNP)和待測(cè)性狀的相關(guān)矩陣(corPhe)，代碼為：

corPhe=estcov(Zs，Zs=True)

plink2-file hapmap3-extract SNP_ id-keep CEU_ hapmap-r2 inter-chr with-freqs-ld-window-r20-make-bed-out corPhe

(2)生成一個(gè)SNP信息矩陣和基因信息矩陣，代碼為：

snp.info<-as.data.frame(snp.info)

gene.info<-as.data.frame(gene.info)

(3)下載aSPU安裝包，執(zhí)行MTaSPUsSetPath命令，代碼為：

library(aSPU)

out<-MTaSPUsSetPath(Zs，corPhe=corPhe，corSNP=corSNP，n.perm=100，snp.info=snp.info，gene.info=gene.info)

out

即可得出不同SNP權(quán)重、不同基因權(quán)重和不同性狀權(quán)重下該通路與多個(gè)性狀間關(guān)聯(lián)測(cè)試的P值，結(jié)果中MTaSPUsSetPath即為該測(cè)試在不同權(quán)重取值下均適用的P值。

SPU測(cè)試的應(yīng)用

II-Youp Kwak等[12]將MTaSPUsSet測(cè)試應(yīng)用于“人體特征基因調(diào)查聯(lián)盟”(genetic investigation of ANthropometric traits，GIANT)的匯總統(tǒng)計(jì)數(shù)據(jù)，分別對(duì)男性和女性的身高、體重、BMI、腰圍、臀圍和腰臀圍這六個(gè)人體測(cè)量學(xué)特征進(jìn)行了基于基因水平的關(guān)聯(lián)測(cè)試。通過(guò)MTaSPUsSet測(cè)試，共有2722976個(gè)SNPs被定位到17562個(gè)基因(每個(gè)基因加上2kb上游和2kb下游區(qū)域)，共鑒定出137個(gè)對(duì)男性或女性人體測(cè)量特征具有全基因組意義的基因：男性為81個(gè)，女性為125個(gè)，兩者共有為69個(gè)。而在相同的參考面板下采用MGAS方法，使用“kgg”軟件僅能識(shí)別出19個(gè)顯著基因。同時(shí)使用MTaSPUsSet和MGAS兩種方法，MTaSPUsSet識(shí)別出27個(gè)對(duì)男性顯著的基因和39個(gè)對(duì)女性顯著的基因，而MGAS分別識(shí)別出7個(gè)和14個(gè)基因，結(jié)果顯示基因RPGRIP1L和RPS10-NUDT3等僅能通過(guò)MTaSPUsSet檢測(cè)得到，這表明MTaSPUsSet有更高的效能。另外，由于metaCCA需要所有單核苷酸多態(tài)性-性狀對(duì)的樣本量相同，而一些單核苷酸多態(tài)性的樣本量在性狀間從大約200到大約70000不等，因此metaCCA不適用于不同疾病之間含有重復(fù)研究對(duì)象的數(shù)據(jù)，如GIANT等[12]。目前還有其他進(jìn)行全基因組關(guān)聯(lián)性研究的相關(guān)應(yīng)用的研究[8，14]，結(jié)果表明相比于一些新的適應(yīng)性測(cè)試，如KBAC(kernel-based adaptive clustering test)，PWST(P-value weighted sum test)和aSSU(adaptive sum ofsquared score test)，SPU測(cè)試具有更高的適應(yīng)性和效能并在模擬實(shí)驗(yàn)中效果更好[8]。

SPU測(cè)試的優(yōu)缺點(diǎn)

首先，與傳統(tǒng)的GWAS相比，SPU測(cè)試在任何參考面板下性能都較優(yōu)，其估計(jì)膨脹因子k接近于1[8]。其次，不同情況下可選擇不同的SPU測(cè)試，由于SPU測(cè)試都是基于一般回歸模型的得分向量，當(dāng)同一SNP具有多個(gè)特征且這些特征的影響范圍很小時(shí)，aSPU測(cè)試效能相比其他方法更強(qiáng)[8，15]。最后，MTaSPUsSet與傳統(tǒng)多基因關(guān)聯(lián)分析相比，在SNP和性狀層面上都是自適應(yīng)的，該方法考慮了SNP和性狀之間可能存在的不同關(guān)聯(lián)模式，例如關(guān)聯(lián)強(qiáng)度和方向，從而在多數(shù)情況下保持高效能。另外，MTaSPUsSet可應(yīng)用于混合類型的性狀，也適用于從單個(gè)GWAS或多個(gè)GWAS的meta分析中獲得的Z統(tǒng)計(jì)量或P值的匯總數(shù)據(jù)。同時(shí)，仿真和實(shí)際數(shù)據(jù)的數(shù)值研究表明，此法具有良好的應(yīng)用前景[12]。此外，除本文介紹的幾種典型SPU測(cè)試外還有許多適用于不同情況的自適應(yīng)關(guān)聯(lián)測(cè)試[11，16]。

然而，SPU測(cè)試是基于匯總數(shù)據(jù)統(tǒng)計(jì)分析得到，而沒(méi)有進(jìn)行生物驗(yàn)證，因此我們無(wú)法推測(cè)任何確定位點(diǎn)的因果影響[17]。可在此基礎(chǔ)上進(jìn)行轉(zhuǎn)錄組廣泛關(guān)聯(lián)研究，以推斷基因表達(dá)狀態(tài)并進(jìn)一步探索候選基因與結(jié)果之間的因果關(guān)系[18-19]。

小 結(jié)

SPU測(cè)試是一種新的基于基因和路徑的自適應(yīng)關(guān)聯(lián)測(cè)試，該測(cè)試可使用GWAS匯總統(tǒng)計(jì)數(shù)據(jù)，其I類錯(cuò)誤率得到了很好控制，克服了目前已知方法不能應(yīng)用于大規(guī)模數(shù)據(jù)以及弱相關(guān)或多重相關(guān)時(shí)不敏感等缺點(diǎn)[20-21]。該方法在基因組學(xué)、蛋白質(zhì)組學(xué)等方面有較好的應(yīng)用前景，為人類了解復(fù)雜性疾病的發(fā)病機(jī)制提供更多的線索，但其理論和方法仍需在應(yīng)用中進(jìn)一步完善。