生信分析16：基于PSMC估计有效群体大小

本次推送是文献分享16的对应内容。

Fig 1

有效群体大小的估计或者种群历史动态的估计经常出现在基因组/重测序相关文章中，可视化形式一般是Fig 1。

Fig 1的解读：由于该分析是以个体来推测群体，所以作者随机选取了4份圆果化香 (Platycarya longipes) 和4份化香树 (Platycarya strobilacea) 的重测序数据，分别比对到各自的参考基因组上。然后基于每一份数据去估计该物种的历史有效群体大小。图中横坐标是距离今天的时间，从左到右是距离今天104年-距离今天107年。纵坐标是有效群体大小。图中浅蓝色线有四条，代表随机选取的4份圆果化香的估计值，一条深蓝色线代表4份圆果化香的估计值拟合成的值。黄色代表化香树。阴影部分的竖线代表两物种的估计分化时间。g是指该物种繁殖一代的时间，单位为年。μ是该物种的突变速率。

通过这个图我们可以看出两物种的分化时间，以及分化之后各自种群大小的变化，可以与各种地质事件结合。

相关概念

有效群体大小 (Effective population size, Ne)

（1）是指与实际群体具有相同基因频率方差或相同杂合度衰减率的理想群体大小,它反映了群体平均近交系数增量的大小以及群体遗传结构中基因的平均纯合度。

（2）有效群体数量并不是根据统计群体数量来定义的，而是根据遗传变异在群体中的表现来定义群体大小的。如果一个群体在遗传漂变的作用下，等位基因频率在一定时间内发生的改变很缓慢，那么我们就认为该群体的有效群体数量很大；反之，如果等位基因频率变化很大，那么我们就认为该群体的有效群体数量很小。

（3）在一个理想群体中，在随机遗传漂变影响下，能够产生相同的等位基因分布或者等量的同系繁殖的个体数量，通常小于绝对的种群大小。

以上是我整理的不同来源的定义方式，看起来还是有些模糊。我自己的理解是在一个群体中，可以交配产生后代的个体数。比如在一个10000人的群体中，有2000人不具有生殖能力，有1000人不具有生殖意愿，只有剩余的7000个可以“自由交配”产生后代。那我们根据后代的基因型只能估计出这7000个自由交配的人，而无法判断出实际人群大小。（纯属个人的理解，可能有错误，欢迎讨论）

物种的突变速率μ

这是估计Ne时必须用到的一个值。Fig 1中标记的突变速率为2.06×10-9，单位是每代每碱基。不同物种的突变速率是有差异的，因此不能套用同一个值。对于常见的物种拟南芥、水稻、玉米等，突变率一般是可以查到。但是大多数物种的μ是查不到的，可以根据公式μ = Ks/2T进行估计，其中Ks是两物种间的synonymous divergence, T是两物种间的divergence time。

PSMC原理介绍

Fig 2

估计Ne的最经典软件是2011年李恒博士开发的PSMC软件，尽管该软件是针对动物模型，但是在植物中也有较多的应用，并且化香树的文章也是用的PSMC，所以本次推送以PSMC为例进行讲解。

PSMC的全称是成对序列马可夫共祖先分析 (Pairwise Sequentially Markovian coalescent)，利用单个个体的重测序数据来推测该个体所属的种群在历史上各阶段的有效群体大小。该流程的核心是计算最近共同祖先时间（the Time since the Most Recent Common Ancestor, TMRCA）。假设TMRCA属于某个时期的片段比例越少，则该时期的Ne越大。

PSMC适用于二倍体物种，且该群体需要随机交配。对于二倍体物种，基因组上是存在广泛的纯合和杂合位点的（Fig 2）。在最近共同祖先中，认为均是纯合位点，随着时间的推移，经历了重组和突变积累了杂合位点。某区段的杂合位点越多说明经过的代数/时间越久，所以该区段最近共祖先离现在越远。

基于以上原理，可以根据单个个体的重测序数据估计该个体所在群体的Ne.

需要注意的是PSMC只适用于2万年-300万年之间的Ne估计，更短或更长均会导致准确性降低，且重测序数据的深度不低于18×.

PSMC实操

Fig 3

软件的安装参考https://github.com/lh3/psmc

输入文件（Fig 3）：

1、重测序数据比对到参考基因组的bam文件

2、参考基因组

第一步 生成二倍体一致性序列（Fig 4）

Fig 4

需要注意两个参数：

-d 指定最小深度，作者建议设置成平均深度的三分之一。

-D 指定最大深度，作者建议是平均深度的两倍。

Fq2psmcfa是psmc自带的一个命令，用于将fq.gz转换为fa格式。

第二步 运行psmc (Fig 5)

Fig 5

Fig 6

-N  迭代的最大次数

-t  2N0的最大世代时间，代表TMRCA的上限

-r  起始θ/ρ率

以上三个参数对于大部分植物使用，可不用更改。

-p  代表PSMC有效群体大小变化的图随时间变化划分的时间阶段单位数量，例如默认的-p “4+25*2+4+6”表示从古至今依次经历了一个4单位+25个2单位+一个4单位+一个6单位的时期。数字越大需要的计算资源越多。可根据经验调整，也可查找近缘物种的值用作参考。

（参考自生信技工）

第三步 画图（Fig 7）

Fig 7

Fig 8

PSMC内置了画图脚本psmc_plot.pl

其中两个最关键的参数是-u指定该物种的突变率，-g指定该物种繁殖一代的时间，单位为年，比如人默认为25年，黄瓜、西瓜等认为是1年。

-p指定输出pdf格式（Fig 8）。

Fig 9

当然除了PSMC软件，近些年也开发出了更多的用来推断种群历史动态/有效群体大小的软件，可能会在后续的文献分享-生信分析中涉及。

最后，推荐一个网站（https://methodspopgen.com/methods-to-infer-populations-history/，Fig9），里面推荐了很多群体遗传学软件。

参考链接

[1] https://www.jianshu.com/p/33ab4d2c81f6

[2]https://zhuanlan.zhihu.com/p/364652712

[3]https://baike.baidu.com/item/%E6%9C%89%E6%95%88%E7%BE%A4%E4%BD%93%E5%A4%A7%E5%B0%8F/5975113

[4]https://yanzhongsino.github.io/2022/10/17/bioinfo_psmc/

[5]https://cloud.tencent.com/developer/article/1911384

[6]Li, H., Durbin, R. Inference of human population history from individualwhole-genome sequences. Nature 475, 493–496 (2011). 

本文使用 文章同步助手 同步