标签： Linkage disequilibrium(LD)

使用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

为什么在PCA或估计GRM时要去除长LD区域 Remove long-LD region

本文内容

1.背景介绍
2.为什么在PCA等时要去除LongLD区域？
3.长LD区域起始位置的列表
4.使用PLINK去除长LD区域里的SNP：

1.背景介绍：

LD ：连锁不平衡 linkage disequilibrium LD

PCA ：群体分层与主成分分析 Population structure & PCA

2.为什么在QC时要去除LongLD区域？

人类基因组中存在若干长LD的区域，这些区域多位于染色体的着丝粒附近，还有一些位于HLA等区域。如下图所示：

这些区域跨度很长（长度超过2Mb），单次LD-pruning无法完全去除互相成LD的SNP，在进行诸如PCA，或是计算GRM，进行基于LMM模型的GWAS分析时，我们应当去除掉这些区域。

长LD区域的形成并不一定是因为选择，其他原因诸如倒位多态性（inversion polymorphism）也可能造成长LD区域的存在。在进行研究时，应当谨慎区分这些区域形成的原因。如果在计算模型中没有对这些长LD区域进行处理，就可能影响群体遗传结构中对于本地群体的估计，造成系统性的偏倚。

3.长LD区域起始位置的列表（hg38,hg19与hg18参考基因组版本）

hg38版本

Chr	Start	Stop
chr1	47761740	51761740
chr1	125169943	125170022
chr1	144106678	144106709
chr1	181955019	181955047
chr2	85919365	100517106
chr2	87416141	87416186
chr2	87417804	87417863
chr2	87418924	87418981
chr2	89917298	89917322
chr2	135275091	135275210
chr2	182427027	189427029
chr2	207609786	207609808
chr3	47483505	49987563
chr3	83368158	86868160
chr5	44464140	51168409
chr5	129636407	132636409
chr6	25391792	33424245
chr6	26726947	26726981
chr6	57788603	58453888
chr6	61109122	61357029
chr6	61424410	61424451
chr6	139637169	142137170
chr7	54964812	66897578
chr7	62182500	62277073
chr8	8105067	12105082
chr8	43025699	48924888
chr8	47303500	47317337
chr8	110918594	113918595
chr9	40365644	40365693
chr9	64198500	64200392
chr9	88958735	88959017
chr10	36671065	43184546
chr10	41693521	41885273
chr11	88127183	91127184
chr12	32955798	41319931
chr12	34639034	34639084
chr14	87391719	87391996
chr14	94658026	94658080
chr17	43159541	43159574
chr20	4031884	4032441
chr20	33948532	36438183
chr22	30060084	30060162
chr22	42980497	42980522

hg19版本

Chr	Start	Stop ID
1 48000000 52000000 1
2 86000000 100500000 2
2 134500000 138000000 3
2 183000000 190000000 4
3 47500000 50000000 5
3 83500000 87000000 6
3 89000000 97500000 7
5 44500000 50500000 8
5 98000000 100500000 9
5 129000000 132000000 10
5 135500000 138500000 11
6 25000000 35000000 12
6 57000000 64000000 13
6 140000000 142500000 14
7 55000000 66000000 15
8 7000000 13000000 16
8 43000000 50000000 17
8 112000000 115000000 18
10 37000000 43000000 19
11 46000000 57000000 20
11 87500000 90500000 21
12 33000000 40000000 22
12 109500000 112000000 23
20 32000000 34500000 24

hg18版本

Chr	Start	Stop	ID
1	48060567	52060567	hild1
2	85941853	100407914	hild
2	134382738	137882738	hild3
2	182882739	189882739	hild4
3	47500000	50000000	hild5
3	83500000	87000000	hild6
3	89000000	97500000	hild7
5	44500000	50500000	hild8
5	98000000	100500000	hild9
5	129000000	132000000	hild10
5	135500000	138500000	hild11
6	25500000	33500000	hild12
6	57000000	64000000	hild13
6	140000000	142500000	hild14
7	55193285	66193285	hild15
8	8000000	12000000	hild16
8	43000000	50000000	hild17
8	112000000	115000000	hild18
10	37000000	43000000	hild19
11	46000000	57000000	hild20
11	87500000	90500000	hild21
12	33000000	40000000	hild22
12	109521663	112021663	hild23
20	32000000	34500000	hild24
X	14150264	16650264	hild25
X	25650264	28650264	hild26
X	33150264	35650264	hild27
X	55133704	60500000	hild28
X	65133704	67633704	hild29
X	71633704	77580511	hild30
X	80080511	86080511	hild31
X	100580511	103080511	hild32
X	125602146	128102146	hild33
X	129102146	131602146	hild34

4.使用PLINK去除长LD区域里的SNP：

我们可以使用PLINK来去除长LD区域里的SNP，主要分为两步：

1.将上一节中的列表拷贝进high-ld.txt文件中（使用时记得去掉header），使用--make-set选项提取区域中的SNP

2.在分析时利用--exclude去除掉上一步所生成列表中的SNP

plink --file mydata --make-set high-ld.txt --write-set --out hild
plink --file mydata --exclude hild.set --recode --out mydatatrimmed

参考：

https://genome.sph.umich.edu/wiki/Regions_of_high_linkage_disequilibrium_(LD)

Price et al. (2008) Long-Range LD Can Confound Genome Scans in Admixed Populations. Am. J. Hum. Genet. 86, 127-147

更新：

20220905 修改表述错误，更新PCA链接，并增加hg38版本

通过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 

连锁不平衡 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.