SPU測試及其在全基因組多元統(tǒng)計分析中的應(yīng)用*

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

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

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

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

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

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

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

aSPU測試的主要思想是在廣義估計模型(generalized estimation equation，GEE)的框架下構(gòu)建不同權(quán)重的測試，從這類加權(quán)測試中選擇最強大的一種，使其能夠在多種情況下保持高效能。aSPU測試旨在分析在個體水平或匯總統(tǒng)計下的單性狀多SNP相關(guān)性。因此，Wei Pan等[11]又將aSPU擴展到適應(yīng)于匯總統(tǒng)計數(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)變得越來越重要，不僅因為它有潛力提高統(tǒng)計效能，也因為其考慮了復(fù)雜疾病不同表型之間的相關(guān)性。2016年II-Youp Kwak等[12]提出了基于基因的自適應(yīng)的測試(MTaSPUsSet)，用GWAS匯總統(tǒng)計對多個性狀進行關(guān)聯(lián)分析并將其擴展到通路水平(MTaSPUsSetPath)。與傳統(tǒng)多基因關(guān)聯(lián)分析相比，該方法在SNP和性狀層面上都是自適應(yīng)的，并考慮了SNP之間和不同性狀之間可能存在的各種關(guān)聯(lián)模式，使該測試在大部分情況下保持高功率。該方法也適用于從單個GWAS或多個GWAS的meta分析中獲得的Z統(tǒng)計量或P值的匯總數(shù)據(jù)[12-13]。

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

1.aSPU模型的建立

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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分別對SNP、基因和性狀進行加權(quán)以自適應(yīng)的選擇(γ1，γ2，γ3)，從而提出：

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

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

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

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

(1)利用Plink計算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ù)計算表型之間的相關(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計算的LCORL基因與身高、體重、BMI、腰圍、臀圍和腰臀圍這六個性狀相關(guān)性的P值，結(jié)果如表1。MTaSPUsSet=0.821782，即在不同權(quán)重取值下均適用的該基因-多性狀關(guān)聯(lián)測試的P值為0.821782。

表1 MTaSPUsSet的模型實現(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)單核苷酸與六種性狀的對數(shù)轉(zhuǎn)換P值

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

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

(1)生成待測SNP的相關(guān)矩陣(corSNP)和待測性狀的相關(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)生成一個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)重下該通路與多個性狀間關(guān)聯(lián)測試的P值，結(jié)果中MTaSPUsSetPath即為該測試在不同權(quán)重取值下均適用的P值。

SPU測試的應(yīng)用

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

SPU測試的優(yōu)缺點

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

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

小 結(jié)

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