If EPS is <10^(-3) or <10^(-6), then pix is NAN, because lineages coalesce within the interval with probability (almost) 1.

I tested out running ngsPSMC using the instructions outlined in the tutorial.

I first ran this script on a single .bam file to generate an input file

./angsd -i S13_2949_trim_cleaned.bam -ref canFam3_withY.fa -out S13_2949_psmcinput -dopsmc 1 -gl 1 -minq 30 -minmapq 30

Looking at the bottom of my .arg file it appears it worked and generated the following output files

abcPSMC.cpp:
-writePSMC 1
1: (still beta, not really working)
->"S13_2949_psmcinput.arg"
->"S13_2949_psmcinput.psmc.gz"
->"S13_2949_psmcinput.psmc.pos.gz"
->"S13_2949_psmcinput.psmc.idx"

    [ALL done] cpu-time used =  4683.03 sec
    [ALL done] walltime used =  4776.00 sec

I then tried running ngsPSMC and am getting the following error:

./ngsPSMC S13_2949_psmcinput.psmc.idx -p "1*4+25*2+1*4+1*6" -dospline 0 -nthreads 8 -nIter 20 -init 1  -theta 0.000233095 -rho 0.005357 > S13_2949.psmc.ml.txt

-> we are in file: main_psmc.cpp function: main_psmc line:187
    -> Version of fname: '/home/sophiepq/GrayFoxWGS/GrayFox1/angsd/ngsPSMC/S13_2949_psmcinput.psmc.idx' is:1
    -> Assuming .psmc.gz file: /home/sophiepq/GrayFoxWGS/GrayFox1/angsd/ngsPSMC/S13_2949_psmcinput.psmc.gz
    -> Assuming .psmc.pos.gz: /home/sophiepq/GrayFoxWGS/GrayFox1/angsd/ngsPSMC/S13_2949_psmcinput.psmc.pos.gz
    -> Problem opening file: '-dospline'

I don't really know what this error means or how to address it. I believe -dospline 0 or -dospline 1 are requirements to run program. I'm not sure whether there is another input file it's looking for that I am missing here.

fpsmc.cpp::306-309
rho and theta are set in the code to particular values.

Hi,
we had problems compiling ngsPSMC on one of our current systems. This was related both to the code in ngsPSMC and htslib. Below is the workaround (or "fix") we tried. For what it's worth.
Cheers
/Johan

(Used g++ (GCC) 4.8.5 20150623 (Red Hat 4.8.5-44), running CentOS Linux release 7.9.2009)

1. Get htslib from the developers release on GitHub

    $ cd /path/to/folder
    $ wget https://github.com/samtools/htslib/releases/download/1.15/htslib-1.15.tar.bz2
    $ tar -xf htslib-1.15.tar.bz2
    $ rm htslib-1.15.tar.bz2
   

2. Make changes in htslib linking, and build (need to generate libhts.a)

    $ cd htslib-1.15
    $ sed -i '/htslib_default_libs/ s/$/ -lssl -lcrypto/' Makefile
    $ ./configure
    $ make
   

3. Clone ngsPSMC from GitHub

    $ cd /path/to/folder
    $ git clone https://github.com/ANGSD/ngsPSMC
    $ cd ngsPSMC
   

4. Make changes in the code (linking, and add include statement), and build (expect and ignore Warnings).

    $ sed -i '/o ngsPSMC/ s/$/ -lssl -lcrypto/' Makefile
    $ sed -i '/header.h/ s/$/\n#include <stdio.h>/' header.cpp
    $ sed -i '/cstring/ s/$/\n#include <stdio.h>/' header.h
    $ make HTSSRC=../htslib-1.15/
   

Final executable file is /path/to/folder/ngsPSMC/ngsPSMC

Hi,

I have cloned the ngsPSMC into my computer and I have ran a test with the chr38 of a dog. As recommended in README.md, I ran:

$ ../../angsd/angsd -i dog.qatar.chr38.bam -dopsmc 1 -out chr38 -gl 1 -minq 20 -minmapq 30
with output filenames chr38.arg; chr38.psmc.gz, chr38.psmc.poz.gz, chr38.psmc.idx

After running ngsPSMC in this file as recommended in README.md:
../ngsPSMC chr38.psmc.idx -p "14+252+14+16"
-dospline 0 -nthreads 8 -nIter 20 -init 1
-theta 0.000233095 -rho 0.005357 > chr38.psmc.ml.txt

The file seems correct:

$ more chr38.psmc.ml.txt MM ../../ngsPSMC/ngsPSMC chr38.psmc.psmc.idx -p 1*4+25*2+1*4+1*6 -dospline 0 -nthreads 8 -nIt er 20 -init 1 -theta 0.000233095 -rho 0.005357 MM -> ngsPSMC version: v0.02-54-gbcda080 (htslib: 1.9) build(Oct 14 2019 15:57:42) MM TR is theta rho, theta in this context is given by the persite theta, and not the per wind ow theta, this is different from the original PSMC MM filereading took: (wall(min),cpu(min)):(0.050000,0.042487) RD 0 LK -30553205.977373 QD 0.000000 -> 0.000000 RI ? TR 0.000233 0.005357 MT 23.861429 MM buildhmm(wall(min),cpu(min)):(3.450000,3.458922) tk_l:64 RS 0 0.000000 1.000000 1000000.0 1000000.0 1000000.0 RS 1 0.009086 1.000000 1000000.0 1000000.0 1000000.0 RS 2 0.018998 1.000000 1000000.0 1000000.0 1000000.0 RS 3 0.029811 1.000000 1000000.0 1000000.0 1000000.0

But there is no specifications as what to do next. I tried to decompress chr38.psmc.gz into chr38.psmc and run psmc_plot.py:

$ python3.7 psmc_plot.py -psmc ../../test_ngspsmc/testchr38/chr38.psmc
Units: mutation rate = 1.25e-08 binsize = 100.0 N0 = 10000.0 generation time = 29.0
Argument should have 4 or 5 parameters.

Am I doing something wrong or is the code still under development?

Maybe I can directly use the psmc_plot.py code of PSMC (/utils/psmc_plot.py) to run in my chr38?

Thanks

Carlos

https://github.com/ANGSD/ngsPSMC/blob/master/hmm_psmc.cpp#L238-L240