这篇文章是对以前啊啊救救我,为什么个人QQ图那么飘(全基因组关联分析)这篇文章的一个补坑。git
LD SCore除了查看显著SNP位点对表型是否为基因多效性外,还额外补充了怎么计算表型的遗传度和遗传相关性。github
git clone https://github.com/bulik/ldsc.git
less
cd ldsc
测试
conda env create --file environment.yml
ui
source activate ldsc
3d
若是安装成功,输入./ldsc.py -h
代码会出现以下图:
code
输入./munge_sumstats.py -h
代码会出现以下图:
blog
summary.txt
为关联分析的summary数据,包含rs编号、染色体编号、位置、A1(效应等位基因)、A2(无效等位基因)、效应值(OR或BETA)、P值,以下图所示:token
munge_sumstats.py --sumstats summary.txt --N 17115 --out scz --merge-alleles w_hm3.snplist
ci
这里的N指的是研究的样本数量;
scz是输出的文件名;
w_hm3.snplist是被归入分析的SNP,包含三列:包含rs编号、位置、A1(效应等位基因)、A2(无效等位基因)# 这一步无关紧要#
若是想把全部的SNP位点归入分析,那么采用这个命令:
munge_sumstats.py --sumstats summary.txt --N 17115 --out scz
这一步会生成scz.sumstats.gz的文件;
for q in $(seq 1 22); do plink --bfile file --chr $q --make-bed --out chr$q done
这个步骤会生成22个plink格式文件(bed,bim,fam),每个文件表明一条染色体。
for q in $(seq 1 22); do ldsc.py --bfile chr$q --l2 --ld-wind-cm 5 --yes-really --out chr/$q done
生成的文件以下所示:
ldsc.py --h2 scz.sumstats.gz --ref-ld-chr chr/ --w-ld-chr chr/ --out scz_h2
scz.sumstats.gz为步骤5生成的文件
chr/ 为步骤7生成的LD文件路径
scz_h2为回归截距和遗传度的输出文件
less scz_h2.log
输出文件最底部:
Intercept: 1.0252 (0.0075)
截距为1.0252
关于回归截距怎么看,请看以前发过的推文:啊啊救救我,为什么个人QQ图那么飘(全基因组关联分析)
less scz_h2.log
输出文件最底部:
Total Observed scale h2: 0.7153 (0.0386)
遗传度为0.7153
ldsc.py --rg trait1.sumstats.gz,trait2.sumstats.gz --ref-ld-chr chr/ --w-ld-chr chr/ --out trait1_trait2
trait1.sumstats.gz为表型1的ldsc格式文件;
trait2.sumstats.gz为表型2的ldsc格式文件;
chr/ 为步骤7生成的LD文件路径
trait1_trait2为表型1和表型2的遗传相关性输出文件;
less trait1_trait2.log
输出文件最底部:
Genetic Correlation: 0.6561 (0.0605)
表型1和表型2的遗传相关性为0.6561