read nexus tree blocks: with translate tables and networks (JuliaPhylo#197)

readNexusTrees changed to readnexus_treeblock.
assumptions: format used by bacter or by SpeciesNetworks BEAST2 packages (which are different from each other)

Author: cecileane

docs/src/lib/public.md

@@ -82,7 +82,7 @@ readTopology
 readTopologyLevel1
 readInputTrees
 readMultiTopology
-readNexusTrees
+readnexus_treeblock
 readSnaqNetwork
 readTrees2CF
 countquartetsintrees

examples/test_reticulatetreeblock.nex

@@ -0,0 +1,32 @@
+#NEXUS
+
+begin taxa;
+dimensions ntax=2;
+taxlabels
+A
+D
+;
+End;
+
+begin trees;
+
+translate
+5 tax5, [to check robustness]
+1 tax1,
+2 tax2,
+3 tax3,
+4 tax4
+;
+
+Tree	bacter=((4,(3,#7[&conv=1, relSize=0.08, affectedBlocks="{bark,dirty,dry}"]:0.001)[&height_95%_HPD={38.01,87.51}]:0.3):0.6,(2,(1:0.1)#7:0.9):10);
+tree gt1 = ((4,(3,#7[&conv=7,relSize=0.58, posterior=0.87]:0.001::0.57):0.3):0.6,(2,(1)#7));
+tree fake = this line won't match a newick string
+tree gt2 = (error will be caught and warning sent;
+tree speciesnetwork = (((1:1.13,((2:0.21)#H1:0.89,(3:1.03,(#H1[&gamma=0.28]:0.30,4:0.51)S3:0.51)S4:0.08)S5:0.2):0.6,5:1.14):0.16);
+
+end;
+
+bacter example: https://raw.githubusercontent.com/taming-the-beast/Bacter-Tutorial/master/precooked_runs/bacter_tutorial.Ecoli_rplA.trees
+SpeciesNetwork examples:
+MCMC sample: https://raw.githubusercontent.com/mmatschiner/tutorials/master/bayesian_analysis_of_species_networks/res/speciesnetwork.trees
+summary: https://raw.githubusercontent.com/mmatschiner/tutorials/master/bayesian_analysis_of_species_networks/res/speciesnetwork_sum.trees

src/PhyloNetworks.jl

@@ -66,7 +66,7 @@ module PhyloNetworks
         writeTableCF,
         mapAllelesCFtable,
         readInputTrees,
-        readNexusTrees,
+        readnexus_treeblock,
         summarizeDataCF,
         snaq!,
         readSnaqNetwork,

src/deprecated.jl

@@ -1,6 +1,5 @@
-@deprecate phyloNetworklm phylolm
-@deprecate sigma2_estim sigma2_phylo
-@deprecate mu_estim mu_phylo
+@deprecate readNexusTrees readnexus_treeblock
+
 @deprecate getChild getchild false
 @deprecate getChildren getchildren false
 @deprecate getChildEdge getchildedge false

src/readData.jl

@@ -963,7 +963,9 @@ function readTrees2CF(treefile::AbstractString; quartetfile="none"::AbstractStri
                       writeTab=true::Bool, CFfile="none"::AbstractString,
                       taxa::AbstractVector=Vector{String}(),
                       writeQ=false::Bool, writeSummary=true::Bool, nexus=false::Bool)
-    trees = (nexus ? readNexusTrees(treefile, readTopologyUpdate, false, false) : readInputTrees(treefile))
+    trees = (nexus ?
+             readnexus_treeblock(treefile, readTopologyUpdate, false, false; reticulate=false) :
+             readInputTrees(treefile))
     if length(taxa)==0        # unionTaxa(trees) NOT default argument:
       taxa = unionTaxa(trees) # otherwise: tree file is read twice
     end
@@ -1275,49 +1277,3 @@ function createQuartet(taxa::Union{Vector{<:AbstractString},Vector{Int}}, qvec::
     end
     return Quartet(num,names,[1.0,0.0,0.0])
 end
-
-"""
-    readNexusTrees(filename::AbstractString, treereader=readTopology::Function [, args...])
-
-Read trees in nexus-formatted file and return a vector of `HybridNetwork`s.
-For the nexus format, see Maddison, Swofford & Maddison (1997)
-https://doi.org/10.1093/sysbio/46.4.590.
-The optional arguments are passed onto the individual tree reader.
-
-Warnings:
-- "translate" tables are not supported yet
-- only the first tree block is read
-"""
-function readNexusTrees(file::AbstractString, treereader=readTopology::Function, args...)
-    vnet = HybridNetwork[]
-    rx_start = r"^\s*begin\s+trees\s*;"i
-    rx_end = r"^\s*end\s*;"i
-    rx_tree = r"^\s*tree\s+[^(]+(\([^;]*;)"i
-    # spaces,"Tree",spaces,any_symbols_other_than_(, then we capture:
-    # ( any_symbols_other_than_; ;
-    treeblock = false # whether we are currently reading the TREE block or not
-    open(file) do s
-        numl = 0
-        for line in eachline(s)
-            numl += 1
-            if treeblock # currently reading trees, check for END signal
-                occursin(rx_end, line) && break # break if end of tree block
-            else # not reading trees: check for the BEGIN signal
-                if occursin(rx_start, line) treeblock=true; end
-                continue # to next line, either way
-            end
-            # if we get there, it's that we are inside the treeblock (true) and no END signal yet
-            m = match(rx_tree, line)
-            m != nothing || continue # continue to next line if no match
-            phy = m.captures[1] # string
-            try
-                push!(vnet, treereader(phy, args...)) # readTopologyUpdate(phy,false)
-            catch err
-                print("skipped phylogeny on line $(numl) of file $file: ")
-                if :msg in fieldnames(typeof(err)) println(err.msg); else println(typeof(err)); end
-            end
-        end
-    end
-    return vnet # consistent output type: HybridNetwork vector. might be of length 0.
-end
-

src/readwrite.jl

@@ -1371,6 +1371,207 @@ function readMultiTopology(file::AbstractString, fast::Bool=true)
     return vnet
 end
 
+
+@doc raw"""
+    readnexus_treeblock(filename, treereader=readTopology, args...;
+                        reticulate=true, stringmodifier=[r"#(\d+)" => s"#H\1"])
+
+Read the *first* "trees" block of a nexus-formatted file, using the translate
+table if present, and return a vector of `HybridNetwork`s.
+Information inside `[&...]` are interpreted as comments and are discarded by the
+default tree reader. Optional arguments `args` are passed to the tree reader.
+
+For the nexus format, see
+[Maddison, Swofford & Maddison (1997)](https://doi.org/10.1093/sysbio/46.4.590).
+
+Unless `reticulate` is false, the following is done to read networks with reticulations.
+
+Prior to reading each phylogeny, each instance of `#number` is replaced by
+`#Hnumber` to fit the standard extended Newick format at hybrid nodes.
+This behavior can be changed with option `stringmodifier`, which should be a
+vector of pairs accepted by `replace`.
+
+Inheritance γ values are assumed to be given within "comment" blocks at *minor*
+hybrid edges (cut as tips to form the extended Newick) like this for example,
+as output by bacter ([Vaughan et al. 2017](http://dx.doi.org/10.1534/genetics.116.193425)):
+
+    #11[&conv=0, relSize=0.08, ...
+
+or like this, as output by SpeciesNetwork
+([Zhang et al. 2018](https://doi.org/10.1093/molbev/msx307)):
+
+    #H11[&gamma=0.08]
+
+In this example, the corresponding edge to hybrid H11 has γ=0.08.
+"""
+function readnexus_treeblock(file::AbstractString, treereader=readTopology::Function, args...;
+            reticulate=true, stringmodifier=[r"#(\d+)\b" => s"#H\1"]) # add H
+    vnet = HybridNetwork[]
+    rx_start = r"^\s*begin\s+trees\s*;"i
+    rx_end = r"^\s*end\s*;"i
+    rx_tree = r"^\s*tree\s+[^(]+(\([^;]*;)"i
+    treeblock = false
+    translate = false
+    id2name = nothing
+    open(file) do s
+        numl = 0
+        for line in eachline(s)
+            numl += 1
+            if treeblock
+                occursin(rx_end, line) && break # exit if end of tree block
+            elseif occursin(rx_start, line)     # start reading tree block
+                    treeblock = true
+                    line, translate, id2name = readnexus_translatetable(s)
+            else continue
+            end
+            # now we are inside the treeblock
+            m = match(rx_tree, line)
+            isnothing(m) && continue
+            phy = m.captures[1]
+            if reticulate # fix #Hn and extract γ from #Hn[&conv=n, relSize=γ]
+                phy = replace(phy, stringmodifier...)
+                id2gamma = readnexus_extractgamma(phy)
+            end
+            net = nothing
+            try
+                net = treereader(phy, args...)
+            catch err
+                warnmsg = "skipped phylogeny on line $(numl) of file\n$file\n" *
+                    (:msg in fieldnames(typeof(err)) ? err.msg : string(typeof(err)))
+                @warn warnmsg
+                continue # don't push to vnet
+            end
+            reticulate && readnexus_assigngammas!(net, id2gamma)
+            if translate
+                for tip in net.leaf
+                    id = parse(Int, tip.name)
+                    tip.name = id2name[id]
+                end
+            end
+            push!(vnet, net)
+        end
+    end
+    return vnet
+end
+
+"""
+    readnexus_translatetable(io)
+
+Read translate table from IO object `io`, whose first non-empty line should contain
+"translate". Then each line should have "number name" and the end of the table
+is indicated by a ;. Output tuple:
+- line that was last read, and is not part of the translate table, taken from `io`
+- translate: boolean, whether a table was successfully read
+- id2name: dictionary mapping number to name.
+"""
+function readnexus_translatetable(io)
+    rx_translate = r"^\s*translate"i
+    rx_emptyline = r"^\s*$"
+    line = readline(io)
+    translate = false
+    id2name = Dict{Int,String}()
+    while true
+        if occursin(rx_translate, line)
+            translate = true
+            break
+        elseif occursin(rx_emptyline, line)
+            line = readline(io)
+        else
+            translate = false
+            break
+        end
+    end
+    if translate # then read the table
+        rx_end = r"^\s*;"
+        rx_idname = r"\s*(\d+)\s+(\w+)\s*([,;]?)"
+        while true
+            line = readline(io)
+            occursin(rx_emptyline, line) && continue
+            if occursin(rx_end, line)
+                line = readline(io)
+                break
+            end
+            m = match(rx_idname, line)
+            if isnothing(m)
+                @warn "problem reading the translate table at line $line.\nnumbers won't be translated to names"
+                translate = false
+                break
+            end
+            push!(id2name, parse(Int,m.captures[1]) => String(m.captures[2]))
+            if m.captures[3] == ";"
+                line = readline(io)
+                break
+            end
+        end
+    end
+    return line, translate, id2name
+end
+
+"""
+    readnexus_extractgamma(nexus_string)
+
+Extract γ from comments and return a dictionary hybrid number ID => γ, from
+one single phylogeny given as a string.
+The output from BEAST2 uses this format for reticulations at *minor* edges,
+as output by bacter ([Vaughan et al. 2017](http://dx.doi.org/10.1534/genetics.116.193425)):
+
+    #11[&conv=0, relSize=0.08, ...
+
+or as output by SpeciesNetwork ([Zhang et al. 2018](https://doi.org/10.1093/molbev/msx307)):
+
+    #H11[&gamma=0.08]
+
+The function below assumes that the "H" was already added back if not present
+already (from bacter), like this:
+
+    #H11[&conv=0, relSize=0.19, ...
+
+The bacter format is tried first. If this format doesn't give any match,
+then the SpeciesNetwork format is tried next.  
+See [`readnexus_assigngammas!`](@ref).
+"""
+function readnexus_extractgamma(nexstring)
+    rx_gamma_v1 = r"#H(\d+)\[&conv=\d+,\s*relSize=(\d+\.\d+)"
+    rx_gamma_v2 = r"#H(\d+)\[&gamma=(\d+\.\d+)"
+    id2gamma = Dict{Int,Float64}()
+    # first: try format v1
+    for m in eachmatch(rx_gamma_v1, nexstring)
+        push!(id2gamma, parse(Int, m.captures[1]) => parse(Float64,m.captures[2]))
+    end
+    if isempty(id2gamma) # then try format v2
+      for m in eachmatch(rx_gamma_v2, nexstring)
+        push!(id2gamma, parse(Int, m.captures[1]) => parse(Float64,m.captures[2]))
+      end
+    end
+    return id2gamma
+end
+
+"""
+    readnexus_assigngammas!(net, d::Dict)
+
+Assign d[i] as the `.gamma` value of the minor parent edge of hybrid "Hi",
+if this hybrid node name is found, and if its minor parent doesn't already
+have a non-missing γ. See [`readnexus_extractgamma`](@ref)
+"""
+function readnexus_assigngammas!(net::HybridNetwork, id2gamma::Dict)
+    for (i,gam) in id2gamma
+        nam = "H$i"
+        j = findfirst(n -> n.name == nam, net.hybrid)
+        if isnothing(j)
+            @warn "didn't find any hybrid node named $nam."
+            continue
+        end
+        hn = net.hybrid[j]
+        he = getparentedgeminor(hn)
+        if he.gamma == -1.0
+            setGamma!(he, gam)
+        else
+            @warn "hybrid edge number $(he.number) has γ=$(he.gamma). won't erase with $gam."
+        end
+    end
+    return net
+end
+
 """
     writeMultiTopology(nets, file_name; append=false)
     writeMultiTopology(nets, IO)

test/test_lm.jl

@@ -387,9 +387,9 @@ phynetlm = phylolm(@formula(trait ~ pred), dfr, net; reml=false)
 @test bic(phynetlm) ≈ bic(fit_mat)
 
 # Deprecated methods
-@test phynetlm.model === phynetlm
-@test phynetlm.mf.f == formula(phynetlm)
-@test phynetlm.mm.m == modelmatrix(phynetlm)
+@test (@test_logs (:warn,r"^accessing") phynetlm.model) === phynetlm
+@test (@test_logs (:warn,r"^accessing") phynetlm.mf.f) == formula(phynetlm)
+@test (@test_logs (:warn,r"^accessing") phynetlm.mm.m) == modelmatrix(phynetlm)
 
 # unordered data
 dfr = dfr[[2,6,10,5,12,7,4,11,1,8,3,9], :]

test/test_readInputData.jl

@@ -27,13 +27,17 @@
 
 
 @testset "test: reading nexus file" begin
+# with translate table, hybrids, failed γ in net 2, bad net 3, v2 format in net 4
+nexusfile = joinpath(@__DIR__, "..", "examples", "test_reticulatetreeblock.nex")
+# nexusfile = joinpath(dirname(pathof(PhyloNetworks)), "..","examples","test_reticulatetreeblock.nex")
+vnet = (@test_logs (:warn, r"^hybrid edge") (:warn,r"^skipped") readnexus_treeblock(nexusfile));
+@test length(vnet) == 3
+@test writeTopology(vnet[1]) == "((tax4,(tax3,#H7:0.001::0.08):0.3):0.6,(tax2,(tax1:0.1)#H7:0.9::0.92):10.0);"
+@test writeTopology(vnet[3]) == "((tax1:1.13,((tax2:0.21)#H1:0.89::0.72,(tax3:1.03,(#H1:0.3::0.28,tax4:0.51)S3:0.51)S4:0.08)S5:0.2):0.6,tax5:1.14);"
+# example without translate table and without reticulations
 nexusfile = joinpath(@__DIR__, "..", "examples", "test.nex")
 # nexusfile = joinpath(dirname(pathof(PhyloNetworks)), "..","examples","test.nex")
-vnet = readNexusTrees(nexusfile);
-@test length(vnet) == 10
-@test length(vnet[10].edge) == 10
-@test vnet[10].edge[7].length ≈ 0.00035
-vnet = readNexusTrees(nexusfile, PhyloNetworks.readTopologyUpdate, false, false);
+vnet = readnexus_treeblock(nexusfile, PhyloNetworks.readTopologyUpdate, false, false; reticulate=false);
 @test length(vnet) == 10
 @test length(vnet[10].edge) == 9
 @test vnet[10].edge[7].length ≈ 0.00035