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

Posted chenwenyan

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了GWAS:拒绝假阳性之case和control数量比例严重失衡的解决方案(SAIGE模型的应用)相关的知识,希望对你有一定的参考价值。

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

 

 

以上是关于GWAS:拒绝假阳性之case和control数量比例严重失衡的解决方案(SAIGE模型的应用)的主要内容,如果未能解决你的问题,请参考以下文章

一行命令学会全基因组关联分析(GWAS)的meta分析

人脸识别中的假阳性

Python深度学习之路-3.1性能评价指标

Python深度学习之路-3.1性能评价指标

ROC评分中概念之阳性预测值/阴性预测值计算方法

Scikit-learn:如何获得真阳性、真阴性、假阳性和假阴性