【记录】PSMC软件分析群体历史有效群体大小步骤

---

更新于2021年8月23日
现在出现了更好的算法，能用多个样本估计一个群体的历史群体大小：MSMC详情可见此链接

@PSMC软件分析群体历史有效群体大小流程
首先是软件下载，需要用到PSMC和samtool
PSMC下载地址

1、文件转换

基因组文件格式为.bam，callsnp之后，再转换为.fq.gz格式，关于callSNP博客写的比较明确，我们在此只展示bwt之后已经生成bam的代码
在软件readme中有代码：

samtools mpileup -C50 -uf ref.fa aln.bam | bcftools view -c  \
      | vcfutils.pl vcf2fq -d 10 -D 100 | gzip > diploid.fq.gz

这段代码重点需要注意两个问题，第一个是参考基因组ref.fa文件一定要用生成bam的参考文件，否则很可能在psmc.fa中跑出空文档
第二个可能是由于samtool版本的改变，bcftools view -c需要改成bcftools call -c才能运行成功。
运行代码

samtools mpileup -C50 -uf pig.fa sample.recal.bam | bcftools call -c | vcfutils.pl vcf2fq -d 10 -D 100 | gzip > sample.fq.gz

2、生成psmc.fa和psmc文件

这一步在README中也有

 utils/fq2psmcfa -q20 diploid.fq.gz > diploid.psmcfa
 psmc -N25 -t15 -r5 -p "4+25*2+4+6" -o diploid.psmc diploid.psmcfa

-t 的意思在文章里面应该对应着Tmax，和时间间隔有关的一个参数
$$t_i=0.1\ast e^{i/n\ast log(1+10\ast T_{max})}-0.1$$
$i=0,\dots ,n$ n是设定的时间间隔（64），和下一步中的-p有关系
最重要的参数-p，文中给出的解释为

option specifies that there are 64 atomic time intervals and 28 (=1+25+1+1) free interval parameters.The first parameter spans the first 4 atomic time intervals, each of the next 25 parameters spans 2 intervals, the 27th spans 4 intervals and the last parameter spans the last 6 time intervals.

我的理解就是64个时间间隔，在计算时需要估计28个参数，将64个时间间隔划分一下，第一个划分了4个时间间隔，在这四个时间间隔估计1个Ne的参数；25*2划分了25组2个时间间隔，每两个时间间隔估计一个和Ne有关的参数，后面的都是这个意思。大家如果有兴趣可以数一下：

两边的参数分别是psmc -N25 -t15 -r5 -p "4+25 * 2+4+6"和psmc -N30 -t5 -r5 -p “4+30 * 2+4+6+10”
其他的两个参数“-N25 -r5”，这个在参数里面我找了文档和原文章都找不到，但是应该不会有很大的影响，如果不放心的话也可设置几个参数组测试一下。

3、画图

生成psmc文件之后就是画图的步骤，

perl utils/psmc_plot.pl -u 2e-09 -g 1 -p sample_plot sample.psmc

-u为突变率，-g为世代间隔，用perl编写的脚本需要补充一些库，可以根据错误提示补充更新库。

我当时就遇到库没有更新的问题，（上图）应该是libtif库

可以使用他给出的perl直接绘制多个个样本的历史群体大小，这个回答中 davidaray 的解决方法很细致。
说实话，这个问题我看了很多回答，我觉得github里面的这段代码应该是可以的，虽然我用这个跑不出来。

seq 100 | xargs -i echo psmc -N25 -t15 -r5 -b -p "4+25*2+4+6" -o round-{}.psmc split.fa | sh
cat diploid.psmc round-*.psmc > combined.psmc
perl utils/psmc_plot.pl -u 2e-09 -g 1 -p sample_plot sample.psmc

不同物种绘制到同一张图中

这个回答要感谢@幻影刺客的私信回复，具体大家可以验证一下，有问题欢迎一起探讨。
下面是如何使用软件自带脚本拼接的过程

psmc_plot.pl -u mutation_rate -g generation -Y years -M 'sample1,sample2,sample3' combine_plot combine.psmc

只需将需要出图的 *.psmc合并在一个psmc下，然后在-M参数后面按照样本顺序添加样本名(要和你合并的数量样本数一致哦) 就可以实现合并出图了

同一物种多个样本推算PSMC

这个解答 回答了用多个样本可以怎么算

有一些坑

psmc_plot.pl -u -g 参数：-u 的说明是 absolute mutation rate per nucleotide [2.5e-08]
注意，这里填的是per generation per site的mutation，不是per year per site 的mutation，否则会导致PSMC曲线平移一个数量级。

这个来自于简书
做的时候也参考了一些别人的回答
种群密度怎么算&PSMC推测种群大小历史动态