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

  1. GWAS中的赢家诅咒 Winner’s curse
  2. 赢家诅咒的校正 WC correction
  3. winnerscurse R包
  4. 参考

GWAS中的赢家诅咒 Winner’s curse

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

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

赢家诅咒的校正 WC correction

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

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

\beta_{Observed}的一个例子
image

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

上面的式子等价于

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

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

在通过阈值筛选的情况下,\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

1.3 1级脊回归 Level 1 ridge regression

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

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

其中

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

该式有一个封闭的解

同第一次脊回归一样,超系数w从一组个数为Q的收缩系数中选择。

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

1.4 K-fold cross-validation

由于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)。解决方案之一便是,在逻辑斯蒂回归中加入一个惩罚项,

此图像的alt属性为空;文件名为image-87.png

但这个方法计算量巨大,所以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