removed a trace

ASCII (issue #89)

python/obitools3/apps/optiongroups/__init__.py

@@ -1,3 +1,9 @@
+import codecs
+
+def unescaped_str(arg_str):
+    return arg_str.encode('latin-1', 'backslashreplace').decode('unicode-escape')
+
+
 def __addInputOption(optionManager):
     
     optionManager.add_argument(
@@ -43,7 +49,13 @@ def __addImportInputOption(optionManager):
                      action="store_const", dest="obi:inputformat",
                      default=None,
                      const=b'silva',
-                     help="Input file is in SILVA fasta format")
+                     help="Input file is in SILVA fasta format. If NCBI taxonomy provided with --taxonomy, taxid and scientific name will be added for each sequence.")
+
+    group.add_argument('--rdp-input',
+                     action="store_const", dest="obi:inputformat",
+                     default=None,
+                     const=b'rdp',
+                     help="Input file is in RDP training set fasta format. If NCBI taxonomy provided with --taxonomy, taxid and scientific name will be added for each sequence.")
 
     group.add_argument('--embl-input',
                      action="store_const", dest="obi:inputformat",
@@ -125,15 +137,15 @@ def __addImportInputOption(optionManager):
 def __addTabularOption(optionManager):
     group = optionManager.add_argument_group("Input and output format options for tabular files")
 
-    group.add_argument('--header',
-                     action="store_true", dest="obi:header",
-                     default=False,
-                     help="First line of tabular file contains column names")
+    group.add_argument('--no-header',
+                     action="store_false", dest="obi:header",
+                     default=True,
+                     help="Don't print the header (first line with column names")
 
     group.add_argument('--sep',
                      action="store", dest="obi:sep",
                      default="\t",
-                     type=str,
+                     type=unescaped_str,
                      help="Column separator")
 
 
@@ -165,6 +177,16 @@ def __addTabularInputOption(optionManager):
                      help="Lines starting by this char are considered as comment")
 
 
+def __addTabularOutputOption(optionManager):
+    group = optionManager.add_argument_group("Output format options for tabular files")
+    
+    __addTabularOption(optionManager)
+
+    group.add_argument('--na-int-stay-na',
+                     action="store_false", dest="obi:na_int_to_0",
+                     help="NA (Non available) integer values should be exported as NA in tabular output (default: they are converted to 0 for tabular output).") # TODO
+
+
 def __addTaxdumpInputOption(optionManager):  # TODO maybe not the best way to do it
     group = optionManager.add_argument_group("Input format options for taxdump")
 
@@ -198,6 +220,10 @@ def addTabularInputOption(optionManager):
     __addTabularInputOption(optionManager)
 
 
+def addTabularOutputOption(optionManager):
+    __addTabularOutputOption(optionManager)
+
+
 def addTaxonomyOption(optionManager):
     __addTaxonomyOption(optionManager)
 
@@ -210,6 +236,7 @@ def addAllInputOption(optionManager):
     __addInputOption(optionManager)
     __addImportInputOption(optionManager)
     __addTabularInputOption(optionManager)
+    __addTabularOutputOption(optionManager)
     __addTaxonomyOption(optionManager)
     __addTaxdumpInputOption(optionManager)    
 
@@ -270,6 +297,35 @@ def __addExportOutputOption(optionManager):
                      const=b'tabular',
                      help="Output file is in tabular format")
 
+    group.add_argument('--metabaR-output',
+                     action="store_const", dest="obi:outputformat",
+                     default=None,
+                     const=b'metabaR',
+                     help="Export the files needed by the obifiles_to_metabarlist function of the metabaR package")
+
+    group.add_argument('--metabaR-prefix',
+                     action="store", dest="obi:metabarprefix",
+                     type=str,
+                     help="Prefix for the files when using --metabaR-output option")
+
+    group.add_argument('--metabaR-ngsfilter',
+                     action="store", dest="obi:metabarngsfilter",
+                     type=str,
+                     default=None,
+                     help="URI to the ngsfilter view when using --metabaR-output option (if not provided, it is not exported)")
+
+    group.add_argument('--metabaR-samples',
+                     action="store", dest="obi:metabarsamples",
+                     type=str,
+                     default=None,
+                     help="URI to the sample metadata view when using --metabaR-output option (if not provided, it is built as just a list of the sample names)")
+
+    group.add_argument('--only-keys',
+                       action="append", dest="obi:only_keys",
+                       type=str,
+                       default=[],
+                       help="Only export the given keys (columns).")
+
     group.add_argument('--print-na',
                      action="store_true", dest="obi:printna",
                      default=False,
@@ -302,14 +358,14 @@ def addTabularOutputOption(optionManager):
 
 def addExportOutputOption(optionManager):
     __addExportOutputOption(optionManager)
-    __addTabularOption(optionManager)
+    __addTabularOutputOption(optionManager)
 
 
 def addAllOutputOption(optionManager):
     __addOutputOption(optionManager)
     __addDMSOutputOption(optionManager)
     __addExportOutputOption(optionManager)
-    __addTabularOption(optionManager)
+    __addTabularOutputOption(optionManager)
 
 
 def addNoProgressBarOption(optionManager):

python/obitools3/commands/addtaxids.pyx

@@ -26,7 +26,7 @@ import sys
 from cpython.exc cimport PyErr_CheckSignals
 
  
-__title__="Annotate sequences with their corresponding NCBI taxid found from the taxon scientific name."
+__title__="Annotate sequences with their corresponding NCBI taxid found from the taxon scientific name"
 
 
  
@@ -45,7 +45,7 @@ def addOptions(parser):
                        metavar="<TAXID_TAG>",
                        default=b"TAXID",
                        help="Name of the tag to store the found taxid "
-                            "(default: 'TAXID'.")
+                            "(default: 'TAXID').")
 
     group.add_argument('-n', '--taxon-name-tag', 
                        action="store", 
@@ -53,7 +53,7 @@ def addOptions(parser):
                        metavar="<SCIENTIFIC_NAME_TAG>",
                        default=b"SCIENTIFIC_NAME",
                        help="Name of the tag giving the scientific name of the taxon "
-                            "(default: 'SCIENTIFIC_NAME'.")
+                            "(default: 'SCIENTIFIC_NAME').")
 
     group.add_argument('-g', '--try-genus-match',
                      action="store_true", dest="addtaxids:try_genus_match",
@@ -174,6 +174,7 @@ def run(config):
                 taxid_column[i] = taxon.taxid
                 found_count+=1
             elif try_genus: # try finding genus or other parent taxon from the first word
+                #print(i, o_view[i].id)
                 taxon_name_sp = taxon_name.split(b" ")
                 taxon = taxo.get_taxon_by_name(taxon_name_sp[0], res_anc)
                 if taxon is not None:

python/obitools3/commands/align.pyx

@@ -19,7 +19,7 @@ import time
 import sys
 
 
-__title__="Aligns one sequence column with itself or two sequence columns"
+__title__="Align one sequence column with itself or two sequence columns"
 
 
 def addOptions(parser):
@@ -158,7 +158,7 @@ def run(config):
     i_view_name = i_uri.split(b"/")[0]
     i_column_name = b""
     i_element_name = b""
-    if len(i_uri.split(b"/")) == 2:
+    if len(i_uri.split(b"/")) >= 2:
         i_column_name = i_uri.split(b"/")[1]
     if len(i_uri.split(b"/")) == 3:
         i_element_name = i_uri.split(b"/")[2]
@@ -181,7 +181,7 @@ def run(config):
         i_dms_name_2 = i_dms_2.name
         i_uri_2 = input_2[1]
         original_i_view_name_2 = i_uri_2.split(b"/")[0]
-        if len(i_uri_2.split(b"/")) == 2:
+        if len(i_uri_2.split(b"/")) >= 2:
             i_column_name_2 = i_uri_2.split(b"/")[1]
         if len(i_uri_2.split(b"/")) == 3:
             i_element_name_2 = i_uri_2.split(b"/")[2]

python/obitools3/commands/alignpairedend.pyx

@@ -23,7 +23,7 @@ import os
 
 from cpython.exc cimport PyErr_CheckSignals
 
-__title__="Aligns paired-ended reads"
+__title__="Align paired-ended reads"
 
 
 

python/obitools3/commands/build_ref_db.pyx

@@ -16,7 +16,7 @@ import sys
 from cpython.exc cimport PyErr_CheckSignals
 
 
-__title__="Tag a set of sequences for PCR and sequencing errors identification"
+__title__="Build a reference database for ecotag"
 
 
 def addOptions(parser):
@@ -31,10 +31,9 @@ def addOptions(parser):
     group.add_argument('--threshold','-t',
                       action="store", dest="build_ref_db:threshold",
                       metavar='<THRESHOLD>',
-                      default=0.0,
+                      default=0.99,
                       type=float,
-                      help="Score threshold as a normalized identity, e.g. 0.95 for an identity of 95%%. Default: 0.00"
-                           " (no threshold).")
+                      help="Score threshold as a normalized identity, e.g. 0.95 for an identity of 95%%. Default: 0.99.")
 
 
 def run(config):

python/obitools3/commands/cat.pyx

@@ -22,7 +22,7 @@ import sys
 from cpython.exc cimport PyErr_CheckSignals
 
 
-__title__="Concatenate views."
+__title__="Concatenate views"
 
  
 def addOptions(parser):
@@ -97,7 +97,7 @@ def run(config):
             Column.new_column(o_view, REVERSE_QUALITY_COLUMN, OBI_QUAL, associated_column_name=REVERSE_SEQUENCE_COLUMN, associated_column_version=o_view[REVERSE_SEQUENCE_COLUMN].version)
         
     # Initialize multiple elements columns
-    if type(output_0)==BufferedWriter: 
+    if type(output_0)!=BufferedWriter: 
         dict_cols = {}
         for v_uri in config["cat"]["views_to_cat"]:
             v = open_uri(v_uri)[1]
@@ -134,7 +134,11 @@ def run(config):
                 rep = repr(entry)
                 output_0.write(str2bytes(rep)+b"\n")
             else:
+                try:
                     o_view[i] = entry
+                except:
+                    print("\nError with entry:", repr(entry))
+                    print(repr(o_view))
             i+=1
         v.close()
 

python/obitools3/commands/clean.pyx

@@ -54,11 +54,11 @@ def addOptions(parser):
                        default=False,
                        help="Only sequences labeled as heads are kept in the output. Default: False")
 
-    group.add_argument('--cluster-tags', '-C',
-                       action="store_true", 
-                       dest="clean:cluster-tags",
-                       default=False,
-                       help="Adds tags for each sequence giving its cluster's head and weight for each sample.")
+#    group.add_argument('--cluster-tags', '-C',
+#                       action="store_true", 
+#                       dest="clean:cluster-tags",
+#                       default=False,
+#                       help="Adds tags for each sequence giving its cluster's head and weight for each sample.")
 
     group.add_argument('--thread-count','-p',   # TODO should probably be in a specific option group
                        action="store", dest="clean:thread-count",
@@ -143,3 +143,4 @@ def run(config):
     i_dms.close(force=True)
 
     logger("info", "Done.")
+    

python/obitools3/commands/count.pyx

@@ -10,7 +10,7 @@ from obitools3.dms.capi.obiview cimport COUNT_COLUMN
 from cpython.exc cimport PyErr_CheckSignals
 
 
-__title__="Counts sequence records"
+__title__="Count sequence records"
 
  
 def addOptions(parser):
@@ -29,6 +29,12 @@ def addOptions(parser):
                         default=False,
                         help="Prints only the total count of sequence records (if a sequence has no `count` attribute, its default count is 1) (default: False).")
 
+    group.add_argument('-c','--count-tag',
+                        action="store", dest="count:countcol",
+                        default='COUNT',
+                        type=str,
+                        help="Name of the tag/column associated with the count information (default: COUNT).")
+
     
 def run(config):
     
@@ -42,17 +48,19 @@ def run(config):
         raise Exception("Could not read input")
     entries = input[1]
     
+    countcol = config['count']['countcol'] 
+    
     count1 = len(entries)
     count2 = 0
     
-    if COUNT_COLUMN in entries and ((config['count']['sequence'] == config['count']['all']) or (config['count']['all'])) :
+    if countcol in entries and ((config['count']['sequence'] == config['count']['all']) or (config['count']['all'])) :
         for e in entries:
             PyErr_CheckSignals()
-            count2+=e[COUNT_COLUMN]
+            count2+=e[countcol]
 
-    if COUNT_COLUMN in entries and (config['count']['sequence'] == config['count']['all']):
+    if countcol in entries and (config['count']['sequence'] == config['count']['all']):
         print(count1,count2)
-    elif COUNT_COLUMN in entries and config['count']['all']:
+    elif countcol in entries and config['count']['all']:
         print(count2)
     else:
         print(count1)

python/obitools3/commands/ecopcr.pyx

@@ -175,6 +175,14 @@ def run(config):
     o_dms_name = output[0].name
     o_view_name = output[1]
 
+    # Open the taxonomy DMS
+    taxdms = open_uri(config['obi']['taxoURI'],
+                     dms_only=True)
+    if taxdms is None:
+        raise Exception("Could not open taxonomy DMS")
+    tax_dms = taxdms[0]
+    tax_dms_name = taxdms[0].name
+    
     # Read taxonomy name    
     taxonomy_name = config['obi']['taxoURI'].split("/")[-1]   # Robust in theory
 
@@ -197,7 +205,8 @@ def run(config):
 
     # TODO: primers in comments?
 
-    if obi_ecopcr(i_dms.name_with_full_path, tobytes(i_view_name), tobytes(taxonomy_name), \
+    if obi_ecopcr(i_dms.name_with_full_path, tobytes(i_view_name), 
+                  tax_dms.name_with_full_path, tobytes(taxonomy_name), \
                   o_dms.name_with_full_path, tobytes(o_view_name), comments, \
                   tobytes(config['ecopcr']['primer1']), tobytes(config['ecopcr']['primer2']), \
                   config['ecopcr']['error'], \

python/obitools3/commands/export.pyx

@@ -6,6 +6,9 @@ from obitools3.apps.config import logger
 from obitools3.dms import DMS
 from obitools3.dms.obiseq import Nuc_Seq
 from obitools3.dms.capi.obiview cimport QUALITY_COLUMN
+from obitools3.writers.tab import TabWriter
+from obitools3.format.tab import TabFormat
+from obitools3.utils cimport tobytes, tostr
 
 from obitools3.apps.optiongroups import addMinimalInputOption, \
                                         addExportOutputOption, \
@@ -76,6 +79,13 @@ def run(config):
     else:
         pb = ProgressBar(withoutskip - skip, config)
 
+    if config['obi']['outputformat'] == b'metabaR':
+        # Check prefix
+        if "metabarprefix" not in config["obi"]:
+            raise Exception("Prefix needed when exporting for metabaR (--metabaR-prefix option)")
+        else:
+            metabaRprefix = config["obi"]["metabarprefix"]
+        
     i=0
     for seq in iview :
         PyErr_CheckSignals()
@@ -91,6 +101,81 @@ def run(config):
         pb(i, force=True)
         print("", file=sys.stderr)
 
+    if config['obi']['outputformat'] == b'metabaR':
+        
+        # Export ngsfilter file if view provided
+        if 'metabarngsfilter' in config['obi']:
+            ngsfilter_input = open_uri(config['obi']['metabarngsfilter'])
+            if ngsfilter_input is None:
+                raise Exception("Could not read ngsfilter view for metabaR output")
+            ngsfilter_view = ngsfilter_input[1]
+    
+            ngsfilter_output = open(config['obi']['metabarprefix']+'.ngsfilter', 'w')
+            
+            for line in ngsfilter_view:
+
+                line_to_print = b""                
+                line_to_print += line[b'experiment']
+                line_to_print += b"\t"
+                line_to_print += line[b'sample']
+                line_to_print += b"\t"
+                line_to_print += line[b'forward_tag']
+                line_to_print += b":"
+                line_to_print += line[b'reverse_tag']
+                line_to_print += b"\t"
+                line_to_print += line[b'forward_primer']
+                line_to_print += b"\t"
+                line_to_print += line[b'reverse_primer']
+                line_to_print += b"\t"
+                line_to_print += line[b'additional_info']
+                
+                print(tostr(line_to_print), file=ngsfilter_output)
+            
+            if ngsfilter_input[0] != input[0]:
+                ngsfilter_input[0].close()
+            ngsfilter_output.close()
+            
+        # Export sample metadata
+        samples_output = open(config['obi']['metabarprefix']+'_samples.csv', 'w')
+
+        # Export sample metadata file if view provided
+        if 'metabarsamples' in config['obi']:
+            samples_input = open_uri(config['obi']['metabarsamples'])
+            if samples_input is None:
+                raise Exception("Could not read sample view for metabaR output")
+            samples_view = samples_input[1]
+        
+            # Export with tab formatter
+            TabWriter(TabFormat(header=True, sep='\t',), 
+                      samples_output,
+                      header=True)
+            
+            if samples_input[0] != input[0]:
+                samples_input[0].close()
+        
+        # Else export just sample names from main view
+        else:
+            
+            sample_list = []
+            if 'MERGED_sample' in iview:
+                sample_list = iview['MERGED_sample'].keys()
+            elif 'sample' not in iview:
+                for seq in iview:
+                    sample = seq['sample']
+                    if sample not in sample_list:
+                        sample_list.append(sample)
+            else:
+                logger("warning", "Can not read sample list from main view for metabaR sample list export")            
+            
+            print("sample_id", file=samples_output)
+            for sample in sample_list:
+                line_to_print = b""                
+                line_to_print += sample
+                line_to_print += b"\t"
+                print(tostr(line_to_print), file=samples_output)
+        
+        samples_output.close()
+
     # TODO save command in input dms?
     
     if not BrokenPipeError and not IOError:

python/obitools3/commands/grep.pyx

@@ -91,7 +91,7 @@ def addOptions(parser):
                        metavar="<ATTRIBUTE_NAME>",
                        help="Select records with the attribute <ATTRIBUTE_NAME> "
                             "defined (not set to NA value). "
-                            "Several -a options can be used on the same "
+                            "Several -A options can be used on the same "
                             "command line.")
     
     group.add_argument("-L", "--lmax",
@@ -258,6 +258,13 @@ def Filter_generator(options, tax_filter, i_view):
 
 
 def Taxonomy_filter_generator(taxo, options):
+    
+    if (("required_ranks" in options and options["required_ranks"]) or \
+       ("required_taxids" in options and options["required_taxids"]) or \
+       ("ignored_taxids" in options and options["ignored_taxids"])) and \
+       (taxo is None):
+        raise RollbackException("obi grep error: can't use taxonomy options without providing a taxonomy. Rollbacking view")
+    
     if taxo is not None:
         def tax_filter(seq):
             good = True

python/obitools3/commands/head.pyx

@@ -16,7 +16,7 @@ from io import BufferedWriter
 from cpython.exc cimport PyErr_CheckSignals
 
 
-__title__="Keep the N first lines of a view."
+__title__="Keep the N first lines of a view"
 
  
 def addOptions(parser):

python/obitools3/commands/import.pyx

@@ -34,14 +34,17 @@ from obitools3.dms.capi.obidms cimport obi_import_view
 from obitools3.dms.capi.obitypes cimport obitype_t, \
                                          OBI_VOID, \
                                          OBI_QUAL, \
-                                         OBI_STR
+                                         OBI_STR, \
+                                         OBI_INT
 
 from obitools3.dms.capi.obierrno cimport obi_errno
 
 from obitools3.apps.optiongroups import addImportInputOption, \
                                         addTabularInputOption, \
                                         addTaxdumpInputOption, \
-                                        addMinimalOutputOption
+                                        addMinimalOutputOption, \
+                                        addNoProgressBarOption, \
+                                        addTaxonomyOption
 
 from obitools3.uri.decode import open_uri
 
@@ -50,9 +53,10 @@ from obitools3.apps.config import logger
 from cpython.exc cimport PyErr_CheckSignals
 
 from io import BufferedWriter
+import ast
 
 
-__title__="Imports sequences from different formats into a DMS"
+__title__="Import sequences from different formats into a DMS"
  
  
 default_config = {   'destview'     : None,
@@ -69,7 +73,9 @@ def addOptions(parser):
     addImportInputOption(parser)
     addTabularInputOption(parser)
     addTaxdumpInputOption(parser)
+    addTaxonomyOption(parser)
     addMinimalOutputOption(parser)
+    addNoProgressBarOption(parser)
 
     group = parser.add_argument_group('obi import specific options')
 
@@ -85,6 +91,10 @@ def addOptions(parser):
                      help="If importing a view into another DMS, do it by importing each line, saving disk space if the original view "
                           "has a line selection associated.")
 
+#    group.add_argument('--only-id',
+#                     action="store", dest="import:onlyid",
+#                     help="only id")
+
 def run(config):
     
     cdef   tuple       input
@@ -97,6 +107,7 @@ def run(config):
     cdef   bint        get_quality
     cdef   bint        NUC_SEQS_view
     cdef   bint        silva
+    cdef   bint        rdp
     cdef   int         nb_elts
     cdef   object      d
     cdef   View        view
@@ -180,6 +191,16 @@ def run(config):
         logger("info", "Done.")
         return
 
+    # Open taxonomy if there is one
+    if 'taxoURI' in config['obi'] and config['obi']['taxoURI'] is not None:
+        taxo_uri = open_uri(config['obi']['taxoURI'])
+        if taxo_uri is None or taxo_uri[2] == bytes:
+            raise Exception("Couldn't open taxonomy")
+        taxo = taxo_uri[1]
+
+    else :
+        taxo = None
+
     # If importing a view between two DMS and not wanting to save space if line selection in original view, use C API
     if isinstance(input[1], View) and not config['import']['space_priority']:
         if obi_import_view(input[0].name_with_full_path, o_dms.name_with_full_path, input[1].name, tobytes((config['obi']['outputURI'].split('/'))[-1])) < 0 :
@@ -192,8 +213,11 @@ def run(config):
         logger("info", "Done.")
         return
 
-    if entry_count >= 0:
+    # Reinitialize the progress bar
+    if entry_count >= 0 and config['obi']['noprogressbar'] == False:
         pb = ProgressBar(entry_count, config)
+    else:
+        pb = None
         
     NUC_SEQS_view = False
     if isinstance(output[1], View) :
@@ -209,12 +233,17 @@ def run(config):
         def_col = view[DEFINITION_COLUMN]
         seq_col = view[NUC_SEQUENCE_COLUMN]
 
-    # Prepare taxon scientific name if SILVA file
-    if 'inputformat' in config['obi'] and config['obi']['inputformat'] == b"silva":
-        silva = True
-        sci_name_col = Column.new_column(view, SCIENTIFIC_NAME_COLUMN, OBI_STR)
-    else:
+    # Prepare taxon scientific name and taxid refs if RDP or SILVA file
     silva = False
+    rdp = False
+    if 'inputformat' in config['obi'] and (config['obi']['inputformat'] == b"silva" or config['obi']['inputformat'] == b"rdp"):
+        #if taxo is None:
+        #    raise Exception("A taxonomy (as built by 'obi import --taxdump') must be provided for SILVA and RDP files")
+        silva = True
+        rdp = True
+        if taxo is not None:
+            sci_name_col = Column.new_column(view, SCIENTIFIC_NAME_COLUMN, OBI_STR)
+            taxid_col = Column.new_column(view, TAXID_COLUMN, OBI_INT)
 
     dcols = {}
 
@@ -265,8 +294,13 @@ def run(config):
         # Reinitialize the input
         if isinstance(input[0], CompressedFile):
             input_is_file = True
-        if entry_count >= 0:
+        
+        # Reinitialize the progress bar
+        if entry_count >= 0 and config['obi']['noprogressbar'] == False:
             pb = ProgressBar(entry_count, config)
+        else:
+            pb = None
+
         try:
             input[0].close()
         except AttributeError:
@@ -275,6 +309,11 @@ def run(config):
         if input is None:
             raise Exception("Could not open input URI")
     
+#    if 'onlyid' in config['import']:
+#        onlyid = tobytes(config['import']['onlyid'])
+#    else:
+#        onlyid = None
+    
     entries = input[1]
     i = 0
     for entry in entries :
@@ -292,6 +331,9 @@ def run(config):
         elif not i%50000:
             logger("info", "Imported %d entries", i)
         
+#        if onlyid is not None and entry.id != onlyid:
+#            continue
+        
         try:
              
             if NUC_SEQS_view: 
@@ -307,10 +349,17 @@ def run(config):
                 if get_quality:
                     qual_col[i] = entry.quality
                     
-                # Parse taxon scientific name if SILVA file
-                if silva:
-                    sci_name = entry.definition.split(b";")[-1]
-                    sci_name_col[i] = sci_name
+                # Parse taxon scientific name if RDP file
+                if (rdp or silva) and taxo is not None:
+                    sci_names = entry.definition.split(b";")
+                    for sci_name in reversed(sci_names):
+                        if sci_name.split()[0] != b'unidentified' and sci_name.split()[0] != b'uncultured' and sci_name.split()[0] != b'metagenome' :
+                            taxon = taxo.get_taxon_by_name(sci_name)
+                            if taxon is not None:
+                                sci_name_col[i] = taxon.name
+                                taxid_col[i] = taxon.taxid
+                                #print(taxid_col[i], sci_name_col[i])
+                                break
 
             for tag in entry :
                 
@@ -324,6 +373,14 @@ def run(config):
                     if tag[:7] == b"merged_":
                         tag = MERGED_PREFIX+tag[7:]
 
+                    if type(value) == bytes and value[:1]==b"[" :
+                        try:
+                            if type(eval(value)) == list:
+                                value = eval(value)
+                                #print(value)
+                        except:
+                            pass
+                    
                     if tag not in dcols :
                                                 
                         value_type = type(value)
@@ -331,19 +388,23 @@ def run(config):
                         value_obitype = OBI_VOID
                         dict_col = False
                         
-                        if value_type == dict or value_type == list :
+                        if value_type == dict :
                             nb_elts = len(value)
                             elt_names = list(value)
-                            if value_type == dict :
                             dict_col = True
                         else :
                             nb_elts = 1
                             elt_names = None
                         
+                        if value_type == list :
+                            tuples = True
+                        else:
+                            tuples = False
+                        
                         value_obitype = get_obitype(value)
                                                 
                         if value_obitype != OBI_VOID :
-                            dcols[tag] = (Column.new_column(view, tag, value_obitype, nb_elements_per_line=nb_elts, elements_names=elt_names, dict_column=dict_col), value_obitype)
+                            dcols[tag] = (Column.new_column(view, tag, value_obitype, nb_elements_per_line=nb_elts, elements_names=elt_names, dict_column=dict_col, tuples=tuples), value_obitype)
                                                      
                             # Fill value
                             dcols[tag][0][i] = value
@@ -371,7 +432,7 @@ def run(config):
                             # Fill value
                             dcols[tag][0][i] = value
                          
-                        except IndexError :
+                        except (IndexError, OverflowError):
                                                                                  
                             value_type = type(value)
                             old_column = dcols[tag][0]
@@ -419,18 +480,19 @@ def run(config):
                                 dcols[tag][0][i] = value
         
         except Exception as e:
-            print("\nCould not import sequence id:", entry.id, "(error raised:", e, ")")
+            print("\nCould not import sequence:", repr(entry), "(error raised:", e, ")")
             if 'skiperror' in config['obi'] and not config['obi']['skiperror']:
                 raise e
             else:
                 pass
+                i-=1  # overwrite problematic entry
         
         i+=1
  
     if pb is not None:
         pb(i, force=True)
         print("", file=sys.stderr)
-    logger("info", "Imported %d entries", i)
+    logger("info", "Imported %d entries", len(view))
     
     # Save command config in View and DMS comments
     command_line = " ".join(sys.argv[1:])

python/obitools3/commands/ngsfilter.pyx

@@ -24,7 +24,7 @@ from cpython.exc cimport PyErr_CheckSignals
 from io import BufferedWriter
 
 
-__title__="Assigns sequence records to the corresponding experiment/sample based on DNA tags and primers"
+__title__="Assign sequence records to the corresponding experiment/sample based on DNA tags and primers"
 
 
 def addOptions(parser):
@@ -271,7 +271,7 @@ cdef tuple annotate(sequences, infos, no_tags, verbose=False):
         sequences[0][REVERSE_QUALITY_COLUMN] = sequences[1].quality     # used by alignpairedend tool
 
     for seq in sequences:
-        if hasattr(seq, "quality_array"): 
+        if hasattr(seq, "quality_array") and seq.quality_array is not None: 
             q = -reduce(lambda x,y:x+y,(math.log10(z) for z in seq.quality_array),0)/len(seq.quality_array)*10
             seq[b'avg_quality']=q
             q = -reduce(lambda x,y:x+y,(math.log10(z) for z in seq.quality_array[0:10]),0)

python/obitools3/commands/rm.pyx

@@ -0,0 +1,44 @@
+#cython: language_level=3
+
+from obitools3.uri.decode import open_uri
+from obitools3.apps.config import logger
+from obitools3.dms import DMS
+from obitools3.apps.optiongroups import addMinimalInputOption
+from obitools3.dms.view.view cimport View
+import os
+ 
+__title__="Delete a view"
+
+
+def addOptions(parser):    
+    addMinimalInputOption(parser)
+
+def run(config):
+
+    DMS.obi_atexit()
+    
+    logger("info", "obi rm")
+
+    # Open the input
+    input = open_uri(config['obi']['inputURI'])
+    if input is None:
+        raise Exception("Could not read input")
+    
+    # Check that it's a view
+    if isinstance(input[1], View) :
+        view = input[1]
+    else: 
+        raise NotImplementedError()
+    
+    # Get the path to the view file to remove
+    path = input[0].full_path  # dms path
+    path+=b"/VIEWS/"
+    path+=view.name
+    path+=b".obiview"
+    
+    # Close the view and the DMS
+    view.close()
+    input[0].close(force=True)
+    
+    # Rm
+    os.remove(path)

python/obitools3/commands/sort.pyx

@@ -36,7 +36,7 @@ NULL_VALUE = {OBI_BOOL: OBIBool_NA,
               OBI_STR: b""}
  
 
-__title__="Sort view lines according to the value of a given attribute."
+__title__="Sort view lines according to the value of a given attribute"
 
  
 def addOptions(parser):

python/obitools3/commands/stats.pyx

@@ -16,7 +16,7 @@ import sys
 from cpython.exc cimport PyErr_CheckSignals
 
 
-__title__="Compute basic statistics for attribute values."
+__title__="Compute basic statistics for attribute values"
 
 '''
 `obi stats` computes basic statistics for attribute values of sequence records.
@@ -121,7 +121,10 @@ def variance(v):
     if len(v)==1: 
         return 0
     s = reduce(lambda x,y:(x[0]+y,x[1]+y**2),v,(0.,0.))
-    return s[1]/(len(v)-1) - s[0]**2/len(v)/(len(v)-1)
+    var = round(s[1]/(len(v)-1) - s[0]**2/len(v)/(len(v)-1), 5) # round to go around shady python rounding stuff when var is actually 0
+    if var == -0.0: # then fix -0 to +0 if was rounded to -0
+        var = 0.0
+    return var
 
 
 def varpop(values, options):
@@ -154,7 +157,7 @@ def run(config):
     else :
         taxo = None
 
-    statistics = set(config['stats']['minimum']) | set(config['stats']['maximum']) | set(config['stats']['mean'])
+    statistics = set(config['stats']['minimum']) | set(config['stats']['maximum']) | set(config['stats']['mean']) | set(config['stats']['var']) | set(config['stats']['sd'])
     total = 0
     catcount={}
     totcount={}
@@ -195,7 +198,7 @@ def run(config):
         except KeyError:
             totcount[category]=totcount.get(category,0)+1
         for var in statistics:
-            if var in line:
+            if var in line and line[var] is not None:
                 v = line[var]
                 if var not in values:
                     values[var]={}
@@ -238,14 +241,34 @@ def run(config):
     else:
         sdvar= "%s"
         
-    hcat = "\t".join([pcat % x for x in config['stats']['categories']]) + "\t" +\
-           "\t".join([minvar % x for x in config['stats']['minimum']])  + "\t" +\
-           "\t".join([maxvar % x for x in config['stats']['maximum']])  + "\t" +\
-           "\t".join([meanvar % x for x in config['stats']['mean']])  + "\t" +\
-           "\t".join([varvar % x for x in config['stats']['var']])  + "\t" +\
-           "\t".join([sdvar % x for x in config['stats']['sd']]) + \
-           "\t   count" + \
-           "\t   total" 
+    hcat = ""
+    
+    for x in config['stats']['categories']:
+        hcat += pcat % x
+        hcat += "\t"
+    
+    for x in config['stats']['minimum']:
+        hcat += minvar % x
+        hcat += "\t"
+
+    for x in config['stats']['maximum']:
+        hcat += maxvar % x
+        hcat += "\t"
+
+    for x in config['stats']['mean']:
+        hcat += meanvar % x
+        hcat += "\t"
+    
+    for x in config['stats']['var']:
+        hcat += varvar % x
+        hcat += "\t"
+
+    for x in config['stats']['sd']:
+        hcat += sdvar % x
+        hcat += "\t"
+
+    hcat += "count\ttotal" 
+    
     print(hcat)
     sorted_stats = sorted(catcount.items(), key = lambda kv:(totcount[kv[0]]), reverse=True)
     for i in range(len(sorted_stats)):
@@ -265,8 +288,8 @@ def run(config):
             print((("%%%df" % lvarp[m]) % varp[m][c])+"\t", end="")
         for m in config['stats']['sd']:
             print((("%%%df" % lsigma[m]) % sigma[m][c])+"\t", end="")
-        print("%7d" %catcount[c], end="")
-        print("%9d" %totcount[c])
+        print("%d" %catcount[c]+"\t", end="")
+        print("%d" %totcount[c]+"\t")
 
     input[0].close(force=True)
     

python/obitools3/commands/tail.pyx

@@ -15,7 +15,7 @@ from cpython.exc cimport PyErr_CheckSignals
 from io import BufferedWriter
 
 
-__title__="Keep the N last lines of a view."
+__title__="Keep the N last lines of a view"
 
  
 def addOptions(parser):

python/obitools3/commands/test.pyx

@@ -39,7 +39,7 @@ COL_COMMENTS_MAX_LEN = 2048
 MAX_INT = 2147483647    # used to generate random float values
  
  
-__title__="Tests if the obitools are working properly"
+__title__="Test if the obitools are working properly"
 
 
 default_config = {

python/obitools3/dms/capi/obiecopcr.pxd

@@ -8,6 +8,7 @@ cdef extern from "obi_ecopcr.h" nogil:
 
     int obi_ecopcr(const char* input_dms_name,
                    const char* i_view_name,
+                   const char* tax_dms_name,
                    const char* taxonomy_name,
                    const char* output_dms_name,
                    const char* o_view_name,

python/obitools3/dms/capi/obitypes.pxd

@@ -54,6 +54,8 @@ cdef extern from "obitypes.h" nogil:
     extern uint8_t*      OBIQual_int_NA
     extern void*         OBITuple_NA
     
+    extern obiint_t      OBI_INT_MAX
+
     const_char_p name_data_type(int data_type)
 
 ctypedef OBIType_t obitype_t

python/obitools3/dms/column/column.pyx

@@ -7,7 +7,8 @@ __OBIDMS_COLUMN_CLASS__ = {}
 from ..capi.obitypes cimport name_data_type, \
                              obitype_t, \
                              obiversion_t, \
-                             OBI_QUAL
+                             OBI_QUAL, \
+                             OBI_STR
 
 from ..capi.obidms cimport obi_import_column
 
@@ -128,6 +129,10 @@ cdef class Column(OBIWrapper) :
         else:
             elements_names_p = NULL
         
+        if column_name_b == b"SAMPLE" or column_name_b == b"sample":
+        # force str type
+            data_type = OBI_STR
+        
         if data_type == OBI_QUAL:
             if associated_column_name_b == b"":
                 if column_name == QUALITY_COLUMN:

python/obitools3/dms/column/typed_column/str.pyx

@@ -74,6 +74,9 @@ cdef class Column_str(Column_idx):
         if value is None :
             value_b = <char*>OBIStr_NA
         else :
+            if self.name == b'sample' or self.name == b'SAMPLE':
+                if type(value) == int:
+                    value = str(value) # force sample ids to be str
             value_bytes = tobytes(value)
             value_b = <char*>value_bytes
 
@@ -137,6 +140,9 @@ cdef class Column_multi_elts_str(Column_multi_elts_idx):
         if value is None :
             value_b = <char*>OBIStr_NA
         else :
+            if self.name == b'sample' or self.name == b'SAMPLE':
+                if type(value) == int:
+                    value = str(value) # force sample ids to be str
             value_bytes = tobytes(value)
             value_b = <char*>value_bytes
                 
@@ -206,6 +212,9 @@ cdef class Column_tuples_str(Column_idx):
             i = 0
             for elt in value :
                 if elt is not None and elt != '':
+                    if self.name == b'sample' or self.name == b'SAMPLE':
+                        if type(elt) == int:
+                            elt = str(elt) # force sample ids to be str
                     elt_b = tobytes(elt)
                     strcpy(array+i, <char*>elt_b)
                     i = i + len(elt_b) + 1

python/obitools3/dms/taxo/taxo.pyx

@@ -1,5 +1,7 @@
 #cython: language_level=3
 
+import sys
+
 from obitools3.utils cimport str2bytes, bytes2str, tobytes, tostr
 from ..capi.obidms cimport OBIDMS_p, obi_dms_get_full_path
                           
@@ -34,7 +36,7 @@ cdef class Taxonomy(OBIWrapper) :
         return <OBIDMS_taxonomy_p>(self._pointer)        
 
     cdef fill_name_dict(self):
-        print("Indexing taxon names...")
+        print("Indexing taxon names...", file=sys.stderr)
         
         cdef OBIDMS_taxonomy_p pointer = self.pointer()
         cdef ecotx_t*     taxon_p
@@ -91,6 +93,8 @@ cdef class Taxonomy(OBIWrapper) :
             raise RuntimeError("Error : Cannot read taxonomy %s" 
                                % tostr(name))
 
+        print("Taxonomy read", file=sys.stderr)
+
         taxo = OBIWrapper.new_wrapper(Taxonomy, pointer)
 
         dms.register(taxo)
@@ -147,6 +151,8 @@ cdef class Taxonomy(OBIWrapper) :
         for r in range((<OBIDMS_taxonomy_p>pointer).ranks.count) :
             taxo._ranks.append(obi_taxo_rank_index_to_label(r, (<OBIDMS_taxonomy_p>pointer).ranks))
         
+        print('Read %d taxa' % len(taxo), file=sys.stderr)
+        
         return taxo
 
     
@@ -304,6 +310,11 @@ cdef class Taxonomy(OBIWrapper) :
     def name(self):
         return self._name
     
+    # ranks property getter
+    @property
+    def ranks(self):
+        return self._ranks
+
     
     def parental_tree_iterator(self, int taxid):
         """

python/obitools3/dms/view/view.pyx

@@ -77,7 +77,7 @@ cdef class View(OBIWrapper) :
  
  
     @staticmethod
-    def import_view(object dms_1, object dms_2, object view_name_1, object view_name_2):
+    def import_view(object dms_1, object dms_2, object view_name_1, object view_name_2): # TODO argument to import line by line to avoid huge AVL copy
         if obi_import_view(tobytes(dms_1), tobytes(dms_2), tobytes(view_name_1), tobytes(view_name_2)) < 0 :
             raise Exception("Error importing a view")
 
@@ -345,7 +345,7 @@ cdef class View(OBIWrapper) :
                                        nb_elements_per_line=new_nb_elements_per_line, elements_names=new_elements_names, 
                                        dict_column=(new_nb_elements_per_line>1), comments=old_column.comments, alias=column_name_b+tobytes('___new___'))
         
-        switch_to_dict = old_column.nb_elements_per_line == 1 and new_nb_elements_per_line > 1
+        switch_to_dict = not old_column.dict_column and new_nb_elements_per_line > 1
         ori_key = old_column._elements_names[0]
         
         for i in range(length) :
@@ -600,7 +600,8 @@ cdef class View(OBIWrapper) :
                 if element is not None:
                     if element.comments[b"input_dms_name"] is not None :
                         for i in range(len(element.comments[b"input_dms_name"])) :
-                            if element.comments[b"input_dms_name"][i] == element.dms.name and b"/" not in element.comments[b"input_view_name"][i]:  # Same DMS and not a special element like a taxonomy
+                            if b"/" not in element.comments[b"input_view_name"][i] and element.comments[b"input_view_name"][i] in element.dms \
+                            and element.comments[b"input_dms_name"][i] == element.dms.name :  # Same DMS and not a special element like a taxonomy and view was not deleted
                                 top_level.append(element.dms[element.comments[b"input_view_name"][i]])
                             else:
                                 top_level.append(None)
@@ -798,7 +799,8 @@ cdef class Line :
  
  
     def keys(self):
-        return self._view.keys()
+        cdef bytes key
+        return [key for key in self._view.keys()]
  
  
     def __contains__(self, object column_name):

python/obitools3/format/fasta.pyx

@@ -7,11 +7,12 @@ from obitools3.utils cimport bytes2str
 
 cdef class FastaFormat:
 	
-	def __init__(self, list tags=[], bint printNAKeys=False, bytes NAString=b"NA"):
+	def __init__(self, list tags=[], bint printNAKeys=False, bytes NAString=b"NA", bint NAIntTo0=False):
 		self.headerFormatter = HeaderFormat("fasta",
 										    tags=tags,
 										    printNAKeys=printNAKeys,
-										    NAString=NAString)
+										    NAString=NAString,
+										    NAIntTo0=NAIntTo0)
 		
 	@cython.boundscheck(False)	
 	def __call__(self, object data):

python/obitools3/format/fastq.pyx

@@ -8,11 +8,12 @@ from obitools3.utils cimport bytes2str, str2bytes, tobytes
 # TODO quality offset option?
 cdef class FastqFormat:
 	
-	def __init__(self, list tags=[], bint printNAKeys=False, bytes NAString=b"NA"):
+	def __init__(self, list tags=[], bint printNAKeys=False, bytes NAString=b"NA", bint NAIntTo0=False):
 		self.headerFormatter = HeaderFormat("fastq",
 										    tags=tags,
 										    printNAKeys=printNAKeys,
-										    NAString=NAString)
+										    NAString=NAString,
+                                            NAIntTo0=NAIntTo0)
 		
 	@cython.boundscheck(False)	
 	def __call__(self, object data):

python/obitools3/format/header.pxd

@@ -4,5 +4,6 @@ cdef class HeaderFormat:
     cdef set    tags
     cdef bint   printNAKeys
     cdef bytes  NAString
+    cdef bint   NAIntTo0
     cdef size_t headerBufferLength
     

python/obitools3/format/header.pyx

@@ -8,13 +8,14 @@ from obitools3.dms.capi.obiview cimport NUC_SEQUENCE_COLUMN, \
 
 from obitools3.utils cimport str2bytes, bytes2str_object
 from obitools3.dms.column.column cimport Column_line
+from obitools3.dms.column.typed_column.int cimport Column_int, Column_multi_elts_int
 
 
 cdef class HeaderFormat:
     
     SPECIAL_KEYS = [NUC_SEQUENCE_COLUMN, ID_COLUMN, DEFINITION_COLUMN, QUALITY_COLUMN]
     
-    def __init__(self, str format="fasta", list tags=[], bint printNAKeys=False, bytes NAString=b"NA"):
+    def __init__(self, str format="fasta", list tags=[], bint printNAKeys=False, bytes NAString=b"NA", bint NAIntTo0=False):
         '''
             @param format: 
             @type  format: `str`
@@ -32,6 +33,7 @@ cdef class HeaderFormat:
         self.tags           = set(tags)
         self.printNAKeys    = printNAKeys
         self.NAString       = NAString
+        self.NAIntTo0       = NAIntTo0
         
         if format=="fasta":
             self.start=b">"
@@ -57,12 +59,20 @@ cdef class HeaderFormat:
             if k in tags:
                 value = data[k]
                 if value is None or (isinstance(value, Column_line) and value.is_NA()):
-                    if self.printNAKeys:
+                    if isinstance(data.view[k], Column_int) and self.NAIntTo0: # people want missing int values to be 0
+                        value = b'0'
+                    elif self.printNAKeys:
                         value = self.NAString
                     else:
                         value = None
                 else:
                     if type(value) == Column_line:
+                        if isinstance(data.view[k], Column_multi_elts_int) and self.NAIntTo0:
+                            value = dict(value)
+                            for key in data.view[k].keys():
+                                if key not in value or value[key]:
+                                    value[key] = 0
+                        else:
                             value = value.bytes()
                     else:
                         if type(value) == tuple:

python/obitools3/format/tab.pxd

@@ -4,5 +4,8 @@ cdef class TabFormat:
     cdef bint header
     cdef bint first_line
     cdef bytes NAString
-    cdef list tags
+    cdef set   tags
     cdef bytes sep
+    cdef bint NAIntTo0
+    cdef bint metabaR
+    cdef bint ngsfilter

python/obitools3/format/tab.pyx

@@ -4,39 +4,59 @@ cimport cython
 from obitools3.dms.view.view cimport Line
 from obitools3.utils cimport bytes2str_object, str2bytes, tobytes
 from obitools3.dms.column.column cimport Column_line, Column_multi_elts
+from obitools3.dms.column.typed_column.int cimport Column_int, Column_multi_elts_int
 
 import sys
 
 cdef class TabFormat:
     
-    def __init__(self, header=True, bytes NAString=b"NA", bytes sep=b"\t"):
+    def __init__(self, list tags=[], header=True, bytes NAString=b"NA", bytes sep=b"\t", bint NAIntTo0=True, metabaR=False, ngsfilter=False):
+        self.tags = set(tags)
         self.header = header
         self.first_line = True
         self.NAString = NAString
         self.sep = sep
+        self.NAIntTo0 = NAIntTo0
+        self.metabaR = metabaR
+        self.ngsfilter = ngsfilter
         
     @cython.boundscheck(False)    
     def __call__(self, object data):
         
+        cdef set ktags
+        cdef list tags = [key for key in data]
+        
         line = []
         
-        if self.first_line:
-            self.tags = [k for k in data.keys()]
+        if self.tags is not None and self.tags:
+            ktags = self.tags
+        else:
+            ktags = set(tags)            
                 
         if self.header and self.first_line:
-            for k in self.tags:
+            for k in ktags:
+                if k in tags:
+                    if self.metabaR:
+                        if k == b'NUC_SEQ':
+                            ktoprint = b'sequence'
+                        else:
+                            ktoprint = k.lower()
+                        ktoprint = ktoprint.replace(b'merged_', b'')
+                    else:
+                        ktoprint = k
                     if isinstance(data.view[k], Column_multi_elts):
                         keys = data.view[k].keys()
                         keys.sort()
                         for k2 in keys:
-                        line.append(tobytes(k)+b':'+tobytes(k2))
+                            line.append(tobytes(ktoprint)+b':'+tobytes(k2))
                     else:
-                    line.append(tobytes(k))
+                        line.append(tobytes(ktoprint))
             r = self.sep.join(value for value in line)
             r += b'\n'
             line = []
                     
-        for k in self.tags:
+        for k in ktags:
+            if k in tags:
                 value = data[k]
                 if isinstance(data.view[k], Column_multi_elts):
                     keys = data.view[k].keys()
@@ -48,10 +68,13 @@ cdef class TabFormat:
                         for k2 in keys: # TODO could be much more efficient
                             if value[k2] is not None:
                                 line.append(str2bytes(str(bytes2str_object(value[k2]))))  # genius programming
+                            else:
+                                if self.NAIntTo0 and isinstance(data.view[k], Column_multi_elts_int):
+                                    line.append(b"0")
                                 else:
                                     line.append(self.NAString)
                 else:
-                if value is not None:
+                    if value is not None or (self.NAIntTo0 and isinstance(data.view[k], Column_int)):
                         line.append(str2bytes(str(bytes2str_object(value))))
                     else:
                         line.append(self.NAString)

python/obitools3/parsers/genbank.pyx

@@ -22,11 +22,11 @@ from libc.stdlib  cimport free, malloc, realloc
 from libc.string cimport strcpy, strlen
 
 
-_featureMatcher = re.compile(b'^FEATURES.+\n(?=ORIGIN )',re.DOTALL + re.M)
+_featureMatcher = re.compile(b'^FEATURES.+\n(?=ORIGIN(\s*))',re.DOTALL + re.M)
 
 _headerMatcher = re.compile(b'^LOCUS.+(?=\nFEATURES)', re.DOTALL + re.M)
-_seqMatcher    = re.compile(b'ORIGIN .+(?=//\n)', re.DOTALL + re.M)
-_cleanSeq1     = re.compile(b'ORIGIN.+\n')
+_seqMatcher    = re.compile(b'^ORIGIN.+(?=//\n)', re.DOTALL + re.M)
+_cleanSeq1     = re.compile(b'ORIGIN(\s*)\n')
 _cleanSeq2     = re.compile(b'[ \n0-9]+')
 _acMatcher     = re.compile(b'(?<=^ACCESSION   ).+',re.M)
 _deMatcher     = re.compile(b'(?<=^DEFINITION  ).+\n( .+\n)*',re.M)
@@ -155,9 +155,9 @@ def genbankIterator_file(lineiterator,
         yield seq
         read+=1
     
-    # Last sequence
+    # Last sequence if not empty lines
+    if entry.strip():
         seq = genbankParser(entry)
-    
         yield seq
     
     free(entry)

python/obitools3/parsers/ngsfilter.pyx

@@ -48,13 +48,13 @@ def ngsfilterIterator(lineiterator,
         all_lines.insert(0, firstline)
         
     # Insert header for column names
-    column_names = [b"experiment", b"sample", b"forward_tag", b"reverse_tag", b"forward_primer", b"reverse_primer"]
+    column_names = [b"experiment", b"sample", b"forward_tag", b"reverse_tag", b"forward_primer", b"reverse_primer",b"additional_info"]
     header = out_sep.join(column_names)
     
     new_lines.append(header)
         
     for line in all_lines:
-        split_line = line.split()
+        split_line = line.split(maxsplit=5)
         tags = split_line.pop(2)
         tags = tags.split(b":")
         for t_idx in range(len(tags)):
@@ -64,7 +64,7 @@ def ngsfilterIterator(lineiterator,
             tags.append(tags[0])
         split_line.insert(2, tags[0])
         split_line.insert(3, tags[1])
-        new_lines.append(out_sep.join(split_line[0:6]))
+        new_lines.append(out_sep.join(split_line[0:7]))
     
     return tabIterator(iter(new_lines),
                        header = True,

python/obitools3/uri/decode.pyx

@@ -174,6 +174,9 @@ def open_uri(uri,
              dms_only=False,
              force_file=False):
 
+    if type(uri) == str and not uri.isascii():
+        raise Exception("Paths must be ASCII characters only")
+            
     cdef bytes urib = tobytes(uri)
     cdef bytes scheme
     cdef tuple dms
@@ -210,11 +213,12 @@ def open_uri(uri,
     
     error = None
         
-    if urib != b"-" and \
+    if b'/taxonomy/' in urib or \
+        (urib != b"-" and \
         (scheme==b"dms" or \
          (scheme==b"" and \
           (((not input) and "outputformat" not in config["obi"]) or \
-          (input and "inputformat" not in config["obi"])))):   # TODO maybe not best way
+          (input and "inputformat" not in config["obi"]))))):   # TODO maybe not best way
                 
         if default_dms is not None and b"/" not in urip.path:   # assuming view to open in default DMS (TODO discuss)
             dms=(default_dms, urip.path)
@@ -276,7 +280,12 @@ def open_uri(uri,
             iseq = urib
             objclass = bytes
     else:  # TODO update uopen to be able to write? 
-        if not urip.path or urip.path == b'-':
+        if config['obi']['outputformat'] == b'metabaR':
+            if 'metabarprefix' not in config['obi']:
+                raise Exception("Prefix needed when exporting for metabaR (--metabaR-prefix option)")
+            else:
+                file = open(config['obi']['metabarprefix']+'.tab', 'wb')
+        elif not urip.path or urip.path == b'-':
             file = sys.stdout.buffer
         else:
             file = open(urip.path, 'wb')
@@ -302,7 +311,7 @@ def open_uri(uri,
         if b'seqtype' in qualifiers:
             seqtype=qualifiers[b'seqtype'][0]
         else:
-            if format == b"ngsfilter" or format == b"tabular":  # TODO discuss
+            if format == b"ngsfilter" or format == b"tabular" or format == b"metabaR":  # TODO discuss
                 seqtype=None
             else:
                 try:
@@ -386,10 +395,10 @@ def open_uri(uri,
             raise MalformedURIException('Malformed header argument in URI')
        
         if b"sep" in qualifiers:
-            sep=tobytes(qualifiers[b"sep"][0][0])
+            sep = tobytes(qualifiers[b"sep"][0][0])
         else:
             try:
-                sep=tobytes(config["obi"]["sep"])
+                sep = tobytes(config["obi"]["sep"])
             except KeyError:
                 sep=None
         
@@ -427,6 +436,20 @@ def open_uri(uri,
             except KeyError:
                 nastring=b'NA'
 
+        if b"na_int_to_0" in qualifiers:
+            try:
+                na_int_to_0=eval(qualifiers[b"na_int_to_0"][0])
+            except Exception as e:
+                raise MalformedURIException("Malformed 'NA_int_to_0' argument in URI")
+        else:
+            try:
+                na_int_to_0=config["obi"]["na_int_to_0"]
+            except KeyError:
+                if format==b"tabular" or format==b"metabaR":
+                    na_int_to_0=True
+                else:
+                    na_int_to_0=False
+                       
         if b"stripwhite" in qualifiers:
             try:
                 stripwhite=eval(qualifiers[b"stripwhite"][0])
@@ -461,17 +484,36 @@ def open_uri(uri,
             except KeyError:
                 commentchar=b'#'
 
+        if b"only_keys" in qualifiers:
+            only_keys=qualifiers[b"only_keys"][0] # not sure that works but no one ever uses qualifiers
+        else:
+            try:
+                only_keys_str=config["obi"]["only_keys"]
+                only_keys=[]
+                for key in only_keys_str:
+                    only_keys.append(tobytes(key))
+            except KeyError:
+                only_keys=[]
+
+        if b"metabaR_prefix" in qualifiers:
+            metabaR_prefix = tobytes(qualifiers[b"metabaR_prefix"][0][0])
+        else:
+            try:
+                metabaR_prefix = tobytes(config["obi"]["metabarprefix"])
+            except KeyError:
+                metabaR_prefix=None
+
         if format is not None:
             if seqtype==b"nuc":
                 objclass = Nuc_Seq    # Nuc_Seq_Stored? TODO
-                if format==b"fasta" or format==b"silva":
+                if format==b"fasta" or format==b"silva" or format==b"rdp":
                     if input:
                         iseq = fastaNucIterator(file, 
                                                 skip=skip, 
                                                 only=only,
                                                 nastring=nastring)
                     else:
-                        iseq = FastaNucWriter(FastaFormat(printNAKeys=printna, NAString=nastring), 
+                        iseq = FastaNucWriter(FastaFormat(tags=only_keys, printNAKeys=printna, NAString=nastring), 
                                               file,
                                               skip=skip,
                                               only=only)
@@ -484,7 +526,7 @@ def open_uri(uri,
                                              noquality=noquality,
                                              nastring=nastring)
                     else:
-                        iseq = FastqWriter(FastqFormat(printNAKeys=printna, NAString=nastring), 
+                        iseq = FastqWriter(FastqFormat(tags=only_keys, printNAKeys=printna, NAString=nastring), 
                                            file,
                                            skip=skip,
                                            only=only)
@@ -520,7 +562,17 @@ def open_uri(uri,
                                        skip = skip,
                                        only = only)
                 else:
-                    iseq = TabWriter(TabFormat(header=header, NAString=nastring, sep=sep), 
+                    iseq = TabWriter(TabFormat(tags=only_keys, header=header, NAString=nastring, sep=sep, NAIntTo0=na_int_to_0), 
+                                               file,
+                                               skip=skip,
+                                               only=only,
+                                               header=header)
+            elif format==b"metabaR":
+                objclass = dict
+                if input:
+                    raise NotImplementedError('Input data file format not implemented')
+                else:
+                    iseq = TabWriter(TabFormat(tags=only_keys, header=header, NAString=nastring, sep=sep, NAIntTo0=na_int_to_0, metabaR=True), 
                                                file,
                                                skip=skip,
                                                only=only,
@@ -538,7 +590,7 @@ def open_uri(uri,
                                              skip = skip,
                                              only = only)
                 else:
-                    raise NotImplementedError('Output sequence file format not implemented')
+                    raise NotImplementedError('Output data file format not implemented')
         else:
             if input:
                 iseq, objclass, format = entryIteratorFactory(file,
@@ -556,7 +608,7 @@ def open_uri(uri,
                                                               commentchar)
             else:    # default export is in fasta? or tab? TODO
                 objclass = Nuc_Seq   # Nuc_Seq_Stored? TODO
-                iseq = FastaNucWriter(FastaFormat(printNAKeys=printna, NAString=nastring), 
+                iseq = FastaNucWriter(FastaFormat(tags=only_keys, printNAKeys=printna, NAString=nastring), 
                                       file,
                                       skip=skip,
                                       only=only)
@@ -565,6 +617,6 @@ def open_uri(uri,
         
         entry_count = -1
         if input:
-            entry_count = count_entries(file, format)
+            entry_count = count_entries(file, format, header)
         
         return (file, iseq, objclass, urib, entry_count)

python/obitools3/utils.pxd

@@ -3,7 +3,7 @@
 from obitools3.dms.capi.obitypes cimport obitype_t, index_t
 
 cpdef bytes format_uniq_pattern(bytes format)
-cpdef int count_entries(file, bytes format)
+cpdef int count_entries(file, bytes format, bint header)
 
 cdef obi_errno_to_exception(index_t line_nb=*, object elt_id=*, str error_message=*)
 
@@ -18,7 +18,7 @@ cdef object clean_empty_values_from_object(object value, exclude=*)
 
 cdef obitype_t get_obitype_single_value(object value)
 cdef obitype_t update_obitype(obitype_t obitype, object new_value)
-cdef obitype_t get_obitype_iterable_value(object value)
+cdef obitype_t get_obitype_iterable_value(object value, type t)
 cdef obitype_t get_obitype(object value)
 
 cdef object __etag__(bytes x, bytes nastring=*)

python/obitools3/utils.pyx

@@ -9,7 +9,8 @@ from obitools3.dms.capi.obitypes cimport is_a_DNA_seq, \
                                          OBI_QUAL, \
                                          OBI_SEQ, \
                                          OBI_STR, \
-                                         index_t
+                                         index_t, \
+                                         OBI_INT_MAX
 
 from obitools3.dms.capi.obierrno cimport OBI_LINE_IDX_ERROR, \
                                          OBI_ELT_IDX_ERROR, \
@@ -39,7 +40,7 @@ cpdef bytes format_uniq_pattern(bytes format):
         return None
 
 
-cpdef int count_entries(file, bytes format):
+cpdef int count_entries(file, bytes format, bint header):
     
     try:
         sep = format_uniq_pattern(format)
@@ -74,6 +75,8 @@ cpdef int count_entries(file, bytes format):
             total_count += len(re.findall(sep, mmapped_file))
             if format != b"ngsfilter" and format != b"tabular" and format != b"embl" and format != b"genbank" and format != b"fastq":
                 total_count += 1 # adding +1 for 1st entry because separators include \n (ngsfilter and tabular already count one more because of last \n)
+            if format == b"tabular" and header: # not counting header as an entry
+                total_count -= 1
             
     except:
         if len(files) > 1:
@@ -257,38 +260,51 @@ cdef obitype_t update_obitype(obitype_t obitype, object new_value) :
     
     new_type = type(new_value)
     
-    if obitype == OBI_INT :
-        if new_type == float :
-            return OBI_FLOAT
-    # TODO BOOL vers INT/FLOAT
-    elif new_type == str or new_type == bytes :
+    #if new_type == NoneType:  # doesn't work because Cython sucks
+    if new_value == None or new_type==list or new_type==dict or new_type==tuple:
+        return obitype
+    
+    # TODO BOOL to INT/FLOAT
+    if new_type == str or new_type == bytes :
         if obitype == OBI_SEQ and is_a_DNA_seq(tobytes(new_value)) :
             pass
         else :
             return OBI_STR
+    elif obitype == OBI_INT :
+        if new_type == float or new_value > OBI_INT_MAX :
+            return OBI_FLOAT
         
     return obitype
 
 
-cdef obitype_t get_obitype_iterable_value(object value) :
+cdef obitype_t get_obitype_iterable_value(object value, type t) :
     
     cdef obitype_t value_obitype
     
     value_obitype = OBI_VOID
     
+    if t == dict:
         for k in value :
             if value_obitype == OBI_VOID :
                 value_obitype = get_obitype_single_value(value[k])
             else :
                 value_obitype = update_obitype(value_obitype, value[k])
     
+    elif t == list or t == tuple:
+        for v in value :
+            if value_obitype == OBI_VOID :
+                value_obitype = get_obitype_single_value(v)
+            else :
+                value_obitype = update_obitype(value_obitype, v)
+    
     return value_obitype
 
 
 cdef obitype_t get_obitype(object value) :
     
-    if type(value) == dict or type(value) == list or type(value) == tuple :
-        return get_obitype_iterable_value(value)
+    t = type(value)
+    if t == dict or t == list or t == tuple :
+        return get_obitype_iterable_value(value, t)
     
     else :
         return get_obitype_single_value(value)

python/obitools3/version.py

@@ -1,5 +1,5 @@
 major = 3
 minor = 0
-serial= '1b4'
+serial= '1b14'
 
 version ="%d.%d.%s" % (major,minor,serial)

python/obitools3/writers/tab.pyx

@@ -20,8 +20,6 @@ cdef class TabWriter:
             self.only = -1
         else:
             self.only = int(only)
-            if header:
-                self.only += 1
 
         self.formatter = formatter
         self.output = output_object

src/build_reference_db.c

@@ -77,6 +77,7 @@ static inline ecotx_t* get_lca_from_merged_taxids(Obiview_p view, OBIDMS_column_
 {
 	ecotx_t* taxon = NULL;
 	ecotx_t* lca = NULL;
+	ecotx_t* lca1 = NULL;
 	int32_t taxid;
 	index_t taxid_idx;
 	int64_t taxid_str_idx;
@@ -108,10 +109,11 @@ static inline ecotx_t* get_lca_from_merged_taxids(Obiview_p view, OBIDMS_column_
 			else
 			{
 				// Compute LCA
+				lca1 = lca;
 				lca = obi_taxo_get_lca(taxon, lca);
 				if (lca == NULL)
 				{
-					obidebug(1, "\nError getting the last common ancestor of two taxa when building a reference database");
+					obidebug(1, "\nError getting the last common ancestor of two taxa when building a reference database, %d %d", taxid, lca1->taxid);
 					return NULL;
 				}
 			}
@@ -185,7 +187,7 @@ int build_reference_db(const char* dms_name,
 	matrix_view_name = strcpy(matrix_view_name, o_view_name);
 	strcat(matrix_view_name, "_matrix");
 
-	fprintf(stderr, "Aligning queries with reference database...\n");
+	fprintf(stderr, "Aligning sequences...\n");
 	if (obi_lcs_align_one_column(dms_name,
 								 refs_view_name,
 								 "",
@@ -863,7 +865,8 @@ int build_reference_db(const char* dms_name,
 	fprintf(stderr,"\rDone : 100 %%           \n");
 
 	// Add information about the threshold used to build the DB
-	snprintf(threshold_str, 5, "%f", threshold);
+#define snprintf_nowarn(...) (snprintf(__VA_ARGS__) < 0 ? abort() : (void)0)
+	snprintf_nowarn(threshold_str, 5, "%f", threshold);
 
 	new_comments = obi_add_comment((o_view->infos)->comments, DB_THRESHOLD_KEY_IN_COMMENTS, threshold_str);
 	if (new_comments == NULL)

src/kmer_similarity.c

@@ -77,7 +77,6 @@ Obi_ali_p kmer_similarity(Obiview_p view1, OBIDMS_column_p column1, index_t idx1
 	int32_t*        shift_count_array;
 	Obi_ali_p       ali = NULL;
 	int 			i, j;
-	bool 			switched_seqs;
 	bool			reversed;
 	int 			score = 0;
 	Obi_blob_p      blob1 = NULL;
@@ -124,6 +123,8 @@ Obi_ali_p kmer_similarity(Obiview_p view1, OBIDMS_column_p column1, index_t idx1
 	bool			keep_seq2_start;
 	bool			keep_seq1_end;
 	bool			keep_seq2_end;
+	bool			left_ali;
+	bool 			rev_quals = false;
 
 	// Check kmer size
 	if ((kmer_size < 1) || (kmer_size > 4))
@@ -148,19 +149,8 @@ Obi_ali_p kmer_similarity(Obiview_p view1, OBIDMS_column_p column1, index_t idx1
 	}
 
 	// Choose the shortest sequence to save kmer positions in array
-	switched_seqs = false;
 	len1 = blob1->length_decoded_value;
 	len2 = blob2->length_decoded_value;
-	if (len2 < len1)
-	{
-		switched_seqs = true;
-		temp_len = len1;
-		len1 = len2;
-		len2 = temp_len;
-		temp_blob = blob1;
-		blob1 = blob2;
-		blob2 = temp_blob;
-	}
 
 	// Force encoding on 2 bits by replacing ambiguous nucleotides by 'a's
 	if (blob1->element_size == 4)
@@ -196,7 +186,47 @@ Obi_ali_p kmer_similarity(Obiview_p view1, OBIDMS_column_p column1, index_t idx1
 	else
 		reversed = false;
 	if (reversed)
-		switched_seqs = !switched_seqs;
+		// unreverse to make cases simpler. Costly but rare (direct match is reverse primer match)
+	{
+		if (seq2 == NULL)
+			seq2 = obi_blob_to_seq(blob2);
+		seq2 = reverse_complement_sequence(seq2);
+		blob2 = obi_seq_to_blob(seq2);
+
+		if (seq1 == NULL)
+			seq1 = obi_blob_to_seq(blob1);
+		seq1 = reverse_complement_sequence(seq1);
+		blob1 = obi_seq_to_blob(seq1);
+		free_blob1 = true;
+
+		// Get the quality arrays
+		qual1 = obi_get_qual_int_with_elt_idx_and_col_p_in_view(view1, qual_col1, idx1, 0, &qual1_len);
+		if (qual1 == NULL)
+		{
+			obidebug(1, "\nError getting the quality of the 1st sequence when computing the kmer similarity between two sequences");
+			return NULL;
+		}
+		qual2 = obi_get_qual_int_with_elt_idx_and_col_p_in_view(view2, qual_col2, idx2, 0, &qual2_len);
+		if (qual2 == NULL)
+		{
+			obidebug(1, "\nError getting the quality of the 2nd sequence when computing the kmer similarity between two sequences");
+			return NULL;
+		}
+
+		uint8_t* newqual1 = malloc(qual1_len*sizeof(uint8_t));
+		uint8_t* newqual2 = malloc(qual2_len*sizeof(uint8_t));
+
+		for (i=0;i<qual1_len;i++)
+			newqual1[i] = qual1[qual1_len-1-i];
+
+		for (i=0;i<qual2_len;i++)
+			newqual2[i] = qual2[qual2_len-1-i];
+
+		qual1 = newqual1;
+		qual2 = newqual2;
+
+		rev_quals = true;
+	}
 
 	// Save total length for the shift counts array
 	total_len = len1 + len2 + 1;   // +1 for shift 0
@@ -237,7 +267,7 @@ Obi_ali_p kmer_similarity(Obiview_p view1, OBIDMS_column_p column1, index_t idx1
 		 	return NULL;
 		}
 	}
-	else if (len1 >= shift_array_height)
+	else if (total_len >= shift_array_height)
 	{
 		shift_array_height = total_len;
 		*shift_array_p = (int32_t*) realloc(*shift_array_p, ARRAY_LENGTH * shift_array_height * sizeof(int32_t));
@@ -291,7 +321,7 @@ Obi_ali_p kmer_similarity(Obiview_p view1, OBIDMS_column_p column1, index_t idx1
 	*shift_array_height_p = shift_array_height;
 	*shift_count_array_length_p = shift_count_array_length;
 
-	// Fill array with positions of kmers in the shortest sequence
+	// Fill array with positions of kmers in the first sequence
 	encoding = blob1->element_size;
 	kmer_count = len1 - kmer_size + 1;
 	for (kmer_idx=0; kmer_idx < kmer_count; kmer_idx++)
@@ -310,7 +340,7 @@ Obi_ali_p kmer_similarity(Obiview_p view1, OBIDMS_column_p column1, index_t idx1
 		kmer_pos_array[(kmer*kmer_pos_array_height)+i] = kmer_idx;
 	}
 
-	// Compare positions of kmers between both sequences and store shifts
+	// Compare positions of kmers between both sequences and store shifts (a shift corresponds to pos2 - pos1)
 	kmer_count = blob2->length_decoded_value - kmer_size + 1;
 	for (kmer_idx=0; kmer_idx < kmer_count; kmer_idx++)
 	{
@@ -374,35 +404,42 @@ Obi_ali_p kmer_similarity(Obiview_p view1, OBIDMS_column_p column1, index_t idx1
 	// The 873863 cases of hell
 	if (best_shift > 0)
 	{
+		left_ali = false;
 		overlap_len = len2 - best_shift;
 		if (len1 <= overlap_len)
 		{
 			overlap_len = len1;
-			if (! switched_seqs)
 			keep_seq2_end = true;
-			else
-				keep_seq2_start = true;
-		}
-		else if (switched_seqs)
-		{
-			keep_seq2_start = true;
-			keep_seq1_end = true;
 		}
 	}
 	else if (best_shift < 0)
 	{
+		left_ali = true;
 	    overlap_len = len1 + best_shift;
-	    if (!switched_seqs)
+		if (len2 <= overlap_len)
+		{
+			overlap_len = len2;
+			keep_seq1_start = true;
+		}
+		else
 		{
 			keep_seq1_start = true;
 			keep_seq2_end = true;
 		}
 	}
-	else
+	else // if (best_shift == 0)
+	{
+		if (len2 >= len1)
 		{
 			overlap_len = len1;
-		if ((!switched_seqs) && (len2 > len1))
 			keep_seq2_end = true;
+			left_ali = true;
+		}
+		else
+		{
+			overlap_len = len2;
+			left_ali = false; // discussable
+		}
 	}
 
 	ali = (Obi_ali_p) malloc(sizeof(Obi_ali_t));
@@ -433,7 +470,7 @@ Obi_ali_p kmer_similarity(Obiview_p view1, OBIDMS_column_p column1, index_t idx1
 		ali->direction[0] = '\0';
 	else
 	{
-		if (((best_shift <= 0) && (!switched_seqs)) || ((best_shift > 0) && switched_seqs))
+		if (left_ali)
 			strcpy(ali->direction, "left");
 		else
 			strcpy(ali->direction, "right");
@@ -441,6 +478,8 @@ Obi_ali_p kmer_similarity(Obiview_p view1, OBIDMS_column_p column1, index_t idx1
 
 	// Build the consensus sequence if asked
 	if (build_consensus)
+	{
+		if (! rev_quals)
 		{
 			// Get the quality arrays
 			qual1 = obi_get_qual_int_with_elt_idx_and_col_p_in_view(view1, qual_col1, idx1, 0, &qual1_len);
@@ -455,15 +494,13 @@ Obi_ali_p kmer_similarity(Obiview_p view1, OBIDMS_column_p column1, index_t idx1
 				obidebug(1, "\nError getting the quality of the 2nd sequence when computing the kmer similarity between two sequences");
 				return NULL;
 			}
+		}
 
 		// Decode the first sequence if not already done
 		if (seq1 == NULL)
 			seq1 = obi_blob_to_seq(blob1);
 
-		if (! switched_seqs)
 		consensus_len = len2 - best_shift;
-		else
-		    consensus_len = len1 + best_shift;
 
 		// Allocate memory for consensus sequence
 		consensus_seq = (char*) malloc(consensus_len + 1 * sizeof(char)); // TODO keep malloced too maybe
@@ -557,6 +594,12 @@ Obi_ali_p kmer_similarity(Obiview_p view1, OBIDMS_column_p column1, index_t idx1
 	free(seq2);
 	free(blob2);
 
+	if (rev_quals)
+	{
+		free(qual1);
+		free(qual2);
+	}
+
 	return ali;
 }
 

src/obi_ecopcr.c

@@ -645,7 +645,8 @@ static int print_seq(Obiview_p i_view, Obiview_p o_view,
 
 int obi_ecopcr(const char* i_dms_name,
 		       const char* i_view_name,
-			   const char* taxonomy_name,     // TODO discuss that input dms assumed
+			   const char* tax_dms_name,
+			   const char* taxonomy_name,
 			   const char* o_dms_name,
 			   const char* o_view_name,
 			   const char* o_view_comments,
@@ -678,6 +679,7 @@ int obi_ecopcr(const char* i_dms_name,
 
 	OBIDMS_p          i_dms = NULL;
 	OBIDMS_p          o_dms = NULL;
+	OBIDMS_p          tax_dms = NULL;
 	OBIDMS_taxonomy_p taxonomy = NULL;
 	Obiview_p         i_view = NULL;
 	Obiview_p         o_view = NULL;
@@ -965,8 +967,16 @@ int obi_ecopcr(const char* i_dms_name,
 		return -1;
 	}
 
+	// Open taxonomy DMS
+	tax_dms = obi_open_dms(tax_dms_name, false);
+	if (tax_dms == NULL)
+	{
+		obidebug(1, "\nError opening the taxonomy DMS");
+		return -1;
+	}
+
 	// Open the taxonomy
-	taxonomy = obi_read_taxonomy(i_dms, taxonomy_name, false);
+	taxonomy = obi_read_taxonomy(tax_dms, taxonomy_name, false);
 	if (taxonomy == NULL)
 	{
 		obidebug(1, "\nError opening the taxonomy");

src/obi_ecopcr.h

@@ -77,7 +77,8 @@
  *
  * @param i_dms_name The path to the input DMS.
  * @param i_view_name The name of the input view.
- * @param taxonomy_name The name of the taxonomy in the input DMS.
+ * @param tax_dms_name The path to the DMS containing the taxonomy.
+ * @param taxonomy_name The name of the taxonomy.
  * @param o_dms_name The path to the output DMS.
  * @param o_view_name The name of the output view.
  * @param o_view_comments The comments to associate with the output view.
@@ -106,6 +107,7 @@
  */
 int obi_ecopcr(const char* i_dms_name,
 		       const char* i_view_name,
+			   const char* tax_dms_name,
 			   const char* taxonomy_name,
 			   const char* o_dms_name,
 			   const char* o_view_name,

src/obi_ecotag.c

@@ -233,7 +233,7 @@ int obi_ecotag(const char* dms_name,
 		// Write result (max score, threshold, LCA assigned, list of the ids of the best matches)
 
 
-	index_t         i, j, k;
+	index_t         i, j, k, t;
 	ecotx_t*        lca;
 	ecotx_t*        lca_in_array;
 	ecotx_t*        best_match;
@@ -259,16 +259,20 @@ int obi_ecotag(const char* dms_name,
 	int32_t*   		best_match_taxids;
 	int32_t*		best_match_taxids_to_store;
 	int				best_match_count;
+	int				best_match_taxid_count;
 	int             buffer_size;
 	int 			best_match_ids_buffer_size;
 	index_t         best_match_idx;
 	int32_t         lca_array_length;
 	int32_t         lca_taxid;
 	int32_t         taxid_best_match;
+	int32_t         taxid;
+	int32_t         taxid_to_store;
 	bool   			assigned;
 	const char*		lca_name;
 	const char*		id;
 	int				id_len;
+	bool			already_in;
 
 	OBIDMS_p          dms = NULL;
 	OBIDMS_p		  ref_dms = NULL;
@@ -488,10 +492,11 @@ int obi_ecotag(const char* dms_name,
 
 	for (i=0; i < query_count; i++)
 	{
-		if (i%1000 == 0)
+		if (i%10 == 0)
 			fprintf(stderr,"\rDone : %f %%       ", (i / (float) query_count)*100);
 
 		best_match_count = 0;
+		best_match_taxid_count = 0;
 		best_match_ids_length = 0;
 		threshold = ecotag_threshold;
 		best_score = 0.0;
@@ -543,6 +548,7 @@ int obi_ecotag(const char* dms_name,
 					// Reset the array with that match
 					best_match_ids_length = 0;
 					best_match_count = 0;
+					best_match_taxid_count = 0;
 				}
 
 				// Store in best match array
@@ -585,8 +591,27 @@ int obi_ecotag(const char* dms_name,
 
 				// Save match
 				best_match_array[best_match_count] = j;
-				best_match_taxids[best_match_count] = obi_get_int_with_elt_idx_and_col_p_in_view(ref_view, ref_taxid_column, j, 0);
 				best_match_count++;
+
+				// Save best match taxid only if not already in array
+				taxid_to_store = obi_get_int_with_elt_idx_and_col_p_in_view(ref_view, ref_taxid_column, j, 0);
+				already_in = false;
+				for (t=0; t<best_match_taxid_count; t++)
+				{
+					taxid = best_match_taxids[t];
+					//fprintf(stderr, "\ntaxid %d, taxid_to_store %d\n", taxid, taxid_to_store);
+					if (taxid == taxid_to_store)
+					{
+						already_in = true;
+						break;
+					}
+				}
+				if (! already_in)
+				{
+					best_match_taxids[best_match_taxid_count] = taxid_to_store;
+					best_match_taxid_count++;
+				}
+
 				strcpy(best_match_ids+best_match_ids_length, id);
 				best_match_ids_length = best_match_ids_length + id_len + 1;
 			}
@@ -693,7 +718,7 @@ int obi_ecotag(const char* dms_name,
 									assigned_name_column, lca_name,
 									assigned_status_column, assigned,
 									best_match_ids_column, best_match_ids_to_store, best_match_ids_length,
-									best_match_taxids_column, best_match_taxids_to_store, best_match_count,
+									best_match_taxids_column, best_match_taxids_to_store, best_match_taxid_count,
 									score_column, best_score
 									) < 0)
 							return -1;

src/obidms.c

@@ -1417,7 +1417,7 @@ char* obi_dms_formatted_infos(OBIDMS_p dms, bool detailed)
 	char*          view_name = NULL;
 	char*          tax_name = NULL;
 	char*          all_tax_dir_path = NULL;
-	int            i;
+	int            i, last_dot_pos;
 	struct dirent* dp;
 	Obiview_p      view;
 
@@ -1439,17 +1439,21 @@ char* obi_dms_formatted_infos(OBIDMS_p dms, bool detailed)
 		if ((dp->d_name)[0] == '.')
 			continue;
 		i=0;
-		while ((dp->d_name)[i] != '.')
+		while (i < strlen(dp->d_name))
+		{
+			if ((dp->d_name)[i] == '.')
+				last_dot_pos = i;
 			i++;
-		view_name = (char*) malloc((i+1) * sizeof(char));
+		}
+		view_name = (char*) malloc((last_dot_pos+1) * sizeof(char));
 		if (view_name == NULL)
 		{
 			obi_set_errno(OBI_MALLOC_ERROR);
 			obidebug(1, "\nError allocating memory for a view name when getting formatted DMS infos: file %s", dp->d_name);
 			return NULL;
 		}
-		strncpy(view_name, dp->d_name, i);
-		view_name[i] = '\0';
+		strncpy(view_name, dp->d_name, last_dot_pos);
+		view_name[last_dot_pos] = '\0';
 		view = obi_open_view(dms, view_name);
 		if (view == NULL)
 		{

src/obidms_taxonomy.c

@@ -873,7 +873,7 @@ static ecotxidx_t* read_taxonomy_idx(const char* taxa_file_name, const char* loc
 	taxa_index->buffer_size = taxa_index->count;
 
 	taxa_index->max_taxid = 0;
-	printf("Reading %d taxa...\n", count_taxa);
+	fprintf(stderr, "Reading %d taxa...\n", count_taxa);
 	for (i=0; i<count_taxa; i++)
 	{
 		readnext_ecotaxon(f_taxa, &(taxa_index->taxon[i]));
@@ -886,9 +886,9 @@ static ecotxidx_t* read_taxonomy_idx(const char* taxa_file_name, const char* loc
 	}
 
 	if (count_local_taxa > 0)
-		printf("Reading %d local taxa...\n", count_local_taxa);
+		fprintf(stderr, "Reading %d local taxa...\n", count_local_taxa);
 	else
-		printf("No local taxa\n");
+		fprintf(stderr, "No local taxa\n");
 
 	count_taxa = taxa_index->count;
 
@@ -2463,6 +2463,32 @@ int read_merged_dmp(const char* taxdump, OBIDMS_taxonomy_p tax, int32_t* delnode
 				return -1;
 			}
 
+			// Write the rest of the taxa from the current taxa list
+			while (nT < (tax->taxa)->count)
+			{
+			// Add element from taxa list
+				// Enlarge structure if needed
+				if (n == buffer_size)
+				{
+					buffer_size = buffer_size * 2;
+					tax->merged_idx = (ecomergedidx_t*) realloc(tax->merged_idx, sizeof(ecomergedidx_t) + sizeof(ecomerged_t) * buffer_size);
+					if (tax->merged_idx == NULL)
+					{
+						obi_set_errno(OBI_MALLOC_ERROR);
+						obidebug(1, "\nError reallocating memory for a taxonomy structure");
+						closedir(tax_dir);
+						fclose(file);
+						return -1;
+					}
+				}
+
+				(tax->merged_idx)->merged[n].taxid = (tax->taxa)->taxon[nT].taxid;
+				(tax->merged_idx)->merged[n].idx = nT;
+
+				nT++;
+				n++;
+			}
+
 			// Store count
 			(tax->merged_idx)->count = n;
 

src/obidms_taxonomy.h

@@ -75,7 +75,7 @@ typedef struct {
  	 	 	 	 	 	 	 */
 	int32_t max_taxid;		/**< Maximum taxid existing in the taxon index.
 	 	 	 	 	 	 	 */
-	int32_t buffer_size;	/**< Number of taxa.	// TODO kept this but not sure of its use
+	int32_t buffer_size;	/**< .	// TODO kept this but not sure of its use
  	 	 	 	 	 	 	 */
 	ecotx_t taxon[];		/**< Taxon array.
  	 	 	 	 	 	 	 */

src/obidmscolumn.h

@@ -36,7 +36,7 @@
 										   	   	 */
 #define COLUMN_GROWTH_FACTOR (2)	 	  		/**< The growth factor when a column is enlarged.
                                 	   	   	   	 */
-#define MAXIMUM_LINE_COUNT (1000000000)   		/**< The maximum line count for the data of a column (1E9). //TODO
+#define MAXIMUM_LINE_COUNT (1000000000000)      /**< The maximum line count for the data of a column (1E12). //TODO
                                 	       	   	 */
 #define COMMENTS_MAX_LENGTH (4096)        		/**< The maximum length for comments.
  	 	 	 	 	 	 	 	 	       	   	 */

src/obitypes.h

@@ -29,6 +29,8 @@
 #define OBIQual_int_NA (NULL)		/**< NA value for the type OBI_QUAL if the quality is in integer format */
 #define OBITuple_NA (NULL)			/**< NA value for tuples of any type */
 
+#define OBI_INT_MAX (INT32_MAX)     /**< Maximum value for the type OBI_INT */
+
 
 /**
  * @brief enum for the boolean OBIType.

src/obiview.c

@@ -2910,7 +2910,7 @@ int obi_clean_unfinished_views(OBIDMS_p dms)
 		if ((dp->d_name)[0] == '.')
 			continue;
 		i=0;
-		while ((dp->d_name)[i] != '.')
+		while (strncmp((dp->d_name)+i, ".obiview", 8))
 			i++;
 		relative_path = (char*) malloc(strlen(VIEW_DIR_NAME) + strlen(dp->d_name) + 2);
 		strcpy(relative_path, VIEW_DIR_NAME);