一、为什么要校正case和control数量比例不平衡情况

试问作为生信届人员,最怕的是什么,当然是统计结果不靠谱。统计结果不靠谱包括两方面:一个是假阴性,一个是假阳性。假阴性可以理解为白天鹅被误当成丑小鸭了,假阳性可以理解为一大堆青蛙,你不知道哪个才是你的真命天子。假阴性就罢了,最多让你错过发现真理的机会,但万一假阳性呢,你拿着一个看似完美的结果吭哧吭哧做实验验证,一年半载的周期下来,什么结果都验证不出来,岂不是坑了做实验的人。因此,我们就要在源头上,把这个不靠谱的统计结果杜绝出去。

上一篇文章什么!GWAS研究中case和control的比例是有讲究的?就讲到GWAS分析中,如果case和control数量比例失衡的话,会产出非常多的假阳性结果,而用SAIGE模型做GWAS分析可以校正这种数量比例不平衡的情况。下面具体讲讲怎么应用SAIGE模型。

二、怎么校正:SAIGE的下载和安装

1、下载SAIGE

此操作在Linux上进行,系统要求R-3.5.1, gcc >= 5.5.0, cmake 3.8.1

wget https://github.com/weizhouUMICH/SAIGE/blob/master/SAIGE_0.35.8.1_R_x86_64-pc-linux-gnu.tar.gz

  

2、安装SAIGE

R CMD INSTALL SAIGE_XX_R_x86_64-pc-linux-gnu.tar.gz

  

3、安装SAIGE所依赖的其他R包:Rcpp, RcppArmadillo, RcppParallel, data.table, SPAtest, RcppEigen, Matrix, methods, optparse

以下两个方法二选一:

如果是用conda的话,则用以下命令:

conda install -n r-env r-Rcpp #r-env是指conda下的R环境
conda install -n r-env r-RcppArmadillo
conda install -n r-env r-RcppParallel
conda install -n r-env r-SPAtest
conda install -n r-env r-RcppEigen
conda install -n r-env r-optparse

  

也可以进入R,用install.packages( )安装:

install.packages("Rcpp")
install.packages("RcppArmadillo")
install.packages("RcppParallel")
install.packages("SPAtest")
install.packages("RcppEigen")
install.packages("optparse")

  

三、怎么校正:SAIGE的分析、解读

1、第一步,计算NULL GLMM

1)计算NULL GLMM的命令:

Rscript step1_fitNULLGLMM.R \
--plinkFile=./input/plinkforGRM_1000samples_10kMarkers \
--phenoFile=./input/pheno_1000samples.txt \
--phenoCol=y \
--covarColList=x1,x2 \
--sampleIDColinphenoFile=IID \
--traitType=binary \
--outputPrefix=./output/example \
--nThreads=4 \
--LOCO=TRUE

  

plinkFile为plink的输入文件(bed, bim, fam格式)

phenoFile文件格式如下,第一列代表研究的表型,第二列及第N-1列代表协变量,最后一列IID为个体的ID号:

GWAS:拒绝假阳性之case和control数量比例严重失衡的解决方案(SAIGE模型的应用)-LMLPHP

--phenoCol=y # 指定你要研究的表型列名,在本次例子中,指定y的表型分析。

--covarColList=x1,x2 #指定加入的协变量

--sampleIDColinphenoFile=IID #指定样本的ID

--traitType=binary  #指定研究的表型的类型,binary指二分类,即case和control

--outputPrefix #生成文件的输出路径

2)输出文件的结果解读:

这个步骤会生成三个文件,分别为:example.rda、example.varianceRatio.txt、example_30markers.SAIGE.results.txt

第一个文件:example.rda,是一个model file

可以用R打开:

load("./output/example.rda")
names(modglmm)
modglmm$theta

  

第二个文件:example_30markers.SAIGE.results.txt,随意选取位点的关联分析结果

GWAS:拒绝假阳性之case和control数量比例严重失衡的解决方案(SAIGE模型的应用)-LMLPHP

第三个文件:example.varianceRatio.txt

2、计算每个marker的SPA得分

1)计算每个marker的SNP得分命令:

Rscript step2_SPAtests.R \
--dosageFile=./input/dosage_10markers.txt \
--dosageFileNrowSkip=1 \
--dosageFileNcolSkip=6 \
--dosageFilecolnamesSkip=CHR,SNP,CM,POS,EFFECT_ALLELE,ALT_ALLELE \
--minMAF=0.0001 \
--sampleFile=./input/sampleIDindosage.txt \
--GMMATmodelFile=./output/example.rda \
--varianceRatioFile=./output/example.varianceRatio.txt \
--SAIGEOutputFile=./output/example.plainDosage.SAIGE.txt \
--numLinesOutput=2 \
--IsOutputAFinCaseCtrl=TRUE

  

dosage_10markers.txt的文件格式如下,类似于plink的ped格式,前六列分别为:CHR, SNP, CM, POS, COUNTED, ALT, 后面为个体的ID号:

GWAS:拒绝假阳性之case和control数量比例严重失衡的解决方案(SAIGE模型的应用)-LMLPHP

sampleIDindosage.txt文件为样本ID名:

GWAS:拒绝假阳性之case和control数量比例严重失衡的解决方案(SAIGE模型的应用)-LMLPHP

example.rda、example.varianceRatio.txt为第一步生成的两个文件。

2)输出文件的结果解读:

生成example.plainDosage.SAIGE.txt文件,其内容如下:

GWAS:拒绝假阳性之case和control数量比例严重失衡的解决方案(SAIGE模型的应用)-LMLPHP

其中,P值(红框)即为我们校正case和control数量比例不平衡以后得到的GWAS结果,p.value.NA为不校正case和control数量不平衡的结果。

参数说明:

CHR: chromosome

POS: genome position

SNPID: variant ID

Allele1: Ref allele

Allele2: Alt allele

AC_Allele2: allele count of Alt allele

AF_Allele2: allele frequency of Alt allele

N: sample size

BETA: effect size

SE: standard error of BETA

Tstat: score statistic

p.value: p value with SPA applied

p.value.NA: p value when SPA is not applied

Is.SPA.converge: whether SPA is converged or not

varT: estimated variance of score statistic with sample related incorporated

varTstar: variance of score statistic without sample related incorporated

AF.Cases: allele frequency of allele 2 in cases (only for binary traits and if --IsOutputAFinCaseCtrl=TRUE)

AF.Controls: allele frequency of allele 2 in controls (only for binary traits and if --IsOutputAFinCaseCtrl=TRUE)

至此,校正GWAS分析中case和control数量比例严重失衡的工作就完成了。导师再也不用担心你给出一堆假阳性结果了。

05-02 08:04