1 軟件準(zhǔn)備
bwa samtools bcftools?picard? psmc
2 數(shù)據(jù)準(zhǔn)備
二倍體物種（該分析是結(jié)合雜合度來(lái)分析的）??！
Ref_genome.fa: 組裝好的核基因參考序列（個(gè)體A）；
A1_R1.fa.gz,A1_R2.fa.gz: 二代重測(cè)序clean_reads（個(gè)體B）【注：必須兩個(gè)個(gè)體A與B】
3 PSMC分析
##1 index reference by bwa對(duì)參考建立索引
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線(xiàn)程數(shù)設(shè)置為8
##3 將SAM格式轉(zhuǎn)換為BAM? ?
for i in *.sam;
do samtools view? -@ 8 -bS ${i%.sam}.sam >${i%.sam}.bam;
# -@ 8 線(xiàn)程數(shù)設(shè)置為8，下同
##4 對(duì)BAM格式進(jìn)行排序
for i in *.bam;
do samtools sort? -@ 8 ${i%.bam}.bam ${i%.bam}.sort;
##5 移除重復(fù)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 用于降低比對(duì)質(zhì)量的系數(shù)，如果reads中含有過(guò)多的錯(cuò)配,默認(rèn)為0，samtools使用書(shū)推薦50，設(shè)不設(shè)定對(duì)i結(jié)果影響還挺大,網(wǎng)上教程大多設(shè)置為50。
##8? 運(yùn)行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 生活史中的世代時(shí)間，比如人設(shè)為為25，一年生草本設(shè)置為1，世代設(shè)置越長(zhǎng)，估計(jì)的有效群體越大。
##10 結(jié)果展示
將會(huì)生成一個(gè).eps文件，用AI或PS打開(kāi)即可。
參考來(lái)源：
PSMC軟件分析群體歷史有效群體大小步驟? https://blog.csdn.net/zaprily/article/details/108764219
PSMC分析流程??http://www.wuchangsong.com/?p=555
mpileup命令簡(jiǎn)介?https://blog.csdn.net/u013553061/article/details/53293302
