Added reroot command and Corrected tests

cmd/reroot.go

@@ -0,0 +1,104 @@
+package cmd
+
+import (
+	"bufio"
+	"compress/gzip"
+	"github.com/fredericlemoine/gotree/io"
+	"github.com/fredericlemoine/gotree/io/utils"
+	"github.com/fredericlemoine/gotree/tree"
+	"github.com/spf13/cobra"
+	"os"
+	"strings"
+)
+
+var reroottipfile string
+var rerootinputfile string
+var rerootoutputfile string
+
+// rerootCmd represents the reroot command
+var rerootCmd = &cobra.Command{
+	Use:   "reroot",
+	Short: "Reroot the tree using an outgroup",
+	Long: `Reroot the tree using an outgroup given in argument or in stdin.
+
+Example:
+
+
+`,
+	Run: func(cmd *cobra.Command, args []string) {
+		tips := parseTipsFile(reroottipfile)
+
+		var err error
+		var nbtrees int
+
+		compareChannel := make(chan tree.Trees, 15)
+
+		go func() {
+			if nbtrees, err = utils.ReadCompTrees(rerootinputfile, compareChannel); err != nil {
+				io.ExitWithMessage(err)
+			}
+		}()
+
+		var f *os.File
+		if rerootoutputfile != "stdout" {
+			f, err = os.Create(rerootoutputfile)
+		} else {
+			f = os.Stdout
+		}
+		if err != nil {
+			io.ExitWithMessage(err)
+		}
+
+		for t2 := range compareChannel {
+			err = t2.Tree.RerootOutGroup(tips...)
+			if err != nil {
+				io.ExitWithMessage(err)
+			}
+
+			f.WriteString(t2.Tree.Newick() + "\n")
+		}
+
+		f.Close()
+	},
+}
+
+func init() {
+	RootCmd.AddCommand(rerootCmd)
+	rerootCmd.PersistentFlags().StringVarP(&reroottipfile, "tip-file", "l", "stdin", "File containing names of tips of the outgroup")
+	rerootCmd.PersistentFlags().StringVarP(&rerootinputfile, "input", "i", "stdin", "Input Tree")
+	rerootCmd.PersistentFlags().StringVarP(&rerootoutputfile, "output", "o", "stdout", "Rerooted output tree file")
+}
+
+func parseTipsFile(file string) []string {
+	var f *os.File
+	var r *bufio.Reader
+	tips := make([]string, 0, 100)
+	var err error
+	if file == "stdin" || file == "-" {
+		f = os.Stdin
+	} else {
+		f, err = os.Open(file)
+		if err != nil {
+			io.ExitWithMessage(err)
+		}
+	}
+
+	if strings.HasSuffix(file, ".gz") {
+		if gr, err := gzip.NewReader(f); err != nil {
+			io.ExitWithMessage(err)
+		} else {
+			r = bufio.NewReader(gr)
+		}
+	} else {
+		r = bufio.NewReader(f)
+	}
+
+	l, e := Readln(r)
+	for e == nil {
+		for _, name := range strings.Split(l, ",") {
+			tips = append(tips, name)
+		}
+		l, e = Readln(r)
+	}
+	return tips
+}

cmd/root.go

@@ -1,6 +1,7 @@
 package cmd
 
 import (
+	"bufio"
 	"fmt"
 	"github.com/fredericlemoine/gotree/io"
 	"github.com/spf13/cobra"
@@ -59,3 +60,20 @@ func openWriteFile(file string) *os.File {
 		return f
 	}
 }
+
+// Readln returns a single line (without the ending \n)
+// from the input buffered reader.
+// An error is returned iff there is an error with the
+// buffered reader.
+func Readln(r *bufio.Reader) (string, error) {
+	var (
+		isPrefix bool  = true
+		err      error = nil
+		line, ln []byte
+	)
+	for isPrefix && err == nil {
+		line, isPrefix, err = r.ReadLine()
+		ln = append(ln, line...)
+	}
+	return string(ln), err
+}

tests/clone_test.go

@@ -0,0 +1,44 @@
+package tests
+
+import (
+	"github.com/fredericlemoine/gotree/tree"
+	"testing"
+)
+
+/*
+Generates 100 random 1000 tip trees, clone them, and compare them to the original trees
+*/
+func TestCloneTree(t *testing.T) {
+	for i := 0; i < 100; i++ {
+		tr, err := tree.RandomYuleBinaryTree(1000, true)
+		clone := tr.Clone()
+
+		// Comparing tip names
+		tips := tr.Tips()
+		copyTips := clone.Tips()
+		for i, _ := range tips {
+			if tips[i].Name() != copyTips[i].Name() {
+				t.Error("A tip is not found in the cloned tree")
+			}
+		}
+
+		// Check wether the 2 trees have the same set of tip names
+		if err = tr.CompareTipIndexes(clone); err != nil {
+			t.Error(err)
+		}
+
+		// Comparing edges
+		edges := tr.Edges()
+		edges2 := clone.Edges()
+		index := tree.NewEdgeIndex(int64(len(edges)*2), 0.75)
+		for i, e := range edges {
+			index.PutEdgeValue(e, i, e.Length())
+		}
+		for _, e2 := range edges2 {
+			_, ok := index.Value(e2)
+			if !ok {
+				t.Error("An edge of the original tree is not found in the cloned tree")
+			}
+		}
+	}
+}

tests/consensus_test.go

@@ -74,13 +74,13 @@ func TestConsensus2(t *testing.T) {
 	var randtree1, randtree2, randtree3 *tree.Tree
 	var err error
 
-	if randtree1, err = tree.RandomBinaryTree(1000); err != nil {
+	if randtree1, err = tree.RandomUniformBinaryTree(1000, false); err != nil {
 		t.Error(err)
 	}
-	if randtree2, err = tree.RandomBinaryTree(1000); err != nil {
+	if randtree2, err = tree.RandomUniformBinaryTree(1000, false); err != nil {
 		t.Error(err)
 	}
-	if randtree3, err = tree.RandomBinaryTree(1000); err != nil {
+	if randtree3, err = tree.RandomUniformBinaryTree(1000, false); err != nil {
 		t.Error(err)
 	}
 

tests/gentree_test.go

@@ -11,7 +11,7 @@ func benchmarkBinaryTreeGeneration(nbtips int, b *testing.B) {
 	var t *tree.Tree
 	for n := 0; n < b.N; n++ {
 		var err error
-		t, err = tree.RandomBinaryTree(nbtips)
+		t, err = tree.RandomUniformBinaryTree(nbtips, false)
 		if err != nil {
 			b.Error(err)
 		}

tests/parsetree_test.go

@@ -14,7 +14,7 @@ func benchmarkTreeParse(nbtips int, b *testing.B) {
 	var t *tree.Tree
 	for n := 0; n < b.N; n++ {
 		var err error
-		t, err = tree.RandomBinaryTree(nbtips)
+		t, err = tree.RandomUniformBinaryTree(nbtips, false)
 		if err != nil {
 			b.Error(err)
 		}

tests/reroot_test.go

@@ -0,0 +1,50 @@
+package tests
+
+import (
+	"github.com/fredericlemoine/gotree/tree"
+	"testing"
+)
+
+/*
+Generates a 1000 tip random tree, then reroot it at each tip
+and compare all bipartitions of the rerooted tree with the original tree
+*/
+func TestRootOutgroup(t *testing.T) {
+	tr, err := tree.RandomYuleBinaryTree(1000, true)
+	clone := tr.Clone()
+	if err != nil {
+		t.Error(err)
+	}
+
+	edges := tr.Edges()
+	index := tree.NewEdgeIndex(int64(len(edges)*2), 0.75)
+	for i, e := range edges {
+		index.PutEdgeValue(e, i, e.Length())
+	}
+	tips := tr.Tips()
+
+	for _, tip := range tips {
+		err = clone.RerootOutGroup(tip.Name())
+		found := false
+		for _, n := range clone.Root().Neigh() {
+			if n.Tip() && n.Name() == tip.Name() {
+				found = true
+			}
+		}
+		if !found {
+			t.Error("Outgroup (tip) not found in the children of the root on the rerooted tree")
+		}
+		edges2 := clone.Edges()
+		// Check wether the 2 trees have the same set of tip names
+		if err = tr.CompareTipIndexes(clone); err != nil {
+			t.Error(err)
+		}
+
+		for _, e2 := range edges2 {
+			_, ok := index.Value(e2)
+			if !ok {
+				t.Error("An edge of the original tree is not found in the rerooted tree")
+			}
+		}
+	}
+}

tree/algo.go

@@ -297,3 +297,62 @@ func Consensus(trees <-chan Trees, cutoff float64) *Tree {
 
 	return startree
 }
+
+/*
+This function first unroot the input tree and reroot it using the outgroup in argument
+if the outgroup is not monophyletic, it will take all the descendant of the LCA.
+An error is triggered if the LCA is multifurcated, and several branches are possible
+for the placement of the root.
+*/
+func (t *Tree) RerootOutGroup(tips ...string) error {
+	t.UnRoot()
+
+	n, edges, _ := t.LeastCommonAncestorUnrooted(nil, tips...)
+	var rootedge *Edge
+
+	if len(n.br) == 1 {
+		rootedge = n.br[0]
+	} else {
+		if len(n.br)-len(edges) != 1 {
+			return errors.New("Reroot error: Several possible branches for root placement (multifurcated node)")
+		}
+		/**
+		We search the unique branch connecting "n" and that is not part of the outgroup
+		If there were several branches, there should have been an error above
+		*/
+		for _, e := range n.br {
+			found := false
+			for _, e2 := range edges {
+				if e == e2 {
+					found = true
+					break
+				}
+			}
+			if !found {
+				rootedge = e
+				break
+			}
+		}
+	}
+
+	root := t.NewNode()
+	length := rootedge.Length()
+	lnode := rootedge.Left()
+	rnode := rootedge.Right()
+	lnode.delNeighbor(rnode)
+	rnode.delNeighbor(lnode)
+
+	ne := t.ConnectNodes(root, lnode)
+	ne2 := t.ConnectNodes(root, rnode)
+
+	if length > 0 {
+		ne.SetLength(length / 2.0)
+		ne2.SetLength(length / 2.0)
+	}
+
+	t.Reroot(root)
+	t.ClearBitSets()
+	t.UpdateBitSet()
+	t.ComputeDepths()
+	return nil
+}

tree/tree.go

@@ -1212,3 +1212,50 @@ func (t *Tree) MedianSupport() float64 {
 	}
 	return result
 }
+
+// Copy attributes of the node into a new node
+func (t *Tree) CopyNode(n *Node) *Node {
+	out := t.NewNode()
+	out.name = n.name
+	out.depth = n.depth
+	out.id = n.id
+	out.comment = make([]string, len(n.comment))
+	for i, c := range n.comment {
+		out.comment[i] = c
+	}
+	return out
+}
+
+// Copy attributes of the given edge to the other given edge
+func (t *Tree) CopyEdge(e *Edge, copy *Edge) {
+	copy.length = e.length
+	copy.support = e.support
+	copy.pvalue = e.pvalue
+	copy.id = e.id
+	copy.bitset = e.bitset.Clone()
+}
+
+// Clone the input tree
+func (t *Tree) Clone() *Tree {
+	copy := NewTree()
+	root := t.CopyNode(t.Root())
+	copy.SetRoot(root)
+	for _, e := range t.Root().br {
+		t.copyTreeRecur(copy, root, t.Root(), e)
+	}
+	copy.UpdateTipIndex()
+	return (copy)
+}
+
+// Recursive function to clone the tree
+func (t *Tree) copyTreeRecur(copytree *Tree, copynode, node *Node, edge *Edge) {
+	child := edge.Right()
+	copychild := t.CopyNode(child)
+	copyedge := copytree.ConnectNodes(copynode, copychild)
+	t.CopyEdge(edge, copyedge)
+	for _, e := range child.br {
+		if e != edge {
+			t.copyTreeRecur(copytree, copychild, child, e)
+		}
+	}
+}