用 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/bashconda
git clone https://github.com/bulik/ldsc.git
cd ldscconda 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.0137615head 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.bz2tar -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/