From 15b211f7a83be0ed069c539ba88ded2921f977c7 Mon Sep 17 00:00:00 2001 From: Cecile Ane Date: Mon, 18 Jun 2018 17:28:43 -0500 Subject: [PATCH] better cooordinates for non-tree-child networks --- src/phylonetworksPlots.jl | 37 +++++++++++++++++++++++++++------- test/test_phylonetworkPlots.jl | 27 +++++++++++++++++++++++++ test/test_rexport.jl | 8 +++++--- 3 files changed, 62 insertions(+), 10 deletions(-) diff --git a/src/phylonetworksPlots.jl b/src/phylonetworksPlots.jl index 0d01043..065e2db 100644 --- a/src/phylonetworksPlots.jl +++ b/src/phylonetworksPlots.jl @@ -28,24 +28,47 @@ function getEdgeNodeCoordinates(net::HybridNetwork, useEdgeLength::Bool) node_yB = zeros(Float64,net.numNodes) # min (B=begin) and max (E=end) node_yE = zeros(Float64,net.numNodes) # of at children's nodes nexty = ymax # first tips at the top, last at bottom - for i=length(net.node):-1:1 # post-order traversal in major tree + # set node_y of leaves: follow cladewise order + for i=length(net.node):-1:1 ni = net.cladewiseorder_nodeIndex[i] if net.node[ni].leaf node_y[ni] = nexty nexty -= 1.0 - else - node_yB[ni]=ymax; node_yE[ni]=ymin; - for e in net.node[ni].edge - if net.node[ni] == (e.isChild1 ? e.node[2] : e.node[1]) # if e = child of node - if (!e.isMajor) continue; end + end + end + # set node_y of internal nodes: follow post-order + for i=length(net.node):-1:1 + nn = net.nodes_changed[i] + !nn.leaf || continue # previous loop took care of leaves + ni = findfirst(net.node, nn) + node_yB[ni]=ymax; node_yE[ni]=ymin; + minor_yB = ymax; minor_yE = ymin; + nomajorchild=true + for e in nn.edge + if nn == PhyloNetworks.getParent(e) # if e = child of node + if e.isMajor + nomajorchild = false # we found a child edge that is a major edge yy = node_y[findfirst(net.node, PhyloNetworks.getChild(e))] yy!=0 || error("oops, child has not been visited and its y value is 0.") node_yB[ni] = min(node_yB[ni], yy) node_yE[ni] = max(node_yE[ni], yy) + elseif nomajorchild # e is minor edge, and no major found so far + yy = node_y[findfirst(net.node, PhyloNetworks.getChild(e))] + yy!=0 || error("oops, child has not been visited and its y value is 0.") + minor_yB = min(minor_yB, yy) + minor_yE = max(minor_yE, yy) end - node_y[ni] = (node_yB[ni]+node_yE[ni])/2 end end + if nomajorchild # if all children edges are minor hybrid edges. If so: not level 1, not tree-child + if minor_yB == minor_yE # one single child. jitter by 0.1 to make the plot readable + minor_yB += (minor_yB < (ymax+ymin)/2 ? 0.1 : -0.1) + minor_yE = minor_yB + end + node_yB[ni] = minor_yB + node_yE[ni] = minor_yE + end + node_y[ni] = (node_yB[ni]+node_yE[ni])/2 end # setting branch lengths for plotting diff --git a/test/test_phylonetworkPlots.jl b/test/test_phylonetworkPlots.jl index c6a1bd1..318b6df 100644 --- a/test/test_phylonetworkPlots.jl +++ b/test/test_phylonetworkPlots.jl @@ -53,4 +53,31 @@ min=[false,false,true,false,false,false,false,false,false], x=collect(6.:14),y=collect(26.:34))) + # example with level-2 network, non-tree child: + # one hybrid node ends up as a leaf in the major tree. + # no major child edge to follow to set coordinates + net = readTopology("((((B)#H1:::0.2)#H2,((D,C,#H2:::0.8)S1,(#H1,A)S2)S3)S4);") + @test_nowarn plot(net, :R, showNodeNumber=true, showGamma=true); + @test PhyloPlots.getEdgeNodeCoordinates(net, false) == ( + [5.,4.,1.,3.,3.,3.,2.,4.,4.,2.,1.], + [6.,5.,4.,6.,6.,4.,3.,5.,6.,4.,2.], + [2.,2.1,2.275,4.,3.,2.1,3.05,2.,1.,1.5,2.275], + [2.,2. ,2.1, 4.,3.,2.1,3.05,2.,1.,1.5,2.275], + [6.,5.,4., 6.,6.,3., 6.,4., 2., 1.], + [2.,2.,2.1,4.,3.,3.05,1.,1.5,2.275,2.275], + [0.,2.,2.1,0.,0.,2.1,0.,1.,1.5, 2.275], + [0.,2.,2.1,0.,0.,4. ,0.,2.,3.05,2.275], + 1.,6.,1.,4.) + net = readTopology("((((B)#H1:::0.2)#H2,((D,C,#H2)S1,(#H1,A)S2)S3)S4);") + @test_nowarn plot(net, :R, showNodeNumber=true, showGamma=true); + @test PhyloPlots.getEdgeNodeCoordinates(net, false) == ( + [5.,4.,1.,3.,3.,3.,2.,4.,4.,2.,1.], + [6.,5.,4.,6.,6.,4.,3.,5.,6.,4.,2.], + [2.,2.1,2.1,4.,3.,3.5,3.5,2.,1.,1.5,2.5], + [2.,2. ,2.1,4.,3.,2.1,3.5,2.,1.,1.5,2.5], + [6.,5.,4., 6.,6.,3., 6.,4., 2., 1.], + [2.,2.,2.1,4.,3.,3.5,1.,1.5,2.5,2.3], + [0.,2.,2.1,0.,0.,3.,0.,1.,1.5,2.1], + [0.,2.,2.1,0.,0.,4.,0.,2.,3.5,2.5], + 1.,6.,1.,4.0) end diff --git a/test/test_rexport.jl b/test/test_rexport.jl index ccf0f7a..bbae4fd 100644 --- a/test/test_rexport.jl +++ b/test/test_rexport.jl @@ -32,8 +32,8 @@ class(tree2r) = "phylo" @rput net2 @test convert(Bool, R"isTRUE(all.equal(net2, tree2r))") -# network, h=1, some missing gamma values -s = "(((A:4.0,(B:1.0)#H1:1.1::0.9):0.5,(C:0.6,#H1:1.0):1.0):3.0,D:5.0);"; +# network, h=1, missing gamma values +s = "(((A:4.0,(B:1.0)#H1:1.1):0.5,(C:0.6,#H1:1.0):1.0):3.0,D:5.0);"; phy1 = PhyloPlots.rexport(readTopology(s)); net2 = readTopology(s); if useape # needs ape version > 4.1 with read.evonet (not in 4.1) @@ -73,7 +73,9 @@ storage.mode(tree1r$edge) <- "integer" class(tree1r) = "phylo" phy2r = list(edge=matrix(c(5,5,6,6,7,8,8,9, 6,4,8,7,3,1,9,2), 8,2), Nnode=5L, edge.length=c(3,5,.5,1,.6,4,1.1,NA), - reticulation=matrix(c(7L,9L),1,2), tip.label=c("A","B","C","D")) + reticulation=matrix(c(7L,9L),1,2), + reticulation.gamma=0.1, + tip.label=c("A","B","C","D")) storage.mode(phy2r$edge) <- "integer" class(phy2r) = c("evonet","phylo") """;