Association test – GWASLab

GWAS中的赢家诅咒与其校正 Winner’s curse correction

GWAS中的赢家诅咒 Winner’s curse

GWAS中的赢家诅咒是指遗传效应的大小由于GWAS中的筛选过程（通过全基因组显著阈值筛选lead SNP）而被系统性地过高估计。

赢家诅咒本用来指代在拍卖中类似的现象。即使一件拍卖品对所有买家来说都有相同的价值（出价是无偏的），最后拍得物品的赢家很可能过高估计了拍卖偏的内在价值。类比于GWAS，lead SNP即为赢家，而它的效应量可能过高估计了真实的遗传效应。

赢家诅咒的校正 WC correction

假设观察到的 $\beta_{Observed}$ 的近似分布为:

$\beta_{Observed} \sim N(\beta_{True},\sigma^2)$

$\beta_{Observed}$ 的一个例子

$c$ : 显著性阈值对应的Z分数

上面的式子等价于

${{\beta_{Observed} - \beta_{True}}\over{\sigma}} \sim N(0,1)$

${{\beta_{Observed} - \beta_{True}}\over{\sigma}}$ 的一个例子

在通过阈值筛选的情况下， $\beta_{Observed}$ 的近似抽样分布（实际上为一个截断正态分布 truncated normal distribution）为：

$f(x,\beta_{True}) ={{1}\over{\sigma}} {{\phi({{{x - \beta_{True}}\over{\sigma}}})} \over {\Phi({{{\beta_{True}}\over{\sigma}}-c}) + \Phi({{{-\beta_{True}}\over{\sigma}}-c})}}$

其中

$|{{x}\over{\sigma}}|\geq c$

$\phi(x)$ : 标准正态分布的概率密度函数
$\Phi(x)$ : 标准正态分布的累积分布函数

从以上的近似抽样分布可以得到，筛选出来的SNP的效应量的期望分布为：

$E(\beta_{Observed}; \beta_{True}) = \beta_{True} + \sigma {{\phi({{{\beta_{True}}\over{\sigma}}-c}) - \phi({{{-\beta_{True}}\over{\sigma}}-c})} \over {\Phi({{{\beta_{True}}\over{\sigma}}-c}) + \Phi({{{-\beta_{True}}\over{\sigma}}-c})}}$

$\beta_{Observed}$ is biased.
偏差的大小由 $\beta_{True}$ , SE $\sigma$ , 以及用于筛选的显著性阈值决定.

公式推导可以参考 Ghosh, A., Zou, F., & Wright, F. A. (2008). Estimating odds ratios in genome scans: an approximate conditional likelihood approach. The American Journal of Human Genetics, 82(5), 1064-1074. 中的Appendix A

用这个式子便可以对效应量进行赢家诅咒的校正。

winnerscurse R包

可以使用这个R包进行赢家诅咒的校正。

https://amandaforde.github.io/winnerscurse/articles/winners_curse_methods.html

参考

Bazerman, M. H., & Samuelson, W. F. (1983). I won the auction but don’t want the prize. Journal of conflict resolution, 27(4), 618-634.
Göring, H. H., Terwilliger, J. D., & Blangero, J. (2001). Large upward bias in estimation of locus-specific effects from genomewide scans. The American Journal of Human Genetics, 69(6), 1357-1369.

Zhong, H., & Prentice, R. L. (2008). Bias-reduced estimators and confidence intervals for odds ratios in genome-wide association studies. Biostatistics, 9(4), 621-634.
Ghosh, A., Zou, F., & Wright, F. A. (2008). Estimating odds ratios in genome scans: an approximate conditional likelihood approach. The American Journal of Human Genetics, 82(5), 1064-1074.

Also see reference: https://amandaforde.github.io/winnerscurse/articles/winners_curse_methods.html

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).

基于LMM的一种快速的关联检验方法 fastGWA

TL；DR

一句话总结：fastGWA是基于LMM的检验方法，该方法应用了一种快速的估算方法（fastGWA-REML，也就是网格搜索），可以避免计算协方差矩阵V的逆，从而大幅提升计算速度。

背景：

目前基于线性混合模型（LMM）的检验方法已经被广泛应用于GWAS研究中，因为LMM可以矫正群体分层以及亲缘关系。基本的原理就是以从所有SNP估计而来的样本结构为条件，检验变异与表型的关联。详见 LMM模型。

然而该模型的运算时间过长，目前的方法时间复杂度约为 O(mn2）到 O(m2n)之间，其中m是变异的数量，n是样本大小。

基于此背景，fastGWA旨在运用一种及其节省运算资源的算法，来进行基于LMM的GWAS分析，该方法内置在GCTA软件中。

fastGWA的LMM模型：

一般情况下LMM模型需要通过REML来估计参数，fastGWAS采用了以下的算法来简化计算，提升速度。

fastGWA-REML 算法：grid search（加速的重点）：

计算 log-likelihood scores 时，将g的方差可能的取值范围，等间距取100个值，

所以每一步的间隔就是

注意，这里可能的取值上限取了大于表型方差的值，原因是要尽可能的包含各种情况，包括当真实的遗传力较大，且环境因素影响也十分显著时， g的方差估计值可能大于表型方差。

接下来，要细化搜索窗口，在前面100个取值点中，log-likelihood scores 取最大值的点周围，20%的范围内，再细分16步：

也就是说如果第一步

那么就要在下面的范围里，再细分出16个取值点，

每步的大小如下，

接下来不断重复这个步骤，直到两次细分后找的最大值只差小于：

这样我们就可以相对快速而精准的估计出随机效应方差的大小。

有了估计的协方差矩阵，我们就能利用一般化的最小二乘法估计效应量beta，算式如下：

参考：

Jiang L, Zheng Z, Qi T, et al. A resource-efficient tool for mixed model association analysis of large-scale data[J]. Nature genetics, 2019, 51(12): 1749-1755.

https://cnsgenomics.com/software/gcta/#Overview

基于高斯混合模型的关联检验 Bolt-LMM – GWAS方法

Key words:　LMM，高斯混合模型，贝叶斯，无穷小模型，长尾分布

TL；DR

Bolt-LMM对SNP的效应拟合高斯混合模型，该算法使用一种快速的方差近似方法，计算近似的表型残差，并通过回顾性的分数检验统计量检验残差与待检验SNP的关联，这样就构建了表型预测的贝叶斯模型与频率学派关联检验的桥梁。同时，基于LD分数回归，还会对统计量进行调整。

背景：

在bolt-lmm论文发表时已出现的LMM方法有以下的不足：

需要大量的计算资源，时间复杂度高。
现有模型由于对遗传结构非最优化（suboptimal modeling assumptions regarding the genetic architectures）的假设，会导致power降低。当前标准的线性混合模型是基于无穷小模型（infinitesimal model），该模型假设所有的变异都是效应量很小的因果变异，各效应量相互独立，服从高斯分布，但实际情况中，对于复杂表型，因果loci的数量大约为1000个左右。

为了解决以上问题，Bolt-LMM采取了贝叶斯的观点修改了混合模型，新模型中SNP效应量服从非高斯的先验分布，以更好地反映效应量大小不同的loci的遗传效应。

方法详解：

BOLT-LMM算法包含四个步骤，每一步都需要若干次时间复杂度为O(MN)的迭代。 (1a) 估计方差系数； (1b) 计算无穷小混合模型下的关联统计量（Bolt-LMM-inf） (2a) 估计高斯混合模型的系数 (2b) 计算高斯混合模型下的关联统计量（Bolt-LMM）
简要推导：

标准的线性混合模型如下所示：

Y是表型，x是待检验的SNP，g是遗传效应，e是环境因素在无穷小模型下，遗传效应g可以表示为：

其中XGRM是一个NxM的矩阵，每一列都是某个SNP标准化后的基因型，βGRM是长度为M的向量，包含了SNP的随机效应，效应都从相同的正态分布中抽取，并且整体上服从协方差矩阵如下所示的多元正态分布，

这里BOlt-LMM为了避免近端污染（proximal contamination），采用了LOCO方法。

这个矩阵在习惯上称为GRM，或是亲缘关系矩阵K，于是有：

σg2是方差系数。环境效应也被认为是独立同分布，服从多元正态分布，

σe2是方差系数，I是单位矩阵。实际上σg2与σe2是未知的，所以我们要先通过REML来估计。然后计算前瞻性的卡方检验统计量：

其中，

使σg2与σe2为空模型：β=0是的估计值，在LOCO下，检验统计量变为：

BOLT-LMM-inf 无穷小混合模型统计量：

cinf是一个常数的校正因子，由下式估计：

使得，

实际操作中选取30个伪随机的SNP来估计cinf。我们可以将BOLT-LMM-inf统计量视为前瞻性统计量（将表型视为随机）的近似，或是回顾性的统计量（将基因型视为随机，基于SNP构建空模型）
BOLT-LMM 高斯混合模型关联统计量：
我们注意到，

是以下最优无偏估计（BLUP）的表型残差向量的标量倍数，

因此，BOLT-LMM-inf统计量就等价于计算（并调整）待检验的SNP xtest与BLUP残差的相关系数的平方。混合模型的power是基于以下事实，SNPs是基于对“去噪声”后的表型残差进行检验，即被混合模型估计的其他SNP的效应已经被矫正。我们可以一般化这个过程，定义：

其中 yresidual-LOCO表示拟合标准LMM的高斯混合扩展（用于待测SNP不在一条染色体上的SNP）后的一般化的表型残差向量，C表示校正因子，通过LD分数回归估计，以使得BOLT-LMM的卡方统计量回归后的截距匹配BOLT-LMM-inf的截距。在无穷小模型下，yresidual-LOCO与Vy成正比，BOLT-LMM的卡方统计量即为BOLT-LMM-inf的卡方统计量。
为了定义高斯混合模型LMM扩展，我们首先构建贝叶斯框架下的标准LMM模型，BOLT-LMM-inf的空模型是

其中，SNP效应βm（m是指除m号染色体之外染色体上的SNP）互相独立且服从以下的高斯先验分布

环境效应也互相独立，服从以下分布：

这里BLUP估计等同于计算遗传效应XLOCOβLOCO的后验均值

为了一般化这个模型（非无穷小模型），对于SNP效应，我们采用一个更一般化的先验分布，在BOLT-LMM中，使用了两个高斯分布的spike and slab的混合分布作为先验分布，如下所示：

这种混合更灵活的表示了遗传效应更为典型的长尾分布（heavier-tailed distributions）。在这个一般化的模型中，后验均值不再与BLUP相一致，但我们仍可以拟合这个贝叶斯模型以或得残差：

最后将残差带入前面的算式就可以得到BOLT-LMM 高斯混合模型关联统计量。

参考：

Loh, P. R. et al. Efficient Bayesian mixed-model analysis increases association power in large cohorts. Nature Genetics 47, 284–290 (2015).

Regenie – 高效的全基因组回归

本方法由Regeneron Genetics Center开发，regenie基于一种机器学习方法（Stacked block ridge regression）。Regenie相比于之前的基于LMM的方法，如bolt-LMM，fastGWA，以及SAIGE等方法，具有以下主要特点：

1 相比其他方法更快，更省内存。
2 可适用于数量表型，或是二分类表性。
3 可以同时对多个表型进行检验。
4 可以矫正在case-control不平衡时，mac过低造成的SNP效应值的过高估计（SAIGE-SPA存在此问题）。

regenie 方法简介：

regenie与saige或fastGWA一样，同样采用了两个步骤的模式来进行检验。

STEP 1，估计空模型（主要利用 Stacked block ridge regression的方法）

1 层叠区块脊回归 Stacked block ridge regression

1.1 分割区块 Block partitions

将转化后的基因型矩阵G分割为B个连续的区块，每个区块所包含的SNP都来自同一条染色体。定义区块的方法有很多种，例如物理距离，遗传距离等等，实际使用中，一般将区块大小固定。设mi为第i个区块的SNP数量，则有：

1.2 0级脊回归 Level 0 ridge regressions

这一步的目的是减小估计值的数量，使估计值能够捕获（1）区块中SNP之间的LD信息（2）SNP对表性的任何效应。

这一步考虑拟合带有L2惩罚项的线性模型（也就是脊回归）。

其中，

我们可以通过以下式子估计SNP的效应值γ

其中λ为收缩系数， λ 越高会使估计的效应值趋近于0. 效应值得估计也可以被视为贝叶斯框架下的最大后验估计值（MAP），这个框架里对于SNP效应值，我们使用一个正态先验概率：

基于以上假设，我们可以将脊回归的收缩系数视为一个SNP遗传力的方程，即

由于我们并不知道SNP对于表型真实的效应值，所以考虑使用一组系数来反应SNP的遗传力可能的值。

更准确的说，就是选择R个介于0到1之间等间隔的值。然后利用该值来计算λ。

对于每个大小为Mi的区块，我们都可以计算出R个估计值，

第一步的结果就是一个Nx（BR）的预测值矩阵。

1.3 1级脊回归 Level 1 ridge regression

为了进一步整合第一步结果W中的预测值，在第一步脊回归的基础上再次脊回归（称为stacking），目的是通过惩罚项获得一个最优化的区块预测值的线性组合。

通过变换使W有单位方差，然后进行第二次脊回归，有

其中

该式的解可以由以下式子表示，

该式有一个封闭的解

一个最适的w*通过最小化预测值误差来获得，

由于W与表型相关，采用K-fold cross-validation来避免过拟合。

1.5 LOCO

为了避免近端污染（proximal contamination），在 cross-validation 的计算中会将一整条染色体的区块的值，设为0。这样最终，我们会有22个预测值。

1.6 对于二分类表型，

第一次回归与数量表型一样，但在第二次回归时，采用逻辑斯蒂回归而不是脊回归。

STEP 2 关联检验 Association testing

第二步，检验，

2.1 对于数量表型，考虑一下模型，

使用分数检验，统计量如下：

2.2 对于二分类表型，采用以下模型

2.3 Firth logistic regression（重点）

当case-control不平衡时，稀有变异的检验通常会有不稳定的效应估计值与标准差，主要是在逻辑斯蒂回归中，case里不存在minor allele时会出现拟完全分离（quasi-complete separation）。解决方案之一便是，在逻辑斯蒂回归中加入一个惩罚项，

但这个方法计算量巨大，所以regenie使用了近似的 Firth logistic regression ，使用了一个调整后的惩罚项：

这样就大大减少了计算量，最后应用likelihood ratio test (LRT) 来检验。

2.4 regnie在第二步中也支持saige所使用的SPA方法。这里就略过。

最后放上文章中的总结图：

参考：

https://rgcgithub.github.io/regenie/

https://www.biorxiv.org/content/10.1101/2020.06.19.162354v2