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分析
代码语言：javascript复制 ##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 "@RGtID:${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打开即可。
