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	ID:${i%_R1.fq.gz}	SM:${i%_R1.fq.gz}	PL: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 结果展示