Repository: smdabdoub/kraken-biom

Issue:
Taxonomic hierarchy incorrect for certain cases
We've found that in some edge cases, the taxonomic hierarchy kraken-biom assigns to a given ID is incorrect.  It looks like the assumption that a given rank in the Kraken report always falls under the most recent higher rank (for example, for a given "S" entry for species, the closest "P" for phylum in the previous lines in the file) is not always true.

This is with kraken-biom 1.0.1 under Python 3.5 (in Anaconda) on Linux.

Here's a chunk from an example Kraken report, where I'm specifically looking at 5693 (Trypanosoma cruzi) near the bottom.

```
...
  0.00  1       0       P       3041            Chlorophyta
  0.00  1       1       C       75966             Trebouxiophyceae
  0.00  2       0       -       556282        Jakobida
  0.00  2       0       G       221723          Seculamonas
  0.00  2       2       S       221724            Seculamonas ecuadoriensis
  0.00  1       0       -       33682         Euglenozoa
  0.00  1       0       O       5653            Kinetoplastida
  0.00  1       0       F       5654              Trypanosomatidae
  0.00  1       0       G       5690                Trypanosoma
  0.00  1       0       -       47570                 Schizotrypanum
  0.00  1       0       S       5693                    Trypanosoma cruzi
  0.00  1       1       -       353153                    Trypanosoma cruzi strain CL Brener
...
```

Looking from ID 5693 on up, in terms of indentation in the last column: Kraken shows taxa up through "O", then an un-ranked taxon (Euglenozoa), and then nothing for a *very* long time until Eukaryota many lines above (not shown).  The phylum Chlorophyta and class Trebouxiophyceae do not actually contain Trypanosoma cruzi; they're just the closest previous phylum and class shown above that species in the file.  But kraken-biom's output gives this Consensus Lineage for ID 5693:

```
k__Eukaryota; p__Chlorophyta; c__Trebouxiophyceae; o__Kinetoplastida; f__Trypanosomatidae; g__Trypanosoma; s__cruzi
```

The [NCBI Taxonomy Browser](https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?mode=Info&id=5693&lvl=3&p=has_linkout&p=blast_url&p=genome_blast&p=mapview&lin=f&keep=1&srchmode=1&unlock) seems to match what I saw in Kraken, with Kingdom=Eukaryota; Unranked=Euglenozoa; Order=Kinetoplastida; Family=Trypanosomatidae; etc. (No explicit phylum or class listed.)

I can't say for sure because [the Kraken documentation](https://ccb.jhu.edu/software/kraken/MANUAL.html) doesn't go into detail, but it looks to me like it's the *indentation for the scientific name* that corresponds to the hierarchy and to what rank sits above a given entry, and not necessarily the *rank of the previous taxa listed*.  So in my case even though the previous "P" in the report file is Chlorophyta, that group doesn't actually include ID 5693, so we shouldn't have a phylum or class assigned.

Discussion:

Patch:
diff --git a/.travis.yml b/.travis.yml
index 8c3a27c..dd55a5f 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -2,9 +2,9 @@ cache: apt
 sudo: false
 language: python
 python:
+  - "3.3"
   - "3.4"
   - "3.5"
-  - "3.6"
 addons:
   apt:
     packages:
diff --git a/kraken_biom.py b/kraken_biom.py
index 728206b..0f7b369 100755
--- a/kraken_biom.py
+++ b/kraken_biom.py
@@ -6,13 +6,13 @@ Kraken output (http://ccb.jhu.edu/software/kraken/).
 """
 from __future__ import absolute_import, division, print_function
 
-from pathlib import Path
 import argparse
 from collections import OrderedDict
 import csv
 from datetime import datetime as dt
 from gzip import open as gzip_open
 import os.path as osp
+import re
 import sys
 from textwrap import dedent as twdd
 
@@ -33,7 +33,7 @@ __license__ = "MIT"
 __url__ = "http://github.com/smdabdoub/kraken-biom"
 __maintainer__ = "Shareef M. Dabdoub"
 __email__ = "dabdoub.2@osu.edu"
-__version__ = '1.1.0'
+__version__ = '1.0.1'
 
 
 field_names = ["pct_reads", "clade_reads", "taxon_reads", 
@@ -81,6 +81,31 @@ def tax_fmt(tax_lvl, end):
     # print(tax)
     return tax
 
+def parse_tax_lvl(entry, tax_lvl_depth=[]):
+    """
+    Parse a single kraken-report entry and return a dictionary of taxa for its
+    named ranks.
+
+    :type entry: dict
+    :param entry: attributes of a single kraken-report row.
+    :type tax_lvl_depth: list
+    :param tax_lvl_depth: running record of taxon levels encountered in
+    previous calls.
+    """
+    # How deep in the hierarchy are we currently?  Each two spaces of
+    # indentation is one level deeper.  Also parse the scientific name at this
+    # level.
+    depth_and_name = re.match('^( *)(.*)', entry['sci_name'])
+    depth = len(depth_and_name.group(1))//2
+    name = depth_and_name.group(2)
+    # Remove the previous levels so we're one higher than the level of the new
+    # taxon.  (This also works if we're just starting out or are going deeper.)
+    del tax_lvl_depth[depth:]
+    # Append the new taxon.
+    tax_lvl_depth.append((entry['rank'], name))
+    # Create a tax_lvl dict for the named ranks.
+    tax_lvl = {x[0]: x[1] for x in tax_lvl_depth if x[0] in ranks}
+    return(tax_lvl)
 
 def parse_kraken_report(kdata, max_rank, min_rank):
     """
@@ -95,8 +120,6 @@ def parse_kraken_report(kdata, max_rank, min_rank):
     taxa = OrderedDict()
     # the master collection of read counts (keyed on NCBI ID)
     counts = OrderedDict()
-    # running record of the current taxonomic hierarchy
-    tax_lvl = {}
     # current rank
     r = 0
     max_rank_idx = ranks.index(max_rank)
@@ -105,18 +128,13 @@ def parse_kraken_report(kdata, max_rank, min_rank):
     for entry in kdata:
         erank = entry['rank'].strip()
         # print("erank: "+erank)
-        
-        # move back up the taxa tree to the current level
-        if erank in ranks and ranks.index(erank) < r:
-            r = ranks.index(erank)
-            tax_lvl = {r: tax_lvl[r] if r in tax_lvl else '' for r in ranks[:r]}
-        
-        # add current rank to running tally of ranks
+
         if erank in ranks:
-            tax_lvl[erank] = entry["sci_name"].strip()
-            # print("Recording tax level ({}): {}".format(erank, tax_lvl[erank]))
             r = ranks.index(erank)
         
+        # update running tally of ranks
+        tax_lvl = parse_tax_lvl(entry)
+
         # record the reads assigned to this taxon level, and record the taxonomy string with the NCBI ID
         if erank in ranks and min_rank_idx >= ranks.index(entry['rank']) >= max_rank_idx:
             taxon_reads = int(entry["taxon_reads"])
@@ -313,10 +331,8 @@ def handle_program_options():
 
     parser = argparse.ArgumentParser(description=twdd(descr),
                            formatter_class=argparse.RawDescriptionHelpFormatter)
-    parser.add_argument('kraken_reports', nargs='*',
+    parser.add_argument('kraken_reports', nargs='+',
                         help="Results files from the kraken-report tool.")
-    parser.add_argument('-k', '--kraken_reports_fp', metavar="REPORTS_FP",
-                        help="Folder containing kraken reports")
     parser.add_argument('--max', default="O", choices=ranks[:-1],
                         help="Assigned reads will be recorded only if \
                               they are at or below max rank. Default: O.")
@@ -366,16 +382,12 @@ def main():
         Defaulting to BIOM 1.0 (JSON)."""
         print(twdd(msg))
 
-    if ranks.index(args.max) > ranks.index(args.min):
+    if ranks.index(args.max) >= ranks.index(args.min):
         msg = "ERROR: Max and Min ranks are out of order: {} < {}"
         sys.exit(msg.format(args.max, args.min))
 
-    reports = args.kraken_reports
-    if args.kraken_reports_fp:
-        reports += [str(p) for p in Path(args.kraken_reports_fp).glob('*')]
-
     # load all kraken-report files and parse them
-    sample_counts, taxa = process_samples(reports, 
+    sample_counts, taxa = process_samples(args.kraken_reports, 
                                           max_rank=args.max, 
                                           min_rank=args.min)
 


Test patch:
diff --git a/test/test_parsing.py b/test/test_parsing.py
index cc092dd..ac91bd5 100644
--- a/test/test_parsing.py
+++ b/test/test_parsing.py
@@ -129,6 +129,16 @@ class kraken_biom_Test(unittest.TestCase):
             0.00	1	0	-	206037	          Human endogenous retroviruses
             0.00	1	0	S	45617	            Human endogenous retrovirus K
             0.00	1	1	-	166122	              Human endogenous retrovirus K113
+            0.00	3	1	D	2759	    Eukaryota
+            0.00	1	0	P	3041	        Chlorophyta
+            0.00	1	1	C	75966	          Trebouxiophyceae
+            0.00	1	0	-	33682	      Euglenozoa
+            0.00	1	0	O	5653	        Kinetoplastida
+            0.00	1	0	F	5654	          Trypanosomatidae
+            0.00	1	0	G	5690	            Trypanosoma
+            0.00	1	0	-	47570	              Schizotrypanum
+            0.00	1	0	S	5693	                Trypanosoma cruzi
+            0.00	1	1	-	353153	                  Trypanosoma cruzi strain CL Brener
             """)))
 
     def run_parse_kraken_report(self, manual, max_rank, min_rank):
@@ -169,7 +179,8 @@ class kraken_biom_Test(unittest.TestCase):
                   '37734': 9,
                   '45617': 1,
                   '543': 6,
-                  '562': 1}
+                  '562': 1,
+                  '5693': 1}
 
         self.run_parse_kraken_report(manual,  max_rank="O", min_rank="S")
 
@@ -198,7 +209,8 @@ class kraken_biom_Test(unittest.TestCase):
                   '37734': 9,
                   '45617': 1,
                   '543': 6,
-                  '562': 1}
+                  '562': 1,
+                  '5693': 1}
 
         self.run_parse_kraken_report(manual,  max_rank="F", min_rank="S")
 
@@ -226,7 +238,8 @@ class kraken_biom_Test(unittest.TestCase):
                   '374840': 2,
                   '37734': 9,
                   '45617': 1,
-                  '562': 1}
+                  '562': 1,
+                  '5693': 1}
 
         self.run_parse_kraken_report(manual,  max_rank="G", min_rank="S")
 
@@ -252,7 +265,8 @@ class kraken_biom_Test(unittest.TestCase):
                   '374840': 2,
                   '37734': 9,
                   '45617': 1,
-                  '562': 1}
+                  '562': 1,
+                  '5693': 1}
 
         self.run_parse_kraken_report(manual,  max_rank="S", min_rank="S")
 
@@ -279,7 +293,9 @@ class kraken_biom_Test(unittest.TestCase):
                   '29465': 2,
                   '48736': 11,
                   '543': 6,
-                  '561': 1}
+                  '561': 1,
+                  '75966': 1,
+                  '5690': 1}
 
         self.run_parse_kraken_report(manual,  max_rank="C", min_rank="G")
 
@@ -307,7 +323,9 @@ class kraken_biom_Test(unittest.TestCase):
                             "s__Hyposoter fugitivus ichnovirus"],
                   "374840":["k__Viruses", "p__", "c__", "o__", 
                             "f__Microviridae", "g__Microvirus", 
-                            "s__Enterobacteria phage phiX174 sensu lato"]
+                            "s__Enterobacteria phage phiX174 sensu lato"],
+                  "5693":["k__Eukaryota", "p__", "c__", "o__Kinetoplastida",
+                          "f__Trypanosomatidae", "g__Trypanosoma", "s__cruzi"]
                  }
 
         _, taxa = kb.parse_kraken_report(self.sample_kraken_rep,
@@ -324,4 +342,4 @@ class kraken_biom_Test(unittest.TestCase):
 
 
 if __name__ == '__main__':
-    unittest.main()
\ No newline at end of file
+    unittest.main()