GWAS：拒绝假阳性之case和control数量比例严重失衡的解决方案（SAIGE模型的应用）

Posted 2021-11-25 chenwenyan

一、为什么要校正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号：

--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，随意选取位点的关联分析结果

第三个文件：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号：

sampleIDindosage.txt文件为样本ID名：

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

 

2）输出文件的结果解读：

生成example.plainDosage.SAIGE.txt文件，其内容如下：

其中，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数量比例严重失衡的工作就完成了。导师再也不用担心你给出一堆假阳性结果了。