Added option -r to reroot outgroup + added command collapse single

cmd/collapsesingle.go

@@ -0,0 +1,46 @@
+package cmd
+
+import (
+	"github.com/fredericlemoine/gotree/io"
+	"github.com/spf13/cobra"
+)
+
+// collapseCmd represents the collapse command
+var collapsesingleCmd = &cobra.Command{
+	Use:   "single",
+	Short: "Collapse branches that connect single nodes",
+	Long: `Collapse branches that connect single nodes.
+
+Single nodes are defined as nodes that:
+- Connect only 2 neighbors
+- Are not the root
+
+* Branch lengths are added
+* Max branch support is assigned to remaining branch
+
+Ex:
+           t1           t1
+           /	       /
+ n0--n1--n2   => n0--n2
+           \	       \
+            t2          t2
+
+`,
+	Run: func(cmd *cobra.Command, args []string) {
+		f := openWriteFile(outtreefile)
+		treefile, treechan := readTrees(intreefile)
+		defer treefile.Close()
+		for t := range treechan {
+			if t.Err != nil {
+				io.ExitWithMessage(t.Err)
+			}
+			t.Tree.RemoveSingleNodes()
+			f.WriteString(t.Tree.Newick() + "\n")
+		}
+		f.Close()
+	},
+}
+
+func init() {
+	collapseCmd.AddCommand(collapsesingleCmd)
+}

cmd/outgroup.go

@@ -29,6 +29,8 @@ If the outgroup includes a tip that is not present in the tree,
 this tip will not be taken into account for the reroot. A warning
 will be issued.
 
+If the option -r|--remove-outgroup is given, then the outgroup is
+removed after reroot.
 `,
 	Run: func(cmd *cobra.Command, args []string) {
 		var tips []string
@@ -48,7 +50,7 @@ will be issued.
 			if t.Err != nil {
 				io.ExitWithMessage(t.Err)
 			}
-			err := t.Tree.RerootOutGroup(tips...)
+			err := t.Tree.RerootOutGroup(removeoutgroup, tips...)
 			if err != nil {
 				io.ExitWithMessage(err)
 			}
@@ -63,4 +65,5 @@ will be issued.
 func init() {
 	rerootCmd.AddCommand(outgroupCmd)
 	outgroupCmd.PersistentFlags().StringVarP(&tipfile, "tip-file", "l", "none", "File containing names of tips of the outgroup")
+	outgroupCmd.PersistentFlags().BoolVarP(&removeoutgroup, "remove-outgroup", "r", false, "Removes the outgroup after reroot")
 }

cmd/root.go

@@ -19,6 +19,7 @@ import (
 )
 
 // Variables used in lots of commands
+var inalignfile string
 var intreefile, intree2file, outtreefile string
 var seed int64
 var inputname string
@@ -36,6 +37,7 @@ var treeformat = utils.FORMAT_NEWICK
 var cfgFile string
 var rootCpus int
 var rootInputFormat string
+var removeoutgroup bool
 
 // RootCmd represents the base command when called without any subcommands
 var RootCmd = &cobra.Command{

support/booster.go

@@ -383,6 +383,7 @@ func (supporter *BoosterSupporter) ComputeValue(refTree *tree.Tree, cpu int, edg
 				}
 			}
 		}
+		treeV.Tree.Delete()
 	}
 	return err
 }

test.sh

@@ -164,6 +164,17 @@ gotree generate yuletree -s 10 | gotree collapse depth -m 2 -M 2 | gotree brlen
 diff result expected
 rm -f expected result
 
+# gotree collapse single
+echo "->gotree collapse single"
+cat > test_input <<EOF
+((((A,B)),((C))),(D,(E)));
+EOF
+cat > expected <<EOF
+(((A,B),C),(D,E));
+EOF
+gotree collapse single -i test_input > result
+diff result expected
+rm -f expected result test_input
 
 # gotree compare trees
 echo "->gotree compare trees"

tests/collapse_test.go

@@ -37,3 +37,40 @@ func TestCollapse(t *testing.T) {
 		t.Error(fmt.Sprintf("The sum of branch lengths after collapse is not 14.127 (%f)", sumlen))
 	}
 }
+
+var treestring2 string = "(A:1,(B:1):1,C:1);"
+
+func TestCollapseSingle(t *testing.T) {
+	tree, err2 := newick.NewParser(strings.NewReader(treestring2)).Parse()
+	if err2 != nil {
+		t.Error(err2)
+	}
+	nbranches := len(tree.Edges())
+	nnodes := len(tree.Nodes())
+	sumlen := tree.SumBranchLengths()
+	if nbranches != 4 {
+		t.Error(fmt.Sprintf("The number of edges before collapse is not 4 (%d)", nbranches))
+	}
+	if sumlen != 4.0 {
+		t.Error(fmt.Sprintf("The sum of branch lengths before collapse is not 4.0  (%f)", sumlen))
+	}
+	if nnodes != 5 {
+		t.Error(fmt.Sprintf("The number of nodes before collapse is not 5  (%d)", nnodes))
+	}
+
+	tree.RemoveSingleNodes()
+	nbranches = len(tree.Edges())
+	sumlen = tree.SumBranchLengths()
+	nnodes = len(tree.Nodes())
+
+	if nbranches != 3 {
+		t.Error(fmt.Sprintf("The number of edges after collapse is not 3 (%d)", nbranches))
+	}
+
+	if sumlen != 4.0 {
+		t.Error(fmt.Sprintf("The sum of branch lengths after collapse is not 3.0 (%f)", sumlen))
+	}
+	if nnodes != 4 {
+		t.Error(fmt.Sprintf("The number of nodes after collapse is not 4  (%d)", nnodes))
+	}
+}

tests/reroot_test.go

@@ -1,6 +1,7 @@
 package tests
 
 import (
+	"fmt"
 	"github.com/fredericlemoine/gotree/tree"
 	"testing"
 )
@@ -24,7 +25,7 @@ func TestRootOutgroup(t *testing.T) {
 	tips := tr.Tips()
 
 	for _, tip := range tips {
-		err = clone.RerootOutGroup(tip.Name())
+		err = clone.RerootOutGroup(false, tip.Name())
 		found := false
 		for _, n := range clone.Root().Neigh() {
 			if n.Tip() && n.Name() == tip.Name() {
@@ -48,3 +49,76 @@ func TestRootOutgroup(t *testing.T) {
 		}
 	}
 }
+
+/*
+Generates a 1000 tip ROOTED random tree, then reroot it at each tip
+and compare all bipartitions of the rerooted tree with the original tree
+*/
+func TestReRootOutgroupRemove(t *testing.T) {
+	tr, err := tree.RandomYuleBinaryTree(1000, true)
+	tr.ReinitIndexes()
+	edges := tr.Edges()
+	tips := tr.Tips()
+	nodes := tr.Nodes()
+
+	for _, tip := range tips {
+		clone := tr.Clone()
+		if err != nil {
+			t.Error(err)
+		}
+
+		err = clone.RerootOutGroup(true, tip.Name())
+
+		if !clone.Rooted() {
+			t.Error("Output tree should be rooted")
+		}
+		if len(clone.Edges()) != len(edges)-2 {
+			t.Error(fmt.Sprintf("Rerooted tree should have %d edges and has %d", len(edges)-1, len(clone.Edges())))
+		}
+		if len(clone.Nodes()) != len(nodes)-2 {
+			t.Error(fmt.Sprintf("Rerooted tree should have %d nodes and has %d", len(nodes)-1, len(clone.Nodes())))
+		}
+
+		for _, t2 := range clone.Tips() {
+			if t2.Name() == tip.Name() {
+				t.Error(fmt.Sprintf("Outgroup Tip %s should not be present in the rerooted tree ", t2.Name()))
+			}
+		}
+	}
+}
+
+/*
+Generates a 1000 tip UNROOTED random tree, then reroot it at each tip
+and compare all bipartitions of the rerooted tree with the original tree
+*/
+func TestReRootOutgroupRemoveUnRooted(t *testing.T) {
+	tr, err := tree.RandomYuleBinaryTree(1000, false)
+	tr.ReinitIndexes()
+	edges := tr.Edges()
+	tips := tr.Tips()
+	nodes := tr.Nodes()
+
+	for _, tip := range tips {
+		clone := tr.Clone()
+		if err != nil {
+			t.Error(err)
+		}
+
+		err = clone.RerootOutGroup(true, tip.Name())
+
+		if !clone.Rooted() {
+			t.Error("Output tree should be rooted")
+		}
+		if len(clone.Edges()) != len(edges)-1 {
+			t.Error(fmt.Sprintf("Rerooted tree should have %d edges and has %d", len(edges)-1, len(clone.Edges())))
+		}
+		if len(clone.Nodes()) != len(nodes)-1 {
+			t.Error(fmt.Sprintf("Rerooted tree should have %d nodes and has %d", len(nodes)-1, len(clone.Nodes())))
+		}
+		for _, t2 := range clone.Tips() {
+			if t2.Name() == tip.Name() {
+				t.Error(fmt.Sprintf("Outgroup Tip %s should not be present in the rerooted tree ", t2.Name()))
+			}
+		}
+	}
+}

tests/tree_test.go

@@ -300,3 +300,12 @@ func TestSackin1(t *testing.T) {
 		t.Error(fmt.Sprintf("Sackin index found : %d instead of %d expected", sidx, expected))
 	}
 }
+
+// Test counting the Sackin Index on balanced rooted and unrooted trees
+func TestDelete(t *testing.T) {
+	tr, err := tree.RandomYuleBinaryTree(701, true)
+	if err != nil {
+		t.Error(err)
+	}
+	tr.Delete()
+}

tree/algo.go

@@ -373,7 +373,9 @@ func Consensus(trees <-chan Trees, cutoff float64) (*Tree, error) {
 //
 // If the outgroup includes a tip that is not present in the tree,
 // this tip will not be taken into account for the rerooting.
-func (t *Tree) RerootOutGroup(tips ...string) error {
+//
+// If removeoutgroup is true, then the outgrouped is removed from the rerooted tree.
+func (t *Tree) RerootOutGroup(removeoutgroup bool, tips ...string) error {
 	t.UnRoot()
 
 	n, edges, _, err := t.LeastCommonAncestorUnrooted(nil, tips...)
@@ -407,26 +409,50 @@ func (t *Tree) RerootOutGroup(tips ...string) error {
 		}
 	}
 
-	root := t.NewNode()
-	length := rootedge.Length()
-	support := rootedge.Support()
-
-	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)
-		ne.SetSupport(support)
-		ne2.SetSupport(support)
+	var root *Node
+	// Here we will remove outgroup nodes
+	if removeoutgroup {
+		// We get the new root as the node on the
+		// other side of the edge
+		root = rootedge.Left()
+		if root == n {
+			root = rootedge.Right()
+		}
+		root.delNeighbor(n)
+		n.delNeighbor(root)
+		// We retrieve all nodes of the
+		// subtree we want to remove
+		nodes := make([]*Node, 0, len(tips))
+		t.nodesRecur(&nodes, n, nil)
+		for _, no := range nodes {
+			t.delNode(no)
+		}
+	} else {
+		// Else we add a new root at the middle of the edge
+		// connecting the outgroup to the other subtree
+		root = t.NewNode()
+		length := rootedge.Length()
+		support := rootedge.Support()
+
+		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)
+			ne.SetSupport(support)
+			ne2.SetSupport(support)
+		}
 	}
-
 	t.Reroot(root)
+	if removeoutgroup {
+		t.UpdateTipIndex()
+	}
 	t.ClearBitSets()
 	t.UpdateBitSet()
 	t.ComputeDepths()