群体基因组学——概述和图表详解

最近看群体遗传学的文章，苦于不能理解研究者做的另人眼花缭乱的图表……哎，该补补基础了。这篇笔记先从樊龙江主编的《植物基因组学》教材开始学习基础概念，有部分摘自知乎和一些文献，加上这个领域论文图表的解读和自己的理解做汇总整理。下一篇笔记自己实操跑一下GWAS流程。

1. 群体基因组学概述

群体基因组学，将基因组数据和技术与群体遗传学理论体系结合，通过覆盖全基因组范围的多态性推测全基因组效应和位点特异性效应。说人话，这门学科的作用就是通过检测SNP等的信息来回答关于物种起源、进化以及环境适应的分子机制。

群体基因组研究的基础是高质量的群体基因型数据，在获得群体原始测序数据后，如何高效准确的进行变异检测和群体基因分型是后续研究的关键（因为要用到大量的测序数据，所以也很烧钱）。

变异检测的流程可分为数据质控——读序联配——变异检测——基因分型

上图的bam文件通过软件GATK检测结构变异或者SNP就是变异检测，这部分内容我前面组装三代基因组测序数据中使用过，因为我只有一个二倍体基因组，所以当时也没用到基因分型，GATK用法可以看我的这篇笔记：

0基础学习三代基因组测序组装（9）——GATK检测植物基因组SNP和INDEL变异 - 我的小破站 (shelven.com)

2. 群体基因组进化分析方法

2.1 群体系统发生树和群体结构

群体发生树反映特定群体中个体间亲缘关系。一般我们做进化树都是用单拷贝基因，而做群体研究构建进化树用了整个基因组数据信息（比如全基因组SNP、Indel、SV等），能一定程度抵消横向基因转移和单个基因速率差异带来的分析误差，比单基因构建的结果更接近真实进化关系。全基因组SNP构建进化树（一般是把SNP merge到一起），个体间两两比较，根据SNP差异计算个体差异距离，构建群体差异矩阵。软件有MEGA、PHYLIP、FastTree等。

还有一个非常重要的概念：群体遗传结构，指基因型或者基因在空间和时间上的分布模式，包括种群内的遗传变异和种群间的遗传分化。种群遗传结构反映物种进化历史中的一些特殊进化事件，软件有Structure、Frappe等。总体流程确定亚群数（K），然后计算各个体归属第K亚群的概率Q值。以及主成分分析（PCA）也常用于群体结构划分。

上图是华农王茂军老师课题组做的不同二倍体棉花的物种发生树和群体结构。

2.2 群体遗传分化

2.2.1 固定系数 Fst

在明确一个物种存在群体结构后，通产需要量化该物种亚群间的分化程度。固定系数（Fst）反映群体等位基因杂合性水平，Fst越大，种群间遗传分化程度越大。

Fst= (Ht-Hs）/Ht
Hs:亚群体中的平均杂合度
Ht:复合群体中的平均杂合度

2.2.2 核酸多态性 𝝅

做群体的文献中经常出现这个指标类似𝝅野生/𝝅栽培这样的系数，核酸多态性（Nucleotide diversity）是用于衡量群体的多态性的，为群体内平均每两个样本每个位点的差异核苷酸数量。一般称为𝝅，自然群体多态性一般高于栽培群体。

需要注意下这个𝝅的单位，一般都在 10^(-3) 这个数量级。

这些群体分化的相关系数计算可以用软件Arlequin、FSTAT、VCFtools等。

2.3 基因流（渐渗）

基因流（也称为基因迁移）指从一个物种的一个种群向另一个种群引入遗传物质，从而改变群体的遗传组成。

当一个种群中的一个个体迁移到另一个群体，就会把基因带到新的群体，也就是产生“基因的流动”。基因流越大，群体间的相似性越大，会导致群体间基因频率和基因型频率呈现哈迪-温伯格平衡（在一个随机交配的大群体中，没有突变、选择等外界环境因素的影响，基因频率和基因型频率世代恒定）。

自然选择和遗传漂变会使群体之间的差异增加，而基因流的作用是“弱化”群体间的遗传差异，群体之间趋于一致。在植物中，种子扩散和花粉传播是植物基因流的最主要方式。从上面的描述也能看出，基因流与Fst成反比；与种群间地理距离成反比。检测方法主要有D-统计（ABBA-BABA）、Treemix、Migrate-N等。

2.4 种群动态和进化历史

种群动态指种群大小或密度随时间或空间的变化。

有效群体大小（Ne）指与实际群体有相同基因频率方差或者相同杂合度衰减率的理想群体含量，通常小于绝对群体大小，决定群体平均近郊系数增量大小，反映群体遗传结构中基因的平均纯和速度。种群动态理论与方法主要包括种群大小及其变化、种群生长模式的量化描述，以及引起种群变化的外在环境因素。基于全基因组的种群历史动态分析方法主要有PSMC、MSMC等。

上图展示了不同地理来源的拟南芥的有效群体大小变化。可以看到拟南芥在9万～12万年前，由祖先群体在非洲开始分离形成亚群。欧亚祖先群体在 8 万年前从非洲分化形成，然后欧亚拟南芥群体在4万年前进一步分化。这一模式与包括人类在内的多种物种分化时间模式非常相似，这意味着间冰期和洪积世时期（即9万～12万年前）的气候事件对物种分布十分重要。

3. 基因组选择信号

自然群体区别于栽培群体的最大特征是丰富的遗传多样性。草本作物祖先自然群体多态性（𝝅）一般在 0.003~0.008，高于相应的驯化作物遗传多态性。

遗传瓶颈效应：由于环境骤变或人类活动，某一生物种群的规模迅速减少，仅有少部分个体能够顺利通过瓶颈事件，之后可能经历一段恢复期并产生大量后代。由于连锁不平衡（LD）的作用，不仅仅是驯化基因的多态性降低，其侧翼区域遗传多样性也会降低。
有害性突变（deleterious SNP，dSNP）指导致生物个体的整体适合度（fitness）减低的突变。植物被驯化的过程中，基因组上有害突变的数量、频率会不断增加和累积，导致驯化后的作物在原本自然环境中适合度降低，称为作物的“驯化成本”。
在作物群体基因组中，一个明显的特征是存在大量来自野生祖先种的遗传渐渗。来自野生群体的基因渗入可以有效增加栽培种的环境适应性。

选择分析可以在宏观和微观两个尺度下进行选择基因和信号鉴定。

3.1 宏观进化尺度

宏观进化尺度

基于基因：Ka/Ks（非同义替换比同义替换）、MKT检验（比较物种内和物种间的Ka/Ks值）

基于进化速率：HKA

宏观尺度检测选择信号的方法，通常在亲缘关系相近的物种之间进行同源基因序列的比较，判断基因是否在某一物种或进化分支上存在加速进化的情况。

判断加速进化的背景指标是Ka/Ks，比较每个位点区域中非同义替换率与每个位点的同义替换率。同义突变总体是中性进化的（中性基准），如果存在过多的非同义突变意味着存在倾向于蛋白序列改变的正向选择（正选择表现为固定某种有利等位基因，负选择表现为清除不利的等位基因）。MKT检验本质是比较物种内和物种间的Ka/Ks值，中性情况下相等，如果物种间的比值显著高于物种内比值，意味着物种中存在正向选择信号，如果是在作物驯化中，可能是人工选择的信号。

3.2 微观进化尺度

微观进化尺度

基于等位基因频谱：Tajima’s D, Fu andLi’sD, Fay and Wu’sH, CLR, Hp

基于连锁不平衡：LRH、iHS、连锁不平衡衰减（LDD）、IBD分析

基于群体分化：LKT、LSBL

组合方法：CLR、XP-CLR（适合SNP大数据集）、DH检验、CMS

微观尺度鉴定方法，在正向选择的作用下，有利等位基因在群体中频率会变高甚至固定。当有利等位基因和周围变异达到较高频率的时候，这段区间在群体水品就会呈现遗传多态性降低的现象，也就是选择性清除（selective sweep）。

比如，栽培物种在选择性清除区域的遗传多样性显著降低，就是驯化区域的典型特征。选择性清除还会使连锁不平衡区块延长、群体间固定指数增大、Tajima‘s D值为负且显著偏离零。中性模型为基础的检测方法很多，Tajima‘s D检验是最经典的。

总体来说，选择信号检测原理和方法可以分为以下5类：

A 基于群体多态性。原理：受选择位点两侧序列多态性因连带效应保持很低水平。在驯化选择区间内，𝝅野生/𝝅栽培的值变高。
B 基于等位基因频谱。原理：有利突变不断在群体中固定，当完全固定时，周围可能由于随机突变出现稀有等位基因
C 基于连锁不平衡。原理：选择性清除会引起包含等位基因的单倍体纯和性提高，与周围变异的连锁不平衡程度变高。当群体处于正向选择作用下时，基因突变及其连锁位点在正选择的作用下，在短时间内会达到较高频率，形成大片段的纯合单倍型。
D 基于群体分化。随着受选择的等位基因频率上升，导致与原群体的Fst变大。Fst的取值范围为0-1，1表示群体间完全分化的位点，0表示在群体间完全没有分化的位点。
E 组合方法。

在检测选择信号时，需要考虑遗传漂移（genetic drift）对选择信号的影响，为了降低该影响，通常在全基因组扫描窗口的基础上，仅考虑极端值（前5%）区域作为选择位点（这点很重要，后面实操的时候会说）。

选择信号检测方法都是基于自下而上的，通过群体基因组扫描鉴定具有选择信号的基因组区域，从而挖掘可能与某种表型或者适应性相关联的基因，假阳性高。

也可以用自上而下的方法进行佐证，比如获得该物种多个时间点或者多个世代的群体基因组数据，就可以直接检测该位点等位基因频率变化加以验证，比如用QTL、GWAS等传统定位方法。

4. 连锁不平衡（LD）

连锁不平衡理论对群体研究非常重要，前面说到基于连锁不平衡的原理检测选择信号，这里再加深一下理解。

连锁不平衡（linkage disequilibrium，LD）：在某一群体中，两个基因同时遗传的频率大于随机组合的频率。

连锁不平衡的三个衡量指标：

D值：D = P(AB) - P(A) X P(B)，也就是实际概率减去独立遗传的理论概率，D值绝对值越大，连锁程度越大。
D‘值：理解维归一化之后的D值，归一化之的值可以用于比较不同基因连锁程度的大小。0表示完全连锁平衡，独立遗传；1表示完全连锁不平衡。
r平方

r值定义如下：

通常文献里也是用r平方来表征连锁不平衡程度，r平方等于0时，表示完全连锁平衡，独立遗传；r平方等于1时, 表示完全连锁不平衡。

这里就要引出另一个很重要的概念——LD衰减，也是做GWAS出现最多的图。

4.1 LD衰减曲线

LD衰减指位点间由连锁不平衡到连锁平衡的演变过程；LD衰减的速度在不同物种间或同物种的不同亚群间差异非常大。所以，通常会使用“LD衰减距离”来描述LD衰减速度的快慢。

衰减距离，是当平均LD系数衰减到一定大小的时候，对应的物理距离。（这个一定大小是没有特别规定的，比如可以到一半大小，可以到0.5，可以到0.1）

上面的图描绘了不同水稻群体的LD衰减曲线，横坐标是物理距离（kb），纵坐标是LD系数（r^2）。

LD衰减曲线就是利用曲线图来呈现基因组上分子标记间的平均LD系数随着标记间距离增加而降低的过程。大概的计算原理就是先统计基因组上两两标记间的LD系数大小，再按照标记间的距离对LD系数进行分类，最终计算出一定距离的分子标记间的平均LD系数大小。

我们看左边的图，Japonic这个亚群（红色的线）在基因组100Kb距离的平均LD系数大小为0.35，到了200kb距离，对应的LD系数降低到了不到0.3。LD衰减速度在不同亚群是不一样，Japonica衰减速度最慢（其他亚群在物理距离很小的时候LD系数降就从0.35衰减到了0.3），衰减距离最大。

知道这个衰减速度快慢和衰减距离有什么用呢？

LD衰减速度越慢，形成的单体型区块越大，关联分析中需要的群体和标记数目越少，但定位越不精准。此外，同一个连锁群上，LD衰减地慢说明该群体受到选择，一般来说，野生群体比驯化改良群体LD衰减快，遗传多样性高。而驯化选择会导致群体遗传多样性下降，位点间的相关性（连锁程度）加强。所以，驯化程度越高，选择强度越大的群体，LD衰减速度越慢。

也就是说左边这个图的4个水稻品种中，Japonic是驯化程度最高的。

LD衰减距离应用：

辅助分析进化和选择，同一个连锁群上，LD衰减地慢说明该群体受到选择。
一般而言，两个位点在基因组上离得越近，相关性就越强，LD系数就越大。反之，LD系数越小。也就是说，随着位点间的距离不断增加，LD系数通常情况下会慢慢下降。
一般而言，LD系数大于0.8就是强相关。如果LD系数小于0.1，则可以认为没有相关性。如果LD衰减到0.1这么大的区间内都没有标记覆盖的话，即使这个区间有一个效应很强的功能突变，也是检测不到关联信号的。所以，通常可以通过比较LD衰减（到0.1）距离和标记间的平均距离，来判断标记是否对全基因组有足够的覆盖度。
如果GWAS检测到显著关联的区间后，则可以通过进一步绘制局部的LD单倍型块图（这个后面讲），来进一步判断显著相关的SNP和目标基因间是否存在强LD关系。
基因组上那些受选择的区域相比普通的区域，LD衰减速度也是更慢的。
判断GWAS所需标记量，决定GWAS的检测效力以及精度； GWAS标记量 = 基因组大小 / LD衰减距离。

4.2 单倍型块（Haplotype Block）

描述LD不平衡的另一个常用的图就是LD单倍型块图。

单倍型块，即连锁不平衡区域，是指同一条染色体上处于连锁不平衡状态的一段连续的区域。单倍型块分析可以用于筛选tag SNP、确定候选基因的范围等。

颜色从白色到红色，代表连锁程度从低到高，方框中的数值为LD系数r^2，为了美观，这里将 r2^2乘以了100。

上面的图中每个正方形是两个相邻SNP的LD值计算结果，如果大于设定的阈值（这里估计是大于94），那么就构成一个block，图中黑框为一个block，共有三个。可以从三个block中分别选择一个SNP作为Tag SNP来分别代表这三个block。

5. 全基因组关联分析（GWAS）

全基因组关联分析（Genome-Wide Association Study, GWAS），以连锁不平衡（linkage disequilibrium, LD）为基础，通过分析数百个或者数千个个体的高密度分子标记的分离特征（一般是上万个或者上百万个SNP标记），筛选与复杂性状表现型变异相关联的分子标记，进而分析分子标记对表现型的遗传效应。

插一句题外话，GWAS中最后一个字母已经是study，在写英文文章的时候不要写成GWAS study……

传统的QTL定位研究通常以2个亲本杂交群体为研究对象, 通过连锁作图定位目标性状位点。局限性是投入大量资源构建数量庞大的重组群体。而关联分析则可以利用个体在全基因组范围的遗传变异标记进行多态性检测，获得更高分辨率的定位结果（单碱基水平）, 遗传变异来源也更为广泛，根据统计量或者显著性P值筛选最有可能影响该性状的遗传变异，往往能定位到比双亲本作图群体中更多的性状关联位点（这可节省太多时间了）。

GWAS的局限性在于可以确定相关位点但不能直接确定基因本身，假阳性也比较高。解决这一局限性有两种策略，一是算法，二是群体的选取。

5.1 GWAS算法模型

目前用的算法模型有以下几种：

1.GLM （Generalized linear model；一般线性模型）
2.MLM （Mixed linear model；混合线性模型）
3.GEMMA、CMLM、SUPER、FarmCPU、Blink

主要介绍前两个，GLM模型以群体结构矩阵 Q或主成分分析矩阵PCs为协变量做回归拟合。

如果两个表型差异很大，但群体本身还含有其他的遗传差异（如地域等），则那些与该表型无关的遗传差异也会影响到相关性。MLM模型可以把群体结构矩阵 Q、亲缘关系矩阵（kinship）或联合利用主成分分析矩阵PCs和亲缘关系矩阵为协变量，把这种位点校正掉（也就是抑制假关联）。

其他模型大多是基于GLM和MLM进行优化，比如GEMMA 计算个体基因型相似性矩阵，排除了LD的干扰。其他暂时就不多说了，感兴趣再深入研究研究。

5.2 GWAS群体选取

群体中丰富的表型变异和充分的遗传重组是GWAS 成功的关键条件，有以下选取策略：

(1) 群体内没有明显的群体结构，样本间没有过近的亲缘关系，同时具有丰富的表型变异
(2) 群体来自具有一定水平遗传分化的不同类群(如不同亚种与亚群)，具有丰富的遗传和表型变异，但同时不同类群之间存在频繁的遗传交流，保证目标性状在不同类群内部也存在一定水平的变异
目前GWAS样本量普遍大于100份

前面也说过，同时要考虑GWAS标记量，计算公式GWAS标记量 = 基因组大小 / LD衰减距离。

在用GATK做基因分型的时候也要注意，结果一般保留maf值大于0.05的 SNP，通常认为这是在驯化选择区间内。

5.3 GWAS结果图解

还是先说明一下，做GWAS需要有三类数据：SNP数据、表型数据和群体结构。这些数据在实操会再次提到。

GWAS的结果通常以曼哈顿图和QQ图来展示。曼哈顿图显示每个SNP在关联分析中的显著性水平; QQ 图反映关联分析的效果。

以下图721份水稻GWAS结果图为例：

5.3.1 曼哈顿图（Manhattan Plot）

图4（A）就是一个曼哈顿图，每个点代表一个SNP，x轴代表SNP在基因组上的遗传位置，y轴为–log10 (P-value)，显著性水平线有红蓝两条（P = 0.01和P= 0.05）。y轴值越大，说明该位点表型的关联程度越大，受到连锁不平衡（LD）的影响，基因组上强关联位点周围的SNP也会呈现出关联性由高到低连续变化的信号强度，从而在P值小的地方出现尖峰（peak）。

上面文章深入分析水稻6号染色体上与抽穗期相关的一个尖峰，将其定位在Hd3a附近, 估计候选区间大概在 2.68–4.62 Mb (图4C)。这个候选区间是如何确定的呢？

我们通常将显著关联SNP在N kb以内的位置确认为相关区间，这个N就是对应的LD衰减距离。确定GWAS鉴定的候选区间后，就可以在候选区间内找到基因功能注释和表型相关的基因，或者有其他组学或者功能研究支持的基因进行验证和深入研究。

5.3.2 QQ图（QQ Plot）

图4（B）是QQ图，本质上就是做两组数据的比较，判断是否一致。每个点代表一个SNP，横坐标是期望的P-value，纵坐标为实际观测到的P-value，虚线代表实际观测和期望的P-value值一致的情况。

我们做关联分析，当然是希望大部分的位点观测值和期望值一样（在对角线，也就是这里的虚线上），也就是这部分位点确定是与性状不关联的。一部分位点在虚线的上方，也就是观测值超过期望值，说明这些位点的效应超过随机效应，也就是这些位点是与性状显著相关的。也就是出现4（B）这个图才是理想的结果。

还有以下三种其他情况：

1 观察值低于期望值（点在对角线下方），可能是模型不合理，P-value被过度矫正。或者是群体中大量SNP位点间存在连锁不平衡，有效位点数（相互间不存在连锁不平衡的位点）明显低于实际位点数，所以P-value的期望值被低估了。
2 观察值和期望值相同，说明没有找到与性状显著关联的位点。
3 观察值显著高于期望值（点全都在对角线上方），也就是所有位点都与某个性状显著相关，说明分析模型不合理，假阳性太多了。

现在基于GWAS分析还衍生出其他诸如TWAS（基因表达量做标记）、PWAS（蛋白做标记）、EWAS（甲基化表观信息做标记）等关联分析方法，做标记的信息不同，本质上也都是和表型做关联分析。

6. 泛基因组（Pan-Genome）

想了想还是把泛基因组也放了进来，虽然泛基因组研究方法和上面的研究方法不一样，但也是群体基因组学的一个重要分支，主要是开展群体基因组重测序和遗传变异的分析工作。

我们知道单一的参考基因组仅能代表物种基因组空间范围的一小部分（同一个物种不同亚群基因组存在差异），以线性参考基因组进行群体变异的分析就会错失一些变异位点（SV&CNV&PAV）。所以这个时候我们就需要获取一个物种的全部遗传信息（只能说是相对的全部），并解析每个个体间的遗传变异，这就诞生了泛基因组（Pan-Genome）研究。