Merge lp:~jon-hill/supertree-toolkit/fill_in_taxa into lp:supertree-toolkit

[Jon Hill]

Adds functionality to fill in a tree using taxonomy

[Jon Hill (jon-hill)]

Missing test for load taxonomy

[lp:~jon-hill/supertree-toolkit/fill_in_taxa]

272. By Jon Hill on 2015-03-26

Adding test for load_taxonomy and fixing broken test

[Jon Hill (jon-hill)]

review: Approve

=== modified file 'stk/p4/Tree_muck.py'
--- stk/p4/Tree_muck.py	2012-01-11 08:57:43 +0000
+++ stk/p4/Tree_muck.py	2015-03-26 09:59:01 +0000
@@ -769,8 +769,8 @@
     else:
         gm.append("The 2 specified nodes should have a parent-child relationship")
         raise Glitch, gm
-
-    self.deleteCStuff()
+    if var.usePfAndNumpy:
+        self.deleteCStuff()
 
     hasBrLens = False
     for n in self.iterNodes():

=== added file 'stk/scripts/fill_in_with_taxonomy.py'
--- stk/scripts/fill_in_with_taxonomy.py	1970-01-01 00:00:00 +0000
+++ stk/scripts/fill_in_with_taxonomy.py	2015-03-26 09:59:01 +0000
@@ -0,0 +1,408 @@
+#!/usr/bin/env python
+#
+#    Supertree Toolkit. Software for managing and manipulating sources
+#    trees ready for supretree construction.
+#    Copyright (C) 2015, Jon Hill, Katie Davis
+#
+#    This program is free software: you can redistribute it and/or modify
+#    it under the terms of the GNU General Public License as published by
+#    the Free Software Foundation, either version 3 of the License, or
+#    (at your option) any later version.
+#
+#    This program is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#    GNU General Public License for more details.
+#
+#    You should have received a copy of the GNU General Public License
+#    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+#    Jon Hill. [e-mail]
+
+import urllib2
+from urllib import quote_plus
+import simplejson as json
+import argparse
+import os
+import sys
+stk_path = os.path.join( os.path.realpath(os.path.dirname(__file__)), os.pardir )
+sys.path.insert(0, stk_path)
+import supertree_toolkit as stk
+import csv
+
+# What we get from EOL
+current_taxonomy_levels = ['species','genus','family','order','class','phylum','kingdom']
+# And the extra ones from ITIS
+extra_taxonomy_levels = ['superfamily','infraorder','suborder','superorder','subclass','subphylum','superphylum','infrakingdom','subkingdom']
+# all of them in order
+taxonomy_levels = ['species','genus','subfamily','family','superfamily','infraorder','suborder','order','superorder','subclass','class','subphylum','phylum','superphylum','infrakingdom','subkingdom','kingdom']
+
+def get_tree_taxa_taxonomy(taxon,wsdlObjectWoRMS):
+
+    taxon_data = wsdlObjectWoRMS.getAphiaRecords(taxon.replace('_',' '))
+    if taxon_data == None:
+        return {}
+
+    taxon_id = taxon_data[0]['valid_AphiaID'] # there might be records that aren't valid - they point to the valid one though
+    # call it again via the ID this time to make sure we've got the right one.
+    taxon_data = wsdlObjectWoRMS.getAphiaRecordByID(taxon_id)
+    # add data to taxonomy dictionary
+    # get the taxonomy of this species
+    classification = wsdlObjectWoRMS.getAphiaClassificationByID(taxon_id)
+    # construct array
+    tax_array = {}
+    # classification is a nested dictionary, so we need to iterate down it
+    current_child = classification.child
+    while True:
+        tax_array[current_child.rank.lower()] = current_child.scientificname
+        current_child = current_child.child
+        if current_child == '': # empty one is a string for some reason
+            break
+    return tax_array
+
+
+
+def get_taxonomy_worms(taxonomy, start_otu):
+    """ Gets and processes a taxon from the queue to get its taxonomy."""
+    from SOAPpy import WSDL    
+
+    wsdlObjectWoRMS = WSDL.Proxy('http://www.marinespecies.org/aphia.php?p=soap&wsdl=1')
+
+    # this is the recursive function
+    def get_children(taxonomy, ID):
+
+        # get data
+        this_item = wsdlObjectWoRMS.getAphiaRecordByID(ID)
+        if this_item == None:
+            return taxonomy
+        if this_item['rank'].lower() == 'species':
+            # add data to taxonomy dictionary
+            # get the taxonomy of this species
+            classification = wsdlObjectWoRMS.getAphiaClassificationByID(ID)
+            taxon = this_item.scientificname
+            if not taxon in taxonomy: # is a new taxon, not previously in the taxonomy
+                # construct array
+                tax_array = {}
+                # classification is a nested dictionary, so we need to iterate down it
+                current_child = classification.child
+                while True:
+                    if taxonomy_levels.index(current_child.rank.lower()) <= taxonomy_levels.index(start_taxonomy_level):
+                        # we need this - we're closer to the tips of the tree than we started                    
+                        tax_array[current_child.rank.lower()] = current_child.scientificname
+                    current_child = current_child.child
+                    if current_child == '': # empty one is a string for some reason
+                        break
+                taxonomy[this_item.scientificname] = tax_array
+                return taxonomy
+            else:
+                return taxonomy
+
+        children = wsdlObjectWoRMS.getAphiaChildrenByID(ID, 1, False)
+        
+        for child in children:
+            taxonomy = get_children(taxonomy, child['valid_AphiaID'])
+
+        return taxonomy
+            
+
+    # main bit of the get_taxonomy_worms function
+    try:
+        start_taxa = wsdlObjectWoRMS.getAphiaRecords(start_otu)
+        start_id = start_taxa[0]['valid_AphiaID'] # there might be records that aren't valid - they point to the valid one though
+        # call it again via the ID this time to make sure we've got the right one.
+        start_taxa = wsdlObjectWoRMS.getAphiaRecordByID(start_id)
+        start_taxonomy_level = start_taxa['rank'].lower()
+    except HTTPError:
+        print "Error"
+        sys.exit(-1)
+
+    taxonomy = get_children(taxonomy,start_id)
+
+    return taxonomy, start_taxonomy_level
+            
+
+def main():
+
+    # do stuff
+    parser = argparse.ArgumentParser(
+         prog="Fill tree in using taxonomy",
+         description="Fills in the taxonomic gaps using polytomies within a tree to increase coverage",
+         )
+    parser.add_argument(
+            '-v', 
+            '--verbose', 
+            action='store_true', 
+            help="Verbose output: mainly progress reports.",
+            default=False
+            )
+    parser.add_argument(
+            '--pref_db',
+            help="Taxonomy database to use. Default is Species 2000/ITIS",
+            choices=['itis', 'worms', 'ncbi'],
+            default = 'worms'
+            )
+    parser.add_argument(
+            '--save_taxonomy',
+            help="Save the taxonomy downloaded. Give a filename"
+            )
+    parser.add_argument(
+            '--taxonomy_from_file',
+            help='Use a downloaded taxonomy database from the chosen database, rather than online. Much quicker for large datasets. Give the filename',
+            )
+    parser.add_argument(
+            '--tree_taxonomy',
+            help="Supply a STK taxonomy file for taxa in the tree. If not, one will be created from the database being used here."
+            )
+    parser.add_argument(
+            'top_level', 
+            nargs=1,
+            help="The top level group to look in, e.g. Arthropoda, Decapoda. Must match the database."
+            )
+    parser.add_argument(
+            'input_file', 
+            metavar='input_file',
+            nargs=1,
+            help="Your tree file"
+            )
+    parser.add_argument(
+            'output_file', 
+            metavar='output_file',
+            nargs=1,
+            help="Your new tree file"
+            )
+
+    args = parser.parse_args()
+    verbose = args.verbose
+    input_file = args.input_file[0]
+    output_file = args.output_file[0]
+    top_level = args.top_level[0]
+    save_taxonomy_file = args.save_taxonomy
+    tree_taxonomy = args.tree_taxonomy
+    pref_db = args.pref_db
+    if (save_taxonomy_file == None):
+        save_taxonomy = False
+    else:
+        save_taxonomy = True
+
+    # grab taxa in tree
+    tree = stk.import_tree(input_file)
+    taxa_list = stk._getTaxaFromNewick(tree)
+
+    taxonomy = {}
+
+    # we're going to add the taxa in the tree to the taxonomy, to stop them
+    # being fetched in first place. We delete them later
+    for taxon in taxa_list:
+        taxon = taxon.replace('_',' ')
+        taxonomy[taxon] = []
+
+
+    if (pref_db == 'itis'):
+        # get taxonomy info from itis
+        print "Sorry, ITIS is not implemented yet"
+        pass
+    elif (pref_db == 'worms'):
+        # get tree taxonomy from worms
+        if (tree_taxonomy == None):
+            tree_taxonomy = {}
+            for t in taxa_list:
+                from SOAPpy import WSDL    
+                wsdlObjectWoRMS = WSDL.Proxy('http://www.marinespecies.org/aphia.php?p=soap&wsdl=1')
+                tree_taxonomy[t] = get_tree_taxa_taxonomy(t,wsdlObjectWoRMS)
+        else:
+            tree_taxonomy = stk.load_taxonomy(tree_taxonomy)
+        # get taxonomy from worms
+        taxonomy, start_level = get_taxonomy_worms(taxonomy,top_level)
+
+    elif (pref_db == 'ncbi'):
+        # get taxonomy from ncbi
+        print "Sorry, NCBI is not implemented yet"        
+        pass
+    else:
+        print "ERROR: Didn't understand you database choice"
+        sys.exit(-1)
+
+    # clean up taxonomy, deleting the ones already in the tree
+    for taxon in taxa_list:
+        taxon = taxon.replace('_',' ')        
+        del taxonomy[taxon]
+
+    # step up the taxonomy levels from genus, adding taxa to the correct node
+    # as a polytomy
+    for level in taxonomy_levels[1::]: # skip species....
+        new_taxa = []
+        for t in taxonomy:
+            # skip odd ones that should be in there
+            if start_level in taxonomy[t] and taxonomy[t][start_level] == top_level:
+                try:
+                    new_taxa.append(taxonomy[t][level])
+                except KeyError:
+                    continue # don't have this info
+        new_taxa = _uniquify(new_taxa)
+        for nt in new_taxa:
+            taxa_to_add = []
+            taxa_in_clade = []
+            for t in taxonomy:
+                if start_level in taxonomy[t] and taxonomy[t][start_level] == top_level:
+                    try:
+                        if taxonomy[t][level] == nt:
+                            taxa_to_add.append(t.replace(' ','_'))
+                    except KeyError:
+                        continue
+            # add to tree
+            for t in taxa_list:
+                if level in tree_taxonomy[t] and tree_taxonomy[t][level] == nt:
+                    taxa_in_clade.append(t)
+            if len(taxa_in_clade) > 0:
+                tree = add_taxa(tree, taxa_to_add, taxa_in_clade)
+                for t in taxa_to_add: # clean up taxonomy
+                    del taxonomy[t.replace('_',' ')]
+
+
+    trees = {}
+    trees['tree_1'] = tree
+    output = stk._amalgamate_trees(trees,format='nexus')
+    f = open(output_file, "w")
+    f.write(output)
+    f.close()
+
+    if not save_taxonomy_file == None:
+        with open(save_taxonomy_file, 'w') as f:
+            writer = csv.writer(f)
+            headers = []
+            headers.append("OTU")
+            headers.extend(taxonomy_levels)
+            headers.append("Data source")
+            writer.writerow(headers)
+            for t in taxonomy:
+                otu = t
+                try:
+                    species = taxonomy[t]['species']
+                except KeyError:
+                    species = "-"
+                try:
+                    genus = taxonomy[t]['genus']
+                except KeyError:
+                    genus = "-"
+                try:
+                    family = taxonomy[t]['family']
+                except KeyError:
+                    family = "-"
+                try:
+                    superfamily = taxonomy[t]['superfamily']
+                except KeyError:
+                    superfamily = "-"
+                try:
+                    infraorder = taxonomy[t]['infraorder']
+                except KeyError:
+                    infraorder = "-"
+                try:
+                    suborder = taxonomy[t]['suborder']
+                except KeyError:
+                    suborder = "-"
+                try:
+                    order = taxonomy[t]['order']
+                except KeyError:
+                    order = "-"
+                try:
+                    superorder = taxonomy[t]['superorder']
+                except KeyError:
+                    superorder = "-"
+                try:
+                    subclass = taxonomy[t]['subclass']
+                except KeyError:
+                    subclass = "-"
+                try:
+                    tclass = taxonomy[t]['class']
+                except KeyError:
+                    tclass = "-"
+                try:
+                    subphylum = taxonomy[t]['subphylum']
+                except KeyError:
+                    subphylum = "-"
+                try:
+                    phylum = taxonomy[t]['phylum']
+                except KeyError:
+                    phylum = "-"
+                try:
+                    superphylum = taxonomy[t]['superphylum']
+                except KeyError:
+                    superphylum = "-"
+                try:
+                    infrakingdom = taxonomy[t]['infrakingdom']
+                except:
+                    infrakingdom = "-"
+                try:
+                    subkingdom = taxonomy[t]['subkingdom']
+                except:
+                    subkingdom = "-"
+                try:
+                    kingdom = taxonomy[t]['kingdom']
+                except KeyError:
+                    kingdom = "-"
+                try:
+                    provider = taxonomy[t]['provider']
+                except KeyError:
+                    provider = "-"
+
+                if (isinstance(species, list)):
+                    species = " ".join(species)
+                this_classification = [
+                        otu.encode('utf-8'),
+                        species.encode('utf-8'),
+                        genus.encode('utf-8'),
+                        family.encode('utf-8'),
+                        superfamily.encode('utf-8'),
+                        infraorder.encode('utf-8'),
+                        suborder.encode('utf-8'),
+                        order.encode('utf-8'),
+                        superorder.encode('utf-8'),
+                        subclass.encode('utf-8'),
+                        tclass.encode('utf-8'),
+                        subphylum.encode('utf-8'),
+                        phylum.encode('utf-8'),
+                        superphylum.encode('utf-8'),
+                        infrakingdom.encode('utf-8'),
+                        subkingdom.encode('utf-8'),
+                        kingdom.encode('utf-8'),
+                        provider.encode('utf-8')]
+                writer.writerow(this_classification)
+
+
+def _uniquify(l):
+    """
+    Make a list, l, contain only unique data
+    """
+    keys = {}
+    for e in l:
+        keys[e] = 1
+
+    return keys.keys()
+
+def add_taxa(tree, new_taxa, taxa_in_clade):
+
+    # create new tree of the new taxa
+    #tree_string = "(" + ",".join(new_taxa) + ");"
+    #additionalTaxa = stk._parse_tree(tree_string) 
+
+    # find mrca parent
+    treeobj = stk._parse_tree(tree)
+    mrca = stk.get_mrca(tree,taxa_in_clade)
+    mrca_parent = treeobj.node(mrca).parent
+
+    # insert a node into the tree between the MRCA and it's parent (p4.addNodeBetweenNodes)
+    newNode = treeobj.addNodeBetweenNodes(mrca, mrca_parent)
+
+    # add the new tree at the new node using p4.addSubTree(self, selfNode, theSubTree, subTreeTaxNames=None)
+    #treeobj.addSubTree(newNode, additionalTaxa)
+    for t in new_taxa:
+        treeobj.addSibLeaf(newNode,t)
+
+    # return new tree
+    return treeobj.writeNewick(fName=None,toString=True).strip()
+
+if __name__ == "__main__":
+    main()
+
+
+

=== modified file 'stk/supertree_toolkit.py'
--- stk/supertree_toolkit.py	2014-12-10 08:55:43 +0000
+++ stk/supertree_toolkit.py	2015-03-26 09:59:01 +0000
@@ -52,6 +52,7 @@
 IDENTICAL = 0
 SUBSET = 1
 PLATFORM = sys.platform
+taxonomy_levels = ['species','genus','family','superfamily','infraorder','suborder','order','superorder','subclass','class','subphylum','phylum','superphylum','infrakingdom','subkingdom','kingdom']
 
 # supertree_toolkit is the backend for the STK. Loaded by both the GUI and
 # CLI, this contains all the functions to actually *do* something
@@ -1991,6 +1992,32 @@
 
     return output_string
 
+
+
+def load_taxonomy(taxonomy_csv):
+    """Load in a taxonomy CSV file and convert to taxonomy Dict"""
+    
+    import csv
+
+    taxonomy = {}
+
+    with open(taxonomy_csv, 'rU') as csvfile:
+        tax_reader = csv.reader(csvfile, delimiter=',')
+        tax_reader.next()
+        for row in tax_reader:
+            current_taxonomy = {}
+            i = 1
+            for t in taxonomy_levels:
+                if not row[i] == '-':
+                    current_taxonomy[t] = row[i]
+                i = i+ 1
+
+            current_taxonomy['provider'] = row[17] # data source
+            taxonomy[row[0]] = current_taxonomy
+    
+    return taxonomy
+
+
 def data_overlap(XML, overlap_amount=2, filename=None, detailed=False, show=False, verbose=False, ignoreWarnings=False):
     """ Calculate the amount of taxonomic overlap between source trees.
     The output is a True/False by default, but you can specify an 
@@ -2852,6 +2879,58 @@
 
     return XML
 
+def get_mrca(tree,taxa_list):
+    """Return the node number for the MRCA of the list of given taxa
+       This node number must be used in conjection with a p4 tree object, along
+       the lines of:
+       treeobj = _parse_tree(tree_string)
+       treeobj.node(mrca).parent 
+    """
+
+    # find MRCA of all taxa within this clade, already in the tree
+    node_ids = []
+    # get the nodes of the taxa in question
+    node_id_for_taxa = []
+    treeobj = _parse_tree(tree)
+    for t in taxa_list:
+        node_id_for_taxa.append(treeobj.node(t).nodeNum)
+    # for each, get all parents to root
+    for n in node_id_for_taxa:
+        nodes = []
+        nodes.append(treeobj.node(n).parent.nodeNum)
+        while 1:
+            nn = treeobj.node(nodes[-1]).parent
+            if nn == None:
+                break
+            else:
+                nodes.append(nn.nodeNum)
+        node_ids.append(nodes)
+    # in the shortest list, loop through the values, check they exist in all lists. If it does, 
+    # that node is your MRCA
+    big = sys.maxsize
+    node_ids
+    shortest = 0
+    for n in node_ids:
+        if len(n) < big:
+            big = len(n)
+            shortest = n
+    mrca = -1
+    for s in shortest:
+        found = True
+        for n in node_ids:
+            if not s in n:
+                found = False
+                break # move to next s
+        # if we get here, we have the MRCA
+        if (found):
+            mrca = s
+            break
+    if mrca == -1:
+        # something went wrong!
+        raise InvalidSTKData("Error finding MRCA of: "+" ".join(taxa_list))
+
+    return mrca
+
 ################ PRIVATE FUNCTIONS ########################
 
 def _uniquify(l):

=== modified file 'stk/test/_supertree_toolkit.py'
--- stk/test/_supertree_toolkit.py	2014-12-10 08:55:43 +0000
+++ stk/test/_supertree_toolkit.py	2015-03-26 09:59:01 +0000
@@ -12,7 +12,7 @@
 from stk.supertree_toolkit import data_overlap, read_matrix, subs_file_from_str, clean_data, obtain_trees, get_all_source_names
 from stk.supertree_toolkit import add_historical_event, _sort_data, _parse_xml, _check_sources, _swap_tree_in_XML, replace_genera
 from stk.supertree_toolkit import get_all_taxa, _get_all_siblings, _parse_tree, get_characters_used, _trees_equal, get_weights
-from stk.supertree_toolkit import get_outgroup, set_all_tree_names, create_tree_name
+from stk.supertree_toolkit import get_outgroup, set_all_tree_names, create_tree_name, load_taxonomy
 from lxml import etree
 from util import *
 from stk.stk_exceptions import *
@@ -558,6 +558,18 @@
             self.assert_(c in expected_characters)
         self.assert_(len(characters) == len(expected_characters))
 
+    def test_load_taxonomy(self):
+        csv_file = "data/input/create_taxonomy.csv"
+        expected = {'Archaeopteryx lithographica': {'subkingdom': 'Metazoa', 'subclass': 'Tetrapodomorpha', 'suborder': 'Coelurosauria', 'provider': 'Paleobiology Database', 'genus': 'Archaeopteryx', 'class': 'Aves'},
+                    'Egretta tricolor': {'kingdom': 'Animalia', 'family': 'Ardeidae', 'subkingdom': 'Bilateria', 'subclass': 'Neoloricata', 'class': 'Aves', 'phylum': 'Chordata', 'superphylum': 'Lophozoa', 'suborder': 'Ischnochitonina', 'provider': 'Species 2000 & ITIS Catalogue of Life: April 2013', 'infrakingdom': 'Protostomia', 'genus': 'Egretta', 'order': 'Pelecaniformes', 'species': 'Egretta tricolor'}, 
+                    'Gallus gallus': {'kingdom': 'Animalia', 'infrakingdom': 'Protostomia', 'family': 'Phasianidae', 'subkingdom': 'Bilateria', 'class': 'Aves', 'phylum': 'Chordata', 'superphylum': 'Lophozoa', 'provider': 'Species 2000 & ITIS Catalogue of Life: April 2013', 'genus': 'Gallus', 'order': 'Galliformes', 'species': 'Gallus gallus'}, 
+                    'Thalassarche melanophris': {'kingdom': 'Animalia', 'family': 'Diomedeidae', 'subkingdom': 'Bilateria', 'class': 'Aves', 'phylum': 'Chordata', 'provider': 'Species 2000 & ITIS Catalogue of Life: April 2013', 'infrakingdom': 'Deuterostomia', 'subphylum': 'Vertebrata', 'genus': 'Thalassarche', 'order': 'Procellariiformes', 'species': 'Thalassarche melanophris'},
+                    'Jeletzkytes criptonodosus': {'kingdom': 'Metazoa', 'family': 'Scaphitidae', 'order': 'Ammonoidea', 'phylum': 'Mollusca', 'provider': 'PBDB', 'species': 'Jeletzkytes criptonodosus', 'class': 'Cephalopoda'}}
+        taxonomy = load_taxonomy(csv_file)
+        self.maxDiff = None
+
+        self.assertDictEqual(taxonomy, expected)
+
     def test_name_tree(self):
         XML = etree.tostring(etree.parse('data/input/single_source_no_names.phyml',parser),pretty_print=True)
         xml_root = _parse_xml(XML)

=== modified file 'stk/test/_trees.py'
--- stk/test/_trees.py	2014-06-09 15:36:19 +0000
+++ stk/test/_trees.py	2015-03-26 09:59:01 +0000
@@ -6,7 +6,7 @@
 from stk.supertree_toolkit import import_tree, obtain_trees, get_all_taxa, _assemble_tree_matrix, create_matrix, _delete_taxon, _sub_taxon,_tree_contains
 from stk.supertree_toolkit import _swap_tree_in_XML, substitute_taxa, get_taxa_from_tree, get_characters_from_tree, amalgamate_trees, _uniquify
 from stk.supertree_toolkit import import_trees, import_tree, _trees_equal, _find_trees_for_permuting, permute_tree, get_all_source_names, _getTaxaFromNewick
-
+from stk.supertree_toolkit import get_mrca
 import os
 from lxml import etree
 from util import *
@@ -33,8 +33,8 @@
         self.assert_(e_tree == tree)
 
     def test_import_tutorial_tree(self):
-        test_file = "../../doc/tutorial/Cebezas_etal_tree1.tre"
-        e_tree = "(Onconida_alaini, ((Paramunida_granulata, ((Paramunida_pronoe, ((Paramunida_thalie, (Paramunida_pictura, Paramunida_labis)), (Paramunida_luminata, (Paramunida_belone, (Paramunida_salai, Paramunida_lophia))))), (Paramunida_stichas, Paramunida_proxima))), (Plesionida_concava, Plesionida_aliena)));"
+        test_file = "../../doc/tutorial/5.3_DataEntry/HallThatje_2009.tre"
+        e_tree = "((Aegla_sp., (Pagurus_bernhardus, Pagurus_hirsutiusculus)), (((Cryptolithodes_sitchensis, Cryptolithodes_typicus), (Phyllolithodes_papillosus, (Lopholithodes_mandtii, (Glyptolithodes_cristatipes, (Paralomis_formosa, Paralomis_spinosissima))), (Neolithodes_brodiei, (Paralithodes_camtschaticus, Paralithodes_brevipes), (Lithodes_confundens, Lithodes_ferox)))), (Oedignathus_inermis, (Hapalogaster_dentata, Hapalogaster_mertensii))));"
         tree = import_tree(test_file)
         self.assert_(e_tree == tree)
 
@@ -209,6 +209,12 @@
 
 class TestTreeManipulation(unittest.TestCase): 
 
+    
+    def test_get_mrca(self):
+        tree = "(B,(C,(D,(E,((A,F),((I,(G,H)),(J,(K,L))))))));"
+        mrca = get_mrca(tree,["A","I", "L"])
+        self.assert_(mrca == 8)
+
     def test_get_all_trees(self):
         XML = etree.tostring(etree.parse(single_source_input,parser),pretty_print=True)
         tree = obtain_trees(XML)

=== added file 'stk/test/data/input/create_taxonomy.csv'
--- stk/test/data/input/create_taxonomy.csv	1970-01-01 00:00:00 +0000
+++ stk/test/data/input/create_taxonomy.csv	2015-03-26 09:59:01 +0000
@@ -0,0 +1,6 @@
+OTU,species,genus,family,superfamily,infraorder,suborder,order,superorder,subclass,class,subphylum,phylum,superphylum,infrakingdom,subkingdom,kingdom,Data source
+Archaeopteryx lithographica,-,Archaeopteryx,-,-,-,Coelurosauria,-,-,Tetrapodomorpha,Aves,-,-,-,-,Metazoa,-,Paleobiology Database
+Thalassarche melanophris,Thalassarche melanophris,Thalassarche,Diomedeidae,-,-,-,Procellariiformes,-,-,Aves,Vertebrata,Chordata,-,Deuterostomia,Bilateria,Animalia,Species 2000 & ITIS Catalogue of Life: April 2013
+Egretta tricolor,Egretta tricolor,Egretta,Ardeidae,-,-,Ischnochitonina,Pelecaniformes,-,Neoloricata,Aves,-,Chordata,Lophozoa,Protostomia,Bilateria,Animalia,Species 2000 & ITIS Catalogue of Life: April 2013
+Gallus gallus,Gallus gallus,Gallus,Phasianidae,-,-,-,Galliformes,-,-,Aves,-,Chordata,Lophozoa,Protostomia,Bilateria,Animalia,Species 2000 & ITIS Catalogue of Life: April 2013
+Jeletzkytes criptonodosus,Jeletzkytes criptonodosus,-,Scaphitidae,-,-,-,Ammonoidea,-,-,Cephalopoda,-,Mollusca,-,-,-,Metazoa,PBDB