群体遗传 — 核苷酸多样性π

**核苷酸多样性（nucleotide diversity），记为π，是分子遗传学中一个重要的概念，用于量化种群内部或不同种群间的遗传多样性。**这一概念由根井正利和李文雄在 1979 年提出。核苷酸多样性的计算基于从一个种群中获得的多个样本 DNA 序列上相同位点碱基差异的平均值，反映了群体内不同个体 DNA 序列间的平均碱基差异比例，从而广泛用于表征种群的遗传多样性水平。这种方法用数学公式表示如下：

xi 、xj分别代表第i个和第j个序列的对应频率[注 1]。πij则为两个序列之间不同位点所占百分比[注 2]。n则为总样本数。

选择压力分析

**通过核苷酸多样性进行选择压力分析的本质是基因组上受选择区域的变异位点杂合率下降。**如下图所示，当一个有利突变发生后，这个突变基因的适合度越高，就越容易被选择固定。当这个位点被快速固定之后，与此基因座紧密连锁的染色体区域，由于搭车效应也被固定下来，该区域及其连锁的区域表现为多态性降低，纯和度增加。

基因组上受选择程度越高，则杂合度程度越低，对应的核苷酸多样性越低。因此对全基因组的核苷酸多样性进行检测，可以推断出这一种群所受到的选择压力及基因组中受到选择的区域。

上图中展示了甘蓝型油菜地方品种的全基因核苷酸多样性水平显著高于改良后的栽培种，表明油菜栽培种受到了较强的人工选择影响。

扫码关注微信公众号【生信F3】获取更多生物信息学最新知识。

计算方法

1 计算全基因组核苷酸多样性

**输入文件：**单个种群的全基因组VCF文件。

**使用工具：**VCFtools v0.1.16 (https://vcftools.github.io/man_latest.html；Danecek et al.，2011)

使用VCFtools对VCF文件中的特定样本进行扫描以计算单个种群的基因组核苷酸多样性，指令如下：

vcftools --vcf Samples.vcf \--window-pi 20000 \--window-pi-step 10000 \--keep PopA.SampleID.txt \--out PopA.Win20k-Step10k

• --vcf Samples.vcf 表示输入文件为Samples.vcf；

• --window-pi 20000表示20Kb的窗口大小扫描该VCF文件；

• --window-pi-step 10000表示窗口以10Kb的步长在基因组上进行滑动和计算（该参数一般以窗口大小的一半进行设定）；

• --keep PopA.SampleID.txt 表示包含种群样本ID的文件名。

  当VCF文件包含来自不同种群的个体时需要以 --keep（保留一个或多个个体）或 --remove（剔除一个或多个个体）选项对样本加以区分和筛选，以此确保分析结果的准确性。样本ID需要与VCF文件中保持一致，每行单独存储一个样本ID。

• --out PopA.Win100k-Step10k 表示输出结果的文件名前缀。

基于上述命令，可获得PopA.Win20k-Step10k.windowed.pi 结果文件（其文件格式如下图所示）。

注：每个结果文件第一行为表头信息，CHROM为基因组中染色体/Scaffold编号；BIN_START 为每个窗口的起始坐标，BIN_END为终止坐标；N_VARIANTS表示该窗口内所包含的SNPs数目；PI为核苷酸多样度。

2 滑动窗口设置

选择清除分析为什么使用滑动窗口？

利用滑动窗口寻找受进化选择的基因组区间主要有以下三个原因：

1. 由于连锁不平衡的存在，受选择压力的位点不可能单独存在，因此我们需要以基因组区间为单位，而非单个位点；
2. 在测序深度有限或采用混池测序的时候，单个位点检测到的SNP多态性信息准确性不高，但如果我们将范围放大到一定程度做整体考量则可尽可能减小SNP分型错误位点的影响；
3. 加入滑动窗口步长是为了探究连续滑动窗口之间核苷酸多样新是否存在剧烈波动，有利于较为精细的定位受到自然选择的核心区域。

因此，在进行滑动窗口分析的时候，窗口大小的选择是非常重要的。**如果窗口过小，窗口内的SNP数可能会不足，杂合率的计算可能会受到极端值的影响，从而造成假阳性；窗口过大的话，如果潜在受选择区间的杂合率本来很小，但是由于受到其他相邻未受选择区域的影响，又可能会稀释掉原本显著的受选择区间。**既然这么重要，那滑动窗口的大小怎样选择呢？

如何确定滑动窗口的最优大小？

利用核苷酸多样性进行选择清除分析时窗口大小的设定一般要考虑三个主要因素：

1. 基因组组装质量

   如果在上游SNPs的获取等数据处理中，其所用参考基因组质量不高且序列拼装得较为碎片化（scaffold水平）时，则建议窗口不要取得太大（容易导致在拼装很短的scaffold上出现扫描区域长度的不均匀），一般选择10或20K大小的窗口进行扫描；

2. LD衰减距离

   不同物种随繁殖力和世代间隔的差异，LD 衰减速度也不同。与全基因组关联分析中通过LD衰减距离确定候选基因的范围相似，在种群规模较大时，我们最好结合LD衰减距离来确定滑动窗口大小，此时推荐选择LD衰减距离的一半左右作为窗口大小。

3. 变异位点密度

   在2012年PNAS的一篇研究猪选择进化的文章里面，作者给我们提供了一个思路（文章标题：Strong signatures of selection in the domestic pig genome）。文章作者在进行滑窗分析的时候，尝试了从10-550K的窗口大小，并且统计了不同窗口大小里面的SNP分布的直方图。

作者发现，当窗口太小的时候，小于20个SNP的窗口有很多，不利于后面的杂合率的分析。但是如果窗口过大的话又像我们刚刚提到的，可能会检测不到较小的选择区间。所以文章作者在这里选择了150Kb的这样一个适中窗口大小来进行后续的分析。

综上所述，在针对不同物种的研究的时候，可以选择尝试不同的滑动窗口来进行分析，选择比较合适的窗口来进行后续的选择分析。如果我们还是拿不太准，则可按以下建议确定窗口大小：

1）先选择一个较小（10K，20K或50K等）的窗口进行分析（窗口设置太大容易导致选择信号被基因组的背景信号掩盖）；

2）查看输出结果文件中 N_VARIANTS 这一列记录的各个窗口的SNP数量，一般而言，只要保证大部分窗口含有足量的SNP（>20个SNP）可用就可以保证分析的可靠性了。

3 定位选择压力区间

得到结果后，我们首先可以从全基因组层面解析该群体的平均核苷酸多样度水平（平均值或者中位值）和分布模式。以PopA.Win20k-Step10k.windowed.pi文件为例，使用R语言进一步对全基因组窗口化的核苷酸多样性计算结果进行解析。

# 在R环境下将结果文件“PopA.Win20k-Step10k.windowed.pi”读入‘Diversity’变量中用于后续统计分析；
Diversity<-read.table(file="PopA.Win20k-Step10k.windowed.pi", header=T)
# 计算该群体全基因组窗口化核苷酸多样度平均值、标准差、中位值
mean(Diversity$PI)
sd(Diversity$PI)
median(Diversity$PI)
# 统计该群体全基因组窗口化核苷酸多样度的分布密度
distr<-density(Diversity$PI)；
plot(distr, main="Distribution of Nucleotide diversity", xlab="Bins of Windowed θπ (20K)", ylab="Frequency", col="skyblue2")；

全基因组范围的核苷酸多样度符合正态分布，其平均值和标准差的大小往往反映了该种群基因组总体的遗传进化特征。

通过计算和比较种群基因组的遗传多样性，可用于揭示基因组中哪些区域经历了快速演化，例如受到了自然或人工选择的影响。因此，对于单个种群，我们可以在基因组不同区域进行扫描，观察基因型多样性（杂合率）的变化，那些杂合率下降的区域可能就是受进化选择的区域，其中的基因就是潜在受选择的基因。

例如：从结果文件‘PopA.Win20k-Step10k.windowed.pi’，中提取chr01染色体上的核苷酸多样性计算结果，并解析窗口化π在染色体上的分布和波动。θπ显著偏低的遗传区域暗示可能经历了自然或人工选择作用。

# 从全基因组的结果文件中提取chr01区域的结果用以分析
grep "^chr01\s" PopA.Win20k-Step10k.windowed.pi >> chr01.PopA.Win20k-Step10k.windowed.pi

利用R语言可视化：

# 读取该染色体区域的结果信息到变量‘chr1’中
chr1<-read.table(file="chr01.PopA.Win20k-Step10k.windowed.pi", header=T)
# 直方图显示该染色体上窗口化θπ的分布和波动
plot(chr1$BIN_START, chr1$PI, ylab="Windowed θπ", xlab="Scans on chr01", pch=20, type="h", col="gray")

由图可见，基因组中某些区域的多样度水平会急剧变小而偏离于基因组平均水平，统计上这些区域一般可以第一或第五分位数分别表征其对基因组偏离的显著性（1%和5%显著性差异）。它们往往反映了该区域可能经历了非随机的演化事件（例如自然选择）；上图中，核苷酸多样度水平低于红色虚线（表征基因组第五分位数阈值）的遗传区域则可能经历了选择作用。

4 选择清除分析ROD

对单个种群的分析表明，核苷酸多样性在染色体上并不是恒定不变，不同区域存在较大波动。**然而，多样性显著低的区域并非全部由选择压力导致，其它显著差异的DNA遗传特性，如重组率和突变速率等也会显著影响基因组杂合度。**此外，选择压力一般发生在两个种群之间。因此，对单个种群分析核苷酸多样性找到的潜在受选择区间会存在大量的假阳性。

在当前的群体分化研究中，我们更多的选择同时计算两个种群（野生和栽培种群间，驯化和未驯化间）的核苷酸多样性并进行比较，根据比较后的显著差异结果得到潜在受选择的区间。

ROD (reduction of diversity) 一般通过比较野生群体和驯化群体（或其他参照群体）的核苷酸多样性计算得来，其计算公式可以简单表示为：

ROD = π(野生群体) / π(驯化群体)

这个公式的结果是一个比值，表明了驯化种群相对于野生种群的核苷酸多样性减少程度。ROD值越大，表明驯化群体的核苷酸多样性相对降低的越多，可能受到的选择压力越大。

我们一般选取最大的1%和5%ROD值所在窗口作为选择清除分析候选区间。

扫码关注微信公众号【生信F3】获取更多生物信息学最新知识。

参考

1. Danecek, P., Auton, A., Abecasis, G., Albers, C. A., Banks, E., DePristo, M. A., Handsaker, R. E., Lunter, G., Marth, G. T., Sherry, S. T., McVean, G., Durbin, R. and Genomes Project Analysis, G. (2011). The variant call format and VCFtools. Bioinformatics 27(15): 2156-2158.
2. https://cn.bio-protocol.org/bio101/e1010601

本文由mdnice多平台发布