翻译整理自:https://github.com/bulik/ldsc/wiki/Heritability-and-Genetic-Correlation
本教程包含用 ldsc
计算:
•精神分裂症的 LD Score 回归截距•精神分裂症的 SNP 遗传度•精神分裂症和躁郁症之间的遗传相关性
使用2013 年发表在柳叶刀上的 PGC 文章[1]数据为例。
GitHub:https://github.com/bulik/ldsc
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
git clone https://github.com/bulik/ldsc.git
cd ldsc
conda env create --file environment.yml
source activate ldsc
进入 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
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 回归。
有了 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/