PSMC估计种群大小历史变化

https://github.com/lh3/psmc

#path
export GNUPLOT_PS_DIR="/usr/share/gnuplot/gnuplot/5.2/PostScript"
export ngsdir=""
export sample=""
export ref=""
export bcftools=""
export psmc=""
export generation=""
export mutation_rate=""

#bwa&rmdup
bwa index $ref
bwa mem -t 40 $ref $ngsdir/${sample}_1.fastq.gz $ngsdir/${sample}_2.fastq.gz | samtools view -@ 40 -bS | samtools sort -@ 40 -m 1G > ${sample}.sort.bam
samtools rmdup ${sample}.sort.bam ${sample}.sort.rmdup.bam
samtools index ${sample}.sort.rmdup.bam

#bam2fq
$bcftools/bin/bcftools mpileup -Ou -I -f $ref ${sample}.sort.rmdup.bam | $bcftools/bin/bcftools call -c -Ov | $bcftools/bin/vcfutils.pl vcf2fq -d 10 -D 100 | gzip > ${sample}.output.fq.gz

#generate psmc input
$psmc/utils/fq2psmcfa -q20 ${sample}.output.fq.gz > ${sample}.psmcfa

#psmc run
$psmc/psmc -N25 -t15 -r5 -p "4+25*2+4+6" -o ${sample}.psmc ${sample}.psmcfa

#plot
$psmc/utils/psmc_plot.pl -X 10000000 -Y 25 -g $generation -u $mutation_rate -p ${sample} ${sample}.psmc

作图时g,u的设置可参考相关文献。

合并作图

#按顺序合并
cat A.psmc B.psmc > combined.psmc
#-M设置与顺序对应
$psmc/utils/psmc_plot.pl -X 10000000 -Y 25 -g $generation -u $mutation_rate -M "A,B" -p combined combined.psmc