Heritability – GWASLab – GWAS实验室

通过Bivariate LD Score regression估计遗传相关性 genetic correlation

背景回顾：

LD score 回归：连锁不平衡分数回归 LD score regression

双变量/跨性状LD分数回归（Bivariate/Cross-trait LD Score regression）

简介与原理：

该方法使用GWAS的概括性数据来估计一对表型的遗传相关。该方法基于以下的事实：GWAS检验中，对一个SNP效应量的估计通常也会包含与该SNP成LD的其他SNP的效应，也就是说一个与其他SNP成高LD的SNP，通常也会有更高的卡方检验量。而当我们把卡方检验量替换为，来自两个相关表型的GWAS的Z分数的乘积时，这个事实依然存在。所以基于此，我们可以将单表型的LD分数回归扩展成如下所示的 双变量/跨性状LD分数回归 ：

详细推导见：

Bulik-Sullivan, B., Finucane, H., Anttila, V. et al. An atlas of genetic correlations across human diseases and traits. Nat Genet 47, 1236–1241 (2015). https://doi.org/10.1038/ng.3406

下面简单介绍使用方法：

数据下载，数据清理等步骤与单表型的LD分数回归基本相同，详见连锁不平衡分数回归 LD score regression

我们只需要将单表型的LD分数中 -h2 的选项改为 -rg，再加上要估计遗传相关的两个表型的数据清理后的文件即可：

ldsc.py \
--rg scz.sumstats.gz,bip.sumstats.gz \
--ref-ld-chr eur_w_ld_chr/ \
--w-ld-chr eur_w_ld_chr/ \
--out scz_bip

`--rg` ：计算遗传相关性（genetic correlation），参数为上一步处理好的数据文件名

`--ref-ld-chr` ：使用的按染色体号分类的LD score文件名，参数为LD score文件所在文件夹的路径。默认情况下ldsc会在文件名末尾添加染色体号，例如 `--ref-ld-chr eur_w_ld_chr/` 的意思就是使用 `eur_w_ld_chr/1.l2.ldscore ... eur_w_ld_chr/22.l2.ldscore 这些文件`。如果你的染色体号在其他位置，也可以使用@来告诉LDSC，例如 `--ref-ld-chr ld/chr@` 。当然也可以用 `--ref-ld` 来指定一个整体的LD score文件。

`--w-ld-chr`：指定ldsc回归权重所用的LDscore文件，理论上对于SNPj的LD分数，应当包含这个SNP与其他所有SNP的R2之和，但实际操作中，LD score回归对于计算权重的SNP的选择并不敏感，所以一般情况下我们可以使用与 `--ref-ld` 与相同的文件。.

`--out` ：输出文件的路径与前缀

运行后我们就能得到如下结果：

Genetic Covariance
------------------
Total Observed scale gencov: 0.3644 (0.0368)
Mean z1*z2: 0.1226
Intercept: 0.0037 (0.0071)

Genetic Correlation
-------------------

Genetic Correlation: 0.6561 (0.0605)
Z-score: 10.8503
P: 1.9880e-27

Summary of Genetic Correlation Results
              p1               p2     rg    se       z          p  h2_obs  h2_obs_se  h2_int  h2_int_se  gcov_int  gcov_int_se
 scz.sumstats.gz  bip.sumstats.gz  0.656  0.06   10.85  1.988e-27   0.522      0.053   1.001      0.009     0.004        0.007

Genetic Correlation: 0.6561 (0.0605) 即为我们需要估计的遗传相关

除此之外，LDSC还会输出一个Summary of Genetic Correlation Results，当有多组遗传相关需要估计时，查看该表格会更加方便。

常见问题：

为什么rg会在[-1,1]的范围以外: （https://github.com/bulik/ldsc/issues/78）

ldsc的算法并没有限定[-1,1]的范围，当rg接近一时，加上估计误差有时就会导致rg在范围外

参考：

Bulik-Sullivan, B., Finucane, H., Anttila, V. et al. An atlas of genetic correlations across human diseases and traits. Nat Genet 47, 1236–1241 (2015). https://doi.org/10.1038/ng.3406

https://github.com/bulik/ldsc/wiki/Heritability-and-Genetic-Correlation

解释复杂疾病的四种主流模型 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甲基化的继承。

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/

易感性尺度遗传力与观测尺度遗传力 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

使用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是标准化的基因型矩阵，其中第ij^th 个元素为：

其中 x_ij 是第j个个体的第i个SNP的参考allele的拷贝数， p_i该 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

遗传力 Heritability 与 Missing heritability

本文关键词： Heritability， h2, family heritability , SNP heritability, GWAS heritability , Missing heritability

遗传力的基础概念

遗传力（Heritability）是我们理解遗传与环境因素对性状影响的基础，定义为遗传方差占性状方差（总方差）的比值，可以理解为遗传因素对性状的影响，数学上以h2表示。

通常根据对遗传方差的定义而分为广义与狭义遗传力：

广义遗传力 broad-sense heritability ：

其中VP = VG + VE, VP为性状方差，VE为环境方差（也包括测量误差等），而分子的VG为遗传方差。

遗传方差VG也可进一步细分为

VA是加性遗传效应的方差，VNA是非加性遗传效应的方差（上位与显性遗传效应）

加性遗传效应是指当两个或多个基因对于某一性状，或是单个基因的不同等位基因对于某一性状的整体作用，等于它们单独作用之和。

非加性遗传效应则包括上位与显性遗传效应。

因为对于绝大多数复杂性状，很少有证据证明有非加性遗传效应存在，所以我们目前聚焦于主要考虑加性遗传效应的狭义遗传力，

狭义遗传力 narrow-sense heritability ：

遗传力的估计

我们有多种方法可以估计遗传力h2的大小，目前主要的方法有三种，通过双胞胎研究，SNP或是GWAS来估计。

h2 family ：双胞胎研究，通过比较同卵与异卵双胞胎的相似性，计算得到h2，通常为这三种中最高。

h2 SNP ：GWAS研究所用chip上所有variants共同解释的方差与性状方差的比值，比 h2 family 低，但会显著高于h2 GWAS。可以使用GCTA的GREML模型来估计。（使用GCTA （GREML）来估计SNP-遗传力 SNP Heritability ）

h2 GWAS ：仅由GWAS所发现的某疾病相关variants解释的方差与性状方差的比值，三者中最低。

一般情况下，三者的关系：

更直观的关系如下图所示：

图一：三种遗传力的关系。（引自：The contribution of genetic variants to disease depends on the ruler，2014）

Missing and hidden heritability

GWAS研究中一个核心问题便是 Missing heritability，定义为 h2 family 与 h2 GWAS 之间的差值。 h2 GWAS 之所以低于 h2 family ，潜在的原因包括：非加性遗传效应（尽管目前证据很少），效应量大的稀有变异（rare variants），或是双胞胎研究中由于共同的环境因素而造成的过高估计。

Missing heritability 又可细分为 still- missing heritability 与 hidden heritability 。 still- missing heritability 为 h2 family 与 h2 SNP 之间的差，Yang 认为可能的原因是在GWAS研究中由于样本数量的限制，大多数效应量较小的遗传效应无法被可靠地检测。

而 hidden heritability 则为 h2 SNP 与 h2 GWAS 的差。对它的理解建立在Fisher最初对于无穷小模型（infinitesimal model），即多数变异都只有很小的效应。在GWAS研究中，由于我们所选显著阈值的高低，遗传力或许并不是消失（ missing ）而是被隐藏（ hidden ）了。另一种可能则是，人群的异质性（heterogeneity.），因为 h2 GWAS 大多来自包含多群体的meta分析，而遗传效应在这些群体中的异质性也可能使 h2 GWAS 偏低。

如何估计SNP遗传力：

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

参考：

An Introduction to Statistical Genetic Data Analysis.

--rg ：计算遗传相关性 （genetic correlation），参数为上一步处理好的数据文件名

--out ： 输出文件的路径与前缀

`--rg` ：计算遗传相关性（genetic correlation），参数为上一步处理好的数据文件名

`--out` ：输出文件的路径与前缀