added compare_tree command to ghost-tree script and did unit tests

ghosttree/tests/test_util.py

@@ -0,0 +1,37 @@
+import unittest
+from StringIO import StringIO
+
+import numpy as np
+
+from ghosttree.util import compare_tip_to_tip_distances
+
+
+# name of class corresponds to the function we are testing
+# we test library code NOT command line so we use the library code function
+# names
+class TestCompareTipToTipDistances(unittest.TestCase):
+    def setUp(self):
+        self.tree1 = StringIO(tree1)
+        self.tree1_copy = StringIO(tree1)
+        self.tree2 = StringIO(tree2)
+
+    def test_same_trees(self):
+        np.random.seed(1234)  # remove randomness for sake of testing p-value
+        result = compare_tip_to_tip_distances(self.tree1, self.tree1_copy)
+        self.assertEqual(result, (1.0, 0.042, 4))
+
+    def test_different_trees(self):
+        np.random.seed(1234)  # remove randomness for sake of testing p-value
+        coeff, p_value, n = compare_tip_to_tip_distances(self.tree1,
+                                                         self.tree2)
+        self.assertAlmostEqual(coeff, 0.59603956067926978)
+        self.assertAlmostEqual(p_value, 0.69599999999999995)
+        self.assertEqual(n, 3)
+
+tree1 = "[example](a:0.1, 'b_b''':0.2, (c:0.3, d_d:0.4)e:0.5)f:0.0;"
+tree2 = "[example](a:0.1, 'b_b''':0.2, (g:0.3, d_d:0.1)e:0.1)f:0.0;"
+
+# lets you run tests in this file individually, if you want to test one
+# piece at a time
+if __name__ == "__main__":
+    unittest.main()

ghosttree/util.py

@@ -0,0 +1,12 @@
+from skbio import TreeNode
+from skbio.stats.distance import mantel
+
+
+def compare_tip_to_tip_distances(tree_fh1, tree_fh2, method="pearson"):
+    tree1 = TreeNode.read(tree_fh1)
+    tree2 = TreeNode.read(tree_fh2)
+
+    dm1 = tree1.tip_tip_distances()
+    dm2 = tree2.tip_tip_distances()
+
+    return mantel(dm1, dm2, strict=False, method=method)

scripts/ghost-tree

@@ -1,15 +1,35 @@
 #!/usr/bin/env python
+# third party dependencies (i.e. packages not included with python)
 import click
-
 import skbio
 
+# library code imports
 from ghosttree.silva.filter import fungi_from_fasta
+from ghosttree.util import compare_tip_to_tip_distances
 
 
 @click.group()
 def cli():
     pass
 
+# click command and args here are for the user
+# cli group is the "invisible" top level group
+@cli.command("compare-trees")
+@click.argument("tree-file1", type=click.File("U"))
+@click.argument("tree-file2", type=click.File("U"))
+@click.option("-m", "--method", default="pearson",
+              type=click.Choice(["pearson", "spearman"]),
+              help="correlation method to use in Mantel test")
+def compare_trees(tree_file1, tree_file2, method):
+    coeff, p_value, n = compare_tip_to_tip_distances(tree_file1,
+                                                     tree_file2,
+                                                     method)
+    click.echo("Correlation coefficient: %f" % coeff)
+    click.echo("p-value: %f" % p_value)
+    click.echo("Number of overlapping tips: %d" % n)
+
+# silva is nested under cli group (silva shows up on command line
+# interface but cli does not
 @cli.group()
 def silva():
     pass