PSMC软件分析群体历史有效群体大小流程

1 文件转换

基因组文件格式为.bam，callsnp之后，再转换为.fq.gz格式才能做为PSMC输入文件

bcftools mpileup 和 bcftools call可以进行SNP calling，

bcftools mpileup -Ou -I -f ref.fa .sorted.mardup.bam | bcftools call -c -Ov | vcfutiles vcf2fq -d 10 -D 100 | gzip > dilpoid.fq.gz

需要注意的是mpileup命令虽然也会输出vcf格式文件，但是并不直接进行snp calling。
下面的命令可以生成vcf格式文件

bcftools mpileup -I -f ref.fa .sorted.markdup.bam > mpileup.vcf 

直接 生成的vcf文件内容如下：

mpileup直接生成的vcf文件

里面的每一条记录并不是一个SNP位点，而是染色体上每个碱基的比对情况汇总，这种信息官方称为genotype likelihoods。

call命令才是真正执行SNP calling的程序，基本用法如下

bcftools call mpileup.vcf -c  -Ov -o variants.vcf

在进行SNP calling 时，必须选择一种算法，有两种calling算法可供选择，分别是-c 和 -m 参数。-c参数对应consensus-caller算法，-m参数对应multiallelic-caller算法，后者更适合多种allel和罕见变异的calling。

2 生成psmc.fa和psmc文件

utils/fq2psmcfa -q20 diploid.fq.gz > diploid.psmcfa

psmc -N25 -t15 -r5 -p "4+25*2+4+6" -o diploid.psmc diploid.psmcfa

-p : 64个时间间隔，在计算时需要估计28个参数，将64个时间间隔划分一下，第一个划分了4个时间间隔，在这四个时间间隔估计1个Ne的参数；25*2划分了25组2个时间间隔，每两个时间间隔估计一个和Ne有关的参数。

3 画图

生成psmc文件之后开始画图

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

-u为突变率，-g为世代间隔
-u 的说明是 absolute mutation rate per nucleotide [2.5e-08]
注意，这里填的是per generation per site的mutation，不是per year per site 的mutation，否则会导致PSMC曲线平移一个数量级。

4 不同物种绘制在同一个图中

cat sample1.psmc sample2.psmc sample3.psmc > combine.psmc 

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

-u 一般根据文献来看用什么，不过突变率一般哺乳动物不会相差很大，-g时代间隔就是从个体出生到他的孩子出生的时间，可以根据经验定