Population split and merge (-es, -ej) indexing

GoogleCodeExporter commented 9 years ago

What steps will reproduce the problem?

./macs 2 100000 -I 2 1 1 -es 0.01 1 0.05 -ej 0.05 3 2

What is the expected output? What do you see instead?

This run terminates because of an indexing error:

------------------------------
INPUT: Sample size is now 2
INPUT: Seq length is now 100000
INPUT: Setting chr sampled for pop 1 to 1
INPUT: Setting chr sampled for pop 2 to 1
INPUT: Global migration rate to 0
INPUT: At time 0.01: Population 1 splits at proportion 0.05
INPUT: At time 0.05: WARNING: The pop IDs used in pop join is greater than the 
number specified in -I.  You must have a split event before this join event.
Population 3 will merge with 2
SEED:   1357694700
Debugging: 0

Graph Builder begin
Simulator caught exception with message:
Infinite coalescent time. No migration
Program is complete
----------------------

What version of the product are you using? On what operating system?
Mac OS 10.6.8, macs-0.4e

Please provide any additional information below.

In the command README that appears when you type ./macs into the terminal, the 
-es command is described as follows:

-es <t> <i> <p> (Split two populations.  At time t, a proportion p of 
chromosomes from pop i will migrate to a population i+1.

This is different from how -es works in Hudson's MS: at time t, if there exist 
n total populations and you split population i into two parts, the new 
population index is n+1, not i+1. I wrote the above command line assuming the 
discrepancy was a typo and that the indexing works the same as in MS, but it 
didn't seem to work that way. However, the following command line assumes that 
-es is correctly described in the macs manual, and it doesn't work either:

-------------------
COMMAND:        ./macs 2 100000 -I 2 1 1 -es 0.01 2 0.05 -ej 0.05 3 1
INPUT: Sample size is now 2
INPUT: Seq length is now 100000
INPUT: Setting chr sampled for pop 1 to 1
INPUT: Setting chr sampled for pop 2 to 1
INPUT: Global migration rate to 0
INPUT: At time 0.01: Population 2 splits at proportion 0.05
INPUT: At time 0.05: WARNING: The pop IDs used in pop join is greater than the 
number specified in -I.  You must have a split event before this join event.
Population 3 will merge with 1
SEED:   1357695229
Debugging: 0

Graph Builder begin
Simulator caught exception with message:
Infinite coalescent time. No migration
Program is complete
-------------------------

Thanks so much for your help with this issue!
-Kelley Harris

Original issue reported on code.google.com by harris.k...@gmail.com on 9 Jan 2013 at 1:47

GoogleCodeExporter commented 9 years ago

[deleted comment]

GoogleCodeExporter commented 9 years ago

You didn't specify a migration rate at the end of the -I flag so it defaulted 
to zero.  Try something like:

./macs 2 100000 -I 2 1 1 1 -es 0.01 1 0.05 -ej 0.05 3 2

Original comment by gche...@gmail.com on 9 Jan 2013 at 6:24

Changed state: Started

GoogleCodeExporter commented 9 years ago

Thanks--I accidentally left out a second population merge event when trying to 
reduce a complicated command line to this one! I think my code is working now, 
but is there documentation about the following error message that was giving me 
trouble?

Graph Builder begin
Attempting to add edge of population 10
Simulator caught exception with message:
edge population is too large

Sorry for the confusion...

Original comment by harris.k...@gmail.com on 9 Jan 2013 at 8:42

GoogleCodeExporter commented 9 years ago

Can you give me the full command line that caused the last error message?

Original comment by gche...@gmail.com on 9 Jan 2013 at 8:45

GoogleCodeExporter commented 9 years ago

Sure, this is a command line that makes it happen:

./macs 2 1e7 -t 0.001 -r 0.0004 -h 1 -I 2 1 1 -es 0.0115 1 0.0466 -es 0.0115 2 
0.0466 -ej 0.0115001 3 2 -ej 0.0115001 4 1 -es 0.0478 1 0.961 -ej 0.0478001 5 2 
-ej  0.431 1 2 

This command doesn't always give the same output. Sometimes it triggers 'edge 
population too large' like this:

-------
COMMAND:    ./macs 2 1e7 -t 0.001 -r 0.0004 -h 1 -I 2 1 1 -es 0.0115 1 0.0466 -es 
0.0115 2 0.0466 -ej 0.0115001 3 2 -ej 0.0115001 4 1 -es 0.0478 1 0.961 -ej 
0.0478001 5 2 -ej 0.431 1 2
INPUT: Sample size is now 2
INPUT: Seq length is now 1e+07
INPUT: Scaled mutation rate is now 10000
INPUT: Scaled recombination rate is now 4000
INPUT: Base pairs to track is 1
INPUT: Setting chr sampled for pop 1 to 1
INPUT: Setting chr sampled for pop 2 to 1
INPUT: Global migration rate to 0
INPUT: At time 0.0115: Population 1 splits at proportion 0.0466
INPUT: At time 0.0115: Population 2 splits at proportion 0.0466
INPUT: At time 0.0115001: WARNING: The pop IDs used in pop join is greater than 
the number specified in -I.  You must have a split event before this join event.
Population 3 will merge with 2
INPUT: At time 0.0115001: WARNING: The pop IDs used in pop join is greater than 
the number specified in -I.  You must have a split event before this join event.
Population 4 will merge with 1
INPUT: At time 0.0478: Population 1 splits at proportion 0.961
INPUT: At time 0.0478001: WARNING: The pop IDs used in pop join is greater than 
the number specified in -I.  You must have a split event before this join event.
Population 5 will merge with 2
INPUT: At time 0.431: Population 1 will merge with 2
SEED:   1357764719
Debugging: 0

Graph Builder begin
Attempting to add edge of population 4
Simulator caught exception with message:
edge population is too large
Program is complete
---------------

Other times, it triggers a segmentation fault:

-------------
COMMAND:    ./macs 2 1e7 -t 0.001 -r 0.0004 -h 1 -I 2 1 1 -es 0.0115 1 0.0466 -es 
0.0115 2 0.0466 -ej 0.0115001 3 2 -ej 0.0115001 4 1 -es 0.0478 1 0.961 -ej 
0.0478001 5 2 -ej 0.431 1 2
INPUT: Sample size is now 2
INPUT: Seq length is now 1e+07
INPUT: Scaled mutation rate is now 10000
INPUT: Scaled recombination rate is now 4000
INPUT: Base pairs to track is 1
INPUT: Setting chr sampled for pop 1 to 1
INPUT: Setting chr sampled for pop 2 to 1
INPUT: Global migration rate to 0
INPUT: At time 0.0115: Population 1 splits at proportion 0.0466
INPUT: At time 0.0115: Population 2 splits at proportion 0.0466
INPUT: At time 0.0115001: WARNING: The pop IDs used in pop join is greater than 
the number specified in -I.  You must have a split event before this join event.
Population 3 will merge with 2
INPUT: At time 0.0115001: WARNING: The pop IDs used in pop join is greater than 
the number specified in -I.  You must have a split event before this join event.
Population 4 will merge with 1
INPUT: At time 0.0478: Population 1 splits at proportion 0.961
INPUT: At time 0.0478001: WARNING: The pop IDs used in pop join is greater than 
the number specified in -I.  You must have a split event before this join event.
Population 5 will merge with 2
INPUT: At time 0.431: Population 1 will merge with 2
SEED:   1357764706
Debugging: 0

Graph Builder begin
Time for build prior tree: 0 seconds
Tree:0,pos:0,len:3.05332,TMRCA:1.52666,ARG:,len:3.05332,TMRCA:1.52666
SITE:   0    1.28106006e-05 10
SITE:   1    1.28598812e-05 10
iHistoryMax: 0
Tree:1,pos:3.39317e-05,len:3.76017,TMRCA:1.88008,ARG:,len:3.76017,TMRCA:1.88008
SITE:   2    4.88701756e-05 01
SITE:   3    5.01489262e-05 01
iHistoryMax: 1
Tree:2,pos:8.2074e-05,len:3.23806,TMRCA:1.61903,ARG:,len:3.23806,TMRCA:1.61903
SITE:   4    0.000130348503 01
SITE:   5    0.000139158248 10
iHistoryMax: 2
Tree:3,pos:0.000175768,len:3.23806,TMRCA:1.61903,ARG:,len:3.23806,TMRCA:1.61903
iHistoryMax: 3
Tree:4,pos:0.000185621,len:0.656081,TMRCA:0.328041,ARG:,len:0.656081,TMRCA:0.328
041
SITE:   6    0.000315345344 01
SITE:   7    0.000471867992 10
SITE:   8     0.00047927398 10
iHistoryMax: 4
Tree:5,pos:0.00056717,len:0.985478,TMRCA:0.492739,ARG:,len:0.985478,TMRCA:0.4927
39
SITE:   9    0.000654781904 01
SITE:   10   0.000693792245 10
SITE:   11   0.000859340054 10
iHistoryMax: 5
Tree:6,pos:0.00088489,len:2.77008,TMRCA:1.38504,ARG:,len:2.77008,TMRCA:1.38504
SITE:   12   0.000887199999 01
iHistoryMax: 6
Tree:7,pos:0.000887541,len:1.49568,TMRCA:0.74784,ARG:,len:1.49568,TMRCA:0.74784
SITE:   13   0.000952663631 10
iHistoryMax: 7
Tree:8,pos:0.0010439,len:1.49568,TMRCA:0.74784,ARG:,len:1.49568,TMRCA:0.74784
SITE:   14    0.00105116557 10
SITE:   15     0.0010629279 01
SITE:   16    0.00108205987 01
SITE:   17    0.00108370083 10
SITE:   18    0.00108418391 01
SITE:   19    0.00109058072 01
SITE:   20    0.00110299058 10
SITE:   21    0.00118907367 01
SITE:   22    0.00125116708 10
iHistoryMax: 8
Tree:9,pos:0.00128794,len:0.971942,TMRCA:0.485971,ARG:,len:0.971942,TMRCA:0.4859
71
SITE:   23    0.00141130235 01
SITE:   24    0.00153928963 01
SITE:   25    0.00166629385 10
SITE:   26    0.00174952866 01
SITE:   27    0.00176435511 10
SITE:   28    0.00191876966 01
SITE:   29    0.00192085395 01
SITE:   30    0.00195752302 01
SITE:   31    0.00196465375 01
SITE:   32    0.00208921013 10
iHistoryMax: 9
Tree:10,pos:0.00216419,len:0.971942,TMRCA:0.485971,ARG:,len:0.971942,TMRCA:0.485
971
SITE:   33     0.0021661666 10
Segmentation fault

----------------------

Original comment by harris.k...@gmail.com on 9 Jan 2013 at 9:06

GoogleCodeExporter commented 9 years ago

Please download the latest version 0.5b or check it out from svn.  It should 
provide friendlier error messages now.

Original comment by gche...@gmail.com on 10 Jan 2013 at 10:51

GoogleCodeExporter commented 9 years ago

I have not heard updates and will close this issue now.

Original comment by gche...@gmail.com on 22 Feb 2013 at 2:30

Changed state: Fixed

gchen98 / macs

Population split and merge (-es, -ej) indexing #11