Merge pull request #106 from McTavishLab/dev

Dev

bin/rf_distance.py

@@ -13,6 +13,7 @@
 
 
 parser = argparse.ArgumentParser()
+parser.add_argument("-d","--results_dir", help="results directory from run")
 parser.add_argument("-t1","--original_tree", help="Original tree")
 parser.add_argument("-t2","--updated_tree", help="Updated Tree")
 parser.add_argument("-otu", "--otu_info", help="File with taxon information JSON")
@@ -24,22 +25,44 @@
 args = parser.parse_args()
 
 tns = dendropy.TaxonNamespace()
-tree1 = dendropy.Tree.get_from_path(args.original_tree,
+
+
+if args.results_dir:
+    workdir = args.results_dir
+    files = [f for f in os.listdir(workdir)]
+    for file in files:
+        if file.startswith('inputs_'):
+            tag = file.split('.')[0].replace('inputs_', '')
+    rundir = "{}/run_{}".format(workdir, tag)
+    outputsdir = "{}/outputs_{}".format(workdir, tag)
+    inputsdir = "{}/inputs_{}".format(workdir, tag)
+    tree1_path = "{}/physcraper_{}.tre".format(inputsdir, tag)
+    otu_json_path = "{}/otu_info_{}.json".format(rundir, tag)
+    tree2_path = "{}/physcraper_{}.tre".format(outputsdir, tag)
+else:
+    tree1_path = args.original_tree
+    tree2_path =- args.updated_tree
+    otu_json_path = args.otu_info
+
+
+tree1 = dendropy.Tree.get_from_path(tree1_path,
                                     schema,
                                     taxon_namespace=tns,
                                     preserve_underscores=True)
-tree2 = dendropy.Tree.get_from_path(args.updated_tree,
+tree2 = dendropy.Tree.get_from_path(tree2_path,
                                     schema,
                                     taxon_namespace=tns,
                                     preserve_underscores=True)
 
+otu_dict = json.load(open(otu_json_path, "r"))
+
+
 leaves_t1 = set([leaf.taxon for leaf in tree1.leaf_nodes()])
 leaves_t2 = set([leaf.taxon for leaf in tree2.leaf_nodes()])
 
 new_tips = len(leaves_t1) - len(leaves_t2)
 sys.stdout.write("{} new tips were added\n".format(new_tips))
 
-otu_dict = json.load(open(args.otu_info, "r"))
 
 old_spp = set()
 new_spp = set()
@@ -70,14 +93,14 @@
 weightedrf = dendropy.calculate.treecompare.weighted_robinson_foulds_distance(tree1, tree2)
 
 tree2.write(path = "pruned_updated.tre", schema="newick")
-print("The RobinsonFoulds distance between these trees is {} and the weighted RF is {}".format(RF, weightedrf))
+sys.stdout.write("The RobinsonFoulds distance between these trees is {} and the weighted RF is {}".format(RF, weightedrf))
 
 
 ##Write out t2 for conflict with opentree
 workdir = os.getcwd()
 
 
-def write_conflict_tree(inputtree, otu_dict):
+def conflict_tree(inputtree, otu_dict):
         tmp_tree = copy.deepcopy(inputtree)
         new_names = set()
         i = 1
@@ -90,15 +113,51 @@ def write_conflict_tree(inputtree, otu_dict):
                 node.taxon.label = new_label
             else:
                 node.label = "_nd{}_".format(i)
-        return tmp_tree.as_string(schema="newick")
-
-treestr_updated = write_conflict_tree(unpruned_tree2, otu_dict)
-treestr_orig = write_conflict_tree(tree1, otu_dict)
-
-
-resp_updated = OT.conflict_str(treestr_updated, 'ott')
-resp_orig = OT.conflict_str(treestr_orig, 'ott')
-
-print(resp_orig.response_dict)
-
-
+        return tmp_tree
+
+tree_updated = conflict_tree(unpruned_tree2, otu_dict)
+tree_orig = conflict_tree(tree1, otu_dict)
+
+
+resp_updated = OT.conflict_str(tree_updated.as_string(schema="newick"), 'ott')
+resp_orig = OT.conflict_str(tree_orig.as_string(schema="newick"), 'ott')
+conflict = resp_updated.response_dict
+
+for node in tree_updated:
+    if node.taxon:
+        node_id = node.taxon.label.split('_')[1]
+        conf_node = conflict.get(node_id, False)
+        if conf_node:
+            new_label = "{} {}".format(conf_node['status'], conf_node['witness_name'])
+            node.taxon.label = new_label
+    else:
+        node_id = node.label.split('_')[1]
+        conf_node = conflict.get(node_id, False)
+        if conf_node:
+            new_label = "{} {}".format(conf_node['status'], conf_node['witness_name'])
+            node.label = new_label
+
+tree_updated.write(path = "conflict_label.tre", schema="newick")
+
+
+'''
+tree_updated_synth = conflict_tree(unpruned_tree2, otu_dict)
+resp_updated_synth = OT.conflict_str(tree_updated_synth.as_string(schema="newick"), 'synth')
+conflict_synth = resp_updated_synth.response_dict
+
+
+for node in tree_updated_synth:
+    if node.taxon:
+        node_id = node.taxon.label.split('_')[1]
+        conf_node = conflict_synth.get(node_id, False)
+        if conf_node:
+            new_label = "{} {}".format(conf_node['status'], conf_node['witness'])
+            node.taxon.label = new_label
+    else:
+        node_id = node.label.split('_')[1]
+        conf_node = conflict_synth.get(node_id, False)
+        if conf_node:
+            new_label = "{} {}".format(conf_node['status'], conf_node['witness'])
+            node.label = new_label
+
+'''

physcraper/aligntreetax.py

@@ -71,7 +71,7 @@ def generate_ATT_from_run(workdir, start_files='output', tag=None, configfile=No
        configfile = "{}/run.config".format(rundir)
     try:
         alnfi = "{}/physcraper_{}.fas".format(outputsdir, tag)
-        treefile = "{}/physcraper_{}.tre".format(inputsdir, tag)
+        treefile = "{}/physcraper_{}.tre".format(outputsdir, tag)
         otu_json = "{}/otu_info_{}.json".format(rundir, tag)
         assert(os.path.exists(alnfi))
         assert(os.path.exists(treefile))

physcraper/scrape.py

@@ -386,9 +386,16 @@ def get_full_seq(self, gb_id, blast_seq):
         seq_path = "{}/{}.fasta".format(self.ids.full_seq_path, gb_id)
         if not os.path.exists(seq_path):
             db_path = "{}/nt".format(self.config.blastdb)
-            cmd1 = "blastdbcmd -db {}  -entry {} -outfmt %f -out {}".format(db_path, gb_id, seq_path)
-            # debug(cmd1)
-            os.system(cmd1)
+            try:
+                subprocess.check_call(["blastdbcmd",
+                                        "-db",  db_path,
+                                        "-entry", gb_id,
+                                        "-outfmt", "%f",
+                                        "-out", seq_path])
+
+            except subprocess.CalledProcessError as grepexc:
+                sys.stderr.write("error code {}, {}".format(grepexc.returncode, grepexc.output))
+                sys.exit()
             # read in file to get full seq
         f = open(seq_path)
         seq = ""