专栏首页生信菜鸟团用 LDSC 计算遗传度以及遗传相关性

用 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 disorderSchizophrenia 这两个数据集,在下拉菜单中分别选择,进行下载。

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.logscz.sumstats.gzbip.logbip.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.gzbip.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_barldsc 会将结果输出到 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/

本文分享自微信公众号 - 生信菜鸟团(bio_123456789),作者:鲍志炜

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-07-26

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • R tips:使用!!来增加dplyr的可操作性

    dplyr包在数据变换方面非常的好用,它有很多易用性的体现:比如书写数据内的变量名时不需要引号包裹,也不需要绝对引用,而这在多数baseR函数中都不是这样的,比...

    生信菜鸟团
  • 一文读懂KEGG数据库

    在进行生物学实验或者生物信息的学习中,都会听说KEGG富集分析,而且该方法在高通量测序分析中已然成为数据分析中必不可少的一环。

    生信菜鸟团
  • 一文极速读懂 Gene Ontology (GO)数据库

    一、介绍1、分子功能(Molecular Function,MF )2、细胞组分(Cellular Component ,CC)3、生物过程(Biologic...

    生信菜鸟团
  • 重温统计学之五——归一化

    标准正态分布(Standard Normal Distribution):标准正态分布式一个特殊的正态分布。其随机变量均值为0,标准偏差为1。普通的随机变量在标...

    统计学家
  • 二、DNS如何工作?

      如前所述,DNS 的核心是一个分层系统。在这个系统的顶部是所谓的 “根服务器”。这些服务器由各种组织控制,并由 ICANN(互联网名称和数字地址分配公司)授...

    咻一咻
  • Linux中Homebrew的正确使用方法

    首先要避免将 Homebrew 的 bin 目录添加到PATH ,而仅仅将你需要使用的几个可执行做软连接放到~/bin 下面(这个目录在PATH 中),以避免环...

    砸漏
  • 【Excel系列】Excel数据分析:相关与回归分析

    相关系数 15.1 相关系数的概念 著名统计学家卡尔·皮尔逊设计了统计指标——相关系数(Correlation coefficient)。相关系数是用以反映变量...

    数据科学社区
  • electron 快速开始

    阅读 https://github.com/electron/electron 要运行一个最小demo,只需

    mafeifan
  • .NET 的文本转语音合成

    我经常飞去芬兰见我的妈妈。每次飞机降落在万塔机场时,我都会对鲜有旅客前往机场出口感到惊讶。绝大多数的旅客会转机到跨越所有中欧及东欧的目的地。所以难怪在飞机开始下...

    Edison.Ma
  • 你离黑客的距离,就差这20个神器了

    在不少电影电视剧中,主角的身边都有这么一位电脑高手:他们分分钟可以黑进反派的网络,攻破安全防线,破解口令密码,拿到重要文件。他们的电脑屏幕上都是一些看不懂的图形...

    轩辕之风

扫码关注云+社区

领取腾讯云代金券