psmc报错及参数问题 2020-09-07

population size history from a diploid sequence using the Pairwise Sequentially Markovian Coalescent (PSMC) model
https://github.com/lh3/psmc

#!/usr/bin/perl -w

=head1 Name

        perl PSMC.pipe.dedup.H.pl -- script for population size evaluation using bwa+samtools and PSMC

=head1 Description

        The script using bwa+samtools and PSMC to calculate the population size:

=head1 Usage

        perl PSMC.pipe.dedup.H.pl <reads.list> <referance>  [-option]
        reads.list         PE reads or Single reads pathway list.
        reference          reference genomics
        --outdir <dir>     output directory, default=./
        --quality <int>    the min average quality of reads [20]
        --score <int>      the min score of align of bwa [20]
        --u FLOAT          absolute mutation rate per nucleotide [3.51e-9]
        --p STR            pattern of parameters [4+10*1+20*2+4+6]
        --gmin INT         minimum number of years per generation [1]
        --gmax INT         maximum number of years per generation [25]
        --sp <str>         set the name of spcies [animal]
        --Y FLOAT          maximum popsize, 0 for auto [0]
        --I                the input is in the Illumina 1.3+ FASTQ-like format
        --step <int>       step as follow, default=123456
                           Step 1: run bwt-builder at 00.index/
                           Step 2: run aln at 01.bwa/
                           Step 3: run map at 01.bwa/
                           Step 4: dedup and merge at 01.bwa/
                           Step 5: call snp at 01.bwa/
                           Step 6: PSMC at 02.psmc/
        --queue <str>  set the queue ,default no
        --pro_code <str>  set the project code ,default no

=head1 Example

        nohup perl PSMC.pipe.dedup.H.pl all.reads.lst reference.fa &

00.index
01.bwa
02.psmc

报错

Use of uninitialized value in division (/) at /zfssz3/NASCT_BACKUP/MS_PMO2017/lijia1/ifs1/bin/PSMC_pipe/PSMC_kit/psmc/utils/psmc_plot.pl line 119, <> line 233209.
Illegal division by zero at /zfssz3/NASCT_BACKUP/MS_PMO2017/lijia1/ifs1/bin/PSMC_pipe/PSMC_kit/psmc/utils/psmc_plot.pl line 119, <> line 233209.
Error: /undefinedfilename in (des.his_g1_u_3.51e-09.eps)
Operand stack:

Execution stack:
   %interp_exit   .runexec2   --nostringval--   --nostringval--   --nostringval--   2   %stopped_push   --nostringval--   --nostringval--   --nostringval--   false   1   %stopped_push
Dictionary stack:
   --dict:1151/1684(ro)(G)--   --dict:0/20(G)--   --dict:70/200(L)--
Current allocation mode is local
Last OS error: 2
GPL Ghostscript 8.70: Unrecoverable error, exit code 1
Use of uninitialized value in division (/) at /zfssz3/NASCT_BACKUP/MS_PMO2017/lijia1/ifs1/bin/PSMC_pipe/PSMC_kit/psmc/utils/psmc_plot.pl line 119, <> line 233209.
Use of uninitialized value in division (/) at /zfssz3/NASCT_BACKUP/MS_PMO2017/lijia1/ifs1/bin/PSMC_pipe/PSMC_kit/psmc/utils/psmc_plot.pl line 119, <> line 233209.
Illegal division by zero at /zfssz3/NASCT_BACKUP/MS_PMO2017/lijia1/ifs1/bin/PSMC_pipe/PSMC_kit/psmc/utils/psmc_plot.pl line 119, <> line 233209.
Error: /undefinedfilename in (des.his_g2_u_7.02e-09.eps)

寻错思路历程：

1. Use of uninitialized value in division (/) at /zfssz3/NASCT_BACKUP/MS_PMO2017/lijia1/ifs1/bin/PSMC_pipe/PSMC_kit/psmc/utils/psmc_plot.pl line 119, <> line 233209.
   "大神说，这里是因为line 119对应的有个除数是0"
   于是去排查psmc_plot.pl line 119后发现很多生成文件不对，包括des.SNP.psmc

MM      Version: 0.6.5-r67
MM      pattern:4+10*1+20*2+4+6, n:63, n_free_lambdas:33
MM      n_iterations:30, skip:1, max_t:15, theta/rho:5
MM      is_decoding:0
MM      n_seqs:2, sum_L:1101790, sum_n:0
RD      0
LK      0.000000
QD      0.000000 -> 0.000000
RI      -nan
TR      -0.000000       -0.000000
MT      15.000000
MM      C_pi: 1.000000, n_recomb: -0.000000
RS      0       0.000000        1.000000        -nan    -nan    0.000000
RS      1       0.008290        1.000000        -nan    -nan    0.000000
RS      2       0.017266        1.000000        -nan    -nan    0.000000
RS      3       0.026987        1.000000        -nan    -nan    0.000000
RS      4       0.037514        1.000000        -nan    -nan    0.000000

2. 于是去逐步查看，发现01.bwa/00.get_sai.sh.5104.qsub的第一步就没跑完，磁盘满额

[bwa_aln_core] write to the disk... 0.03 sec
[bwa_aln_core] 260046848 sequences have been processed.
[bwa_aln_core] calculate SA coordinate... 205.79 sec
[bwa_aln_core] write to the disk... 0.03 sec
[bwa_aln_core] 260308992 sequences have been processed.
[bwa_aln_core] calculate SA coordinate... 202.33 sec
[bwa_aln_core] write to the disk... [fwrite] Disk quota exceeded

第二步也没跑完01.bwa/01.run_mapping.sh.28381.qsub/run_map_001.sh

Java HotSpot(TM) 64-Bit Server VM warning: ignoring option MaxPermSize=2g; support was removed in 8.0
[Mon Sep 07 15:38:17 CST 2020] picard.sam.SortSam INPUT=/ldfssz1/MS_OP/USER/lifan/05.Desis/00.genome/09.Evolution/04.population/./01.bwa/L500/.sam OUTPUT=/ldfssz1/MS_OP/USER/lifan/05.Desis/00.genome/09.Evolution/04.population/./01.bwa/L500/.sort.bam SORT_ORDER=coordinate TMP_DIR=[/ldfssz1/MS_OP/USER/lifan/05.Desis/00.genome/09.Evolution/04.population/./01.bwa/L500/TMP1] VALIDATION_STRINGENCY=SILENT MAX_RECORDS_IN_RAM=1000000    VERBOSITY=INFO QUIET=false COMPRESSION_LEVEL=5 CREATE_INDEX=false CREATE_MD5_FILE=false
[Mon Sep 07 15:38:17 CST 2020] Executing as lifan@cngb-compute-e06-2.cngb.sz.hpc on Linux 2.6.32-696.30.1.el6.x86_64 amd64; Java HotSpot(TM) 64-Bit Server VM 1.8.0_121-b13; Picard version: 1.117(107391d3f3e72b31589868c250262ca79659f577_1405353489) JdkDeflater
[Mon Sep 07 15:39:29 CST 2020] picard.sam.SortSam done. Elapsed time: 1.20 minutes.
Runtime.totalMemory()=4584898560
To get help, see http://picard.sourceforge.net/index.shtml#GettingHelp
Exception in thread "main" htsjdk.samtools.util.RuntimeIOException: java.io.IOException: 超出磁盘限额
        at htsjdk.samtools.util.SortingCollection.spillToDisk(SortingCollection.java:245)
        at htsjdk.samtools.util.SortingCollection.add(SortingCollection.java:165)
        at htsjdk.samtools.SAMFileWriterImpl.addAlignment(SAMFileWriterImpl.java:179)
        at picard.sam.SortSam.doWork(SortSam.java:75)
        at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:183)
        at picard.cmdline.CommandLineProgram.instanceMainWithExit(CommandLineProgram.java:124)
        at picard.sam.SortSam.main(SortSam.java:61)
Caused by: java.io.IOException: 超出磁盘限额

解决方案：清理磁盘，重新跑

以下是参数

per-generation mutation rate -u
the generation time in years -g

$ /PSMC_pipe/PSMC_kit/psmc/utils/psmc_plot.pl

Usage:   psmc_plot.pl [options] <out.prefix> <in.psmc>

Options: -u FLOAT   absolute mutation rate per nucleotide [2.5e-08]
         -s INT     skip used in data preparation [100]
         -X FLOAT   maximum generations, 0 for auto [0]
         -x FLOAT   minimum generations, 0 for auto [10000]
         -Y FLOAT   maximum popsize, 0 for auto [0]
         -m INT     minimum number of iteration [5]
         -n INT     take n-th iteration (suppress GOF) [20]
         -M titles  multiline mode [null]
         -f STR     font for title, labels and tics [Helvetica,16]
         -g INT     number of years per generation [25]
         -w INT     line width [4]
         -P STR     position of the keys [right top]
         -T STR     figure title [null]
         -N FLOAT   false negative rate [0]
         -S         no scaling
         -L         show the last bin
         -p         convert to PDF (with epstopdf)
         -R         do not remove temporary files
         -G         plot grid