bug with reading network, else question about network format

[cecileane]

Hi @shajoezhu and @jamdeg, I gave a time-consistent (ultrametric) network as input, but hybrid-lambda complains that it's not ultrametric, and the simulated gene trees are not ultrametric. So here is the example:

hybrid-Lambda -spcu '((A:1.0)H19#0.6:4.0,(B:3.0,(C:1.0,H19#0.6:0.0)I1:2.0)I2:2.0)I3;' -num 10 -seed 1234 -o net8trial

Default Kingman coalescent on all branches.
Default population size of 10000 on all branches. 
WARNING! NOT ULTRAMETRIC!!!
Random seed: 1234
WARNING: Gene tree is not ultrametric
WARNING: Gene tree is not ultrametric
...
Produced gene tree files: 
net8trial_coal_unit

Here is a plot of the network, with edges labelled by their lengths: it is ultrametric (contrary to hybrid-lambda's warning), so it seems that my network string description is wrong, or that it's bug:

This network image, the standard Newick string for the network and the network string used as input to hybrid-lambda were obtained in julia like this:

using PhyloNetworks
net8 = readTopology("((A:1)#H19:4::0.6,(B:3,(C:1,#H19:0.0::0.4):2):2);")
plot(net8,:R; showEdgeLength=true, useEdgeLength=true);
netHL = hybridlambdaformat(net8) # "((A:1.0)H19#0.6:4.0,(B:3.0,(C:1.0,H19#0.6:0.0)I1:2.0)I2:2.0)I3;"

The first 6 gene trees look like this, where we can see a pattern of species A being too "far" by 4 units, or too close to the root by 4 units, depending on how it was inherited at the reticulation.

If I misunderstand the format required by hybrid-lambda to write the network, I would love to learn more! Many thanks. @RC0515

[cecileane]

It looks like hybrid-lambda is actually using the following network to simulate gene trees (with edges not drawn proportional to their lengths, because it's impossible on a time-inconsistent network):

So the 2 hybrid edges would have had their branch lengths switched... Here is the standard extended Newick description for this alternate network:

net8wrong = readTopology("((A:1)#H19:0::0.6,(B:3,(C:1,#H19:4::0.4):2):2);")
hybridlambdaformat(net8wrong)
"((A:1.0)H19#0.6:0.0,(B:3.0,(C:1.0,H19#0.6:4.0)I1:2.0)I2:2.0)I3;"

[cecileane]

It looks like the problem goes away if I "rotate" the 2 edges at the root. In the parenthetical description, it means that the "BC clade" comes first, and major reticulation edge + A comes later:

rotate!(net8, -2)
writeTopology(net8)
"((B:3.0,(C:1.0,#H19:0.0::0.4):2.0):2.0,(A:1.0)#H19:4.0::0.6);"
netHL = hybridlambdaformat(net8)
"((B:3.0,(C:1.0,H19#0.4:0.0)I1:2.0)I2:2.0,(A:1.0)H19#0.4:4.0)I3;"

When calling hybrid-lambda with this alternative parenthetical description, there aren't any warnings, and the simulations seem to come from the intended network:

hybrid-Lambda -spcu '((B:3.0,(C:1.0,H19#0.4:0.0)I1:2.0)I2:2.0,(A:1.0)H19#0.4:4.0)I3;' -num 10 -seed 1234 -o net8trial

with these first 6 simulated gene trees:

Is there a way to formalize the assumption made by hybrid-lambda, with regards to the ordering of clades, to make sure the network parenthetical description is correctly parsed & interpreted?

[jamdeg]

HI Cécile, It's been a while since I've tried to set these up. I'm sorry that it isn't clear how to formalize this. Hopefully Joe can shed some light on it.

James

On Mon, Sep 7, 2020 at 10:17 AM Cécile Ané notifications@github.com wrote:

It looks like hybrid-lambda is actually using the following network to simulate gene trees (with edges not drawn proportional to their lengths, because it's impossible on a time-inconsistent network): [image: unwanted species network] https://user-images.githubusercontent.com/9857294/92405486-4ac36b80-f0fb-11ea-8425-4ab284483786.png

So the 2 hybrid edges would have had their branch lengths switched... Here is the standard extended Newick description for this alternate network:

net8wrong = readTopology("((A:1)#H19:0::0.6,(B:3,(C:1,#H19:4::0.4):2):2);")hybridlambdaformat(net8wrong)"((A:1.0)H19#0.6:0.0,(B:3.0,(C:1.0,H19#0.6:4.0)I1:2.0)I2:2.0)I3;"

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/hybridLambda/hybrid-Lambda/issues/36#issuecomment-688419040, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACHBVJ5C6W7QVRC6TSQNOV3SEUBTLANCNFSM4Q6N23LA .

[shajoezhu]

HI Cécile

Thanks for reporting this. Definitely a bug while parsing. I will look into it.

Usually, I use the hybridlambda plot flag -plot or -dot to check the input network if it makes sense before I run my simulation.

Best Joe

On Mon, 7 Sep 2020 at 19:04, jamdeg notifications@github.com wrote:

HI Cécile,

It's been a while since I've tried to set these up. I'm sorry that it isn't

clear how to formalize this. Hopefully Joe can shed some light on it.

James

On Mon, Sep 7, 2020 at 10:17 AM Cécile Ané notifications@github.com wrote:

It looks like hybrid-lambda is actually using the following network to

simulate gene trees (with edges not drawn proportional to their

lengths, because it's impossible on a time-inconsistent network):

[image: unwanted species network]

< https://user-images.githubusercontent.com/9857294/92405486-4ac36b80-f0fb-11ea-8425-4ab284483786.png

So the 2 hybrid edges would have had their branch lengths switched... Here

is the standard extended Newick description for this alternate network:

net8wrong = readTopology("((A:1)#H19:0::0.6,(B:3,(C:1,#H19:4::0.4):2):2);")hybridlambdaformat(net8wrong)"((A:1.0)H19#0.6:0.0,(B:3.0,(C:1.0,H19#0.6:4.0)I1:2.0)I2:2.0)I3;"

—

You are receiving this because you were mentioned.

Reply to this email directly, view it on GitHub

< https://github.com/hybridLambda/hybrid-Lambda/issues/36#issuecomment-688419040 ,

or unsubscribe

< https://github.com/notifications/unsubscribe-auth/ACHBVJ5C6W7QVRC6TSQNOV3SEUBTLANCNFSM4Q6N23LA

.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/hybridLambda/hybrid-Lambda/issues/36#issuecomment-688459937, or unsubscribe https://github.com/notifications/unsubscribe-auth/AA4FP7JAHOFAIA54XEE7ZA3SEUOBHANCNFSM4Q6N23LA .

[shajoezhu]

@cecileane I think I understand the parsing issue now. For '((A:1.0)H19#0.6:4.0,(B:3.0,(C:1.0,H19#0.6:0.0)I1:2.0)I2:2.0)I3;'

H19 has first branch length as 4.0, second branch length as 0.0. However the first parent node should be defined as I3, but it is defined after I1, because hybrid-Lambda read from left to right, and create node, as it sees it appears. So the 4.0 branch was wrongly assigned with parent node of I1 because it is defined before I3.

hybrid-Lambda -spcu '((A:1.0)H19#0.6:4.0,(B:3.0,(C:1.0,H19#0.6:0.0)I1:2.0)I2:2.0)I3;' -dot -o net1 -branch

net1.pdf

Whereas, '((B:3.0,(C:1.0,H19#0.4:0.0)I1:2.0)I2:2.0,(A:1.0)H19#0.4:4.0)I3;' branch 0.0 is the first branch which is correctly connected to the first parent node I1.

hybrid-Lambda -spcu '((B:3.0,(C:1.0,H19#0.4:0.0)I1:2.0)I2:2.0,(A:1.0)H19#0.4:4.0)I3;' -dot -o net2 -branch

net2.pdf

[jamdeg]

Thanks for looking into it, Joe.

On Thu, Sep 10, 2020 at 2:56 PM Joe Zhu notifications@github.com wrote:

@cecileane https://github.com/cecileane I think I understand the parsing issue now. For '((A:1.0)H19#0.6:4.0,(B:3.0,(C:1.0,H19#0.6:0.0)I1:2.0)I2:2.0)I3;'

H19 has first branch length as 4.0, second branch length as 0.0. However the first parent node should be defined as I3, but it is defined after I1, because hybrid-Lambda read from left to right, and create node, as it sees it appears. So the 4.0 branch was wrongly assigned with parent node of I1 because it is defined before I3.

hybrid-Lambda -spcu '((A:1.0)H19#0.6:4.0,(B:3.0,(C:1.0,H19#0.6:0.0)I1:2.0)I2:2.0)I3;' -dot -o net1 -branch

net1.pdf https://github.com/hybridLambda/hybrid-Lambda/files/5204611/net1.pdf

Whereas, '((B:3.0,(C:1.0,H19#0.4:0.0)I1:2.0)I2:2.0,(A:1.0)H19#0.4:4.0)I3;' branch 0.0 is the first branch which is correctly connected to the first parent node I1.

hybrid-Lambda -spcu '((B:3.0,(C:1.0,H19#0.4:0.0)I1:2.0)I2:2.0,(A:1.0)H19#0.4:4.0)I3;' -dot -o net2 -branch

net2.pdf https://github.com/hybridLambda/hybrid-Lambda/files/5204614/net2.pdf

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/hybridLambda/hybrid-Lambda/issues/36#issuecomment-690725110, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACHBVJ56ACLK6EYHDHHG6E3SFE4QBANCNFSM4Q6N23LA .

[cecileane]

Do you think you might have the bandwidth to fix this parsing issue in hybrid-lambda?

In any case, this diagnostic is very helpful. To work around the issue, we would need to write a function to check in advance the order in which the 2 partner hybrid edges will appear, the order in which the 2 parent node names will appear, and then "rotate" sibling edges at appropriate nodes to ensure that the "first" hybrid edge will correspond to the "first" parent node. For now, I don't see what algorithm can find an appropriate rotation, however.

Please let me know if my understanding of this order assumption if incorrect or incomplete.

Parsing in PhyloNetworks also goes from left to right, but each node is created when its first opening parenthesis is read (not when its name is read at the closing parenthesis). Edges can be attached to both of their nodes when edges are "read", since the nodes already exist. When the "first" edge above H19 is read, it can be created with correct parameters (length 4.0 and gamma inheritance 0.6). The root node I3 is created after reading the very first parenthesis, although without a name at first, since the name comes at the very end. A node name may stay empty, or be updated after reading the closing parenthesis, if the name exists.

[shajoezhu]

Hi Cecile

Do you think you might have the bandwidth to fix this parsing issue in hybrid-lambda?

I will try and find some time. Apologies if I can not prioritize this. I think plot and check is the best option for now. I will do my best and find some time to work on the algorithm.

In any case, this diagnostic is very helpful. To work around the issue, we would need to write a function to check in advance the order in which the 2 partner hybrid edges will appear, the order in which the 2 parent node names will appear, and then "rotate" sibling edges at appropriate nodes to ensure that the "first" hybrid edge will correspond to the "first" parent node. For now, I don't see what algorithm can find an appropriate rotation, however.

Yes, I always prefer to write my trees and networks with brackets ( together or close by, I.e. ((a, b), c) rather than (c,(a, b)).

Please let me know if my understanding of this order assumption if incorrect or incomplete.

Parsing in PhyloNetworks also goes from left to right, but each node is created when its first opening parenthesis is read (not when its name is read at the closing parenthesis). Edges can be attached to both of their nodes when edges are "read", since the nodes already exist. When the "first" edge above H19 is read, it can be created with correct parameters (length 4.0 and gamma inheritance 0.6). The root node I3 is created after reading the very first parenthesis, although without a name at first, since the name comes at the very end. A node name may stay empty, or be updated after reading the closing parenthesis, if the name exists.

Spot on! I might have to create a cache space for the entire string, so when moving the cursor I can also scan what's on the right of the cursor. So I will be able to extract the correct parent node, not sure how to make this process efficient at the moment.

Any idea is welcome!

Best Joe

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/hybridLambda/hybrid-Lambda/issues/36#issuecomment-691602697, or unsubscribe https://github.com/notifications/unsubscribe-auth/AA4FP7LPTIVWFYZ2XRSBLN3SFQ36XANCNFSM4Q6N23LA .

[shajoezhu]

Okay, the above approach didn't seem to work.

Trying new approach, for each node, Mark the index of node name, for hybrid nodes, Mark the first index as it appears the first time, and second index, as it appears the second time. When linking the hybrid node, for first index, scan to the right, til the first ), and get the node name. Note that, for the first ) that is not paired with ( that appears after the first hybrid node index.

Do the same for the second node.

I will try this now

[jamdeg]

Thanks for working on this, Joe.

On Sun, Sep 20, 2020 at 2:28 PM Joe Zhu notifications@github.com wrote:

Okay, the above approach didn't seem to work.

Trying new approach, for each node, Mark the index of node name, for hybrid nodes, Mark the first index as it appears the first time, and second index, as it appears the second time. When linking the hybrid node, for first index, scan to the right, til the first ), and get the node name. Note that, for the first ) that is not paired with ( that appears after the first hybrid node index.

Do the same for the second node.

I will try this now

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/hybridLambda/hybrid-Lambda/issues/36#issuecomment-695832264, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACHBVJ2RPC67XFFSTRUFIJLSGZQVXANCNFSM4Q6N23LA .

[shajoezhu]

Hi @jamdeg and @cecileane

Hope all is well. Sorry this took a bit longer. Just made some quick changes on branch 36_network_reading,

The test case seem to work now,

hybrid-Lambda -spcu '((A:1.0)H19#0.6:4.0,(B:3.0,(C:1.0,H19#0.6:0.0)I1:2.0)I2:2.0)I3;' -dot -o net1 -branch

hybrid-Lambda -spcu '((B:3.0,(C:1.0,H19#0.4:0.0)I1:2.0)I2:2.0,(A:1.0)H19#0.4:4.0)I3;' -dot -o net2 -branch

net1.pdf net2.pdf

Please let me know if there is problem.

Thanks! Joe

[cecileane]

Thank you so much!!! I will find time to try out the new version on various examples that gave us trouble. will let you know!

[cecileane]

All the examples above run well with the version on your new branch. So I tried on 9-taxon network that had lead me to the 3-taxon example above. Unfortunately I hit the same bug again... and a new one. (sorry for being the messenger of bad news!). Here is a reproducible example.

I used 4 networks: all with the same topology: 9 taxa, 2 reticulations, level 1, ultrametric.

• n2: everything works fine
• n3: same as n2, but with daughter edges "rotated" around a particular node
• n2s: same as n2, but with species names modified to be "s1" through "s9"
• n3s: same as n3, but with species names modified to be "s1" through "s9"

Simulations look good on n2 but the bug appears on n3: see warnings below. (this is using your bug fix on commit 781eb35)

% hybrid-Lambda -spcu '(((HS03:6.12,TMS01:6.12)I1:0.71,(NBLZG:5.98,H18#0.21:0.0)I2:0.85)I3:0.47,((AYEM:7.13,(XYBX:5.98)H18#0.21:1.15)I4:0.1,ZGQCS:7.23)I5:0.07,((YNEM:5.64,(DLYQCS:0.0,H19#0.11:0.0)I6:5.64)I7:0.56,(JYYJFS:0.0)H19#0.11:6.2)I8:1.1)I9;' -num 10 -seed 42 -o n2

Default Kingman coalescent on all branches.
Default population size of 10000 on all branches. 
Random seed: 42
Produced gene tree files: 
n2_coal_unit

% hybrid-Lambda -spcu '(((HS03:6.12,TMS01:6.12)I1:0.71,(NBLZG:5.98,H18#0.21:0.0)I2:0.85)I3:0.47,((AYEM:7.13,(XYBX:5.98)H18#0.21:1.15)I4:0.1,ZGQCS:7.23)I5:0.07,((JYYJFS:0.0)H19#0.89:6.2,(YNEM:5.64,(DLYQCS:0.0,H19#0.89:0.0)I6:5.64)I7:0.56)I8:1.1)I9;' -num 10 -seed 42 -o n3

Default Kingman coalescent on all branches.
Default population size of 10000 on all branches. 
WARNING! NOT ULTRAMETRIC!!!
Random seed: 42
WARNING: Gene tree is not ultrametric
WARNING: Gene tree is not ultrametric
WARNING: Gene tree is not ultrametric
WARNING: Gene tree is not ultrametric
WARNING: Gene tree is not ultrametric
WARNING: Gene tree is not ultrametric
WARNING: Gene tree is not ultrametric
WARNING: Gene tree is not ultrametric
WARNING: Gene tree is not ultrametric
WARNING: Gene tree is not ultrametric
Produced gene tree files: 
n3_coal_unit

The taxon JYYJFS (just below the "second" reticulation) is not aligned with the other taxa: the external branch to this taxon is either too short or too long by 6.2 exactly. This looks like the same bug as before, as if the parents edges at that reticulations had their branch lengths switched (6.2 and 0).

I then wanted to modify the tip labels to write this bug report. And then surprise: the behavior changed after changing the tip labels. On n2s, no warnings (as expected). But on n3s, the warnings disappeared and the gene trees were ultrametric!

%hybrid-Lambda -spcu '(((s1:6.12,s2:6.12)I1:0.71,(s3:5.98,H18#0.21:0.0)I2:0.85)I3:0.47,((s4:7.13,(s5:5.98)H18#0.21:1.15)I4:0.1,s6:7.23)I5:0.07,((s7:5.64,(s8:0.0,H19#0.11:0.0)I6:5.64)I7:0.56,(s9:0.0)H19#0.11:6.2)I8:1.1)I9;' -num 10 -seed 42 -o n2s

Default Kingman coalescent on all branches.
Default population size of 10000 on all branches. 
Random seed: 42
Produced gene tree files: 
n2s_coal_unit

% hybrid-Lambda -spcu '(((s1:6.12,s2:6.12)I1:0.71,(s3:5.98,H18#0.21:0.0)I2:0.85)I3:0.47,((s4:7.13,(s5:5.98)H18#0.21:1.15)I4:0.1,s6:7.23)I5:0.07,((s9:0.0)H19#0.89:6.2,(s7:5.64,(s8:0.0,H19#0.89:0.0)I6:5.64)I7:0.56)I8:1.1)I9;' -num 10 -seed 42 -o n3s

Default Kingman coalescent on all branches.
Default population size of 10000 on all branches. 
Random seed: 42
Produced gene tree files: 
n3s_coal_unit

Why the taxon labels affect the simulations?? It looks like whether the branch lengths are switched, or not, is affected by the branch labels.

For completeness, here is the julia code that I used to make the plots, with networks in extended newick format. (I also used julia to modify taxon names, rotate edges around a node, and get the network string for hybrid-Lambda).

using PhyloNetworks # v0.12.0
using RCall         # v0.13.10
using PhyloPlots    # add PhyloPlots#master

n2 = readTopology("(((HS03:6.12,TMS01:6.12):0.71,(NBLZG:5.98,#H18:0.0::0.21):0.85):0.47,((AYEM:7.13,(XYBX:5.98)#H18:1.15::0.79):0.1,ZGQCS:7.23):0.07,((YNEM:5.64,(DLYQCS:0.0,#H19:0.0::0.11):5.64):0.56,(JYYJFS:0.0)#H19:6.2::0.89):1.1);")
hybridlambdaformat(n2) # to get format required by hybrid-Lambda: "(((HS03:6.12,TMS01:6.12)I1:0.71,(NBLZG:5.98,H18#0.21:0.0)I2:0.85)I3:0.47,((AYEM:7.13,(XYBX:5.98)H18#0.21:1.15)I4:0.1,ZGQCS:7.23)I5:0.07,((YNEM:5.64,(DLYQCS:0.0,H19#0.11:0.0)I6:5.64)I7:0.56,(JYYJFS:0.0)H19#0.11:6.2)I8:1.1)I9;"
n3 = readTopology("(((HS03:6.12,TMS01:6.12):0.71,(NBLZG:5.98,#H18:0.0::0.21):0.85):0.47,((AYEM:7.13,(XYBX:5.98)#H18:1.15::0.79):0.1,ZGQCS:7.23):0.07,((JYYJFS:0.0)#H19:6.2::0.89,(YNEM:5.64,(DLYQCS:0.0,#H19:0.0::0.11):5.64):0.56):1.1);")
n2s = readTopology("(((s1:6.12,s2:6.12):0.71,(s3:5.98,#H18:0.0::0.21):0.85):0.47,((s4:7.13,(s5:5.98)#H18:1.15::0.79):0.1,s6:7.23):0.07,((s7:5.64,(s8:0.0,#H19:0.0::0.11):5.64):0.56,(s9:0.0)#H19:6.2::0.89):1.1);")
n3s = readTopology("(((s1:6.12,s2:6.12):0.71,(s3:5.98,#H18:0.0::0.21):0.85):0.47,((s4:7.13,(s5:5.98)#H18:1.15::0.79):0.1,s6:7.23):0.07,((s9:0.0)#H19:6.2::0.89,(s7:5.64,(s8:0.0,#H19:0.0::0.11):5.64):0.56):1.1);")

R"layout"([1 2; 3 4])
plot(n2,:R; showEdgeLength=true, useEdgeLength=true);
R"mtext"("n2", line=-2, adj=0.1);
plot(n3,:R; showEdgeLength=true, useEdgeLength=true);
R"mtext"("n3", line=-2, adj=0.1);
plot(n2s,:R; showEdgeLength=true, useEdgeLength=true);
R"mtext"("n2s", line=-2, adj=0.1);
plot(n3s,:R; showEdgeLength=true, useEdgeLength=true);
R"mtext"("n3s", line=-2, adj=0.1);

R"layout"([1 2]); R"par"(mar=[0,0,0,1]);
for gt in readMultiTopology("n2_coal_unit")[1:2]
    plot(gt,:R; showEdgeLength=true, useEdgeLength=true)
end
R"mtext"("first 2 gene trees from n2", line=-1, adj=0.05, outer=true);
# and similar for gene trees from n3, n2s and n3s.

[shajoezhu]

Hi Cecile,

Many thanks for this! Let me take another look at it.

Best Joe

On Sun, 29 Nov 2020 at 02:51, Cécile Ané notifications@github.com wrote:

All the examples above run well with the version on your new branch. So I tried on 9-taxon network that had lead me to the 3-taxon example above. Unfortunately I hit the same bug again... and a new one. (sorry for being the messenger of bad news!). Here is a reproducible example.

I used 4 networks: all with the same topology: 9 taxa, 2 reticulations, level 1, ultrametric.

• n2: everything works fine
• n3: same as n2, but with daughter edges "rotated" around a particular node
• n2s: same as n2, but with species names modified to be "s1" through "s9"
• n3s: same as n3, but with species names modified to be "s1" through "s9"

[image: nets_n2n3n2sn3s] https://user-images.githubusercontent.com/9857294/100530415-7bc1b280-31b7-11eb-8122-1177e5360830.png

Simulations look good on n2 but the bug appears on n3: see warnings below. (this is using your bug fix on commit 781eb35 https://github.com/hybridLambda/hybrid-Lambda/commit/781eb352bc7eacaec0800fb2b97f8c01910af708 )

% hybrid-Lambda -spcu '(((HS03:6.12,TMS01:6.12)I1:0.71,(NBLZG:5.98,H18#0.21:0.0)I2:0.85)I3:0.47,((AYEM:7.13,(XYBX:5.98)H18#0.21:1.15)I4:0.1,ZGQCS:7.23)I5:0.07,((YNEM:5.64,(DLYQCS:0.0,H19#0.11:0.0)I6:5.64)I7:0.56,(JYYJFS:0.0)H19#0.11:6.2)I8:1.1)I9;' -num 10 -seed 42 -o n2

Default Kingman coalescent on all branches. Default population size of 10000 on all branches. Random seed: 42 Produced gene tree files: n2_coal_unit

[image: gt_n2] https://user-images.githubusercontent.com/9857294/100530560-3bfbca80-31b9-11eb-9135-bed5745bb4f7.png

% hybrid-Lambda -spcu '(((HS03:6.12,TMS01:6.12)I1:0.71,(NBLZG:5.98,H18#0.21:0.0)I2:0.85)I3:0.47,((AYEM:7.13,(XYBX:5.98)H18#0.21:1.15)I4:0.1,ZGQCS:7.23)I5:0.07,((JYYJFS:0.0)H19#0.89:6.2,(YNEM:5.64,(DLYQCS:0.0,H19#0.89:0.0)I6:5.64)I7:0.56)I8:1.1)I9;' -num 10 -seed 42 -o n3

Default Kingman coalescent on all branches. Default population size of 10000 on all branches. WARNING! NOT ULTRAMETRIC!!! Random seed: 42 WARNING: Gene tree is not ultrametric WARNING: Gene tree is not ultrametric WARNING: Gene tree is not ultrametric WARNING: Gene tree is not ultrametric WARNING: Gene tree is not ultrametric WARNING: Gene tree is not ultrametric WARNING: Gene tree is not ultrametric WARNING: Gene tree is not ultrametric WARNING: Gene tree is not ultrametric WARNING: Gene tree is not ultrametric Produced gene tree files: n3_coal_unit

[image: gt_n3] https://user-images.githubusercontent.com/9857294/100530570-52098b00-31b9-11eb-9490-30da8db1b1f3.png

The taxon JYYJFS (just below the "second" reticulation) is not aligned with the other taxa: the external branch to this taxon is either too short or too long by 6.2 exactly. This looks like the same bug as before, as if the parents edges at that reticulations had their branch lengths switched (6.2 and 0).

I then wanted to modify the tip labels to write this bug report. And then surprise: the behavior changed after changing the tip labels. On n2s, no warnings (as expected). But on n3s, the warnings disappeared and the gene trees were ultrametric!

%hybrid-Lambda -spcu '(((s1:6.12,s2:6.12)I1:0.71,(s3:5.98,H18#0.21:0.0)I2:0.85)I3:0.47,((s4:7.13,(s5:5.98)H18#0.21:1.15)I4:0.1,s6:7.23)I5:0.07,((s7:5.64,(s8:0.0,H19#0.11:0.0)I6:5.64)I7:0.56,(s9:0.0)H19#0.11:6.2)I8:1.1)I9;' -num 10 -seed 42 -o n2s

Default Kingman coalescent on all branches. Default population size of 10000 on all branches. Random seed: 42 Produced gene tree files: n2s_coal_unit

[image: gt_n2s] https://user-images.githubusercontent.com/9857294/100530576-5a61c600-31b9-11eb-8778-0e6ea4d533a0.png

% hybrid-Lambda -spcu '(((s1:6.12,s2:6.12)I1:0.71,(s3:5.98,H18#0.21:0.0)I2:0.85)I3:0.47,((s4:7.13,(s5:5.98)H18#0.21:1.15)I4:0.1,s6:7.23)I5:0.07,((s9:0.0)H19#0.89:6.2,(s7:5.64,(s8:0.0,H19#0.89:0.0)I6:5.64)I7:0.56)I8:1.1)I9;' -num 10 -seed 42 -o n3s

Default Kingman coalescent on all branches. Default population size of 10000 on all branches. Random seed: 42 Produced gene tree files: n3s_coal_unit

[image: gt_n3s] https://user-images.githubusercontent.com/9857294/100530587-75343a80-31b9-11eb-8322-27ce9387d7b8.png

Why the taxon labels affect the simulations?? It looks like whether the branch lengths are switched, or not, is affected by the branch labels.

For completeness, here is the julia code that I used to make the plots, with networks in extended newick format. (I also used julia to modify taxon names, rotate edges around a node, and get the network string for hybrid-Lambda).

using PhyloNetworks # v0.12.0using RCall # v0.13.10using PhyloPlots # add PhyloPlots#master

n2 = readTopology("(((HS03:6.12,TMS01:6.12):0.71,(NBLZG:5.98,#H18:0.0::0.21):0.85):0.47,((AYEM:7.13,(XYBX:5.98)#H18:1.15::0.79):0.1,ZGQCS:7.23):0.07,((YNEM:5.64,(DLYQCS:0.0,#H19:0.0::0.11):5.64):0.56,(JYYJFS:0.0)#H19:6.2::0.89):1.1);")hybridlambdaformat(n2) # to get format required by hybrid-Lambda: "(((HS03:6.12,TMS01:6.12)I1:0.71,(NBLZG:5.98,H18#0.21:0.0)I2:0.85)I3:0.47,((AYEM:7.13,(XYBX:5.98)H18#0.21:1.15)I4:0.1,ZGQCS:7.23)I5:0.07,((YNEM:5.64,(DLYQCS:0.0,H19#0.11:0.0)I6:5.64)I7:0.56,(JYYJFS:0.0)H19#0.11:6.2)I8:1.1)I9;" n3 = readTopology("(((HS03:6.12,TMS01:6.12):0.71,(NBLZG:5.98,#H18:0.0::0.21):0.85):0.47,((AYEM:7.13,(XYBX:5.98)#H18:1.15::0.79):0.1,ZGQCS:7.23):0.07,((JYYJFS:0.0)#H19:6.2::0.89,(YNEM:5.64,(DLYQCS:0.0,#H19:0.0::0.11):5.64):0.56):1.1);") n2s = readTopology("(((s1:6.12,s2:6.12):0.71,(s3:5.98,#H18:0.0::0.21):0.85):0.47,((s4:7.13,(s5:5.98)#H18:1.15::0.79):0.1,s6:7.23):0.07,((s7:5.64,(s8:0.0,#H19:0.0::0.11):5.64):0.56,(s9:0.0)#H19:6.2::0.89):1.1);") n3s = readTopology("(((s1:6.12,s2:6.12):0.71,(s3:5.98,#H18:0.0::0.21):0.85):0.47,((s4:7.13,(s5:5.98)#H18:1.15::0.79):0.1,s6:7.23):0.07,((s9:0.0)#H19:6.2::0.89,(s7:5.64,(s8:0.0,#H19:0.0::0.11):5.64):0.56):1.1);") R"layout"([1 2; 3 4])plot(n2,:R; showEdgeLength=true, useEdgeLength=true);R"mtext"("n2", line=-2, adj=0.1);plot(n3,:R; showEdgeLength=true, useEdgeLength=true);R"mtext"("n3", line=-2, adj=0.1);plot(n2s,:R; showEdgeLength=true, useEdgeLength=true);R"mtext"("n2s", line=-2, adj=0.1);plot(n3s,:R; showEdgeLength=true, useEdgeLength=true);R"mtext"("n3s", line=-2, adj=0.1); R"layout"([1 2]); R"par"(mar=[0,0,0,1]);for gt in readMultiTopology("n2_coal_unit")[1:2] plot(gt,:R; showEdgeLength=true, useEdgeLength=true)endR"mtext"("first 2 gene trees from n2", line=-1, adj=0.05, outer=true);# and similar for gene trees from n3, n2s and n3s.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/hybridLambda/hybrid-Lambda/issues/36#issuecomment-735321005, or unsubscribe https://github.com/notifications/unsubscribe-auth/AA4FP7OS4AQSKOBSRDRMP2DSSGZKNANCNFSM4Q6N23LA .