Corrected annotate for rooted trees

cmd/annotate.go

@@ -11,6 +11,7 @@ import (
 
 	"github.com/fredericlemoine/gotree/io"
 	"github.com/fredericlemoine/gotree/io/fileutils"
+	"github.com/fredericlemoine/gotree/mutils"
 	"github.com/fredericlemoine/gotree/support"
 	"github.com/fredericlemoine/gotree/tree"
 
@@ -28,7 +29,8 @@ var annotateCmd = &cobra.Command{
 Annotations may be (in order of priority):
 - A tree with labels on internal nodes (-c). in that case, it will label each branch of 
    the input tree with label of the closest branch of the given compared tree (-c) in terms
-   of transfer distance. The labels are of the form: "label_distance_depth";
+   of transfer distance. The labels are of the form: "label_distance_depth"; Only internal branches
+   are annotated, and no internal branch is annotated with a terminal branch.
 - A file with one line per internal node to annotate (-m), and with the following format:
    <name of internal branch/node n1>:<name of taxon n2>,<name of taxon n3>,...,<name of taxon ni>
 	=> If 0 name is given after ':' an error is returned
@@ -127,10 +129,43 @@ If neither -c nor -m are given, gotree annotate will wait for a reference tree o
 						if dist > uint16(len(tips))/2 {
 							dist = uint16(len(tips)) - dist
 						}
-						if annotateComment {
-							e1.Right().AddComment(fmt.Sprintf("%s_%d_%d", e2.Name(true), dist, depth))
-						} else {
-							e1.Right().SetName(fmt.Sprintf("%s_%d_%d", e2.Name(true), dist, depth))
+
+						// If root edge and rooted tree, we take the closest branch
+						if e2.Left().Nneigh() == 2 {
+							var t3, t2, t1 int
+
+							e3 := e2.Left().Edges()[0]
+							if e3 == e2 {
+								e3 = e2.Left().Edges()[1]
+							}
+
+							if t3, err = e3.NumTipsRight(); err != nil {
+								io.LogError(err)
+								return
+							}
+							if t2, err = e2.NumTipsRight(); err != nil {
+								io.LogError(err)
+								return
+							}
+							if t1, err = e1.NumTipsRight(); err != nil {
+								io.LogError(err)
+								return
+							}
+							fmt.Println(t1)
+							fmt.Println(t2)
+							fmt.Println(t3)
+
+							if mutils.Abs(t3-t1) < mutils.Abs(t2-t1) {
+								e2 = e3
+							}
+						}
+
+						if !e2.Right().Tip() {
+							if annotateComment {
+								e1.Right().AddComment(fmt.Sprintf("%s_%d_%d", e2.Name(true), dist, depth))
+							} else {
+								e1.Right().SetName(fmt.Sprintf("%s_%d_%d", e2.Name(true), dist, depth))
+							}
 						}
 					}
 				}

docs/commands/annotate.md

@@ -8,7 +8,8 @@ This command annotates internal nodes/branches of a tree with given information.
 Annotations may be (in order of priority):
 - A tree with labels on internal nodes (-c). in that case, it will label each branch of 
    the input tree with label of the closest branch of the given compared tree (-c) in terms
-   of transfer distance. The labels are of the form: "`label_distance_depth`";
+   of transfer distance. The labels are of the form: "`label_distance_depth`"; Only internal branches
+   are annotated, and no internal branch is annotated with a terminal branch.
 - A file with one line per internal node to annotate (-m), and with the following format:
    `<name of internal branch/node n1>:<name of taxon n2>,<name of taxon n3>,...,<name of taxon ni>`
 	=> If 0 name is given after ':' an error is returned

mutils/utils.go

@@ -0,0 +1,25 @@
+/**
+Package with simple math util functions such as min, max, and abs
+*/
+package mutils
+
+func Min(a, b int) int {
+	if a < b {
+		return a
+	}
+	return b
+}
+
+func Max(a, b int) int {
+	if a < b {
+		return b
+	}
+	return a
+}
+
+func Abs(x int) int {
+	if x < 0 {
+		return -x
+	}
+	return x
+}

test.sh

@@ -68,7 +68,7 @@ cat > intree <<EOF
 ((Tip4,(Tip7,Tip2)n1)n2,Tip0,((Tip8,(Tip9,Tip3)n3)n4,((Tip6,Tip1)n5,Tip5)n6)n7);
 EOF
 cat > expected <<EOF
-((Tip4,(Tip7,Tip2)n1_0_2)n2_0_3,Tip0,((Tip8,(Tip9,Tip3)n3_0_2)n4_0_3,((Tip6,Tip5)Tip6_1_2,Tip1)n6_0_3)n7_0_6);
+((Tip4,(Tip7,Tip2)n1_0_2)n2_0_3,Tip0,((Tip8,(Tip9,Tip3)n3_0_2)n4_0_3,((Tip6,Tip5),Tip1)n6_0_3)n7_0_6);
 EOF
 gotree generate yuletree --seed 10 | gotree brlen clear | gotree annotate -c intree > result
 diff -q -b result expected
@@ -80,7 +80,7 @@ cat > intree <<EOF
 ((Tip4,(Tip7,Tip2)n1)n2,Tip0,((Tip8,(Tip9,Tip3)n3)n4,((Tip6,Tip1)n5,Tip5)n6)n7);
 EOF
 cat > expected <<EOF
-((Tip4,(Tip7,Tip2)[n1_0_2])[n2_0_3],Tip0,((Tip8,(Tip9,Tip3)[n3_0_2])[n4_0_3],((Tip6,Tip5)[Tip6_1_2],Tip1)[n6_0_3])[n7_0_6]);
+((Tip4,(Tip7,Tip2)[n1_0_2])[n2_0_3],Tip0,((Tip8,(Tip9,Tip3)[n3_0_2])[n4_0_3],((Tip6,Tip5),Tip1)[n6_0_3])[n7_0_6]);
 EOF
 gotree generate yuletree --seed 10 | gotree brlen clear | gotree annotate --comment -c intree > result
 diff -q -b result expected
@@ -853,7 +853,7 @@ cat > annotated <<EOF
 ((((((((Gorilla_gorilla_gorilla[subspecies],Pan_troglodytes[species],Homo_sapiens[species])Homo/Pan/Gorilla_group[subfamily],Pongo_pygmaeus_abelii[species])Pongidae[family],Hylobates_pileatus[species])Hominoidea[superfamily],(Macaca_sylvanus[species],Macaca_fascicularis_fascicularis[subspecies],Macaca_fuscata[species],Macaca_mulatta[species])Macaca[genus])Catarrhini[parvorder],Saimiri_sciureus[species])Simiiformes[infraorder],Tarsius_tarsier[species])Haplorrhini[suborder],Lemur_sp.[species])Primates[order],Mus_musculus[species],Bos_taurus[species])Euarchontoglires[superorder];
 EOF
 cat > expected <<EOF
-(((((Hylobates_pileatus:0.23988592,(Pongo_pygmaeus_abelii:0.11809071,(Gorilla_gorilla_gorilla:0.13596645,(Homo_sapiens:0.11344407,Pan_troglodytes:0.11665038)Pan_troglodytes_1_2:0.02364476)Homo/Pan/Gorilla_group_0_3:0.04257513)Pongidae_0_4:0.15711475)Hominoidea_0_5:0.03966791,(Macaca_sylvanus:0.06332916,(Macaca_fascicularis_fascicularis:0.07605049,(Macaca_mulatta:0.06998962,Macaca_fuscata:0)Macaca_fuscata_1_2:0.08492791)Macaca_1_3:0.02236558)Macaca_0_4:0.11208218)Catarrhini_0_9:0.0477543,Saimiri_sciureus:0.25824985)Simiiformes_0_10:0.14311537,(Tarsius_tarsier:0.62272677,Lemur_sp.:0.40249393)Tarsius_tarsier_1_2:0)Primates_0_12:0.077084225,(Mus_musculus:0.4057381,Bos_taurus:0.65776307)Primates_0_2:0.077084225);
+(((((Hylobates_pileatus:0.23988592,(Pongo_pygmaeus_abelii:0.11809071,(Gorilla_gorilla_gorilla:0.13596645,(Homo_sapiens:0.11344407,Pan_troglodytes:0.11665038)0.62:0.02364476)Homo/Pan/Gorilla_group_0_3:0.04257513)Pongidae_0_4:0.15711475)Hominoidea_0_5:0.03966791,(Macaca_sylvanus:0.06332916,(Macaca_fascicularis_fascicularis:0.07605049,(Macaca_mulatta:0.06998962,Macaca_fuscata:0)0.98:0.08492791)Macaca_1_3:0.02236558)Macaca_0_4:0.11208218)Catarrhini_0_9:0.0477543,Saimiri_sciureus:0.25824985)Simiiformes_0_10:0.14311537,(Tarsius_tarsier:0.62272677,Lemur_sp.:0.40249393)0.35:0)Primates_0_12:0.077084225,(Mus_musculus:0.4057381,Bos_taurus:0.65776307)Primates_0_2:0.077084225);
 EOF
 gotree annotate -i inferred -c annotated -o result
 diff -q -b expected result

tree/edge.go

@@ -9,6 +9,7 @@ import (
 
 	"github.com/fredericlemoine/bitset"
 	"github.com/fredericlemoine/gotree/io"
+	"github.com/fredericlemoine/gotree/mutils"
 )
 
 // Structure of an edge
@@ -165,7 +166,7 @@ func (e *Edge) TopoDepth() (int, error) {
 	}
 	count := int(e.bitset.Count())
 	total := int(e.bitset.Len())
-	return min(count, total-count), nil
+	return mutils.Min(count, total-count), nil
 }
 
 // Returns a string representing the bitset (bipartition)
@@ -205,7 +206,7 @@ func (e *Edge) ToStatsString(withedgecomments bool) string {
 	if rightdepth, err = e.Right().Depth(); err != nil {
 		io.ExitWithMessage(err)
 	}
-	depth = min(leftdepth, rightdepth)
+	depth = mutils.Min(leftdepth, rightdepth)
 	var topodepth int
 	topodepth, err = e.TopoDepth()
 	if err != nil {

tree/stats.go

@@ -4,6 +4,8 @@ import (
 	"errors"
 	"math"
 	"sort"
+
+	"github.com/fredericlemoine/gotree/mutils"
 )
 
 const MaxInt = int(^uint(0) >> 1)
@@ -136,7 +138,7 @@ func (t *Tree) CollessIndex() (colless int) {
 		leftindex, lefttips := collessIndexRecur(edge.Left(), edge.Right())
 		rightindex, righttips := collessIndexRecur(edge.Right(), edge.Left())
 		colless += (leftindex + rightindex)
-		colless += max(lefttips, righttips) - min(lefttips, righttips)
+		colless += mutils.Max(lefttips, righttips) - mutils.Min(lefttips, righttips)
 	} else {
 		colless, _ = collessIndexRecur(t.Root(), nil)
 	}
@@ -154,8 +156,8 @@ func collessIndexRecur(n *Node, prev *Node) (colless, tips int) {
 	for _, c := range n.Neigh() {
 		if c != prev {
 			childindex, childtips := collessIndexRecur(c, n)
-			mintips = min(mintips, childtips)
-			maxtips = max(maxtips, childtips)
+			mintips = mutils.Min(mintips, childtips)
+			maxtips = mutils.Max(maxtips, childtips)
 			tips += childtips
 			colless += childindex
 		}

tree/tree.go

@@ -16,6 +16,7 @@ import (
 
 	"github.com/fredericlemoine/bitset"
 	"github.com/fredericlemoine/gotree/io"
+	"github.com/fredericlemoine/gotree/mutils"
 )
 
 // Tree structure having a root and a tip index, that maps tip names to their index
@@ -1580,8 +1581,8 @@ func (t *Tree) deepestEdgeRecur(node, prev *Node, edge *Edge, numtips int) (maxe
 			curtips += t
 		}
 	}
-	if min(numtips-curtips, curtips) > maxdepth {
-		maxdepth = min(numtips-curtips, curtips)
+	if mutils.Min(numtips-curtips, curtips) > maxdepth {
+		maxdepth = mutils.Min(numtips-curtips, curtips)
 		maxedge = edge
 		lefttips = numtips - curtips
 		righttips = curtips