用 LDSC 计算遗传度以及遗传相关性 |
您所在的位置：网站首页 › h2指什么 › 用 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 安装 LDSCDocker代码语言：javascript复制docker pull rticode/ldsc:7618f4943d8f31a37cbd207f867ba5742d03373f # 改 tag docker tag rticode/ldsc:7618f4943d8f31a37cbd207f867ba5742d03373f ldsc:latest创建容器： 代码语言：javascript复制docker run --rm -v $(pwd):/data --name=ldsc -it ldsc:latest /bin/bashconda代码语言：javascript复制git clone https://github.com/bulik/ldsc.git cd ldsc代码语言：javascript复制conda env create --file environment.yml source activate ldsc下载示例数据下载 GWAS 结果数据进入 PGC 官网：https://www.med.unc.edu/pgc/download-results/ 选择 Cross Disorder ，点击下载：  我们需要用到 biopolar disorder 和 Schizophrenia 这两个数据集，在下拉菜单中分别选择，进行下载。  代码语言：javascript复制unzip -o pgc.cross.bip.zip
unzip -o pgc.cross.scz.zip示例数据内容如下： 代码语言：javascript复制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代码语言：javascript复制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代码语言：javascript复制wget https://data.broadinstitute.org/alkesgroup/LDSCORE/eur_w_ld_chr.tar.bz2 wget https://data.broadinstitute.org/alkesgroup/LDSCORE/w_hm3.snplist.bz2代码语言：javascript复制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 参数）。 根据上面的步骤，我们进行格式化的命令如下： 代码语言：javascript复制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 的基本信息： 代码语言：javascript复制********************************************************************** * 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 **********************************************************************下一部分为运行的命令参数： 代码语言：javascript复制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 正确识别了列名。如果不确定脚本是否能够正确识别列名，最简单的方法就是直接运行脚本，如果不能识别命令就会报错。 代码语言：javascript复制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 |
【本文地址】