首先肯定是要call SNPs的，最好用bcftools pileup，避免后面报错
可以用一个管道一次性搞定：
bcftools mpileup -Ou -I -f reference.fasta NL.marked.bam| \
bcftools call -c -Ov| vcfutils.pl vcf2fq -d 10 -D 100| gzip > diploid.fq.gz
## -Ou : 输出未压缩的bcf文件
## -I : skip indels
## -f reference fasta
## -c 使用原始版本的的calling算法(对应于新的-m)
## -Ov ：输出未压缩vcf文件
## -d -D ：depth的上下限
也可以先call出VCF文件：
bcftools mpileup -Ou -I -f reference.fasta \
NL.marked.bam|bcftools call -c -Ov -o vcf_bcftools.vcf
- 一定把这类管道的代码放在脚本里跑，不要nohup，否则很有可能报错。
- 建议用bcftools mpileup，而不是samtools，因为两者有细微差别。
- 用 bcftools call的时候选择参数 -c 不要用-m，两者call variance的方法不同得到的文件格式似乎也有差别。避免pl脚本和psmc识别不了
- 在高版本(我的是1.11)的bcftools中，call和view为两个命令。PSMC软件github上给出的命令为用bcftools view得到vcf，实际在高版本中要用call。
- bcftools call的时候不要使用-v选项，不要只输出variance。因为PSMC的原理是要考虑全基因组片段的，只看SNPs肯定有问题。
- 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曲线平移一个数量级。
