关键词：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模型：
其中
为给定协变量与基因型时，第i个样本为病例的概率，bi为随机效应，服从N(0，τψ)的分布。其中ψ是N X N的遗传关联矩阵（GRM），τ是加性遗传方差。Xi为协变量，Gi表示等位基因的个数，取值为0,1，2，α为 （1+p）x 1 的表示固定效应的系数向量，β为遗传效应的系数。
具体流程如下：
第一步（构建null模型），估计方差与其他参数：
利用 array genotype ，估计 null 逻辑斯蒂混合模型（GLMM）。
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整体架构如下所示：
③ 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).
GWASLab - GWAS实验室 