【记录】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