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文件内容如下：
里面的每一条记录并不是一个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时代间隔就是从个体出生到他的孩子出生的时间，可以根据经验定
