查看: 68105|回复: 32

[动植物重测序] 重测序SCI标配——LD衰减图解读

  [复制链接]

管理员

Rank: 15Rank: 15Rank: 15Rank: 15Rank: 15

主题
422
注册时间
2015.11.23
在线时间
952 小时

宣传达人优秀版主


发表于 2016.5.31 10:41:09 | 显示全部楼层 |阅读模式
即使经历过华丽,你也要接受平凡。看完了华丽的STRUCTURE堆叠图,也还是要读懂一些朴实的图形嘛。例如,今天要介绍的LD衰减图,非常朴实但在重测序类的文章中会经常出现(群体遗传、GWAS等的文章里面)。
LD衰减图
颜值:☆☆
实用性:☆☆☆☆


图1 野生大豆和栽培大豆的LD衰减图

1.基础概念的解读

要理解LD衰减图,我们就必须先理解连锁不平衡(Linkage disequilibrium,LD)的概念。连锁不平衡是由两个名词构成,连锁+不平衡。前者,很容易让我们产生概念混淆;后者,让这个概念变得愈加晦涩。因此从一个类似的概念入手,大家可能更容易理解LD的概念,那就是基因的共表达。

基因的共表达,通常指的是两个基因的表达量呈现相关性。比较常见的例子就是:转录组因子和靶基因间的关系。因为转录因子对它的靶基因有正调控作用,所以转录因子的表达量提高会导致靶基因的表达量也上调,两者往往存在正相关关系。这个正相关关系,可以使用相关系数r^2来度量,这个数值在-1~1之间。总而言之,相关性可以理解为两个元素共同变化,步调一致。

类似的,连锁不平衡(LD)就是度量两个分子标记的基因型变化是否步调一致,存在相关性的指标。如果两个SNP标记位置相邻,那么在群体中也会呈现基因型步调一致的情况。比如有两个基因座,分别对应A/a和B/b两种等位基因。如果两个基因座是相关的,我们将会看到某些基因型往往共同遗传,即某些单倍型的频率会高于期望值。

例如在下图2中,在群体中(A,a,B,b)各个基因型的频率已知的情况下,各种单倍型的期望频率(AB、Ab、aB、ab)都是可以计算出来。例如,AB的频率=(A的频率)X(B的频率)。但我们实际统计群体中各个单倍型的频率的时候,会观察到某些单倍型的频率会大于期望值,例如下图中的单倍型AB的理论频率是0.12,但观察到的实际频率是0.29。那么说明,基因型A更倾向于基因型B共同遗传。

这一般往往是由于在祖先的基因组中,A和B就是位于同一条染色体上,在传代过程中,这种共同遗传的关系被保留了下来。位点间的这种相关性,在杂交家系中一般被称为连锁(孟德尔老师豌豆实验中的发现),在自然群体中则一般被称为连锁不平衡。所以连锁不平衡中的“不平衡”,我认为可以理解为单倍型的频率分布偏离期望值,偏离了平衡。


图2 期望单倍型频率和观测到的单倍型频率

这种不同基因座间的相关性,用一个数值来衡量就是D值(图2中有计算公式)。类似相关系数是标准化后的协方差,LD系数(r^2)则是标准化后的D值(图2中有计算公式),这个数值在0~1波动。r^2=0就是两个位点完全不相关,群体中单倍型分布是随机的(观测值=期望值)。r^2=1就是两个位点完全相关,某些基因型(A)只与特定的基因型(B)共同出现。

一般而言,两个位点在基因组上离得越近,相关性就越强,LD系数就越大。反之,LD系数越小。也就是说,随着位点间的距离不断增加,LD系数通常情况下会慢慢下降。这个规律,通常就会使用LD衰减图来呈现。

2.图形解读和应用

图形的解读


图3 黄瓜群体遗传分析文章中各个亚群的LD衰减图[2]

LD衰减图就是利用曲线图来呈现基因组上分子标记间的平均LD系数随着标记间距离增加而降低的过程。大概的计算原理就是先统计基因组上两两标记间的LD系数大小,再按照标记间的距离对LD系数进行分类,最终可以计算出一定距离的分子标记间的平均LD系数大小。如图3是黄瓜重测序文章中统计各个亚群体的LD衰减速度的图形。横坐标是物理距离(kb),纵坐标是LD系数(r^2)。

从图中我们可以看出,西双版纳这个亚群体(紫色线)在基因组上50kb距离的平均LD系数大小约为0.4,但到了100kb的距离,对应的平均LD系数大小则降低到了不到0.3。而且,我们从图中也可以观察到LD系数的衰减速度在不同的亚群体快慢不同,衰减速度是 india > East Asian& Eurasian > Xishuanbanna。那说明india群体的LD衰减距离最小,可能是india这个群体遗传多样性最高导致。这句话该如何理解呢?

LD衰减距离

实际上,LD衰减的速度在不同物种间或同物种的不同亚群体间,往往差异非常巨大。所以,通常会使用1个标准——“LD衰减距离”来描述LD衰减速度的快慢。

LD衰减距离通常指的是:当平均LD系数衰减到一定大小的时候,对应的物理距离。

“一定大小”是这个定义的关键点,但没有特别统一的标准,在不同文章中标准不同。常见的标准包括:a)LD系数降低到最大值的一半;b)LD系数降低到0.5以下;c)LD系数降低到0.1以下;d)LD系数降低到基线水平(但注意,不同材料的基线值是不同的。比如图3黄瓜群体的基线大概是0.1)。

所以,下次你在文章中看到“LDdecay distance is XXkb”的时候,别忘了看看作者使用的标准是什么。

LD衰减会受什么的影响?


如图3所示, LD系数衰退速度会受到不同因素的影响而有所不同。常见的因素包括:

1)物种类型LD存在的本质是两个位点的连锁遗传导致的相关性。但这种相关性理论上会随着世代的增加、重组次数的增加而不断下降。所以,那些繁殖力强、时代间隔短的物种(例如,昆虫),其LD衰减的速度是非常快的。例如在家蚕和野蚕群体中,LD系数下降到最大值的1/2仅仅需要46bp和7bp的距离[3]。

2)群体类型相同物种的不同群体,由于其遗传背景不同,LD衰减速度也存在很大的差异。驯化选择,会导致群体遗传多样性下降,位点间的相关性(连锁程度)加强。所以,通常驯化程度越高,选择强度越大的群体,LD衰减速度是最慢的。例如,栽培稻比野生稻通常更大的LD衰减距离。类似的,自然选择、遗传漂变导致的群体遗传多样性下降,也会减慢LD衰减的速度。

3)在染色体的位置染色体不同区域的LD衰减距离而是不同的。通常着丝粒区更难重组,所以LD衰减更慢。而基因组上那些受选择的区域相比普通的区域,LD衰减速度也是更慢的[3]。

LD衰减距离的应用

LD衰减速度,在群体遗传分析中本身是对群体特性的评估,与群体类型的特性(自然群体还是驯化群体,选择强度大小)是相关的。但在其他研究中还有更多的应用价值。

基于分子标记(例如,SNP芯片,GBS测序)的GWAS分析,其实并没有检测到功能突变,本质就是利用标记和功能突变的相关性(LD关系),来检测与性状相关的功能突变的位置。一般而言,LD系数大于0.8就是强相关。如果LD系数小于0.1,则可以认为没有相关性。如果LD衰减到0.1这么大的区间内都没有标记覆盖的话,即使这个区间有一个效应很强的功能突变,也是检测不到关联信号的。所以,通常可以通过比较LD衰减(到0.1)距离和标记间的平均距离,来判断标记是否对全基因组有足够的覆盖度。

而如果GWAS检测到显著关联的区间后,则可以通过进一步绘制局部的LD单体型块图,来进一步判断显著相关的SNP和目标基因间是否存在强LD关系。这个图形我们下一篇文章会介绍。

再提一个应用的例子。在之前的文章中我们提到过,在进行STRUCTURE分析的时候理论上必须输入不相关的位点。那么,就可以通过预估LD衰减到0.1的距离,来判断标记间的距离必须大于多少才能保证标记间不具相关性(LD<0.1)。

3.绘制方法


LD衰减图的绘制,实际上有两个步骤:

1)计算marker间两两的LD系数大小

这个可以使用haploview软件完成。计算的时候,只要设定一个关键的参数:区间大小。例如设定为5Mb,那么软件就会计算基因组上所有距离<5Mb的两两位点间的LD系数。实际上这个参数设定更大也没有意义,一般情况下位点间的相关性不会延伸到大于5Mb这么远的距离。

2)绘图

将LD系数按照对应的两个marker间的距离进行分类,例如:距离按照区间大小0~5k,5k~10k,10k~15k…..分别分类。如果重测序的数据,SNP标记密度较大,这个分类区间可以设置小一些;如果是简化基因组数据,SNP标记较为稀疏,则分类区间可以适当加大。然后计算每种距离分类的LD系数的均值。最后在利用均值绘制曲线图就ok了。这一步的绘图,使用excel或R语言都可以轻松完成。


参考文献:
[1]Lam H M, Xu X, Liu X, et al. Resequencing of31 wild and cultivated soybean genomes identifies patterns of genetic diversityand selection[J]. Nature genetics, 2010, 42(12): 1053-1059.
[2] Qi J, Liu X, Shen D, et al. A genomicvariation map provides insights into the genetic basis of cucumberdomestication and diversity[J]. Nature genetics, 2013, 45(12): 1510-1515.
[3]Xia Q, Guo Y, Zhang Z, et al. Complete resequencingof 40 genomes reveals domestication events and genes in silkworm (Bombyx)[J].Science, 2009, 326(5951): 433-436.

PS:论坛里有大牛分享了LDdecay分析流程,大家可以去试试喔!



往期重测序图形专题文章(戳题目查看~):
SCI图形中的白富美——环形图
群体结构图形三剑客——PCA图
技术贴!系统发生树解读
我是structure酱,一个集美貌与智慧于一身的图形

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
有问题请发贴提问
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
3
注册时间
2016.4.22
在线时间
35 小时

发表于 2016.5.31 14:57:32 | 显示全部楼层
回复

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2016.5.30
在线时间
4 小时

发表于 2016.6.2 09:54:21 | 显示全部楼层
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
13
注册时间
2016.4.28
在线时间
134 小时

发表于 2016.6.2 12:50:47 | 显示全部楼层
学习了 感谢
回复 支持 反对

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2016.6.4
在线时间
22 小时

发表于 2016.6.5 20:18:31 | 显示全部楼层
有用!感谢!
回复

使用道具 举报

帝王蝶

Rank: 4

主题
4
注册时间
2016.5.11
在线时间
70 小时

发表于 2016.7.15 11:16:31 | 显示全部楼层
学习一下
回复

使用道具 举报

草履虫

Rank: 2

主题
0
注册时间
2016.9.20
在线时间
1 小时

发表于 2016.9.20 11:36:37 | 显示全部楼层
受益良多,往往数据输出了却不知如何解读。学习了
回复 支持 反对

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2016.3.30
在线时间
64 小时

发表于 2016.10.26 17:06:24 | 显示全部楼层
上文说到“衰减速度是 india > East Asian& Eurasian > Xishuanbanna。那说明india群体的LD衰减距离最大,可能是india这个群体遗传多样性最高导致”,根据图3所示,India的LD衰减距离应该是最小的呀?我有点晕了@小瑶@小瑶
回复 支持 1 反对 0

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
1
注册时间
2017.2.18
在线时间
19 小时

发表于 2017.3.22 15:43:43 | 显示全部楼层
午言 发表于 2016.10.26 17:06
上文说到“衰减速度是 india > East Asian& Eurasian > Xishuanbanna。那说明india群体的LD衰减距离最大, ...

同问?是不是搞错了

点评

是笔误。 india群体的LD衰减距离最小。  发表于 2017.3.24 09:55
回复 支持 反对

使用道具 举报

管理员

Rank: 15Rank: 15Rank: 15Rank: 15Rank: 15

主题
422
注册时间
2015.11.23
在线时间
952 小时

宣传达人优秀版主


 楼主| 发表于 2017.3.24 09:55:22 | 显示全部楼层
午言 发表于 2016.10.26 17:06
上文说到“衰减速度是 india > East Asian& Eurasian > Xishuanbanna。那说明india群体的LD衰减距离最大, ...

是笔误。 india群体的LD衰减距离最小。
话说,你没@到我……
有问题请发贴提问
回复 支持 反对

使用道具 举报

草履虫

Rank: 2

主题
0
注册时间
2017.2.15
在线时间
10 小时

发表于 2017.4.16 17:09:13 | 显示全部楼层
各位好啊,向大家请教一个问题:RADseq后得到的数据能够用来做LD衰减的分析吗,会不会因为RAD本身就是采用酶切的原因方法,call出来的snp位置差异较大,而最终导致r2计算不准确呢
回复 支持 反对

使用道具 举报

帝王蝶

Rank: 4

主题
0
注册时间
2017.5.15
在线时间
74 小时

发表于 2017.5.16 00:26:32 来自手机 | 显示全部楼层
好详细,好文
回复 支持 反对

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2017.6.6
在线时间
183 小时

灌水之王


发表于 2017.6.12 17:35:23 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

帝王蝶

Rank: 4

主题
13
注册时间
2016.1.14
在线时间
99 小时

发表于 2018.4.22 21:13:54 | 显示全部楼层
图3中East Asian和 Eurasian曲线有交叉,怎么解读呢?谢谢!
回复 支持 反对

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
5
注册时间
2017.3.7
在线时间
18 小时

发表于 2018.5.31 12:09:33 | 显示全部楼层
文中说到:基于分子标记(例如,SNP芯片,GBS测序)的GWAS分析,其实并没有检测到功能突变,本质就是利用标记和功能突变的相关性(LD关系),来检测与性状相关的功能突变的位置。一般而言,LD系数大于0.8就是强相关。如果LD系数小于0.1,则可以认为没有相关性。如果LD衰减到0.1这么大的区间内都没有标记覆盖的话,即使这个区间有一个效应很强的功能突变,也是检测不到关联信号的。所以,通常可以通过比较LD衰减(到0.1)距离和标记间的平均距离,来判断标记是否对全基因组有足够的覆盖度。

而如果GWAS检测到显著关联的区间后,则可以通过进一步绘制局部的LD单体型块图,来进一步判断显著相关的SNP和目标基因间是否存在强LD关系。

关于这一段能不能再帮我说仔细一点呢?
我是做了LD衰减,r2=0.1时,10条染色体平均衰减距离300kb,然后做了GWAS分析得到一些显著性SNP标记,请问我该怎样将这两者联系起来再去找候选基因啊??? 谢谢了!!!
回复 支持 反对

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
5
注册时间
2017.3.7
在线时间
18 小时

发表于 2019.2.28 15:06:11 | 显示全部楼层
老师,我用tassel做出来的LD文件,用Excel进行dist-bp区间求了R2的平均值,但是在做出来的图衰减不明显,请问这个tassel出来的数据是可以直接用的还是需要筛选?PS:计算LD 能不能单条染色体做出来再合并画图?
回复 支持 反对

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2017.10.10
在线时间
3 小时

发表于 2019.6.4 10:05:48 | 显示全部楼层
很详细,学习中,感谢分享
回复 支持 反对

使用道具 举报

草履虫

Rank: 2

主题
0
注册时间
2019.2.26
在线时间
3 小时

发表于 2019.10.4 11:41:03 | 显示全部楼层
明白了 一点
回复 支持 反对

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2019.9.29
在线时间
10 小时

发表于 2019.10.13 12:57:07 | 显示全部楼层
留名,收藏
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
0
注册时间
2019.3.19
在线时间
11 小时

发表于 2019.11.6 20:51:28 | 显示全部楼层
感谢楼主分享
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

快速回复 返回顶部 返回列表