Negative branch lengths

GoogleCodeExporter commented 9 years ago

What steps will reproduce the problem?
1./macs 100 3000000 -t 1.0 -T -r .0002 -h 1 >test100Const.tree
2.
3.

What is the expected output? What do you see instead?
I expect to see a sequence of genealogies simulated from SMC'. I see that 
sometimes from one tree to the next one, all the branch lengths change by some 
very small number, so I added a tolerance on my programs. What I am having now 
is some negative branch lengths
?
What version of the product are you using? On what operating system?
I downloaded macs on January this year and my computer is a Mac OS X 10.7.5

Please provide any additional information below.

Original issue reported on code.google.com by julia.pa...@gmail.com on 19 May 2014 at 7:38

GoogleCodeExporter commented 9 years ago

Hi there

Can you please add a -s flag with a desired seed and copy and paste the first 
100 or so lines, just to make sure there's no platform specific differences in 
the RNG?

Once I can reproduce your output, point me to the line numbers you see the 
issue at.

Gary

Original comment by gche...@gmail.com on 19 May 2014 at 8:12

Changed state: Accepted

GoogleCodeExporter commented 9 years ago

Hi Gary,

Thank you so much for looking at this and I am sorry It took me long to reply. 
I couldn't replicate the issue until now and I saved the file this time. The 
seed is 1403270054

>>./macs 100 1000000 -t 1 -T -r .0002 -h 1 -s 1403270054

The header of the output file is:
COMMAND:        ./macs 100 1000000 -t 1 -T -r .0002 -h 1
SEED:   1403270054
NEWICK_TREE:    
[43](68:2.48263,((((29:0.0635055,((((80:0.000159123,23:0.000159123):0.0189092,((
84:0.00115495,89:0.00115495):0.00215123,8:0.00330618):0.0157622):0.0011627,((64:
0.0052641,3:0.0052641):0.0075611,39:0.0128252):0.00740586):0.00244112,(5:0.02166
51,14:0.0216651):0.00100711):0.0408333):0.0584093,(((25:0.018254,((36:0.00140987
,72:0.00140987):0.00806156,(((71:0.0020156,77:0.0020156):0.000956306,90:0.002971
9):0.00545664,((74:0.0068498......

The first lines of output are the following:

Julias-MacBook-Pro-4:macs julie_palacios$ ./macs 100 1000000 -t  1 -T -r .0002 
-h 1 >Const100Idv_5.tree 
INPUT: Sample size is now 100
INPUT: Seq length is now 1e+06
INPUT: Scaled mutation rate is now 1e+06
INPUT: Scaled recombination rate is now 200
INPUT: Base pairs to track is 1
Debugging: 0

Graph Builder begin
Time for build prior tree: 1535 seconds
Tree:0,pos:0,len:7.87279,TMRCA:2.48263,ARG:,len:7.87279,TMRCA:2.48263
iHistoryMax: 0
Tree:1,pos:4.34071e-05,len:7.83686,TMRCA:2.48263,ARG:,len:7.83686,TMRCA:2.48263
iHistoryMax: 1
Tree:2,pos:0.000443469,len:7.89644,TMRCA:2.48263,ARG:,len:7.89644,TMRCA:2.48263
iHistoryMax: 2
Tree:3,pos:0.00045678,len:7.2346,TMRCA:2.15171,ARG:,len:7.2346,TMRCA:2.15171
iHistoryMax: 3
Tree:4,pos:0.00204571,len:5.83316,TMRCA:1.45099,ARG:,len:5.83316,TMRCA:1.45099
iHistoryMax: 4
Tree:5,pos:0.00281437,len:8.02983,TMRCA:2.54932,ARG:,len:8.02983,TMRCA:2.54932
iHistoryMax: 5
Tree:6,pos:0.00285434,len:8.02983,TMRCA:2.54932,ARG:,len:8.02983,TMRCA:2.54932
iHistoryMax: 6
Tree:7,pos:0.00323309,len:8.02983,TMRCA:2.54932,ARG:,len:8.02983,TMRCA:2.54932
iHistoryMax: 7
Tree:8,pos:0.00367939,len:8.02983,TMRCA:2.54932,ARG:,len:8.02983,TMRCA:2.54932
iHistoryMax: 8
Tree:9,pos:0.00562264,len:7.99857,TMRCA:2.54932,ARG:,len:7.99857,TMRCA:2.54932
iHistoryMax: 9
Tree:10,pos:0.00572892,len:6.51183,TMRCA:1.80595,ARG:,len:6.51183,TMRCA:1.80595
iHistoryMax: 10
Tree:11,pos:0.00595076,len:6.50821,TMRCA:1.80595,ARG:,len:6.50821,TMRCA:1.80595
iHistoryMax: 11
Tree:12,pos:0.00617452,len:6.48829,TMRCA:1.80595,ARG:,len:6.48829,TMRCA:1.80595
iHistoryMax: 12
Tree:13,pos:0.0062737,len:6.48829,TMRCA:1.80595,ARG:,len:6.48829,TMRCA:1.80595
iHistoryMax: 13
Tree:14,pos:0.0068805,len:7.43953,TMRCA:2.28157,ARG:,len:7.43953,TMRCA:2.28157
iHistoryMax: 14
Tree:15,pos:0.00693862,len:7.43953,TMRCA:2.28157,ARG:,len:7.43953,TMRCA:2.28157
iHistoryMax: 15
Tree:16,pos:0.00722873,len:7.52084,TMRCA:2.28157,ARG:,len:7.52084,TMRCA:2.28157
iHistoryMax: 16
Tree:17,pos:0.00810503,len:6.60129,TMRCA:1.8218,ARG:,len:6.60129,TMRCA:1.8218
iHistoryMax: 17
Tree:18,pos:0.0086867,len:6.62047,TMRCA:1.8218,ARG:,len:6.62047,TMRCA:1.8218
iHistoryMax: 18
Tree:19,pos:0.00893962,len:4.8401,TMRCA:0.931611,ARG:,len:4.8401,TMRCA:0.931611
iHistoryMax: 19
Tree:20,pos:0.0108458,len:4.86177,TMRCA:0.931611,ARG:,len:4.86177,TMRCA:0.931611
iHistoryMax: 20
Tree:21,pos:0.0141077,len:4.8297,TMRCA:0.931611,ARG:,len:4.8297,TMRCA:0.931611
iHistoryMax: 21
Tree:22,pos:0.0150825,len:4.8383,TMRCA:0.931611,ARG:,len:4.8383,TMRCA:0.931611
iHistoryMax: 22
Tree:23,pos:0.0151725,len:4.93476,TMRCA:0.931611,ARG:,len:4.93476,TMRCA:0.931611
iHistoryMax: 23
Tree:24,pos:0.0160541,len:4.87574,TMRCA:0.931611,ARG:,len:4.87574,TMRCA:0.931611
iHistoryMax: 24
Tree:25,pos:0.0161213,len:4.88225,TMRCA:0.931611,ARG:,len:4.88225,TMRCA:0.931611
iHistoryMax: 25
Tree:26,pos:0.0162664,len:4.86557,TMRCA:0.931611,ARG:,len:4.86557,TMRCA:0.931611
iHistoryMax: 26
Tree:27,pos:0.0176229,len:4.82675,TMRCA:0.931611,ARG:,len:4.82675,TMRCA:0.931611
iHistoryMax: 27
Tree:28,pos:0.0177889,len:4.79918,TMRCA:0.931611,ARG:,len:4.79918,TMRCA:0.931611
iHistoryMax: 28
Tree:29,pos:0.0178136,len:4.80446,TMRCA:0.931611,ARG:,len:4.80446,TMRCA:0.931611
iHistoryMax: 29
Tree:30,pos:0.0183046,len:4.71801,TMRCA:0.888386,ARG:,len:4.71801,TMRCA:0.888386
iHistoryMax: 30
Tree:31,pos:0.0188178,len:4.8188,TMRCA:0.888386,ARG:,len:4.8188,TMRCA:0.888386
iHistoryMax: 31
Tree:32,pos:0.0190226,len:5.07545,TMRCA:0.888386,ARG:,len:5.07545,TMRCA:0.888386
iHistoryMax: 32
Tree:33,pos:0.0200874,len:5.04438,TMRCA:0.888386,ARG:,len:5.04438,TMRCA:0.888386
iHistoryMax: 33
Tree:34,pos:0.0206707,len:5.10663,TMRCA:0.888386,ARG:,len:5.10663,TMRCA:0.888386
iHistoryMax: 34
Tree:35,pos:0.0224353,len:5.12557,TMRCA:0.888386,ARG:,len:5.12557,TMRCA:0.888386
iHistoryMax: 35
Tree:36,pos:0.0245595,len:5.21311,TMRCA:0.888386,ARG:,len:5.21311,TMRCA:0.888386
iHistoryMax: 36
Tree:37,pos:0.0247924,len:5.22454,TMRCA:0.888386,ARG:,len:5.22454,TMRCA:0.888386
iHistoryMax: 37
Tree:38,pos:0.0249896,len:5.27006,TMRCA:0.888386,ARG:,len:5.27006,TMRCA:0.888386
iHistoryMax: 38
Tree:39,pos:0.0256026,len:5.23164,TMRCA:0.888386,ARG:,len:5.23164,TMRCA:0.888386
iHistoryMax: 39
Tree:40,pos:0.0267346,len:5.23164,TMRCA:0.888386,ARG:,len:5.23164,TMRCA:0.888386
iHistoryMax: 40
Tree:41,pos:0.0277982,len:5.12379,TMRCA:0.888386,ARG:,len:5.12379,TMRCA:0.888386
iHistoryMax: 41
Tree:42,pos:0.0287967,len:5.07932,TMRCA:0.888386,ARG:,len:5.07932,TMRCA:0.888386
iHistoryMax: 42
Tree:43,pos:0.030091,len:5.10252,TMRCA:0.888386,ARG:,len:5.10252,TMRCA:0.888386
iHistoryMax: 43
Tree:44,pos:0.0328631,len:5.08487,TMRCA:0.888386,ARG:,len:5.08487,TMRCA:0.888386
iHistoryMax: 44
Tree:45,pos:0.0335227,len:5.11512,TMRCA:0.888386,ARG:,len:5.11512,TMRCA:0.888386
iHistoryMax: 45
Tree:46,pos:0.0335704,len:5.22972,TMRCA:0.888386,ARG:,len:5.22972,TMRCA:0.888386
iHistoryMax: 46
Tree:47,pos:0.034135,len:5.22972,TMRCA:0.888386,ARG:,len:5.22972,TMRCA:0.888386
iHistoryMax: 47
Tree:48,pos:0.034512,len:5.21659,TMRCA:0.888386,ARG:,len:5.21659,TMRCA:0.888386
iHistoryMax: 48
Tree:49,pos:0.0347851,len:5.2289,TMRCA:0.888386,ARG:,len:5.2289,TMRCA:0.888386
iHistoryMax: 49
Tree:50,pos:0.0357164,len:5.06987,TMRCA:0.888386,ARG:,len:5.06987,TMRCA:0.888386
iHistoryMax: 50
Tree:51,pos:0.036649,len:5.0612,TMRCA:0.888386,ARG:,len:5.0612,TMRCA:0.888386
iHistoryMax: 51
Tree:52,pos:0.0367186,len:5.03766,TMRCA:0.888386,ARG:,len:5.03766,TMRCA:0.888386
iHistoryMax: 52
Tree:53,pos:0.0374882,len:5.03766,TMRCA:0.888386,ARG:,len:5.03766,TMRCA:0.888386
iHistoryMax: 53
Tree:54,pos:0.0390489,len:5.03766,TMRCA:0.888386,ARG:,len:5.03766,TMRCA:0.888386
iHistoryMax: 54
Tree:55,pos:0.0393241,len:4.2483,TMRCA:0.493706,ARG:,len:4.2483,TMRCA:0.493706
iHistoryMax: 55
Tree:56,pos:0.0419356,len:3.99931,TMRCA:0.493706,ARG:,len:3.99931,TMRCA:0.493706
iHistoryMax: 56
Tree:57,pos:0.042461,len:3.98027,TMRCA:0.493706,ARG:,len:3.98027,TMRCA:0.493706
iHistoryMax: 57
Tree:58,pos:0.0425497,len:3.92064,TMRCA:0.493706,ARG:,len:3.92064,TMRCA:0.493706
iHistoryMax: 58
Tree:59,pos:0.0429177,len:3.90525,TMRCA:0.493706,ARG:,len:3.90525,TMRCA:0.493706
iHistoryMax: 59
Tree:60,pos:0.0436361,len:4.33299,TMRCA:0.586152,ARG:,len:4.33299,TMRCA:0.586152
iHistoryMax: 60
Tree:61,pos:0.044554,len:4.26819,TMRCA:0.586152,ARG:,len:4.26819,TMRCA:0.586152
iHistoryMax: 61
Tree:62,pos:0.0469577,len:4.24554,TMRCA:0.586152,ARG:,len:4.24554,TMRCA:0.586152
iHistoryMax: 62
Tree:63,pos:0.0472498,len:3.72244,TMRCA:0.428913,ARG:,len:3.72244,TMRCA:0.428913
iHistoryMax: 63
Tree:64,pos:0.0484561,len:3.72087,TMRCA:0.428913,ARG:,len:3.72087,TMRCA:0.428913
iHistoryMax: 64
Tree:65,pos:0.0542318,len:3.71479,TMRCA:0.428913,ARG:,len:3.71479,TMRCA:0.428913
iHistoryMax: 65
Tree:66,pos:0.0550825,len:3.71902,TMRCA:0.428913,ARG:,len:3.71902,TMRCA:0.428913
iHistoryMax: 66
Tree:67,pos:0.0553246,len:3.55166,TMRCA:0.428913,ARG:,len:3.55166,TMRCA:0.428913
iHistoryMax: 67
Tree:68,pos:0.0598516,len:3.54924,TMRCA:0.428913,ARG:,len:3.54924,TMRCA:0.428913
iHistoryMax: 68
Tree:69,pos:0.0599233,len:3.60749,TMRCA:0.428913,ARG:,len:3.60749,TMRCA:0.428913
iHistoryMax: 69
Tree:70,pos:0.0599688,len:3.77489,TMRCA:0.428913,ARG:,len:3.77489,TMRCA:0.428913
iHistoryMax: 70
Tree:71,pos:0.0601909,len:3.77105,TMRCA:0.428913,ARG:,len:3.77105,TMRCA:0.428913
iHistoryMax: 71
Tree:72,pos:0.0602825,len:3.76185,TMRCA:0.428913,ARG:,len:3.76185,TMRCA:0.428913
iHistoryMax: 72
Tree:73,pos:0.0608065,len:3.77507,TMRCA:0.428913,ARG:,len:3.77507,TMRCA:0.428913
iHistoryMax: 73
Tree:74,pos:0.0611727,len:3.78667,TMRCA:0.428913,ARG:,len:3.78667,TMRCA:0.428913
iHistoryMax: 74
Tree:75,pos:0.0619994,len:3.77876,TMRCA:0.428913,ARG:,len:3.77876,TMRCA:0.428913
iHistoryMax: 75
Tree:76,pos:0.06241,len:3.58294,TMRCA:0.401441,ARG:,len:3.58294,TMRCA:0.401441
iHistoryMax: 76
Tree:77,pos:0.0635676,len:3.58945,TMRCA:0.401441,ARG:,len:3.58945,TMRCA:0.401441
iHistoryMax: 77
Tree:78,pos:0.0651993,len:3.61721,TMRCA:0.401441,ARG:,len:3.61721,TMRCA:0.401441

Original comment by julia_pa...@brown.edu on 20 Jun 2014 at 2:37

GoogleCodeExporter commented 9 years ago

Hi Julia

I was able to reproduce these lines you cut and pasted. But what anomaly am I 
specifically supposed to be looking for?

Original comment by gche...@gmail.com on 20 Jun 2014 at 4:48

Changed state: Started

GoogleCodeExporter commented 9 years ago

Hi Gary,
Right! Some genealogies have negative branch lengths. For example genealogy 953 
(See file attached)
Thank you so much,

Original comment by julia_pa...@brown.edu on 20 Jun 2014 at 5:04

GoogleCodeExporter commented 9 years ago

Original comment by julia_pa...@brown.edu on 20 Jun 2014 at 5:04

Attachments:

del2_nbl.pdf

GoogleCodeExporter commented 9 years ago

Can you upload a gzipped file of the stderr and stdout and point me to the line 
where you see this negative value? It's hard to me to map what R is identifying 
as the 953rd tree.

Original comment by gche...@gmail.com on 20 Jun 2014 at 5:09

GoogleCodeExporter commented 9 years ago

Thank you Gary!! 

I looked at the tree file: 
awk '$1~"NEWICK_TREE"' Const100Idv_5.tree | awk -F\] '{print $2}' 
>Const100Idv_5.tree.txt

and then the tree(s) causing me trouble:
awk 'NR==953 {print $0}' ./Const100Idv_5.txt 

and it turns out that the problem is in R - ape with the function 
coalescent.intervals(), so I just fixed that part. Thank you so much! Problem 
solved!

Original comment by julia_pa...@brown.edu on 20 Jun 2014 at 7:45

GoogleCodeExporter commented 9 years ago

Sounds good. Glad it was relatively painless on both out parts. Thanks for 
investigating.

Original comment by gche...@gmail.com on 20 Jun 2014 at 8:34

Changed state: Fixed

gchen98 / macs

Negative branch lengths #24