用 LDSC 计算遗传度以及遗传相关性

翻译整理自：https://github.com/bulik/ldsc/wiki/Heritability-and-Genetic-Correlation

本教程包含用 ldsc 计算：

•精神分裂症的 LD Score 回归截距•精神分裂症的 SNP 遗传度•精神分裂症和躁郁症之间的遗传相关性

使用2013 年发表在柳叶刀上的 PGC 文章[1]数据为例。

GitHub：https://github.com/bulik/ldsc

安装 LDSC

Docker

docker pull rticode/ldsc:7618f4943d8f31a37cbd207f867ba5742d03373f

# 改 tag
docker tag rticode/ldsc:7618f4943d8f31a37cbd207f867ba5742d03373f ldsc:latest

创建容器：

docker run --rm -v $(pwd):/data --name=ldsc -it  ldsc:latest /bin/bash

conda

git clone https://github.com/bulik/ldsc.git
cd ldsc

conda env create --file environment.yml
source activate ldsc

下载示例数据

下载 GWAS 结果数据

进入 PGC 官网：https://www.med.unc.edu/pgc/download-results/

选择 Cross Disorder ，点击下载：

我们需要用到 biopolar disorder 和 Schizophrenia 这两个数据集，在下拉菜单中分别选择，进行下载。

unzip -o pgc.cross.bip.zip
unzip -o pgc.cross.scz.zip

示例数据内容如下：

head pgc.cross.BIP11.2013-05.txt

snpid hg18chr bp a1 a2 or se pval info ngt CEUaf
rs3131972    1    742584    A    G    1.092    0.0817    0.2819    0.694    0    0.16055
rs3131969    1    744045    A    G    1.087    0.0781    0.2855    0.939    0    0.133028
rs3131967    1    744197    T    C    1.093    0.0835    0.2859    0.869    0    .
rs1048488    1    750775    T    C    0.9158    0.0817    0.2817    0.694    0    0.836449
rs12562034    1    758311    A    G    0.9391    0.0807    0.4362    0.977    0    0.0925926
rs4040617    1    769185    A    G    0.9205    0.0777    0.2864    0.98    0    0.87156
rs28576697    1    860508    T    C    1.079    0.2305    0.7423    0.123    0    0.74537
rs1110052    1    863421    T    G    1.088    0.2209    0.702    0.137    0    0.752294
rs7523549    1    869180    T    C    1.823    0.8756    0.4929    0.13    0    0.0137615

head pgc.cross.SCZ17.2013-05.txt

snpid hg18chr bp a1 a2 or se pval info ngt CEUaf
rs3131972    1    742584    A    G    1    0.0966    0.9991    0.702    0    0.16055
rs3131969    1    744045    A    G    1    0.0925    0.9974    0.938    0    0.133028
rs3131967    1    744197    T    C    1.001    0.0991    0.9928    0.866    0    .
rs1048488    1    750775    T    C    0.9999    0.0966    0.9991    0.702    0    0.836449
rs12562034    1    758311    A    G    1.025    0.0843    0.7716    0.988    0    0.0925926
rs4040617    1    769185    A    G    0.9993    0.092    0.994    0.979    0    0.87156
rs4970383    1    828418    A    C    1.096    0.1664    0.5806    0.439    0    0.201835
rs4475691    1    836671    T    C    1.059    0.1181    0.6257    1.02    0    0.146789
rs1806509    1    843817    A    C    0.9462    0.1539    0.7193    0.383    0    0.600917

下载欧洲人群 LD Score 以及 hapmap3 snplist

wget https://data.broadinstitute.org/alkesgroup/LDSCORE/eur_w_ld_chr.tar.bz2
wget https://data.broadinstitute.org/alkesgroup/LDSCORE/w_hm3.snplist.bz2

tar -jxvf eur_w_ld_chr.tar.bz2
bunzip2 w_hm3.snplist.bz2

重新格式化摘要统计数据

在执行 ldsc 前，首先需要使用 munge_sumstats.py 脚本将 GWAS 结果转换为 ldsc 格式，ldsc 的 .sumstats 格式包含六列：

•SNP ID•A1•A2•样本量•P 值•SNP 对表型的效应值，beta，OR，log odds，z-score 等等

Imputation 质量与 LD Score 相关，若质量较差则会导致统计值较低。因此，我们通常根据 INFO > 0.9 进行过滤。本教程的 scz 和 bip 示例数据都包含 INFO 列，munge _ sumstats.py 脚本会自动根据这一列的值进行过滤。若 GWAS 结果中没有 INFO 列，我们建议根据 HapMap3 snplist 进行过滤（用 --merge 或 --merge-alleles 参数，见示例代码）。

由于示例数据中没有描述样本量的列，所以在格式化时需要加上 --N 参数定义样本大小。本例中，scz 研究的样本量为 17115，bip 研究的样本量为 11810。

同时，我们最好检查一下 GWAS 结果文件中列出的等位基因是否与用于计算 LD Score 的数据中列出的等位基因相匹配（使用 --merge-alleles 参数）。

根据上面的步骤，我们进行格式化的命令如下：

munge_sumstats.py 
--sumstats pgc.cross.SCZ17.2013-05.txt 
--N 17115 
--out scz 
--merge-alleles w_hm3.snplist

munge_sumstats.py 
--sumstats pgc.cross.BIP11.2013-05.txt 
--N 11810 
--out bip 
--merge-alleles w_hm3.snplist

每个命令大约会运行 20 秒。命令输出包括： scz.log，scz.sumstats.gz，bip.log 和 bip.sumstats.gz。

格式化步骤的日志文件内容

第一部分为 ldsc 的基本信息：

**********************************************************************
* LD Score Regression (LDSC)
* Version 1.0.0
* (C) 2014-2015 Brendan Bulik-Sullivan and Hilary Finucane
* Broad Institute of MIT and Harvard / MIT Department of Mathematics
* GNU General Public License v3
**********************************************************************

下一部分为运行的命令参数：

Call: 
./munge_sumstats.py 
--out scz 
--merge-alleles w_hm3.snplist 
--N 17115.0 
--sumstats pgc.cross.SCZ17.2013-05.txt

在下一部分中将描述脚本识别列名的情况。默认情况下，munge_sumstats.py 可以识别绝大多数列名，但如果 GWAS 结果中包含一些特殊列名，你可能需要用参数进行指定。比如，输入中的 foobar 列包含 INFO 得分，则命令应改为 munge_sumstats.py -- INFO foobar。所以我们应检查日志文件的这一部分内容，以确保 munge_sumstats.py 正确识别了列名。如果不确定脚本是否能够正确识别列名，最简单的方法就是直接运行脚本，如果不能识别命令就会报错。

Interpreting column names as follows:
info:   INFO score (imputation quality; higher --> better imputation)
snpid:  Variant ID (e.g., rs number)
a1:     Allele 1
pval:   p-Value
a2:     Allele 2
or:     Odds ratio (1 --> no effect; above 1 --> A1 is risk increasing)

下一节中描述了质控的过程。默认情况下，munge_sumstats.py 会根据 INFO > 0.9，MAF > 0.01 以及 0 < P <= 1 进行过滤。它还会删除那些不是 SNPs 的突变（例如 indels）、去除不是 2 个等位位点的 SNPs 以及重复 rs 号的 SNPs。最后，脚本还会检查 SNP 的效应值中位数是否接近无意义的中位数（例如，OR 接近1），以确保该列没有被标记错误。

Reading list of SNPs for allele merge from w_hm3.snplist
Read 1217311 SNPs for allele merge.
Reading sumstats from pgc.cross.SCZ17.2013-05.txt into memory 5000000.0 SNPs at a time.
Read 1237958 SNPs from --sumstats file.
Removed 137131 SNPs not in --merge-alleles.
Removed 0 SNPs with missing values.
Removed 256286 SNPs with INFO <= 0.9.
Removed 0 SNPs with MAF <= 0.01.
Removed 0 SNPs with out-of-bounds p-values.
Removed 2 variants that were not SNPs or were strand-ambiguous.
844539 SNPs remain.
Removed 0 SNPs with duplicated rs numbers (844539 SNPs remain).
Using N = 17115.0
Median value of or was 1.0, which seems sensible.
Removed 39 SNPs whose alleles did not match --merge-alleles (844500 SNPs remain).
Writing summary statistics for 1217311 SNPs (844500 with nonmissing beta) to scz.sumstats.gz.

最后一部分描述了有关 GWAS 结果的一些基本数据。如果平均卡方值小于 1.02，脚本将警告该数据可能不适合进行 LD Score 回归。

计算遗传度，遗传相关性以及 LD Score 回归截距

有了 ldsc 格式的输入文件后，我们就可以用下面的命令进行计算遗传度和遗传相关性：

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

若只是计算遗传度可用下面的命令：

ldsc.py 
--h2 scz.sumstats.gz 
--ref-ld-chr eur_w_ld_chr/ 
--w-ld-chr eur_w_ld_chr/ 
--out scz_h2

--h2 和 --rg 的输出结果非常相似，所以在下面的教程中我们将以 --rg 的输出结果为例。注意，这里使用 --h2 计算的遗传度和 LD Score 回归截距与使用 --rg 得到的遗传度和截距会略有不同。这是因为 --rg scz.sumstats.gz,bip.sumstats.gz 使用的 SNP 是scz.sumstats.gz 和 bip.stats.gz 的交集进行 LD Score 回归，而 --h2 scz.sumstats.gz 则用的是 scz.sumstats.gz 中的所有 SNPs。

运行上面的命令大约需要一分钟，下面为该命令的主要参数：

•--rg：指定 ldsc 计算遗传相关性。后面跟的输入文件应为 .sumstats 格式。在本例中我们只输入了两个文件，但若我们输入更多文件，ldsc.py 将计算第一个文件和之后所有文件之间的遗传相关性（例如--rg a,b,c 将计算 rg(a,b) 和 rg(a,c)）。•--ref-ld-chr：指定每个染色体的 ld score 目录。默认情况下 ldsc 会自动识别文件夹下的染色体文件，例如，输入 --ref-ld-chr eur_w_ld_chr/ ，ldsc 会识别 eur_w_ld_chr/1.l2.ldscore, ... , eur_w_ld_chr/22.l2.ldscore 。若染色体号在文件名中间，可以在染色体号插入的地方加上 @ 占位符，例如 --ref-ld-chr ld/chr@ 。

如果你的研究也是用欧洲人的 GWAS 数据，那就可以直接用 eur_w_ld_chr/ 的 LD Score，不必自己重新计算。计算 LD Score：

python ldsc.py
    --bfile 22
    --l2
    --ld-wind-cm 1
    --out 22

•--w-ld-chr：指定回归分析中每个 SNP 位点的权重。理想情况下，SNP j 的 --w-ld LD Score 应为 r^2_jk 回归中包括的所有 SNP k 的总和。在上面的例子中，我们直接对 --ref-ld 和 --w-ld 使用相同的 LD Score，因为算法对这个权重不敏感。•--out：指定输出文件前缀。例如，设置 --out foo_bar，ldsc 会将结果输出到 foo_bar.log。若不设置这个参数，ldsc 会默认将结果输出到 ldsc.log。

结果文件内容

输出的结果包含以下六个方面：

•输入文件的日志信息；•第一个表型的遗传度（在本例中为 scz）；•第二个表型的遗传度（在本例中为 bip）；•遗传协方差；•遗传相关性；•遗传相关性表。

第一部分为 ldsc 的基本信息以及运行的命令参数：

*********************************************************************
* LD Score Regression (LDSC)
* Version 1.0.0
* (C) 2014-2015 Brendan Bulik-Sullivan and Hilary Finucane
* Broad Institute of MIT and Harvard / MIT Department of Mathematics
* GNU General Public License v3
*********************************************************************

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

下一节的内容包含基本的日志信息以及汇总统计信息。这里需要检查 SNPs 数量是否发生了异常下降。

Beginning analysis at Mon Apr  4 14:24:29 2016
Reading summary statistics from scz.sumstats.gz ...
Read summary statistics for 844500 SNPs.
Reading reference panel LD Score from data/[1-22] ...
Read reference panel LD Scores for 1293150 SNPs.
Removing partitioned LD Scores with zero variance.
Reading regression weight LD Score from data/[1-22] ...
Read regression weight LD Scores for 1293150 SNPs.
After merging with reference panel LD, 840450 SNPs remain.
After merging with regression SNP LD, 840450 SNPs remain.
Computing rg for phenotype 2/2
Reading summary statistics from bip.sumstats.gz ...
Read summary statistics for 1217311 SNPs.
After merging with summary statistics, 840450 SNPs remain.
803373 SNPs with valid alleles.

接下来的两个部分描述了每个表型的遗传度：

Heritability of phenotype 1
---------------------------
Total Observed scale h2: 0.5907 (0.0484)
Lambda GC: 1.2038
Mean Chi^2: 1.2336
Intercept: 1.0014 (0.0113)
Ratio: 0.0059 (0.0482)

Heritability of phenotype 2/2
-----------------------------
Total Observed scale h2: 0.5221 (0.0531)
Lambda GC: 1.1364
Mean Chi^2: 1.1436
Intercept: 1.0013 (0.0094)
Ratio: 0.0093 (0.0652)

下一节显示了遗传协方差：

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

遗传相关性、 z 值和 p 值：

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

引用链接

[1] 2013 年发表在柳叶刀上的 PGC 文章: https://pubmed.ncbi.nlm.nih.gov/23453885/