GWAS方法 – SAIGE 校正样本对照失衡的广义线性混合模型

关键词:saige,GLMM,case-control imbalance,binary phenotype, SPA

本文内容:

SAIGE开发背景
SAIGE原理简介
	第一步
	第二步
	注意事项
SAIGE使用方法
	安装
	第一步
	第二步
	结果解读
参考

一句话解释:SAIGE是一种基于广义线性混合模型的,针对二分类表型(binary phenotype)能够调整样本对照不平衡(case-control imbalance)与隐性关联(cryptic relatedness)的GWAS检验方法。

① SAIGE开发背景:

在大规模生物银行,对上千种表型进行的GWAS中,大多数二分类变量的病例都远小于对照,线性混合模型,或逻辑斯蒂混合模型对于此种病例对照严重失衡(1:10,甚至到1:100)的表型,都不能很好地控制1类错误,容易造成假阳性出现。

SAIGE便是针对此问题开发的方法,目前UKBB等GWAS概括性数据所使用的方法均为SAIGE。SAIGE的特点是利用 SPA(saddlepoint approximation , 鞍点近似法) 在病例对照比例严重失衡时,仍然获得准确的p值。

② SAIGE原理简介:

SAIGE全称是 Scalable and Accurate Implementation of GEneralized mixed model 广义混合模型的可扩展且准确的实现。

SAIGE 与 bolt-lmm 或 regenie 类似,都需要两个步骤。saige基于广义线性混合模型GLMM模型:

https://gwaslab.files.wordpress.com/2021/04/image-26.png?w=301

其中

为给定协变量与基因型时,第i个样本为病例的概率,bi为随机效应,服从N(0,τψ)的分布。其中ψ是N X N的遗传关联矩阵(GRM),τ是加性遗传方差。Xi为协变量,Gi表示等位基因的个数,取值为0,1,2,α为 (1+p)x 1 的表示固定效应的系数向量,β为遗传效应的系数。

具体流程如下:

第一步(构建null模型),估计方差与其他参数:

利用 array genotype ,估计 null 逻辑斯蒂混合模型(GLMM)。

https://gwaslab.files.wordpress.com/2021/04/image-27.png?w=236

SAIGE通过PQL方法 和 AI-REML 算法简化计算提高效率。

迭代估计模型中参数:

每一次迭代中,让

为y的估计的均值,

则权重矩阵为

为N X N的如下工作向量的方差矩阵

目前诸如GMMAT等方法在这一步时需要计算该方差矩阵的逆,

当N较大时,该步骤计算量巨大,为了简化计算,SAIGE采用了PCG算法(PCG是一种找到线性系统解的数值方法,BOLT_LMM也用到此算法)

在如下零假设(null hypothesis)的情况下:

分数检验统计量由下式得到

其方差Var(T)为

其中

因为我们需要对每一个SNP都计算P,计算量巨大,

SAIGE假设真随机效应b已经给定,并计算两种情况下T的方差的比值,来进行简化计算

假设b给定时:

则方差比值r

saige中会取随机取30个snp来估计r。

第二步(检验):分数检验,并通过SPA调整。

采用上述比值r对分数检验的T检验量进行调整,即可得到SAIGE的T检验量。

一般情况下,传统方法假设T渐进服从一个高斯分布,但是当病例对照(case-control)的比例不平衡,以及被检验变异有较低的MAC(Minor allele count)时,T统计量就会与高斯分布有显著偏离。saige方法之所以能调节病例对照不平衡 ,关键点就在于利用SPA对检验的调整,以获得一个较为准确的P值。(saige采用了fastSPA)

SPA:saddlepoint approximation (SPA) 鞍点近似法

SAIGE整体架构如下所示:

https://github.com/weizhouUMICH/img/raw/master/SAIGE2step.png

③ SAIGE的注意事项:

注意:该比值只有在MAC大于30时,才保持恒定,当MAC小于30时,要谨慎使用此方法

④ SAIGE使用方法:

saige的github页面:https://github.com/weizhouUMICH/SAIGE

安装:

saige是基于R的软件,有多种安装方法:

可以通过R安装,也可以通过CONDA安装,但安装过程复杂,个人推荐使用docker进行安装,方便快捷。

docker pull wzhou88/saige:0.44.2

第一步,估计空模型

输入文件:

1 基因型文件PLINK格式

2 表性文件

3 协变量文件

#For Binary traits:
Rscript step1_fitNULLGLMM.R     \
        --plinkFile=./input/nfam_100_nindep_0_step1_includeMoreRareVariants_poly \
        --phenoFile=./input/pheno_1000samples.txt_withdosages_withBothTraitTypes.txt \
        --phenoCol=y_binary \
        --covarColList=x1,x2 \
        --sampleIDColinphenoFile=IID \
        --traitType=binary        \
        --outputPrefix=./output/example_binary \
        --nThreads=4 \
        --LOCO=FALSE \
        --IsOverwriteVarianceRatioFile ## v0.38. Whether to overwrite the variance ratio file if the file already exists

#For Quantitative traits, if not normally distributed, inverse normalization needs to be specified to be TRUE --invNormalize=TRUE
Rscript step1_fitNULLGLMM.R     \
        --plinkFile=./input/nfam_100_nindep_0_step1_includeMoreRareVariants_poly \
        --phenoFile=./input/pheno_1000samples.txt_withdosages_withBothTraitTypes.txt \
        --phenoCol=y_quantitative \
        --covarColList=x1,x2 \
        --sampleIDColinphenoFile=IID \
        --traitType=quantitative       \
		--invNormalize=TRUE	\
        --outputPrefix=./output/example_quantitative \
        --nThreads=4 \\
        --LOCO=FALSE	\\
	--tauInit=1,0

共有三个输出文件:

  • .rda文件,包含空模型
  • 30个随机选择的SNP的关联结果
  • 包含估计的方差比值的文本文件 (回忆上述SAIGE算法中的r

.rda文件,包含空模型

./output/example_quantitative.rda

#open R
R
#load the model file in R
load("./output/example_quantitative.rda")
names(modglmm)
modglmm$theta

#theta: a vector of length 2. The first element is the dispersion parameter estimate and the second one is the variance component parameter estimate
#coefficients: fixed effect parameter estimates
#linear.predictors: a vector of length N (N is the sample size) containing linear predictors
#fitted.values: a vector of length N (N is the sample size) containing fitted mean values on the original scale
#Y: a vector of length N (N is the sample size) containing final working vector
#residuals: a vector of length N (N is the sample size) containing residuals on the original scale
#sampleID: a vector of length N (N is the sample size) containing sample IDs used to fit the null model

30个随机选择的SNP的关联结果

less -S ./output/example_quantitative_30markers.SAIGE.results.txt

方差比值文件 variance ratio file

less -S ./output/example_quantitative.varianceRatio.txt

第二步

输入文件为:

1 插补后的剂量文件,如VCF,BGEN,BCF或SAV,以及

2 对应的索引文件

3 样本文件,无header,包含一列对应剂量文件的样本ID。

4 第一步所得的rda模型文件

5 第一步所得的方差比值文件

这里以常用的bgen与VCF文件为例

#bgen for dosages
Rscript step2_SPAtests.R \
        --bgenFile=./input/genotype_100markers.bgen \
        --bgenFileIndex=./input/genotype_100markers.bgen.bgi \
        --minMAF=0.0001 \ #最小MAF不能低于0.0001
        --minMAC=1 \ #最小MAC不能低于1
        --sampleFile=./input/samplefileforbgen_10000samples.txt \
        --GMMATmodelFile=./output/example.rda \
        --varianceRatioFile=./output/example.varianceRatio.txt \
        --SAIGEOutputFile=./output/example.SAIGE.bgen.txt \
        --numLinesOutput=2 \
        --IsOutputAFinCaseCtrl=TRUE

#VCF containing dosages (--vcfField=DS)
Rscript step2_SPAtests.R \
        --vcfFile=./input/dosage_10markers.vcf.gz \
        --vcfFileIndex=./input/dosage_10markers.vcf.gz.tbi \
        --vcfField=DS \
        --chrom=1 \
        --minMAF=0.0001 \   #最小MAF不能低于0.0001
        --minMAC=1 \   #最小MAC不能低于1
        --sampleFile=./input/sampleIDindosage.txt \
        --GMMATmodelFile=./output/example.rda \
        --varianceRatioFile=./output/example.varianceRatio.txt \
        --SAIGEOutputFile=./output/example.SAIGE.vcf.dosage.txt \
        --numLinesOutput=2 \
        --IsOutputAFinCaseCtrl=TRUE

输出结果各列的解读如下:

CHR: chromosome
POS: genome position 
SNPID: variant ID
Allele1: allele 1
Allele2: allele 2
AC_Allele2: allele count of allele 2
AF_Allele2: allele frequency of allele 2
imputationInfo: imputation info. If not in dosage/genotype input file, will output 1
N: sample size
BETA: effect size of allele 2
SE: standard error of BETA
Tstat: score statistic of allele 2

#SPA后的p值
p.value: p value (with SPA applied for binary traits) 

p.value.NA: p value when SPA is not applied (only for binary traits)
Is.SPA.converge: whether SPA is converged or not (only for binary traits)
varT: estimated variance of score statistic with sample relatedness incorporated
varTstar: variance of score statistic without sample relatedness 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)
Tstat_cond, p.value_cond, varT_cond, BETA_cond, SE_cond: summary stats for conditional analysis

除此之外SAIGE还支持条件分析,男女分层分析,X染色体分析等多种功能:

详细文档见:https://github.com/weizhouUMICH/SAIGE/wiki/Genetic-association-tests-using-SAIGE#installing-saige

参考:

Zhou, W. et al. (SAIGE)Efficiently controlling for case-control imbalance and sample relatedness in large-scale genetic association studies. Nat. Genet. 50, 1335–1341 (2018).

发表评论

Fill in your details below or click an icon to log in:

WordPress.com 徽标

您正在使用您的 WordPress.com 账号评论。 注销 /  更改 )

Twitter picture

您正在使用您的 Twitter 账号评论。 注销 /  更改 )

Facebook photo

您正在使用您的 Facebook 账号评论。 注销 /  更改 )

Connecting to %s