PSMC模型主要利用二代重测序的比对结果进行构建，看到第一篇用三代构建群体历史的文章，发表在MBE上，名字叫“Genome Assembly of a Living Fossil, the Atlantic Horseshoe Crab Limulus polyphemus, Reveals Lineage-Specific Whole-Genome Duplications, Transposable Element-Based Centromeres, and a ZW Sex Chromosome System”。

主要原理如下：

首先，对参考基因组的fai索引进行切片，并用bedtools make window。

代码如下：

samtools faidx ref.fa
cat ref.fa.fai | cut -f 1,2 > ref.gfile
makewindows -g ref.gfile -w 1000 > windows.bed

将三代序列用minimap2，比对到参考基因组上；比对根据数据类型选择参数，如果是HIFI数据，参数为：

minimap2 -ax map-hifi ref.fa  $species.gz   | samtools sort > $species.sort.bam

如果是ONT数据，参数为-ax map-ont，这个根据自己数据类型看minimap2的介绍，不过多赘述。

在得到比对的bam之后，我们使用clair3来CALL变异。clair3的包比较复杂，推荐使用容器运行。这里我们使用singualrity：

singularity execclair3_latest.sif   /opt/bin/run_clair3.sh  \
--bam_fn=$species.sort.bam \
--ref_fn=$ref.fasta --threads=20 \
 --platform="hifi" \
 --model_path="/opt/models/hifi_revio"  \
--output=./  --include_all_ctgs

如果数据是ONT的，需要调整 --platform参数和模型，具体参考clair3官方的说明。这步会得到一个 merge_output.vcf.gz结果文件，就是我们需要的含有变异信息的文件。这里，我们需要对低质量的SNP位点进行过滤：

计算每个窗口的测序深度：

bedtools coverage -a windows.bed -b $species.sort.bam > $species.depth.txt

提取大于20X且小于40X的BED：

cat $species.depth.txt | perl -F'\t' -alne 'print if $F[3] > 20 && $F[3] < 40' | awk '{OFS="\t";print $1,$2,$3}' > target.bed

对clair3 call出来的VCF文件进行过滤，过滤杂合位点，得到双等位位点

bcftools view -m2 -M2 -v snpsmerge_output.vcf.gz -f PASS | \
bcftools filter -i 'GT="het"' | \
bedtools intersect -a stdin -b target.bed -wa -header | \
bgzip >targets.vcf.gz
tabix -p vcf targets.vcf.gz

提取共识序列

bcftools consensus -I -f ref.fa  targets.vcf.gz >consensus.fasta

屏蔽掉没有call出变异的区域

bedtools complement -i targets.bed -g ref.gfile | bgzip >exclude.bed.gz
bedtools maskfasta -fi consensus.fasta -fo maskedconsensus.fasta -bed exclude.bed.gz
samtools faidx maskedconsensus.fasta

将fasta序列转为PSMC所需要的格式：

  ~/bin/psmc/utils/fq2psmcfa maskedconsensus.fasta > maskedconsensus.psmcfa

最后，创建split文件用于进行多次模拟

  ~/bin/psmc/utils/splitfa maskedconsensus.psmcfa>/SPLIT_maskedconsensus.psmcfa

后续用该psmcfa文件进行PSMC分析画图即可。