使用LDSC计算LD分数 – LD score

目录

基于PLINK文件计算LD分数
LD分数文件格式
用于LDSC-SEG的分层LD分数的计算

有很多时候没有现成的LD分数参考,而需要我们自己根据手头基因型数据而计算LD分数(例如一个常用的例子是GWAS数据与scRNA数据的结合 – 使用自己的scRNA数据提取特征基因后计算LD分数并进行partitioned LD 分数回归以检测细胞特异性)。

本文将简要介绍如何计算 单变量的LD分数 以及 分层的LD分数

回顾:

  1. GWASLab:连锁不平衡分数回归 LD score regression -LDSC
  2. GWASLab:通过Bivariate LD Score regression估计遗传相关性 genetic correlation
  3. GWASLab:基于功能分类分割遗传力 – 分层LD分数回归 Stratified LD score regression
  4. GWASLab:使用GWAS数据进行组织细胞特异性分析 – LDSC-SEG

使用LDSC计算单变量LD分数

单变量LD分数为最基本的LD分数的形式,常用在ldsc -h2选项的参考LD分数中。

输入文件: PLINK的 bed / bim / fam文件

使用软件: ldsc

这里推荐使用千人基因组相应族裔的数据(vcf格式),可通过plink -vcf转化为plink格式。

计算LD分数时需要指定窗口的大小,作者推荐 1 centiMorgan (cM) window

一般情况下plink格式中 表示cM的列都为零,可以使用plink —cm-map以及额外的遗传图谱文件来填补bim文件中的cM列

(从ldsc网站下载的plink文件中bim文件的cM列已经被填补,可以直接使用)

下面简单演示如何计算:

示例数据为千人基因组欧洲人群的第22号染色体:

wget <https://data.broadinstitute.org/alkesgroup/LDSCORE/1kg_eur.tar.bz2>
tar -jxvf 1kg_eur.tar.bz2

22.bim
22.fam
22.bed
22.hm3.daf.gz

可以通过以下的代码来计算ld分数 (-l2 选项)

python ldsc.py \
	--bfile ${plinkfile} \           #指定plink文件路径及前缀
	--l2 \                         #计算ld分数
	--ld-wind-cm 1 \                 #指定窗口的大小,单位为cM
	--out 22                         #输出文件的前缀

执行上述代码后得到如下四个文件:

22.log
22.l2.M
22.l2.M_5_50
22.l2.ldscore.gz

LD分数相关文件格式

下面简要介绍LD分数文件的相关格式,主要的输出包括以下五种文件:

.log
.l2.ldscore
.annot
.l2.M_5_50
.l2.M

.log LD分数计算后得到的日志文件,其中会包含LD分数的总结性数据,例如:

Summary of LD Scores in 22.l2.ldscore.gz
      MAF     L2
mean  0.232   18.535
std   0.145   16.104
min   0.001    0.066
25%   0.104    7.839
50%   0.224   13.484
75%   0.355   22.972
max   0.500  109.716

MAF/LD Score Correlation Matrix
     MAF    L2
MAF  1.000  0.275
L2   0.275  1.000
Analysis finished at Fri Jan 30 10:58:46 2015
Total time elapsed: 2.72s

存储LD分数最主要的文件:

.l2.ldscore 文件

包含3+N列,N≥1, 前3列分别为

  1. CHR — chromosome 染色体号
  2. SNP — SNP identifier (rs number) rsID
  3. BP — physical position (base pairs) 物理位置

第三列以后为LD分数的列,可以有多列,对应不同的注释组

  1. all additional columns — LD Scores
CHR SNP	        BP	        L2
22	rs9617528	16061016	1.271
22	rs4911642	16504399	1.805
22	rs140378	16877135	3.849
22	rs131560	16877230	3.769
22	rs7287144	16886873	7.226
22	rs5748616	16888900	7.379
22	rs5748662	16892858	7.195
22	rs5994034	16894090	2.898
22	rs4010554	16894264	6.975

.annot注释文件

上面单变量LD分数的计算中我们使用的是plink文件里所有的SNP,未使用到注释文件,但后续分层LD分数计算时,我们需要用到annot文件,来注释计算LD分数时使用的SNP

该文件每行对应一个SNP,空格分隔,有表头,需要按以下顺序排列:

  1. CHR — chromosome 染色体号
  2. BP — physical position (base pairs) 物理位置
  3. SNP — SNP identifier (rs number) rsID
  4. CM — genetic position (centimorgans) 遗传距离
  5. all additional columns — Annotations 一列或多列注释

注释列既可以为0/1二分类,也可是连续的。

例如 Example:

CHR BP SNP CM AN1 AN2
1   1  rs1  0  1  0
1   2  rs2  0  1  0
1   3  rs3  0  0  1

最后则是两个辅助文件:

.l2.M_5_50

该文件只有一行,列数对应着.l2.ldscore文件中annotation中列的个数,并与其顺序一致。

每一列的数字则表示对应的注释分类中MAF>5%的SNP的个数。

(没有附加注释文件时就是该染色体上全部MAF>5%的SNP数量)

.l2.M

l2.M_5_50相同,只是没有MAF的限制,即表示对应注释总的SNP数

例如:(对应两个annotation分类)

100	1000

计算分层LD分数

通过某种方式,选取注释组中包含的基因

例如使用单细胞测序中某个细胞特异表达的前10%的基因等,

注意:目前ldsc只支持常染色体。

为了计算某种注释特有的LD分数,需要准备一个.annot的注释文件,格式如上所述。

这个注释里的位点的顺序应该与.bim文件中的SNP完全相同。ldsc也允许—thin-annot 格式,也就是上述annot文件忽略掉CHR,BP,SNP以及CM列。

下面介绍如何构建annot文件并以此计算相应注释的LD分数:

Step 1: Creating an annot file 第一步,构建注释文件

输入文件:

—gene-set-file:注释包含的基因集

head GTEx_Cortex.GeneSet
ENSG00000105605
ENSG00000169862
ENSG00000074317
ENSG00000196361
ENSG00000221946

–gene-coord-file:基因坐标文件

head ENSG_coord.txt
GENE	CHR	START	END
ENSG00000243485	chr1	29554	31109
ENSG00000237613	chr1	34554	36081
ENSG00000186092	chr1	69091	70008
ENSG00000238009	chr1	89295	133566
ENSG00000239945	chr1	89551	91105

–windowsize:基因上下游LD计算时所包含窗口大小

–bimfile:计算LD所用PLINK文件中的bim文件

head 1000G.EUR.QC.22.bim 
22	rs587616822	-0.0048996974	16050840	G	C
22	rs62224609	-0.00094708154	16051249	C	T
22	rs587646183	0.010833912	16052463	C	T
22	rs139918843	0.012979739	16052684	C	A
22	rs587743102	0.014465964	16052837	T	C

–annot-file:输出文件

所使用的代码如下所示:使用lsmake_annot.py

python make_annot.py \
		--gene-set-file GTEx_Cortex.GeneSet \
		--gene-coord-file ENSG_coord.txt \
		--windowsize 100000 \
		--bimfile 1000G.EUR.QC.22.bim \
		--annot-file GTEx_Cortex.annot.gz

也可以直接使用bed文件计算

python make_annot.py \
		--bed-file Brain_DPC_H3K27ac.bed \
		--bimfile 1000G.EUR.QC.22.bim \
		--annot-file Brain_DPC_H3K27ac.annot.gz

结果如下:

zcat GTEx_Cortex.annot.gz | head
ANNOT
0
0
0
0
0

Step 2: Computing LD scores with an annot file. 第二步,使用注释文件计算LD分数

使用第一步计算所得的.annot文件(忽略掉了CHR,BP,SNP以及CM列,所以使用–thin-annot选项)

–print-snps hm.22.snp:只打印包含在hm.22.snp的SNP

python ldsc.py \
		--l2 \ 
		--bfile 1000G.EUR.QC.22 \
		--ld-wind-cm 1\
		--annot Brain_DPC_H3K27ac.annot.gz \
		--thin-annot \
		--out Brain_DPC_H3K27ac\
		--print-snps hm.22.snp

如果要与baseline model一起使用的话,要注意使用的SNP应该完全相同。

正式计算时只需要对所有染色体循环即可,确保输出都在同一个文件夹里,且有着相同的前缀,例如:

${prefix}.${chr}.annot.gz , ${prefix}.${chr}.l2.ldscore , and ${prefix}.${chr}.l2.M_5_50

这些文件可以用于后续分层LDSC回归分析中。

GWASLab:使用GWAS数据进行组织细胞特异性分析 – LDSC-SEG

参考

https://github.com/bulik/ldsc

使用GWAS数据估计跨族裔遗传相关 popcorn – Cross-ancestry genetic correlation

popcorn 背景介绍

回顾(在单一族裔中估计遗传相关的方法):通过Bivariate LD Score regression估计遗传相关性 genetic correlation

随着GWAS等遗传变异关联研究在不同人群中不断开展,相关研究数量不断增加,这为我们提供了宝贵的机会以研究遗传结构在不同人种中的差异,有助于解决这个在医疗研究与群体遗传学中十分重要的问题。

该文章基于Brielin等的研究,简单介绍一种估计跨族裔遗传相关的方法,popcorn,用于估计在不同群体中常见SNP的因果效应量的相关性。该方法纳入了研究中所有SNP的关联信息(区别于其他需要事先进行LD-pruning的方法),而且仅使用GWAS的概括数据,减少了计算负担,并避免了无法获取基因型数据的问题。

通常情况下,在单一群体中,两个表性之间的遗传相关的定义为对于两个表型的SNP效应量大小之间的相关系数。而在多群体的情况时,因为不同群体内等位基因频率有显著差异,定义并不统一。因为某个变异可能在某个特定群体中有较高的效应量但较低的频率,所以本方法不仅考虑了效应量大小的相关,同时也考虑遗传冲击(allelic impact)的相关性

该文中,跨族裔遗传效应相关的定义为SNP效应量的相关系数,跨族裔遗传冲击相关的定义为特定群体的等位方差正态化(allele-variance-normalized)的SNP效应量之间的相关系数

直觉上,遗传效应相关性(genetic-effect correlation)衡量的是同一个SNP造成表型改变的程度差异,而遗传冲击相关性(genetic-impact correlation)则是分别在各个群体中,相比于的稀有allele,给予common allele更高的权重。举个例子,如果只看效应相关性,计算得出的相关系数是不完全的,没有反应出两个群体间的效应的完整信息,因为即使效应量的遗传相关是1,但其在两个群体中,对其各自表型分布的冲击也可能由于allele frequency的不同而不同。

所以跨族裔遗传相关的两个核心问题就被引出:

1 同一个变异在两个群体中效应量大小的差异?

2 同一个变异的效应量在两个群体中的差异是否被不同群体的allele frequency所稀释或加重?

为了估计遗传相关,popcorn采用了贝叶斯方法,假设基因型是分别从各个群体内抽取的,并且效应量服从正态的先验分布(无穷小模型)。无穷小模型的假设会产生一个观测统计量(Z scores)的多元正态分布,其中该分布的协方差是遗传力与遗传相关的一个方程。这使得我们可以直接将Z score的膨胀纳入模型中,而不需要进行LD pruning。然后popcorn会最大化近似加权似然方程,来估计遗传力与遗传相关性。

方程推导可参考原文章,这里我们要估计的是以下两个值:

遗传效应相关:pge=Cor(β1,β2)

遗传冲击相关:pgi=Cor(δ1β1,δ2β2)

其中δ为allele的标准差

使用方法:

该文章作者的github仓库: https://github.com/brielin/Popcorn

要求:

popcorn只能使用python2,建议用conda新建一个环境
需要以下python包:
numpy 1.14.2
scipy 1.0.1
pandas 0.22.0
pysnptools 0.3.9
bottleneck 1.0.0
statsmodels 0.8.0
(to use --plot_likelihood) matplotlib 1.5.1

从github clone之后,安装:

cd Popcorn
python setup.py install

可以下载使用千人基因组项目提前计算好的LD分数 (用于EUR-EAS):

https://www.dropbox.com/sh/37n7drt7q4sjrzn/AAAa1HFeeRAE5M3YWG9Ac2Bta .

popcorn 方法分为两步,

第一步为计算LD分数(使用下载好的数据可以跳过这一步,作者提供了EUR-EAS的LD分数)

第二步为计算遗传相关性

第一步,计算LD分数(可选)

python -m Popcorn compute -v 2 --bfile1 Popcorn/test/EUR_ref --bfile2 Popcorn/test/EAS_ref --gen_effect Popcorn/test/EUR_EAS_test_ge.cscore

-v 表示 输出log的详细程度

—bfile表示 输入的plink bfile文件前缀

—gen_effec 表示计算遗传效应相关 (计算遗传冲击相关时不加这个flag)

最后是输出的路径和文件名

第一步log文件

Popcorn version 0.9.6
(C) 2015-2016 Brielin C Brown
University of California, Berkeley
GNU General Public License v3

Invoking command: python /home/brielin/CoHeritability/Popcorn/__main__.py compute -v 2 --bfile1 Popcorn/test/EUR_ref --bfile2 Popcorn/test/EAS_ref Popcorn/test/EUR_EAS_test.cscore
Beginning analysis at DATE
50000 SNPs in file 1
50000 SNPs in file 2
49855 SNPs in file 1 after MAF filter
44021 SNPs in file 2 after MAF filter
43898 SNPs common in both populations
T1
T2
TC
Analysis finished at DATE

结果的格式如下所示:

1	14930	rs75454623	A	G	0.520874751491	0.58630952381	1.4252539451	1.78489131664	1.49214493091
1	15211	rs78601809	T	G	0.731610337972	0.503968253968	1.20243142615	2.3876817626	1.29403313474
1	15820	rs200482301	T	G	0.728628230616	0.605158730159	1.56556914093	1.42934198599	1.40055911225
1	55545	rs28396308	T	C	0.813121272366	0.804563492063	2.59950434496	2.1809181648	2.06610287872

各列分别代表Chromosome, Base Position, SNP ID, Allele 1, Allele 2, Frequency of A1 in POP1, Frequency of A1 in POP2, LD score in POP1, LD score in POP2, and Cross-covariance score

第二步,计算相关

python -m Popcorn fit -v 1 --use_mle --cfile Popcorn/test/EUR_EAS_test_ge.cscore --gen_effect --sfile1 Popcorn/test/EUR_test.txt --sfile2 Popcorn/test/EAS_test.txt Popcorn/test/EAS_EUR_test_ge_result.txt

–sfile1与–sfile2为输入的两个群体的GWAS文件,注意要与LD分数中POP1和POP2的顺序相一致

输入的GWAS数据至少包含一下列(列名应与下面的介绍保持一致):

1 rsid 或者 SNP (与LD分数文件中的id保持一致)

注: 可以使用 chr 以及 pos 来替代

2 a1 和 a2

3 N 每个SNP的样本大小

4 beta和SE 或者 OR和p-value 或者 Z分数 (不要求对应a1为effect allele,但要保持所有snp一致即可)

5 AF (可选)即a2的allele frequency,用于A/T , G/C的SNP的flip

输入文件例如:

rsid	a1	a2	af	N	beta	SE
rs10458597	C	A	0.59421	50000	-0.0152782557183	0.00912151085515
rs3131972	A	C	0.83102	50000	-0.000638379272685	0.0119663325893
rs3131969	A	C	0.8654	50000	-0.000737906302463	0.0131325861467
rs3131967	A	C	0.87044	50000	-0.000189441590945	0.0133506462369

执行代码后的,第二步log文件:

Popcorn version 0.9.6
(C) 2015-2016 Brielin C Brown
University of California, Berkeley
GNU General Public License v3

Invoking command: python /home/brielin/CoHeritability/Popcorn/__main__.py fit -v 1 --cfile Popcorn/test/EUR_EAS_test_ge.cscore --gen_effect --sfile1 Popcorn/test/EUR_test.txt --sfile2 Popcorn/test/EAS_test.txt Popcorn/test/EAS_EUR_test_ge_result.txt
Beginning analysis at DATE
Analyzing a pair of traits in two populations
49999 49999 SNPs detected in input.
43897 SNPs remaining after merging with score file.
43897 SNPs remaining after filtering on AF and self-compliment SNPs.
Initial estimate:
       Val (obs)  SE
h1^2   0.172133 NaN
h2^2   0.086708 NaN
pg     0.302372 NaN
Performing jackknife estimation of the standard error using 200 blocks of size 219 .
Jackknife iter: 200
      Val (obs)        SE          Z     P (Z)
h1^2   0.172133  0.012870  13.374388      0
h2^2   0.086708  0.008555  10.135759      0
pge    0.302372  0.053021  13.157553      0
Analysis finished at DATE
Total time elapsed: T

pge 则为跨族裔遗传相关的估计值 (加–gen_effect时的结果)

pgi 为跨族裔遗传冲击的估计值(不加–gen_effect时的结果)

这里的p值的原假设是估计值等于1,p<0.05则意味着估计值显著小于1

参考:

https://pubmed.ncbi.nlm.nih.gov/27321947/

使用GWAS数据进行组织细胞特异性分析 – LDSC-SEG

目录

LDSC-SEG简介
LDSC-SEG原理
LDSC-SEG教程
参考

LDSC-SEG简介

本文主要介绍LD分数回归的一项功能扩展 ,LDSC-SEG。该方法主要通过整合基因表达数据与GWAS数据,来发现疾病相关的组织或细胞种类,其底层原理是使用分层LD分数回归(stratified LD score regression,回顾:GWASLab:基于功能分类分割遗传力 – 分层LD分数回归 Stratified LD score regression)来检测表型的遗传性是否富集于某个组织或细胞高特异表达基因的周围。

LDSC-SEG原理

LDSC-SEG的主要步骤 (如图所示)

  1. 从一个标准化后的基因表达矩阵(包含不同组织与细胞)开始
  2. 对于每个gene,首先计算其在目标组织中特异表达的t统计量
  3. 根据t统计量从大到小对所有基因进行排序,并将位于t统计量最高的前10%的基因定义为该组织特异表达的基因集。
  4. 在这些特异表达的基因的转录区域两侧加上100kb的窗口,并以此为构建组织特异的基因组注释。
  5. 最后对于以上注释后的不同组别进行分层LD回归来估计不同组织或细胞对于所研究表型遗传力的贡献

联合回归模型中同时包括

1)上述目标组织或细胞特异的基因集,

2)一个包含基因组上所有基因的基因集,

3)包含52种注释的基线模型(这些注释包括增强子区域,保守性区域等,可参考 stratified LD score regression)。

特异基因集的回归系数为正数,则表示矫正其他基因集后,该特异基因集对表型遗传性的贡献为正。

基因表达数据来源(原文中处理了五组公开的数据集,如下所示)

GTExHuman53 tissues/cell typesRNA-seq
Franke labHuman/mouse/rat152 tissues/cell typesArray
CahoyMouse3 brain cell typesArray
PsychENCODEHuman2 neuronal cell typesRNA-seq
ImmGenMouse292 immune cell typesArray

实战操作

ldsc 下载安装可参考: GWASLab:连锁不平衡分数回归 LD score regression -LDSC

教程原文:

https://github.com/bulik/ldsc/wiki/Cell-type-specific-analyses

#以上述5组数据中的Cahoy数据为例:
cts_name=Cahoy 

#下载对应的LD分数文件 (EUR)

#组织、细胞特异LD分数
wget https://data.broadinstitute.org/alkesgroup/LDSCORE/LDSC_SEG_ldscores/${cts_name}_1000Gv3_ldscores.tgz
#基线模型文件
wget https://data.broadinstitute.org/alkesgroup/LDSCORE/1000G_Phase3_baseline_ldscores.tgz
#权重文件
wget https://data.broadinstitute.org/alkesgroup/LDSCORE/weights_hm3_no_hla.tgz

#解压
tar -xvzf ${cts_name}_1000Gv3_ldscores.tgz
tar -xvzf 1000G_Phase3_baseline_ldscores.tgz
tar -xvzf weights_hm3_no_hla.tgz

#下载GWAS数据
wget https://data.broadinstitute.org/alkesgroup/UKBB/body_BMIz.sumstats.gz
#下载SNP list
wget https://data.broadinstitute.org/alkesgroup/LDSCORE/w_hm3.snplist.bz2
#解压
bunzip2 w_hm3.snplist.bz2

#第一步清洗数据 
python munge_sumstats.py \
--sumstats body_BMIz.sumstats.gz \
--merge-alleles w_hm3.snplist \
--out UKBB_BMI

#第二步 进行LDSC-SEG
cts_name=Cahoy
ldsc.py \
    --h2-cts UKBB_BMI.sumstats.gz \
    --ref-ld-chr 1000G_EUR_Phase3_baseline/baseline. \
    --out BMI_${cts_name} \
    --ref-ld-chr-cts $cts_name.ldcts \
    --w-ld-chr weights_hm3_no_hla/weights.

注意:

**--ref-ld-chr-cts 这里输入一个ldcts文件,共两列,第一列为一个便于识别的label,第二列则是逗号分隔的两个路径,第一个为被检测组织细胞的ld,第二个则为control的ld,路径格式与--ref-ld-chr的格式类似。注意这里的都是相对与解压文件夹的相对路径,使用前最好修改为绝对路径。

如下图所示:

#组织细胞标签 检验组LD分数路径,对照组LD分数路径
Astrocyte       Cahoy_1000Gv3_ldscores/Cahoy.1.,Cahoy_1000Gv3_ldscores/Cahoy.control.
Oligodendrocyte Cahoy_1000Gv3_ldscores/Cahoy.2.,Cahoy_1000Gv3_ldscores/Cahoy.control.
Neuron  Cahoy_1000Gv3_ldscores/Cahoy.3.,Cahoy_1000Gv3_ldscores/Cahoy.control.

执行完成后可以查看结果:

.cell_type_results.txt

Name    Coefficient     Coefficient_std_error   Coefficient_P_value
Neuron  7.93194788527e-09       3.02894244784e-09       0.00441303625204
Oligodendrocyte 7.32019970874e-10       3.51868270994e-09       0.417599619801
Astrocyte       -5.76220451287e-09      2.60400594455e-09       0.98654507806

第一列即为dict文件里的label,第二,三列为ldcts文件里第二例第一个路径对应的ld分数的系数与se(对应cell type specific annotation),第四列为该系数的p值。

参考

https://github.com/bulik/ldsc/wiki/Cell-type-specific-analyses

https://www.ncbi.nlm.nih.gov/pmc/ar

通过Bivariate LD Score regression估计遗传相关性 genetic correlation

背景回顾:

LD score 回归: 连锁不平衡分数回归 LD score regression

双变量/跨性状LD分数回归(Bivariate/Cross-trait LD Score regression)

简介与原理:

该方法使用GWAS的概括性数据来估计一对表型的遗传相关。该方法基于以下的事实:GWAS检验中,对一个SNP效应量的估计通常也会包含与该SNP成LD的其他SNP的效应,也就是说一个与其他SNP成高LD的SNP,通常也会有更高的卡方检验量。而当我们把卡方检验量替换为,来自两个相关表型的GWAS的Z分数的乘积时,这个事实依然存在。所以基于此,我们可以将单表型的LD分数回归扩展成如下所示的 双变量/跨性状LD分数回归

详细推导见:

Bulik-Sullivan, B., Finucane, H., Anttila, V. et al. An atlas of genetic correlations across human diseases and traits. Nat Genet 47, 1236–1241 (2015). https://doi.org/10.1038/ng.3406


下面简单介绍使用方法:

数据下载,数据清理等步骤与 单表型的LD分数回归基本相同,详见 连锁不平衡分数回归 LD score regression

我们只需要将 单表型的LD分数 中 -h2 的选项改为 -rg, 再加上要估计遗传相关的两个表型的数据清理后的文件即可:

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
--rg :计算遗传相关性 (genetic correlation),参数为上一步处理好的数据文件名
--ref-ld-chr :使用的按染色体号分类的LD score文件名,参数为LD score文件所在文件夹的路径。默认情况下ldsc会在文件名末尾添加染色体号,例如 --ref-ld-chr eur_w_ld_chr/ 的意思就是使用 eur_w_ld_chr/1.l2.ldscore ... eur_w_ld_chr/22.l2.ldscore 这些文件。如果你的染色体号在其他位置,也可以使用@来告诉LDSC,例如 --ref-ld-chr ld/chr@ 。当然也可以用 --ref-ld 来指定一个整体的LD score文件。
--w-ld-chr:指定ldsc回归权重所用的LDscore文件,理论上对于SNPj的LD分数,应当包含这个SNP与其他所有SNP的R2之和,但实际操作中,LD score回归对于计算权重的SNP的选择并不敏感,所以一般情况下我们可以使用与 --ref-ld 与相同的文件。.
--out : 输出文件的路径与前缀

运行后我们就能得到如下结果:

Genetic Covariance
------------------
Total Observed scale gencov: 0.3644 (0.0368)
Mean z1*z2: 0.1226
Intercept: 0.0037 (0.0071)

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

Genetic Correlation: 0.6561 (0.0605) 即为我们需要估计的遗传相关

除此之外,LDSC还会输出一个Summary of Genetic Correlation Results,当有多组遗传相关需要估计时,查看该表格会更加方便。

常见问题:

为什么rg会在[-1,1]的范围以外: (https://github.com/bulik/ldsc/issues/78)

ldsc的算法并没有限定[-1,1]的范围,当rg接近一时,加上估计误差有时就会导致rg在范围外

参考:

Bulik-Sullivan, B., Finucane, H., Anttila, V. et al. An atlas of genetic correlations across human diseases and traits. Nat Genet 47, 1236–1241 (2015). https://doi.org/10.1038/ng.3406

https://github.com/bulik/ldsc/wiki/Heritability-and-Genetic-Correlation

相关链接:

基于功能分类分割遗传力 – 分层LD分数回归 Stratified LD score regression 

遗传结构 Genetic architecture

1 什么是遗传结构?

人类群体遗传学研究中,“遗传结构”( ‘genetic architecture )是指 对于某一表型,影响广义表型遗传力( broad-sense phenotypic heritability )的遗传变异的总和特征,这是一个整体上的概念。

具体来说,遗传结构包括了 1. 影响某一表性的变异的数量(#),2. 对表性影响效应的大小(Beta / OR),3. 群体中这些变异的频率(MAF),以及 4. 这些变异之间或是与环境的相互作用上位效应G X E)。

2 遗传结构的种类?

遗传结构通常被描述为 单基因的(monogenic)、寡基因的(oligogenic)、多基因的(polygenic)、或是全基因的(omnigenic),表示一个,若干,大量,几乎全部遗传变异都对表型变异性有贡献。

例如,与二型糖尿病相关的变异均为效应量较小的常见变异,但与一型糖尿病相关的变异则既包括效应量较小的常见变异(common variants),还包括部分有着较大效应量的低平率变异和稀有变异(low-frequency and rare variants)。

但无论某个性状或疾病,有单基因的、寡基因的、多基因的、还是全基因的遗传结构,这些遗传结构对表型的遗传贡献还是存在本质上的变异性( variability )。这些变异性可能是表型测量上的不同或不足真实生物学异质性的一个方程。因此,对于不同疾病,遗传变异的数量或是其他相关属性都会有很大差异。(注意遗传结构的定义。)

3 了解疾病遗传构建的重要性?

遗传结构在疾病的筛查与诊断增进生物学理解,药物开发,孟德尔随机化以及基因作图方面都有重要的作用。

4 如何估计某一表型的遗传结构:

目前主要的方法包括GWASgene-based tests

5 什么因素会影响遗传结构?

表型(Phenotype.):不同的表型与遗传变异的相关联的方式不同,与环境的相互作用不同,对不同表型测量的质量的也不尽相同。这些因素都会导致观测到的遗传结构出现差异。

选择(Selection):选择是指遗传变异的频率在局部环境中向着有利的方向变化的进化过程。 遗传结构可能会受到以下影响:性状本身的性质,解释方差的变异的相对年龄与效应大小,以及一些群体的特征。

去渠道化(Decanalization.)Canalization 渠道化是指,虽然某一物种保有大量的遗传变异,并且所处环境也有各种变化,但即使在这种情况下,表型也能保持稳定不变的现象,也就是说表型对于遗传与环境的变化是稳健的。去渠道化则是这个渠道化的逆过程,在某种环境的变化,或者某些效应量极大的遗传变异影响下,上述稳健的系统被破坏,这种长期稳态被破坏时会导致疾病的产生。

参考:

1.Timpson, N. J., Greenwood, C. M. T., Soranzo, N., Lawson, D. J. & Richards, J. B. Genetic architecture: The shape of the genetic contribution to human traits and disease. Nature Reviews Genetics 19, 110–124 (2018).

使用GCTA (GREML)来估计SNP-遗传力 SNP Heritability

GWAS研究中发现的显著SNP只能解释人类群体中复杂性状很小一部分的遗传变异。那么剩余的遗传力在哪里? 很多时候这部分遗传力并没有丢失,而是由于部分snp效应太小以至于无法达到显著水平而没有被检测到。

与单SNP关联检验相对,GCTA中的GREML(genome-based restricted maximum likelihood)方法使用线性混合模型( linear mixed models , LMMs),将全部SNP的效应作为随机效应拟合。

模型如下:

y是 一个n × 1 的表型向量,n是样本大小; β 是固定效应的向量,诸如性别、年龄、以及一个或多个主成分(PC); u 是SNP效应的向量,服从如下分布:

I 是一个 n × n 的单位矩阵, ɛ 是残余效应的向量,服从以下分布:

W是标准化的基因型矩阵,其中第ijth 个元素为: 

其中 xij 是第j个个体的第i个SNP的参考allele的拷贝数, pi该 allele 的频率。如果我们定义  A = WW/N 并且定义 σ2g 为全部SNP解释的方差 , 即 σ2g=Nσ2u, (N为SNP的数量),那么方程 1 等价于:

其中g是一个n x 1的向量,表示各个个体的全部的遗传效应,服从 g∼N(0,Aσ2g) 。A 可以理解为个体两两间遗传关系组成的遗传关系矩阵 genetic relationship matrix (GRM) 。于是基于这个线性混合效应模型我们就可以通过估计的GRM,利用REML( restricted maximum likelihood )算法来无偏地估计全部SNP解释的方差 σ2g (继而估计SNP遗传力)。

基于以上定义,个体j与个体k的遗传关联可以由以下方程估计:

实际操作中我们首先需要排除掉有亲缘关系的个体,主要原因是该模型的目的是估计所有SNP解释的方差,如果纳入有亲缘关系的个体,会由于表型关联( phenotypic correlations )而造成偏差,例如由于共同的环境因素。即使没有上述的偏差,那么对它的解读也会有别于 “无关联的”个体:一个基于家系的估计值捕捉的是所有因果变异的贡献,而该方法捕捉的则是与被基因分型的SNP处于LD的因果变异的贡献


下面简单介绍如何利用GCTA软件估计snp-遗传力:

GCTA下载链接:https://cnsgenomics.com/software/gcta/#Download

GCTA语法上类似PLINK,文件格式也几乎通用,使用起来十分便捷:

第一步: 估计GRM

gcta64 --bfile test --autosome --maf 0.01 --make-grm --out test --thread-num 10

输入为plink格式的bed/bim/fam文件,

  • –autosome 只是用常染色体上的SNP
  • –maf 0.01 过滤掉maf小于0.01的snp
  • –make-grm 计算GRM
  • –out 输出文件的前缀
  • –thread-num 10 要使用的线程数

计算完成后GRM会存储在以下文件中: test.grm.bin, test.grm.N.bin and test.grm.id

我们也可以去除掉有亲缘关系的个体:

gcta64 --grm test --grm-cutoff 0.025 --make-grm --out test_rm025

第二步:利用GCTA-GREML估计SNP遗传力

gcta64 --grm test --pheno test.phen --reml --out test --thread-num 10
  • –pheno 表型文件
  • –reml 利用reml算法计算snp遗传力

计算结果储存在 test.hsq的文件中

当然我们以可以加入协变量,例如前10个主成分 PC1-10:

gcta64 --grm test --pheno test.phen --reml --qcovar test_10PCs.txt --out test --thread-num 10

对于较大的数据量,我们还可以分染色体进行计算,详见:https://cnsgenomics.com/software/gcta/#Tutorial

参考:

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3014363/

Yang, J. et al. Common SNPs explain a large proportion of the heritability for human height. Nature Genetics 42, 565–569 (2010).

https://cnsgenomics.com/software/gcta/#Overview

https://cnsgenomics.com/software/gcta/#Tutorial

基于功能分类分割遗传力 – 分层LD分数回归 Stratified LD score regression

在复杂性状的GWAS研究中,大部分遗传力由并没有达到全基因组显著性水平的SNP关联贡献。然而当前许多利用功能信息与GWAS数据来研究疾病的方法仅仅使用达到显著性水品的loci里的SNP,并且只假设每个loci里只有一个因果SNP,或是完全不考虑LD。基于这些不足,本方法的作者们期望使用所有SNP的信息,并明确的将LD纳入模型,来估计每个功能分类的SNP遗传力,以提升power。

本方法是在基础的LD分数回归上的延伸,可以参考:

在此之前分割SNP遗传力的方法借助REML( restricted maximum likelihood ),例如GCTA,但需要个体的基因型的原始数据,并且需要很大的计算资源。所以作者们开发了分层LD分数回归,只需要GWAS的 summary statistics ,以及从与目标群体相对应的参考群体中计算得到的LD信息。

 该方法的核心基于以下的一个事实,  GWAS中某个SNP的χ2  检验统计量包含了所有被该SNP标记的SNP的效应。因此对于一个多基因性状,一个有高LD分数的SNP,相比于低 LD分数的SNP ,总体上也会有较高的 χ2  检验统计量 。原因主要是这些SNP可能标记了单个有很大效应量的SNP,或是多个较弱效果的SNP。(也就是 “以一敌百” 或者 “人多力量大”)

所以我们将所有SNP分为具有不同遗传力的功能分类,那么与某个高遗传力分类存在LD的SNP就会有更高的 χ2  检验统计量 。如果与某个功能分类存在高LD的SNP有较高的 χ2  检验统计量 ,那么就定义这个功能分类有遗传力聚集。

在多基因模型下,某个SNP j 的 χ2  检验统计量 的期望值为

  • N为样本大小
  • C则是功能分类的索引
  • l(j,C)则是SNP j 对于分类C的LD分数
  • a则衡量了混淆因素的大小
  • τC则表示了每个SNP对于功能分类C遗传力的贡献。

该方程使我们可以估计 τC 的大小(也就是所谓的分割的遗传力)。定义某个分类的聚集为 该分类SNP遗传力的比例除以该分类SNP总数的比例。


该方法的作者们基于多个公开的注释数据库,构建了不针对任何细胞类型的全基线模型 ‘full baseline model’ , 包括 coding, UTR, promoter and intronic regions the histone marks monomethylation (H3K4me1) and trimethylation (H3K4me3) of histone H3 at lysine 等等。除此之外,还基于 全基线模型 ,构建了多个针对特定细胞类型的模型,包含针对细胞类型的注释等。

尽管Stratified LD score regression提供了一种便捷有效的分析 GWAS的 summary statistics 的方法,但我们也要同时注意该方法的不足之处:

  • 为了达到足够的power,需要较大的样本量,或是较大的SNP遗传力,而且性状必须是多基因的
  • 该方法需要针对研究群体的LD参考数据
  • 该方法目前不支持自定义的array
  • 该方法基于加性模型,没考虑上位或非加性的效应
  • 该方法依赖于可用的功能注释,如果没有相应注释则无法检测
  • 等等

使用方法:

下载,安装,配置环境,数据清理详见

  • 连锁不平衡分数回归 LD score regression
  • 除了以上步骤外,我们还需要下载相应的baseline模型:

    https://alkesgroup.broadinstitute.org/LDSCORE/

    从以上链接中我们需要下载以下内容(以欧洲群体为例):

    • 基线模型LD分数 :baseline.* in 1000G_Phase1_baseline_ldscores.tgz
    • 频率 :1000G.mac5eur.* in 1000G_Phase1_frq.tgz
    • 权重 :weights.* in weights_hm3_no_hla.tgz

    解压后即可使用:

    python ldsc.py 
    	--h2 BMI.sumstats.gz\
    	--ref-ld-chr baseline.\ 
    	--w-ld-chr weights.\
    	--overlap-annot\
    	--frqfile-chr 1000G.mac5eur.\
    	--out BMI_baseline
    
    • –h2 : 计算分割的遗传力,参数为之前处理好的gwas summary statistics
    • --ref-ld-chr :下载的参考LD分数文件
    • --w-ld-chr: 权重文件
    • --frqfile-chr: SNP频率文件
    • --overlap-annot: 表示基线模型中功能分类有重叠
    • –out:指定出输出文件的前缀

    参考:

    http://www.github.com/bulik/ldsc

    https://github.com/bulik/ldsc/wiki/Partitioned-Heritability

    Partitioning heritability by functional annotation using genome-wide association summary statistics https://www.nature.com/articles/ng.3404

    GWAS的线性混合模型LMM Linear Mixed Model

    关键词population structurecryptic relatedness,spurious association,LMM,GRM

    早期的GWAS研究使用的模型为固定效应的线性模型,但通常会受两方面混淆因素的影响:

    1. 群体结构/分层 (population structure/ stratification):研究群体中存在有不同祖先(ancestry)的亚群体(subgroup)
    2. 隐性关联(cryptic relatedness):研究样本之间存在未知的亲缘关系

    如果在模型中故意忽略掉这些混淆变量(confounding variable),就很可能导致结果出现假阳性(false positive)或是虚假关联(spurious association),所以这是我们在进行GWAS研究时,必须要考虑的问题。


    Population structure

    首先来看群体结构/分层

    考虑一个case-control研究,如下图所示,红色的群体在整体样本中占了case的大多数,那么一些对疾病并没有影响,但在红蓝两群体之间等位基因频率相差很大的genetic marker,就有可能会造成虚假关联(spurious association)

    在数量形状也是存在类似的情况,

    为了解决这个问题,我们通常会采用genomic control的方法,或是通过PCA来矫正。详见:


    Cryptic relatedness

    由于将主成分作为协变量纳入模型以矫正群体分层的方法,适用于样本中不存在亲缘关系的情况下使用,而样本中存在亲缘关系时,就不一定有效,但通常我们的研究中都会存在有亲缘关系的样本。

    当样本中存在隐性关联 / 错误认定的关联 Cryptic and / or misspecified relatedness时,也可能会造成上述的虚假关联。由于我们通常不能掌握所有研究对象的家谱,所以无法完全去相关个体,例如下面的情况,样本中看似不相关的个体间实际上存在隐性关联:


    Linear Mixed Model

    为了解决以上问题,相比早期的线性模型,我们就采用一种更为灵活的模型,线性混合模型来进行检验。

    Y是 n x 1的向量,表型数据

    固定效应 fixed effects:

    • W为 n x (w + 1)的矩阵,包含了截距,以及协变量
    • β 则是 ( w + 1 ) × 1 的向量,表示协变量的效应量
    • Gs为 n × 1 的向量,某个位点的基因型,每一项的值通常为 0,1,2(等位基因allele的拷贝数)
    • γ则是标量,目标位点基因型的效应量

    随机效应 random effects

    • g为长度为n的随机向量,表示多基因效应
    • δ2g为加性遗传方差
    • Ψ是成对遗传相关的矩阵
    • e是长度为n的随机向量,服从以下的分布,其中δ2e为非遗传效应造成的方差,我们认为这项在个体间是相互独立的
    这样我们就构建了包含群体结构和隐形相关的模型。

    通常成对遗传相关的矩阵Ψ是未知的,在计算中我们使用empirical GRM(genetic relatedness matrix,通过纳入研究的SNP计算得出)


    目前已有多种GWAS检验方法基于LMM模型,包括:

    • EMMAX
    • GEMMA
    • GMMAT
    • Bolt-LMM
    • fastGWAS
    • SAIGE
    • 等等

    参考:

    Lecture 6: GWAS in Samples with Structure ,Summer Institute in Statistical Genetics 2015

    连锁不平衡 linkage disequilibrium LD

    连锁不平衡(linkage disequilibrium)是进化生物学与人类遗传学中一个十分重要的概念,因为遗传过程中很多因素能够影响它,而它又会作用于很多因素,包括选择,重组频率,突变率,遗传漂变,交配模式,群体结构等等。反过来看,连锁不平衡就是反应群体遗传过程的一个强有力的信号。

    连锁不平衡 是指 不同基因座(loci)等位基因(allele)之间非随机(nonrandom)的关联

    首先考虑简单的两基因座情况,设有A, B两个基因座,每个基因做各有两个等位基因,分别用1,2表示。假设每个单倍体型的频率如下所示:

    由上 单倍体型的频率 ,我们也可以简单计算得到各个等位基因的频率:

    如果这两个基因座互相独立不相关(也就是连锁平衡 linkage equilibrium 的状态),那么各个单倍型的频率就可以直接算出,为p1q1 ,p1,q2 , p2q1, p2q2

    而实际情况中单倍型的频率对于不相关情况下的理论值会产生偏离(deviation),这个偏离原因即为连锁不平衡( linkage disequilibrium ),偏离的程度通常记为 D (连锁不平衡系数,coefficient of linkage disequilibrium

    D=x_{11}-p_{1}q_{1}

    下图表示了各单倍型频率,各等位基因频率与D之间的关系。

    但要注意的是,D值并不是一个用来衡量LD的很好的指标,因为D值会受等位基因频率影响,这使得我们无法比较不同频率的等位基因对之间连锁不平衡的大小。

    Lewontin提出通过标准化D值来解决该问题,即用D值除以理论上D可能的最大绝对值:

    {\displaystyle D'={\frac {D}{D_{\max }}}}
    其中

    但更多的时候我们使用相关系数(correlation coefficient)r2来衡量LD:

    在众多LD的数据库中我们查询到的数值也是基于r2的。

    参考:

    https://en.wikipedia.org/wiki/Linkage_disequilibrium

    Montgomery Slatkin. Linkage disequilibrium — understanding the evolutionary past and mapping the medical future.

    GWAS入门文章与书籍推荐

    对于统计遗传学(全基因组关联分析)的初学者来说,最难得莫过于入门,

    好在GWAS问世以来已经过了十多年,该研究领域也形成了一定规模,基础内容逐渐充实,2020年出版的 An Introduction to Statistical Genetic Data Analysis 网罗该领域内研究背景,基础知识,常用工具介绍,代码实操等。

    本书第一部分主要介绍相关基础内容,涵盖群体遗传学的基本概念,统计学基础,人类进化(Human evolution),GWAS,风险评分(PRS)等等基础内容。

    第二部分介绍遗传数据的处理与使用,质控(QC),人群分层,PCA,等等GWAS的操作步骤。

    第三部分则是介绍GWAS下游的分析方法,包括基因环境相互作用,PRS,数据可视化,MTAG,孟德尔随机化等。

    图1: An Introduction to Statistical Genetic Data Analysis 封面

    原书链接: https://mitpress.mit.edu/books/introduction-statistical-genetic-data-analysis

    这本书难度适中,适合初学者入门,对于理解基础概念会有很大帮助,但深度不够,个人建议可以快速阅读这本书,掌握群体遗传学研究的大致框架后,针对感兴趣的领域查找文献原文来了解细节。也可以阅读Nature Reviews Genetics上关于GWAS的综述文章,对本领域入门很有帮助。


    对于希望尝试更硬核内容的同学们,我还强烈推荐这本书,Handbook of statistical genomics,本书更偏向数理原理,有大量的公式推导等,非常刺激,对于无数理基础的初学者则不太推荐此书。

    图2: Handbook of statistical genomics 封面

    除了书籍以外,其他适合入门的资料还包括:

    密西根大学的开设的biostats 666这门课的slides(听着这课的名字是不是就很6)(请记住这个这位大佬的网站Genome Analysis Wiki 以后你会用到无数次的),

    Biostatistics 666: Main Page

    该课程也涵盖了统计遗传学中重要的概念知识点,包括了基础概念以及部分统计原理,适合快速浏览:

    图3 :Biostatistics 666: slides 截图

    华盛顿大学遗传统计学夏季学校 SISG 讲义:

    除了理论基础,还有代码实操等。可以跟着代码一起练习,非常适合初学者。

    GWAS and Sequencing Data

    图4 : SISG slides 截图

    另外就是大阪大学遗传统计学夏季学校的讲义分享, 我最喜欢的零基础入门其实是这个,图文并茂, 概念细节比较到位(但是是日语的,有机会了给大家翻译):

    https://www.slideshare.net/YukinoriOkada/presentations

    最后,当然,还有GWASLab的主页:GWASLab

    博客连接:GWAS文章索引 – Article index

    以及GWAS相关文章汇总:GWASLab:全基因组关联分析GWAS文章汇总(持续更新)

    我会持续更新统计遗传学,生物信息学等相关中文内容,感谢大家的关注!