10群体遗传进化详解
三.PSMC
3.1 公司单次 psmc.cmds
perl /share/nas1/yuj/pipline/genetic_diversity_pip/v1.0/script/06psmc/vcf2psmcGroup.pl -i /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/vcf_filter/samples.pop.snp.recode.vcf -o /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis -g /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/config_dir/group.list
## 会生成两个sh文件，一个分析，一个画图
perl /share/nas1/yuj/pipline/genetic_diversity_pip/v1.0/script/06psmc/cnvfmt_pdf2png.pl -i /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis
# 转换图片格式
3.1.1 分析work1.sh
单次
psmc -N25 -t15 -r5 -p "4+25*2+4+6" -o /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/BX.psmc /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/BX.fasta
psmc -N25 -t15 -r5 -p "4+25*2+4+6" -o /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/HN.psmc /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/HN.fasta
psmc -N25 -t15 -r5 -p "4+25*2+4+6" -o /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/XKH.psmc /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/XKH.fasta
多次
1.运行分析
for i in {1..100};do nohup psmc -b -N25 -t15 -r5 -p "4+25*2+4+6" -o /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/BX.round${i}.psmc /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/BX.fasta & done
for i in {1..100};do nohup psmc -b -N25 -t15 -r5 -p "4+25*2+4+6" -o /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/HN.round${i}.psmc /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/HN.fasta & done
for i in {1..100};do nohup psmc -b -N25 -t15 -r5 -p "4+25*2+4+6" -o /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/XKH.round${i}.psmc /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/XKH.fasta & done
直接运行吧
#-b参数表示bootstrap抽样做psmc，重复100次
#-t：maximum 2N0 coalescent time [15]; 2N0的最大世代时间，代表TMRCA的上限
#-r：initial theta/rho ratio [4]; 起始θ/ρ率，默认4。
#-N：maximum number of iterations [30]；迭代的最大次数，默认30。
#-p：代表PSMC有效群体大小变化的图随时间变化划分的时间阶段单位数量。
例如默认的-p “4+25*2+4+6”表示从古至今依次经历了一个4单位+25个2单位+一个4单位+一个6单位的时期。
数字越大需要的计算资源越多，可根据经验调整，也可查找近缘物种的值用作参考。
2.合并结果
cat BX.psmc BX.round*.psmc > BX.combine.psmc
cat HN.psmc HN.round*.psmc > HN.combine.psmc
cat XKH.psmc XKH.round*.psmc > XKH.combine.psmc
3.1.2 绘图work2.sh
单次
psmc_plot.pl -R -u 1.91e-09 -g 1 -p /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/BX /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/BX.psmc
psmc_plot.pl -R -u 1.91e-09 -g 1 -p /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/HN /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/HN.psmc
psmc_plot.pl -R -u 1.91e-09 -g 1 -p /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/XKH /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/XKH.psmc
psmc_plot.pl -R -u 1.91e-09 -g 1 -M "BX,HN,XKH" -p /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/all.combined.psmc /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/BX.psmc /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/HN.psmc /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/XKH.psmc # 多分组放一张图
多次
psmc_plot.pl -R -u 1.91e-09 -g 1 -p /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/BX /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/BX.combine.psmc
psmc_plot.pl -R -u 1.91e-09 -g 1 -p /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/HN /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/HN.combine.psmc
psmc_plot.pl -R -u 1.91e-09 -g 1 -p /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/XKH /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/XKH.combine.psmc
psmc_plot.pl -R -u 1.91e-09 -g 1 -M "BX,HN,XKH" -p /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/all.combined.psmc /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/BX.combine.psmc /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/HN.combine.psmc /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis/XKH.combine.psmc
放入work2.sh，直接sh就行
#-g：number of years per generation，该物种的平均世代时间（即繁殖一代的时间）
#-u：absolute mutation rate per nucleotide，该物种碱基替换率
3.1.3 转换
perl /share/nas1/yuj/pipline/genetic_diversity_pip/v1.0/script/06psmc/cnvfmt_pdf2png.pl -i /share/nas2/yuj/project/2023/ddRAD-seq/GP-20230823-6676_20231023/analysis/psmc_analysis
多次和单次一样
3.2 运行bootstrap的PSMC【推荐】
如果要运行多次抽样的PSMC，则在上面的格式转换后，先分割长序列，再用for循环语句运行。
1. 分割长序列
splitfa diploid.psmcfa > split.psmcfa
若需要做bootstrap需先把太长序列分割
2. bootstrap抽样运行PSMC
for i in {1..100};do nohup psmc -b -p "6+25*2+6+8" -t15 -N30 -r4 -o round${i}.psmc split.psmcfa &
#-b参数表示bootstrap抽样做psmc，重复100次。估计服务器运算能力，可以10个重复作为一批。
3. 合并结果
cat diploid.psmc round*.psmc >combine.psmc #psmc结果合并
3.3 画图
1. 运行
psmc_plot.pl -p -g 1 -u 4.79e-9 combine combine.psmc
- 画图，生成combine.eps和combine.pdf
2. 参数
-p：指用epstopdf把esp转化成pdf文件，相当于运行完多加一个`epstopdf combined.eps`一步。
-g：number of years per generation [25]，指定该物种的平均世代时间（即繁殖一代的时间），比如人的默认为25年，该处值为-g 25；许多草本植物为1年。
-u：absolute mutation rate per nucleotide [2.5e-08]，指定该物种碱基替换率，单位是subst./syn. site/year，比如默认水稻的2.5e-08，文献里提到的蕨类4.79e-9(Barker 2009)。参考已有文献对研究类群的替换率的推断，也可由进化树的枝长除以r8s评估该物种的分歧时间得到。