1 软件准备
bwa samtools bcftools picard psmc
2 数据准备
二倍体物种（该分析是结合杂合度来分析的）！！
Ref_genome.fa: 组装好的核基因参考序列（个体A）；
A1_R1.fa.gz,A1_R2.fa.gz: 二代重测序clean_reads（个体B）【注：必须两个个体A与B】
3 PSMC分析
##1 index reference by bwa对参考建立索引
bwa index -a is Ref_genome.fa
##2 mapping 将readsmapping到参考序列
for i in *_R1.fq.gz;
do bwa mem -t 8 -R "@RG\tID:${i%_R1.fq.gz}\tSM:${i%_R1.fq.gz}\tPL:ILLUMINA" 202001_sin_ref.fa ${i%_R1.fq.gz}_R1.fq.gz ${i%_R1.fq.gz}_R2.fq.gz > ${i%_R1.fq.gz}.sam;
# -t 8 ，threads线程数设置为8
##3 将SAM格式转换为BAM
for i in *.sam;
do samtools view -@ 8 -bS ${i%.sam}.sam >${i%.sam}.bam;
# -@ 8 线程数设置为8，下同
##4 对BAM格式进行排序
for i in *.bam;
do samtools sort -@ 8 ${i%.bam}.bam ${i%.bam}.sort;
##5 移除重复reads
for i in *.sort.bam;
do samtools rmdup ${i%.sort.bam}.sort.bam ${i%.sort.bam}.sort.rmdup.bam;
##6 建立索引
for i in *.sort.rmdup.bam;
do samtools index ${i%.sort.rmdup.bam}.sort.rmdup.bam ${i%.sort.rmdup.bam}.sort.rmdup.bam.bai;
##7 生成psmc的输入文件
/usr/bin/samtools mpileup -C50 -uf 202001_sin_ref.fa S.sin.sort.rmdup.bam | /usr/bin/bcftools view -c - | vcfutils.pl vcf2fq -d 30 -D 200 | gzip > output.fq.gz
# -C50 用于降低比对质量的系数，如果reads中含有过多的错配,默认为0，samtools使用书推荐50，设不设定对i结果影响还挺大,网上教程大多设置为50。
##8 运行PSMC分析
../psmc-master/utils/fq2psmcfa -q20 output.fq.gz > output.psmcfa
../psmc-master/psmc -N25 -t15 -r5 -p "4+25*2+4+6" -p out -o output.psmc output.psmcfa
../psmc-master/utils/psmc2history.pl output.psmc |../psmc-master/utils/history2ms.pl > output.ms-cmd.sh
##9 作图
perl ../psmc-master/utils/psmc_plot.pl -u 2e-09 -g 1 out_plot output.psmc
#-u 物种碱基替换速率
#-g 生活史中的世代时间，比如人设为为25，一年生草本设置为1，世代设置越长，估计的有效群体越大。
##10 结果展示
将会生成一个.eps文件，用AI或PS打开即可。
参考来源：
PSMC软件分析群体历史有效群体大小步骤 https://blog.csdn.net/zaprily/article/details/108764219
PSMC分析流程 http://www.wuchangsong.com/?p=555
mpileup命令简介 https://blog.csdn.net/u013553061/article/details/53293302
