作者： Yunye

基于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

解释复杂疾病的四种主流模型 CDCV/RAME/infinitesimal/Broad-sense-heritability

一，常见疾病-常见变异假说 Common Disease Common Variant Hypothesis

该模型主要是指常见病的易感性主要由中等数量的常见变异引起。此假设有两个关键点，一是常见变异的效应相比于稀有变异应该较小，二是常见变异只有较小效应而常见病又有遗传力的话，那么一定有多个常见变异共同影响疾病易感性。早期的GWAS都是基于这一简单地假设。但这一假设有明显的不足，CDCV无法解释消失的遗传力的问题（ missing heritability ），即基于CDCV的GWAS所发现的常见变异只能解释很小一部分推测的遗传力。所以其他模型也开始得到重视以解决这一问题。

二，主要效应的稀有变异模型 The rare alleles of major effect (RAME) model

RAME 模型主要是指常见病病因其实异质性非常高，也就是说，少量MAF<0.01的稀有变异可以促进疾病的发展。主要效应的原位变异可以解释个位数百分比的遗传力。该模型主要聚焦于由于单倍剂量不足或是功能获得性突变而引起的显性效应，这类效应能使得风险上升两倍或者更多。

三，无穷小模型 The infinitesimal model

近年来GWAS研究中无穷小模型逐渐流行。该模型认为复杂疾病的遗传变异是由于大量的，效应很弱（相对风险低于1.2）的变异引起。该模型解释了丢失的遗传力其实大部分是被隐藏了，由于大量对疾病有较弱效应的变异无法在检验中达到预设的显著阈值。目前很多GWAS关联检验方法都基于这一模型。

四，广义遗传力模型：非加性 GxG 与 GxE 相互作用，以及表观遗传效应

Broad sense heritability model: non-additive G×G and G×E interactions and epigenetic effects

广义遗传力模型认为常见变异的效应与稀有变异的效应不足以解释丢失的遗传力。该模型的主要支持依据就是目前在模式生物数量遗传研究中发现的基因型-基因型相互作用（G×G interactions; 也叫上位效应epistasis），以及基因型-环境相互作用（G×E interactions）。除此以外还包括研究越来越多的表观遗传效应，最为显著的就是亲源效应的遗传贡献，与DNA甲基化的继承。

图1，基于4种不同疾病模型的GWAS的特征。

X轴是SNP在染色体上的位置，Y轴是单个SNP所解释的疾病易感性方差的百分数（注意：一般曼哈顿图中y轴是-log10（P））。1，在CDCV模型中，少量的中等效应的loci会产生很强的信号。2，在RAME模型中，少数的稀有变异的因果效应在个别个体中有较强的效应，但不足以解释大部分方差。3，无穷小模型有少数显著的loci，若干在LD区块中的SNP也会显著。4，如果关联只出现在某种环境中（绿色与橙色的信号），那么在一个同时有两种环境混合的群体中，整体效应就会减弱（如箭头所示），只有很少的关联能够被检测，降低了所能解释的方差。

参考：

https://www.nature.com/articles/nrg3118

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5635617/

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2914559/

基于高斯混合模型的关联检验 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).

易感性尺度遗传力与观测尺度遗传力 Liability scale heritability & observed scale heritability

什么是易感性阈值模型？

参考：易感性-阈值模型 Genetic liability, Threshold model

什么是易感性尺度遗传力，以及观测尺度遗传力？

遗传力通常可以基于不能直接观测的易感性尺度（hl2），或是基于可观测的二分型表型（ho2）来计算。

通常遗传学家更为关注基于易感性尺度的遗传力（Liability scale heritability），主要原因是如果使用观测尺度遗传力（observed scale heritability），容易产生较大的偏差，原因有以下几点：

1 对于数量性状（ quantitative traits）来说，测量的尺度与遗传力的表现尺度是一致的（都是连续的），对于二分类性状来说（binary traits）, 病例或对照的状态 (case-control status) 是基于0 -1 尺度的，但遗传力在易感性尺度上才是容易解释的。

2 病例的确定。 在病例对照研究中，通常病例所占的比例都要（远）高于人群中的流行率。但病例的确定通常会对遗传力的估计造成偏差。

如下图所示，A为数量变量的易感性分布，B二分类的形状在人群中（即随机）的易感性分布 ，C为病例对照试验中的 易感性分布 ，刻意地过多纳入病例，会引起遗传力估计上的偏差。

3 SNP的质量控制（QC）。相比于数量性状的研究，病例对照实验的QC更为重要。对于数量性状来说，试验或测量上的误差一般不会与表型值相关，但 病例对照实验 中 病例组与对照组一般都是独立采集的，所以试验误差容易造成病例与其他病例更为相似，对照也与其它对照更为相似。这样人工造成的误差会影响基于基因组相似性计算的遗传力。

如何转换？

详细推导过程可以参考下面这篇文章：

Lee SH, Wray NR, Goddard ME, Visscher PM. Estimating missing heritability for disease from genome-wide association studies. Am J Hum Genet. 2011;88(3):294-305. doi:10.1016/j.ajhg.2011.02.002

Lee SH 等人的文章介绍了转换的方法，当case和control并不是从人群中随机抽取的时候，转换公式如下：

hl2为 易感性尺度的遗传力 ，ho2为 观测尺度遗传力

K为人群中病例的比例（流行率）

P为抽取的样本中病例的比例

z则为正态分布的密度函数的在所取阈值t处的值

使用R进行转换

#K = pop prevalence 人群中流行率
#P = proportion of cases in study 样本中病例比例
#h2 = Heritability estimate (on observed scale) 观测尺度遗传力
#T = liability threshold 易感性阈值
#zv = 正态分布的密度函数的在所取阈值t处的值

K=0.0659
P=0.0659
h2=0.0365
zv <- dnorm(qnorm(K))

h2_liab <- h2 * K^2 * ( 1 - K)^2 / P / (1-P) / zv^2
h2_liab

参考文献：

https://www.pnas.org/content/111/49/E5272

Estimating Missing Heritability for Disease from Genome-wide Association Studies

易感性-阈值模型 Genetic liability, Threshold model

易感性-阈值模型是遗传流行病学中重要的理论模型之一。

Liability 易感性

易感性是 遗传因素 与 环境因素 对某多因子疾病的效应 的总称。因为易感性是一个隐变量（latent variable），实际上很难直接测量某个特定个体的易感性，但我们可以通过某群体中发病个体数量来估计该群体对于疾病的易感性。

Threshold Model 阈值模型

易感性-阈值模型（liability -threshold model）是我们用来分析非孟德尔遗传的（non-Mendelian）分类（categorical）表型（例如二分类表型 binary traits）的基本模型之一，如下图所示，通常我们认为易感性服从一个N(0,1)的正态分布。当某个个体的集聚的易感性超过了阈值T时就会发病，阴影区域的面积表述了在这个群体中该疾病的流行率（prevalence）。

对于某个疾病可以有多个阈值，分别对应疾病的不同严重程度。

参考：

https://onlinelibrary.wiley.com/doi/full/10.1002/0470011815.b2a05036

使用GCTA （GREML）来估计SNP-遗传力 SNP Heritability

GWAS研究中发现的显著SNP只能解释人类群体中复杂性状很小一部分的遗传变异。那么剩余的遗传力在哪里？ 很多时候这部分遗传力并没有丢失，而是由于部分snp效应太小以至于无法达到显著水平而没有被检测到。

与单SNP关联检验相对，GCTA中的GREML（genome-based restricted maximum likelihood）方法使用线性混合模型（ linear mixed models ， LMMs），将全部SNP的效应作为随机效应拟合。

模型如下：

y是 一个n × 1 的表型向量，n是样本大小； β 是固定效应的向量，诸如性别、年龄、以及一个或多个主成分（PC）； u 是SNP效应的向量，服从如下分布：

I 是一个 n × n 的单位矩阵, ɛ 是残余效应的向量，服从以下分布：

W是标准化的基因型矩阵，其中第ijth 个元素为： 

其中 xij 是第j个个体的第i个SNP的参考allele的拷贝数， pi该 allele 的频率。如果我们定义  A = WW′/N 并且定义 σ2g 为全部SNP解释的方差 , 即 σ2g=Nσ2u, （N为SNP的数量），那么方程 1 等价于：

其中g是一个n x 1的向量，表示各个个体的全部的遗传效应，服从 g∼N(0,Aσ2g) 。A 可以理解为个体两两间遗传关系组成的遗传关系矩阵 genetic relationship matrix (GRM) 。于是基于这个线性混合效应模型我们就可以通过估计的GRM，利用REML（ restricted maximum likelihood ）算法来无偏地估计全部SNP解释的方差 σ2g （继而估计SNP遗传力）。

基于以上定义，个体j与个体k的遗传关联可以由以下方程估计：

实际操作中我们首先需要排除掉有亲缘关系的个体，主要原因是该模型的目的是估计所有SNP解释的方差，如果纳入有亲缘关系的个体，会由于表型关联（ phenotypic correlations ）而造成偏差，例如由于共同的环境因素。即使没有上述的偏差，那么对它的解读也会有别于 “无关联的”个体：一个基于家系的估计值捕捉的是所有因果变异的贡献，而该方法捕捉的则是与被基因分型的SNP处于LD的因果变异的贡献。

下面简单介绍如何利用GCTA软件估计snp-遗传力：

GCTA下载链接：https://cnsgenomics.com/software/gcta/#Download

GCTA语法上类似PLINK，文件格式也几乎通用，使用起来十分便捷：

第一步： 估计GRM

gcta64 --bfile test --autosome --maf 0.01 --make-grm --out test --thread-num 10

输入为plink格式的bed/bim/fam文件，

  • –autosome 只是用常染色体上的SNP
  • –maf 0.01 过滤掉maf小于0.01的snp
  • –make-grm 计算GRM
  • –out 输出文件的前缀
  • –thread-num 10 要使用的线程数

计算完成后GRM会存储在以下文件中： test.grm.bin, test.grm.N.bin and test.grm.id

我们也可以去除掉有亲缘关系的个体：

gcta64 --grm test --grm-cutoff 0.025 --make-grm --out test_rm025

第二步：利用GCTA-GREML估计SNP遗传力

gcta64 --grm test --pheno test.phen --reml --out test --thread-num 10
  • –pheno 表型文件
  • –reml 利用reml算法计算snp遗传力

计算结果储存在 test.hsq的文件中

当然我们以可以加入协变量，例如前10个主成分 PC1-10：

gcta64 --grm test --pheno test.phen --reml --qcovar test_10PCs.txt --out test --thread-num 10

对于较大的数据量，我们还可以分染色体进行计算，详见：https://cnsgenomics.com/software/gcta/#Tutorial

参考：

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3014363/

Yang, J. et al. Common SNPs explain a large proportion of the heritability for human height. Nature Genetics 42, 565–569 (2010).

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

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

基于功能分类分割遗传力 – 分层LD分数回归 Stratified LD score regression

在复杂性状的GWAS研究中，大部分遗传力由并没有达到全基因组显著性水平的SNP关联贡献。然而当前许多利用功能信息与GWAS数据来研究疾病的方法仅仅使用达到显著性水品的loci里的SNP，并且只假设每个loci里只有一个因果SNP，或是完全不考虑LD。基于这些不足，本方法的作者们期望使用所有SNP的信息，并明确的将LD纳入模型，来估计每个功能分类的SNP遗传力，以提升power。

本方法是在基础的LD分数回归上的延伸，可以参考：

在此之前分割SNP遗传力的方法借助REML（ restricted maximum likelihood ），例如GCTA，但需要个体的基因型的原始数据，并且需要很大的计算资源。所以作者们开发了分层LD分数回归，只需要GWAS的 summary statistics ，以及从与目标群体相对应的参考群体中计算得到的LD信息。

 该方法的核心基于以下的一个事实，  GWAS中某个SNP的χ2  检验统计量包含了所有被该SNP标记的SNP的效应。因此对于一个多基因性状，一个有高LD分数的SNP，相比于低 LD分数的SNP ，总体上也会有较高的 χ2  检验统计量 。原因主要是这些SNP可能标记了单个有很大效应量的SNP，或是多个较弱效果的SNP。（也就是 “以一敌百” 或者 “人多力量大”）

所以我们将所有SNP分为具有不同遗传力的功能分类，那么与某个高遗传力分类存在LD的SNP就会有更高的 χ2  检验统计量 。如果与某个功能分类存在高LD的SNP有较高的 χ2  检验统计量 ，那么就定义这个功能分类有遗传力聚集。

在多基因模型下，某个SNP j 的 χ2  检验统计量 的期望值为

  • N为样本大小
  • C则是功能分类的索引
  • l(j,C)则是SNP j 对于分类C的LD分数
  • a则衡量了混淆因素的大小
  • τC则表示了每个SNP对于功能分类C遗传力的贡献。

该方程使我们可以估计 τC 的大小（也就是所谓的分割的遗传力）。定义某个分类的聚集为 该分类SNP遗传力的比例除以该分类SNP总数的比例。

该方法的作者们基于多个公开的注释数据库，构建了不针对任何细胞类型的全基线模型 ‘full baseline model’ ， 包括 coding, UTR, promoter and intronic regions the histone marks monomethylation (H3K4me1) and trimethylation (H3K4me3) of histone H3 at lysine 等等。除此之外，还基于 全基线模型 ，构建了多个针对特定细胞类型的模型，包含针对细胞类型的注释等。

尽管Stratified LD score regression提供了一种便捷有效的分析 GWAS的 summary statistics 的方法，但我们也要同时注意该方法的不足之处：

  • 为了达到足够的power，需要较大的样本量，或是较大的SNP遗传力，而且性状必须是多基因的
  • 该方法需要针对研究群体的LD参考数据
  • 该方法目前不支持自定义的array
  • 该方法基于加性模型，没考虑上位或非加性的效应
  • 该方法依赖于可用的功能注释，如果没有相应注释则无法检测
  • 等等

使用方法：

下载，安装，配置环境，数据清理详见

  • 连锁不平衡分数回归 LD score regression
    • 除了以上步骤外，我们还需要下载相应的baseline模型：

    https://alkesgroup.broadinstitute.org/LDSCORE/

    从以上链接中我们需要下载以下内容（以欧洲群体为例）：

    • 基线模型LD分数 ：baseline.* in 1000G_Phase1_baseline_ldscores.tgz
    • 频率 ：1000G.mac5eur.* in 1000G_Phase1_frq.tgz
    • 权重 ：weights.* in weights_hm3_no_hla.tgz

    解压后即可使用：

    python ldsc.py 
	--h2 BMI.sumstats.gz\
	--ref-ld-chr baseline.\ 
	--w-ld-chr weights.\
	--overlap-annot\
	--frqfile-chr 1000G.mac5eur.\
	--out BMI_baseline
    • –h2 ： 计算分割的遗传力，参数为之前处理好的gwas summary statistics
    • --ref-ld-chr ：下载的参考LD分数文件
    • --w-ld-chr： 权重文件
    • --frqfile-chr： SNP频率文件
    • --overlap-annot： 表示基线模型中功能分类有重叠
    • –out：指定出输出文件的前缀

    参考：

    http://www.github.com/bulik/ldsc

    https://github.com/bulik/ldsc/wiki/Partitioned-Heritability

    Partitioning heritability by functional annotation using genome-wide association summary statistics https://www.nature.com/articles/ng.3404

    GWAS的线性混合模型LMM Linear Mixed Model

    关键词：population structure，cryptic relatedness，spurious association，LMM，GRM

    早期的GWAS研究使用的模型为固定效应的线性模型，但通常会受两方面混淆因素的影响：

    1. 群体结构/分层 （population structure/ stratification）：研究群体中存在有不同祖先（ancestry）的亚群体（subgroup）
    2. 隐性关联（cryptic relatedness）：研究样本之间存在未知的亲缘关系

    如果在模型中故意忽略掉这些混淆变量（confounding variable），就很可能导致结果出现假阳性（false positive）或是虚假关联（spurious association），所以这是我们在进行GWAS研究时，必须要考虑的问题。

    Population structure

    首先来看群体结构/分层

    考虑一个case-control研究，如下图所示，红色的群体在整体样本中占了case的大多数，那么一些对疾病并没有影响，但在红蓝两群体之间等位基因频率相差很大的genetic marker，就有可能会造成虚假关联（spurious association）。

    在数量形状也是存在类似的情况，

    为了解决这个问题，我们通常会采用genomic control的方法，或是通过PCA来矫正。详见：

    Cryptic relatedness

    由于将主成分作为协变量纳入模型以矫正群体分层的方法，适用于样本中不存在亲缘关系的情况下使用，而样本中存在亲缘关系时，就不一定有效，但通常我们的研究中都会存在有亲缘关系的样本。

    当样本中存在隐性关联 / 错误认定的关联 Cryptic and / or misspecified relatedness时，也可能会造成上述的虚假关联。由于我们通常不能掌握所有研究对象的家谱，所以无法完全去相关个体，例如下面的情况，样本中看似不相关的个体间实际上存在隐性关联：

    Linear Mixed Model

    为了解决以上问题，相比早期的线性模型，我们就采用一种更为灵活的模型，线性混合模型来进行检验。

    Y是 n x 1的向量，表型数据

    固定效应 fixed effects：

    • W为 n x （w + 1）的矩阵，包含了截距，以及协变量
    • β 则是 ( w + 1 ) × 1 的向量，表示协变量的效应量
    • Gs为 n × 1 的向量，某个位点的基因型，每一项的值通常为 0，1，2（等位基因allele的拷贝数）
    • γ则是标量，目标位点基因型的效应量

    随机效应 random effects：

    • g为长度为n的随机向量，表示多基因效应
    • δ2g为加性遗传方差
    • Ψ是成对遗传相关的矩阵
    • e是长度为n的随机向量，服从以下的分布，其中δ2e为非遗传效应造成的方差，我们认为这项在个体间是相互独立的
    这样我们就构建了包含群体结构和隐形相关的模型。

    通常成对遗传相关的矩阵Ψ是未知的，在计算中我们使用empirical GRM（genetic relatedness matrix，通过纳入研究的SNP计算得出）

    目前已有多种GWAS检验方法基于LMM模型，包括：

    • EMMAX
    • GEMMA
    • GMMAT
    • Bolt-LMM
    • fastGWAS
    • SAIGE
    • 等等

    参考：

    Lecture 6: GWAS in Samples with Structure ，Summer Institute in Statistical Genetics 2015

    GWAS线性混合模型中的LOCO Leave-one-chromosome-out

    关键词：LMM，proximal contaminal, LOCO

    目前的GWAS已经开始逐渐使用线性混合模型（ linear mixed models ，LMM）来代替早期的线性模型，主要原因是线性混合模型能够校正多种原因造成的混淆，例如遗传关联（ genetic relatedness ），家庭关联（ familial relatedness ），群体分层（ population structure ）等，LMM模型也因此能够控制假阳性，并提高检验power。

    但在使用混合线性模型中一个重要的问题就是，当我们在GRM中纳入了被检验的SNP时，反而会导致power降低。原因是在模型中我们对待检测SNP进行了二重拟合（ double-fitting ），即:

    • 1 . 作为检验关联时的固定效应（fixed effect）
    • 2. 在GRM中作为随机效应 （random effect）

    这种现象就被称为 临近污染 “proximal contamination”。

    为了避免此现象造成的power损失，理论上在构建null模型中排除掉待检验SNP是正确的做法，但这样太占运算资源，所以在实践中，我们会采用 LOCO Leave-one-chromosome-out ，即使用排除掉待检验SNP所在的染色体的所有SNP，再进行检验（也就是说我们有对应22个常染色体的loco null模型）。

    目前主流的软件都已支持loco，只需要–loco 指定即可。

    参考：

    https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3597090/

    Advantages and pitfalls in the application of mixed-model association methods

    GWAS的条件分析 Conditional analysis

    关键字：conditional analysis，leading SNP, secondary casual variants , fine-mapping

    GWAS研究中，对于某个复杂表型，我们会发现很多显著关联的基因座（loci），每个基因座里有若干显著的SNP，这些SNP通常处于LD。这时我们就会面临一个问题，这个基因座里显著的SNP单单因为与leading SNP（该loci里P值最低的SNP）连锁不平衡而显著，还是因为这个SNP本身就与就与表型相关联。这时就应该进行条件分析（conditional analysis），实际上是一种fine-mapping（对casual SNP 致病SNP的精确定位）的方法，目的是确认是否存在次要的因果SNP secondary causal variants。

    一般来说进行条件分析的方法很简单，

    1.从GWAS结果中抽出每个关联基因座的leading SNP，

    2.将该leading SNP作为协变量加入检验模型

    3.再次进行关联检验，确认 关联基因座里除 leading SNP 以外还是否有次要的致病SNP secondary causal variants。

    例如下图所示的情况：

    图1：针对橙色圈中的leading SNP （该loci里p值最低），进行条件分析后的曼哈顿图。A：该loci只有一个信号，表示没有 secondary causal variants ，之所以有多个显著snp是因为与leading SNP处于连锁不平衡。B：除了leading SNP 还有其他信号，表示这个loci还存在 secondary causal variants 。

    通常conditional analysis的功能已经集成在关联检验的软件中（如PLINK，SAIGE），我们只需要提供lead snp的 id 或是 位置，就可以进行条件分析了。

    例如：

    PLINK中的 –condition 选项，也可以使用 –condition–list(同时将多个SNP作为covariate纳入检验模型)

    SAIGE 第二步中，也提供了 –condition 选项 ，可以对单个或多个snp进行条件分析

    参考：

    Strategies for fine-mapping complex traits.

    https://www.cog-genomics.org/plink/1.9/assoc

    https://github.com/weizhouUMICH/SAIGE/wiki/Genetic-association-tests-using-SAIGE