LD衰减图

Python07

LD衰减图,第1张

LD衰减距离指的是,当平均LD系数衰减到一定大小(最大值的一半/0.5以下)的时候,对应的物理距离。通常用LD衰减距离来描述LD衰减速度。LD衰减速度越快,即衰减距离越小,说明该群体遗传多样性越高;LD衰减速度越慢,通常驯化程度越高,选择强度越大,导致遗传多样性下降。

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

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

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

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

发现不止1000kb

用poplddecay画图时会返回有以下内容的文件,前两列是画图所需的的数据

其中列名#Dist、Mean_r^2在R语言中无法识别,所以首先要改列名

基础知识参考文章:

https://links.jianshu.com/go?to=http%3A%2F%2Fwww.omicshare.com%2Fforum%2Fthread-878-1-1.html

[连锁不平衡粗俗的说就是:这几个基因耍流氓,喜欢抱团遗传,不再随机。而连锁不平衡衰减是指在基因组上,随着物理距离的增大,两个连锁的的等位基因的连锁程度不断减小。]

LD衰减 图,在重测序类的文章中会经常出现群体遗传、GWAS等的文章里面 。

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

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

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

图形的解读

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

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

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

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

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

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

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

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

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语言都可以轻松完成。

分子层面对生物的研究,在个体水平上主要是看单个基因的变化以及全转录本的变化(RNA-seq);在对个体的研究的基础上,开始了群体水平的研究。如果说常规的遗传学主要的研究对象是个体或者个体家系的话,那么群体遗传学则是主要研究由不同个体组成的群体的遗传规律。

在测序技术大力发展之前,对群体主要是依靠表型进行研究,如加拉巴哥群岛的13中鸟雀有着不同的喙,达尔文认为这是自然选择造成的后果 。达尔文的进化论对应的观点可以简单概括为“物竞天择,适者生存”,这也是最为大众所接受的一种进化学说。直到1968年,日本遗传学家提出了中性进化理论[2],也叫中性演化理论。中性理论的提出很大程度上是基于分子生物化学的发展。可以这样理解中性理论:一群人抽奖,在没有内幕的情况下,每个人抽到一等奖的概率是相等的,这个可能性和参与抽奖的人的身高、年龄、爱好等因素都没有关系。中性理论常作为群体遗传研究中的假设理论(CK)来计算其他各种统计指标。

群体遗传学,研究的单位是群体,比如粳稻、籼稻、野生稻,就能够构成不同的群体;我们国内的各省份的水稻也可以作为一个个群体。 群体遗传学大概可以分为群体内的研究和群体间的研究。比如研究云南元阳的水稻的遗传多样性;如果研究是的云南元阳的水稻和东北的水稻,那就可以算成是群体间的研究。群体间和群体内的研究是相互的。

测序价格的急剧下降[3]使得大规模的群体测序得以实现。

常见的变异类型有SNP、IdDel、SV、CNV等。重测序中最关注的是SNP,其次是InDel。其他的几种结构变异的研究不是太多。

有参考基因组的物种的全基因组测序叫做重测序,没有参考基因组的物种的全基因组测序则需要从头组装。随着测序价格的降低,越来越多物种的参考基因组都已经测序组装完成。 plant genomes [4]网站实时显示全基因组测序已经完成的植物,其中2012年以后爆发式增长。在群体遗传学研究中更多的是有参考基因组的物种,尤其是模式物种,植物中常见的是拟南芥、水稻和玉米。

主要的分析流程见下图。现在的测序公司基本上都会帮客户完成整个的分析流程,因为主要耗费的资源是计算资源。我认为在整个分析的流程中最重要的是Linux目录的构建,混乱的目录会导致后续的分析频频出问题,重测序分析会生成很多的中间文件,良好的目录管理会使得项目分析流程井然有序。

该部分涉及到的软件的安装和基础的Linux基础知识就不详细说明了。

正选择似乎可以更好地用自然选择来解释。就是一个基因or位点能够使个体有着更强的生存力或者是育性,这样就会使得这个个体的后代更多,如此一来,这个基因or位点在群体中就越来越多。

正选择能够使有利的突变基因or位点在群体中得到传播,但是与此同时却降低了群体的多态性水平。也就是说原先该位点周围的核苷酸组成是多样性的,在经过正选择之后,这个位点周围核苷酸的多样性就渐渐的趋于同质化了。这就好比一块田,里面本来有水稻和稗草及其他杂草,由于稗草的适应性增强,稗草在逐渐增多,水稻慢慢变少,最后甚至是只剩下了稗草。

我们将这种选择之后多态性降低的情况叫做选择扫荡(Selective Sweep)。检测选择扫荡的软件有SweeD[7]。选择扫荡有可能是人工选择的结果,如2014年 Nature Genetics关于非洲栽培稻的文章就使用了SweeD来检测非洲栽培稻基因组上受人工选择的区域[8]。

负选择和正选择刚好是相反的。简单理解成群体中的某个个体出现了一个致命的突变,从而自己或者是后代从群体中被淘汰。这也导致群体中该位点的多态性的降低。就好比我有10株水稻,其中一株在成长过程中突然不见了,那么对我的这个小的水稻群体来说,这个消失的水稻的独有的位点在群体中就不见了,整体的多态性就降低了。

平衡选择指多个等位基因在一个群体的基因库中以高于遗传漂变预期的频率被保留,如杂合子优势。

平衡选择检测的算法有BetaScan2[10],这是个Python脚本,输入文件只需要过滤好的SNP数据即可。

计算公式为:

其中 是有效群体大小, 是每个位点的突变速率。 但是群体大小往往是无法精确知道的,需要对其进行估计。

分离位点数 是 的估计值,表示相关基因在多序列比对中表现出多态性的位置。计算公式为:

其中 为分离位点数量,比如SNP数量。

为个体数量的倒数和:

指的是核苷酸多样性,值越大说明核苷酸多样性越高。通常用于衡量群体内的核苷酸多样性,也可以用来推演进化关系[11]。计算公式为:

可以理解成现在群体内两两求 ,再计算群体的均值。计算的软件最常见的是 vcftools ,也有对应的R包 PopGenome 。通常是选定有一定的基因组区域,设定好窗口大小,然后滑动窗口进行计算。

3KRGP文章就计算了水稻不同亚群间4号染色体部分区域上的 值[12],能够看出控制水稻籽粒落粒性的基因 Sh4 位置多态性在所有的亚群中都降低了。说明这个基因在所有的亚群中都是受到选择的,这可能是人工选择的结果。

Tajima's D是日本学者Tajima Fumio 1989年提出的一种统计检验方法,用于检验DNA序列在演化过程中是否遵循中性演化模型[14]。计算公式为:

D值大小有如下三种生物学意义:

叫固定分化指数,用于估计亚群间平均多态性大小与整个种群平均多态性大小的差异,反映的是群体结构的变化。其简单估计的计算公式为:

的取值范围是[0,1]。当 时,表明亚群间有着明显的种群分化。

在中性进化条件下, 的大小主要取决于遗传漂变和迁移等因素的影响。假设种群中的某个等位基因因为对特定的生境的适应度较高而经历适应性选择,那该基因的频率在种群中会升高,种群的分化水平增大,使得种群有着较高的 值。

值可以和GWAS的结果一起进行分析, 超过一定阈值的区域往往和GWAS筛选到的位点是一致的,如2018年棉花重测序的文章[15]:

ROD可以基于野生群体和驯化群体间核苷酸多态性参数 的差异识别选择型号,也可以测量驯化群体和野生型群体相比损失的多态性。计算公式为:

和 一样,ROD也可以和GWAS结合起来:

群体结构分析可以简单理解成采样测序的这些个体可以分成几个小组,以及给每个个体之间的远近关系是怎么样的。群体结构分析三剑客, 分别是 进化树 、 PCA 和 群体结构图 。

进化树就是将个体按照远近关系分别连接起来的图。

常用的绘图软件是 Phylip 和 Snpphylo 。进化树修饰的软件有 MEGA , ggtree 等,推荐网页版工具 iTOL ,无比强大。

外群定根法:当群体的个体的差异很小时,可以引入其他物种作为根。如在对三叶草建树时可以引入水稻的序列作为根进行建树。

PCA是很常见的降维方法,如微生物研究中常用来检验样品分群情况。PCA计算的软件很多,plink可以直接用vcf文件计算PCA,R语言也可以进行PCA计算。

PCA图在群体重测序中有如下几种作用:

进化树和PCA能够看出来群体是不是分层的,但是无法知道群体分成几个群合适,也无法看出群体间的基因交流,更无法看出个体的混血程度。这时候就需要群体分层图了。

可以将进化树和群体分层图结合进行展示,如下图:

先了解下概念,此处借鉴基迪奥生物网站的解释[22]。

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

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

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

参照王荣焕等[23]的方法进行LD参数计算:

随着标记间的距离增加,平均的LD程度将降低,呈现出衰减状态,这种情况叫LD衰减。LD衰减分析的作用:

GWAS(genome-wide association study),全基因组关联分析,常用在医学和农学领域。简单理解成将SNP等遗传标记和表型数据进行关联分析,检测和表型相关的位点,然后再倒回去找到对应的基因,研究其对表型的影响。这些被研究的表型在医学上常常是疾病的表型;在农学上常常是受关注的农艺性状,比如水稻的株高、产量、穗粒数等。GWAS思想首次提出是在心肌梗塞的治疗上[24],首次应用是在2005年的文章上[25]。

目前使用最广泛的模型是混合线性模型[26]:

所有的参数软件(如Emmax)会自动完成计算。

GWAS结果文件通常只有两个图,一个是曼哈顿图,另外一个是Q-Q图。一般是先看Q-Q图,如果Q-Q正常,曼哈顿图的结果才有意义。

MSMC(multiple sequentially Markovian coalescent)[27],底层算法很复杂,类似于PSMC。MSMC的主要功能是推断有效群体大小和群体分离历史。

这样看起来更直观:

LAMP(Local Ancestry in Admixed Populations,混杂群体的局部族源推断),用于推断采用聚类的方法假设同时检测的位点间不存在重组情况,对每组相邻的 SNP 进行检测分析[28],在运算速度和推断准确度上都有了质的飞跃。

用于推断群体分离和混合[29]。图是这样的:

测序方案关系到后续的分析,不同的样本量对应不同的测序方法和分析方法。

[1]. 自然选择(维基百科)

[2]. Kimura, Motoo. "Evolutionary rate at the molecular level." Nature . 217.5129 (1968): 624-626 .

[3]. 测序价格变化趋势

[4]. plant genomes

[5]. DePristo, Mark A., et al. "A framework for variation discovery and genotyping using next-generation DNA sequencing data." Nature Genetics . 43.5 (2011): 491.

[6]. Biswas, Shameek, and Joshua M. Akey. "Genomic insights into positive selection." ** TRENDS in Genetics. 22.8 (2006): 437-446.

[7]. Pavlidis, Pavlos, et al. "Sweed: likelihood-based detection of selective sweeps in thousands of genomes." Molecular biology and evolution 30.9 (2013): 2224-2234.

[8]. Wang, Muhua, et al. "The genome sequence of African rice (Oryza glaberrima) and evidence for independent domestication." Nature Genetics 46.9 (2014): 982.

[9]. Bamshad, Michael, and Stephen P. Wooding. "Signatures of natural selection in the human genome." Nature Reviews Genetics 4.2 (2003): 99.

[10]. Siewert, Katherine M., and Benjamin F. Voight. "BetaScan2: Standardized statistics to detect balancing selection utilizing substitution data." BioRxiv (2018): 497255.

[11]. Yu, N.Jensen-Seaman MIChemnick LRyder OLi WH (March 2004). Genetics . 166 (3): 1375–83.

[12]. Wang, Wensheng, et al. "Genomic variation in 3,010 diverse accessions of Asian cultivated rice." Nature 557.7703 (2018): 43.

[13]. Li, C., Zhou, A. &Sang, T. Rice domestication by reducing shattering. Science 311, 1936–1939 (2006).

[14]. Tajima, Fumio. "Statistical method for testing the neutral mutation hypothesis by DNA polymorphism." Genetics 123.3 (1989): 585-595.

[15]. Du, Xiongming, et al. "Resequencing of 243 diploid cotton accessions based on an updated A genome identifies the genetic basis of key agronomic traits." Nature Genetics 50.6 (2018): 796.

[16]. Lu, Kun, et al. "Whole-genome resequencing reveals Brassica napus origin and genetic loci involved in its improvement." Nature communications . 10.1 (2019): 1154.

[17]. Zhou, Z., Jiang, Y., Wang, Z. et al. Resequencing 302 wild and cultivated accessions identifies genes related to domestication and improvement in soybean. Nat Biotechnol 33, 408–414 (2015).

[18]. Liang, Z., Duan, S., Sheng, J. et al. Whole-genome resequencing of 472 Vitis accessions for grapevine diversity and demographic history analyses. Nat Commun 10, 1190 (2019).

[19]. Alexander, D.H., Lange, K. Enhancements to the ADMIXTURE algorithm for individual ancestry estimation. BMC Bioinformatics 12, 246 (2011).

[20]. Francis, Roy M. "pophelper: an R package and web app to analyse and visualize population structure." Molecular ecology resources 17.1 (2017): 27-32.

[21]. http://www.royfrancis.com/pophelper/articles/index.html .

[22]. https://www.omicshare.com/forum/thread-878-1-1.html .

[23]. WANG Rong-Huan, WANG Tian-Yu, LI Yu. Linkage disequilibrium in plant genomes[J]. HEREDITAS , 2007, 29(11): 1317-1323.

[24]. Ozaki, K., Ohnishi, Y., Iida, A. et al. Functional SNPs in the lymphotoxin-α gene that are associated with susceptibility to myocardial infarction. Nat Genet 32, 650–654 (2002).

[25]. Klein, Robert J., et al. "Complement factor H polymorphism in age-related macular degeneration." Science 308.5720 (2005): 385-389.

[26]. Yu, Jianming, et al. "A unified mixed-model method for association mapping that accounts for multiple levels of relatedness." Nature genetics 38.2 (2006): 203.

[27]. Schiffels, Stephan, and Richard Durbin. "Inferring human population size and separation history from multiple genome sequences." Nature genetics 46.8 (2014): 919.

[28]. Sankararaman, Sriram, et al. "Estimating local ancestry in admixed populations." The American Journal of Human Genetics 82.2 (2008): 290-303.

[29]. Pickrell, Joseph K., and Jonathan K. Pritchard. "Inference of population splits and mixtures from genome-wide allele frequency data." PLoS genetics 8.11 (2012): e1002967.

[30]. Chen, Jia-hui, et al. "Genome-wide analysis of Cushion willow provides insights into alpine plant divergence in a biodiversity hotspot." Nature communications 10.1 (2019): 1-12.

[31]. 孙宽,侯一平。法医族源推断的分子生物学进展 [J]. 法医学杂志 ,2018,34 (03):286-293.

[32]. genek.tv