Fixed stdout output

versions of python

python/obitools3/apps/optiongroups/__init__.py

@@ -57,6 +57,18 @@ def __addImportInputOption(optionManager):
                      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('--unite-input',
+                     action="store_const", dest="obi:inputformat",
+                     default=None,
+                     const=b'unite',
+                     help="Input file is in UNITE fasta format. If NCBI taxonomy provided with --taxonomy, taxid and scientific name will be added for each sequence.")
+
+    group.add_argument('--sintax-input',
+                     action="store_const", dest="obi:inputformat",
+                     default=None,
+                     const=b'sintax',
+                     help="Input file is in SINTAX 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",
                      default=None,
@@ -137,10 +149,10 @@ def __addImportInputOption(optionManager):
 def __addTabularOption(optionManager):
     group = optionManager.add_argument_group("Input and output format options for tabular files")
 
-    group.add_argument('--header',
+    group.add_argument('--no-header',
-                     action="store_true", dest="obi:header",
+                     action="store_false", dest="obi:header",
-                     default=False,
+                     default=True,
-                     help="First line of tabular file contains column names")
+                     help="Don't print the header (first line with column names")
 
     group.add_argument('--sep',
                      action="store", dest="obi:sep",
@@ -177,6 +189,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")
 
@@ -210,6 +232,10 @@ def addTabularInputOption(optionManager):
     __addTabularInputOption(optionManager)
 
 
+def addTabularOutputOption(optionManager):
+    __addTabularOutputOption(optionManager)
+
+
 def addTaxonomyOption(optionManager):
     __addTaxonomyOption(optionManager)
 
@@ -222,6 +248,7 @@ def addAllInputOption(optionManager):
     __addInputOption(optionManager)
     __addImportInputOption(optionManager)
     __addTabularInputOption(optionManager)
+    __addTabularOutputOption(optionManager)
     __addTaxonomyOption(optionManager)
     __addTaxdumpInputOption(optionManager)    
 
@@ -282,6 +309,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,
@@ -314,14 +370,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/alignpairedend.pyx

@@ -233,6 +233,10 @@ def run(config):
         PyErr_CheckSignals()
 
         consensus = o_view[i]
+        
+        if two_views:
+            consensus[b"R1_parent"] = forward[i].id
+            consensus[b"R2_parent"] = reverse[i].id
                
         if not two_views:
             seqF = entries[i]

python/obitools3/commands/annotate.pyx

@@ -16,6 +16,8 @@ from obitools3.dms.capi.obiview cimport NUC_SEQUENCE_COLUMN, \
                                         QUALITY_COLUMN, \
                                         COUNT_COLUMN, \
                                         TAXID_COLUMN
+from obitools3.dms.capi.obitypes cimport OBI_STR
+from obitools3.dms.column.column cimport Column
 
 import time
 import math 
@@ -187,6 +189,8 @@ def sequenceTaggerGenerator(config, taxo=None):
                     else:
                         scn=None
                     seq[rank]=rtaxid
+                    if "%s_name"%rank not in seq.view:
+                        Column.new_column(seq.view, "%s_name"%rank, OBI_STR)
                     seq["%s_name"%rank]=scn
              
         if add_rank:

python/obitools3/commands/cat.pyx

@@ -134,7 +134,11 @@ def run(config):
                 rep = repr(entry)
                 output_0.write(str2bytes(rep)+b"\n")
             else:
-                o_view[i] = entry
+                try:
+                    o_view[i] = entry
+                except:
+                    print("\nError with entry:", repr(entry))
+                    print(repr(o_view))
             i+=1
         v.close()
 

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",

python/obitools3/commands/head.pyx

@@ -8,6 +8,7 @@ from obitools3.apps.optiongroups import addMinimalInputOption, addMinimalOutputO
 from obitools3.dms.view import RollbackException
 from obitools3.apps.config import logger
 from obitools3.utils cimport str2bytes
+from obitools3.apps.optiongroups import addExportOutputOption
 
 import time
 import sys
@@ -23,6 +24,7 @@ def addOptions(parser):
     
     addMinimalInputOption(parser)
     addMinimalOutputOption(parser)
+    addExportOutputOption(parser)
     addNoProgressBarOption(parser)
 
     group=parser.add_argument_group('obi head specific options')

python/obitools3/commands/import.pyx

@@ -2,6 +2,7 @@
 
 import sys
 import os
+import re
 
 from obitools3.apps.progress cimport ProgressBar  # @UnresolvedImport
 from obitools3.dms.view.view cimport View
@@ -108,6 +109,8 @@ def run(config):
     cdef   bint        NUC_SEQS_view
     cdef   bint        silva
     cdef   bint        rdp
+    cdef   bint        unite
+    cdef   bint        sintax
     cdef   int         nb_elts
     cdef   object      d
     cdef   View        view
@@ -233,16 +236,27 @@ def run(config):
         def_col = view[DEFINITION_COLUMN]
         seq_col = view[NUC_SEQUENCE_COLUMN]
 
-    # Prepare taxon scientific name and taxid refs if RDP or SILVA file
+    # Prepare taxon scientific name and taxid refs if RDP/SILVA/UNITE/SINTAX formats
     silva = False
     rdp = False
-    if 'inputformat' in config['obi'] and (config['obi']['inputformat'] == b"silva" or config['obi']['inputformat'] == b"rdp"):
+    unite = False
+    sintax=False
+    if 'inputformat' in config['obi'] and (config['obi']['inputformat'] == b"silva" or \
+                                           config['obi']['inputformat'] == b"rdp" or \
+                                           config['obi']['inputformat'] == b"unite" or \
+                                           config['obi']['inputformat'] == b"sintax"):
         #if taxo is None:
         #    raise Exception("A taxonomy (as built by 'obi import --taxdump') must be provided for SILVA and RDP files")
-        silva = True
+        if config['obi']['inputformat'] == b"silva":
-        rdp = True
+            silva = True
+        elif config['obi']['inputformat'] == b"rdp":
+            rdp = True
+        elif config['obi']['inputformat'] == b"unite":
+            unite = True
+        elif config['obi']['inputformat'] == b"sintax":
+            sintax = True
+        sci_name_col = Column.new_column(view, SCIENTIFIC_NAME_COLUMN, OBI_STR)
         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 = {}
@@ -349,27 +363,51 @@ def run(config):
                 if get_quality:
                     qual_col[i] = entry.quality
                     
-                # Parse taxon scientific name if RDP file
+                # Parse taxon scientific name if RDP or Silva or Unite file
-                if (rdp or silva) and taxo is not None:
+                if (rdp or silva or unite or sintax):
-                    sci_names = entry.definition.split(b";")
+                    if rdp or silva:
-                    for sci_name in reversed(sci_names):
+                        sci_names = entry.definition.split(b";")
-                        if sci_name.split()[0] != b'unidentified' and sci_name.split()[0] != b'uncultured' and sci_name.split()[0] != b'metagenome' :
+                        sci_name_col[i] = sci_names[-1]
-                            taxon = taxo.get_taxon_by_name(sci_name)
+                    elif unite:
-                            if taxon is not None:
+                        sci_names = entry.id.split(b'|')[-1].split(b';')
-                                sci_name_col[i] = taxon.name
+                        sci_name_col[i] = re.sub(b'[a-zA-Z]__', b'', sci_names[-1])
-                                taxid_col[i] = taxon.taxid
+                    elif sintax:
-                                #print(taxid_col[i], sci_name_col[i])
+                        reconstructed_line = entry.id+b' '+entry.definition[:-1]
-                                break
+                        splitted_reconstructed_line = reconstructed_line.split(b';')
+                        taxa = splitted_reconstructed_line[1].split(b'=')[1]
+                        taxa = splitted_reconstructed_line[1].split(b',')
+                        sci_names = []
+                        for t in taxa:
+                            tf = t.split(b':')[1]
+                            sci_names.append(tf)
+                        sci_name_col[i] = sci_names[-1]
+                        id_col[i] = reconstructed_line.split(b';')[0]    
+                        def_col[i] = reconstructed_line
+                        
+                    # Fond taxid if taxonomy provided
+                    if taxo is not None :
+                        for sci_name in reversed(sci_names):
+                            if unite:
+                                sci_name = re.sub(b'[a-zA-Z]__', b'', sci_name)
+                            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 :
                 
                 if tag != ID_COLUMN and tag != DEFINITION_COLUMN and tag != NUC_SEQUENCE_COLUMN and tag != QUALITY_COLUMN :  # TODO dirty 
-                                    
+                
                     value = entry[tag]
                     if tag == b"taxid":
                         tag = TAXID_COLUMN
                     if tag == b"count":
                         tag = COUNT_COLUMN
+                    if tag == b"scientific_name":
+                        tag = SCIENTIFIC_NAME_COLUMN
                     if tag[:7] == b"merged_":
                         tag = MERGED_PREFIX+tag[7:]
 
@@ -382,7 +420,7 @@ def run(config):
                             pass
                     
                     if tag not in dcols :
-                                                
+                    
                         value_type = type(value)
                         nb_elts = 1
                         value_obitype = OBI_VOID
@@ -480,7 +518,7 @@ def run(config):
                                 dcols[tag][0][i] = value
         
         except Exception as e:
-            print("\nCould not import sequence:", repr(entry), "(error raised:", e, ")")
+            print("\nCould not import sequence:\n", repr(entry), "\nError raised:", e, "\n/!\ Check if '--input-na-string' option needs to be set")
             if 'skiperror' in config['obi'] and not config['obi']['skiperror']:
                 raise e
             else:

python/obitools3/commands/ngsfilter.pyx

@@ -24,6 +24,9 @@ from cpython.exc cimport PyErr_CheckSignals
 from io import BufferedWriter
 
 
+MAX_PAT_LEN = 31
+
+
 __title__="Assign sequence records to the corresponding experiment/sample based on DNA tags and primers"
 
 
@@ -84,6 +87,8 @@ class Primer:
         @type direct:
         '''
         
+        assert len(sequence) <= MAX_PAT_LEN, "Primer %s is too long, 31 bp max" % sequence
+        
         assert sequence not in Primer.collection        \
             or Primer.collection[sequence]==taglength,  \
             "Primer %s must always be used with tags of the same length" % sequence
@@ -264,14 +269,18 @@ cdef tuple annotate(sequences, infos, no_tags, verbose=False):
     
     not_aligned = len(sequences) > 1
     sequences[0] = sequences[0].clone()
+    
+    if not_aligned:
+        sequences[0][b"R1_parent"] = sequences[0].id
+        sequences[0][b"R2_parent"] = sequences[1].id
 
     if not_aligned:
         sequences[1] = sequences[1].clone()
-        sequences[0][REVERSE_SEQUENCE_COLUMN] = sequences[1].seq             # used by alignpairedend tool
+        sequences[0][REVERSE_SEQUENCE_COLUMN] = sequences[1].seq        # used by alignpairedend tool
         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)
@@ -295,7 +304,7 @@ cdef tuple annotate(sequences, infos, no_tags, verbose=False):
             directmatch.append((p, p(seq, same_sequence=not new_seq, pattern=pattern), seq, p))
             new_seq = False
             pattern+=1
-    
+        
     # Choose match closer to the start of (one of the) sequence(s)
     directmatch = sorted(directmatch, key=sortMatch)
     all_direct_matches = directmatch

python/obitools3/commands/rm.pyx

@@ -5,7 +5,9 @@ from obitools3.apps.config import logger
 from obitools3.dms import DMS
 from obitools3.apps.optiongroups import addMinimalInputOption
 from obitools3.dms.view.view cimport View
+from obitools3.utils cimport tostr
 import os
+import shutil
  
 __title__="Delete a view"
 
@@ -30,15 +32,56 @@ def run(config):
     else: 
         raise NotImplementedError()
     
+    dms = input[0]
+    
     # Get the path to the view file to remove
-    path = input[0].full_path  # dms path
+    path = dms.full_path  # dms path
-    path+=b"/VIEWS/"
+    view_path=path+b"/VIEWS/"
-    path+=view.name
+    view_path+=view.name
-    path+=b".obiview"
+    view_path+=b".obiview"
+    
+    to_remove = {}
+    # For each column:
+    for col_alias in view.keys():
+        col = view[col_alias]
+        col_name = col.original_name
+        col_version = col.version
+        col_type = col.data_type
+        col_ref = (col_name, col_version)
+        # build file name and AVL file names
+        col_file_name = f"{tostr(path)}/{tostr(col.original_name)}.obicol/{tostr(col.original_name)}@{col.version}.odc"
+        if col_type in [b'OBI_CHAR', b'OBI_QUAL', b'OBI_STR', b'OBI_SEQ']:
+            avl_file_name = f"{tostr(path)}/OBIBLOB_INDEXERS/{tostr(col.original_name)}_{col.version}_indexer"
+        else:
+            avl_file_name = None
+        to_remove[col_ref] = [col_file_name, avl_file_name]
+    
+    # For each view:
+    do_not_remove = []
+    for vn in dms:
+        v = dms[vn]
+        # ignore the one being deleted
+        if v.name != view.name:
+            # check that none of the column is referenced, if referenced, remove from list to remove
+            cols = [(v[c].original_name, v[c].version) for c in v.keys()]
+            for col_ref in to_remove:
+                if col_ref in cols:
+                    do_not_remove.append(col_ref)
+    
+    for nr in do_not_remove:
+        to_remove.pop(nr)
     
     # Close the view and the DMS
     view.close()
     input[0].close(force=True)
     
-    # Rm
+    #print(to_remove)
-    os.remove(path)
+    
+    # rm AFTER view and DMS close
+    os.remove(view_path)
+    for col in to_remove:
+        os.remove(to_remove[col][0])
+        if to_remove[col][1] is not None:
+            shutil.rmtree(to_remove[col][1])
+        
+        

python/obitools3/commands/split.pyx

@@ -0,0 +1,105 @@
+#cython: language_level=3
+
+from obitools3.apps.progress cimport ProgressBar  # @UnresolvedImport
+from obitools3.dms import DMS
+from obitools3.dms.view.view cimport View, Line_selection
+from obitools3.uri.decode import open_uri
+from obitools3.apps.optiongroups import addMinimalInputOption, addTaxonomyOption, addMinimalOutputOption, addNoProgressBarOption
+from obitools3.dms.view import RollbackException
+from obitools3.apps.config import logger
+from obitools3.utils cimport tobytes
+
+import sys
+from cpython.exc cimport PyErr_CheckSignals
+
+ 
+__title__="Split"
+
+ 
+def addOptions(parser):
+    
+    addMinimalInputOption(parser)
+    addNoProgressBarOption(parser)
+
+    group=parser.add_argument_group("obi split specific options")
+ 
+    group.add_argument('-p','--prefix',
+                       action="store", dest="split:prefix",
+                       metavar="<PREFIX>",
+                       help="Prefix added to each subview name (included undefined)")
+
+    group.add_argument('-t','--tag-name',
+                       action="store", dest="split:tagname",
+                       metavar="<TAG_NAME>",
+                       help="Attribute/tag used to split the input")
+      
+    group.add_argument('-u','--undefined',
+                       action="store", dest="split:undefined",
+                       default=b'UNDEFINED',
+                       metavar="<VIEW_NAME>",
+                       help="Name of the view where undefined sequenced are stored (will be PREFIX_VIEW_NAME)")
+
+
+def run(config):
+     
+    DMS.obi_atexit()
+    
+    logger("info", "obi split")
+
+    # Open the input
+    input = open_uri(config["obi"]["inputURI"])
+    if input is None:
+        raise Exception("Could not read input view")
+    i_dms = input[0]
+    i_view = input[1]
+
+    # Initialize the progress bar
+    if config['obi']['noprogressbar'] == False:
+        pb = ProgressBar(len(i_view), config)
+    else:
+        pb = None
+    
+    tag_to_split = config["split"]["tagname"]
+    undefined = tobytes(config["split"]["undefined"])
+    selections = {}
+    
+    # Go through input view and split
+    for i in range(len(i_view)):
+        PyErr_CheckSignals()
+        if pb is not None:
+            pb(i)
+        line = i_view[i]
+        if tag_to_split not in line or line[tag_to_split] is None or len(line[tag_to_split])==0:
+            value = undefined
+        else:
+            value = line[tag_to_split]
+        if value not in selections:
+            selections[value] = Line_selection(i_view)
+        selections[value].append(i)
+        
+    if pb is not None:
+        pb(len(i_view), force=True)
+        print("", file=sys.stderr)
+
+    # Create output views with the line selection
+    try:
+        for cat in selections:
+            o_view_name = config["split"]["prefix"].encode()+cat
+            o_view = selections[cat].materialize(o_view_name)
+            # Save command config in View and DMS comments
+            command_line = " ".join(sys.argv[1:])
+            input_dms_name=[input[0].name]
+            input_view_name=[input[1].name]
+            if 'taxoURI' in config['obi'] and config['obi']['taxoURI'] is not None:
+                input_dms_name.append(config['obi']['taxoURI'].split("/")[-3])
+                input_view_name.append("taxonomy/"+config['obi']['taxoURI'].split("/")[-1])
+            o_view.write_config(config, "split", command_line, input_dms_name=input_dms_name, input_view_name=input_view_name)
+            o_view.close()
+    except Exception, e:
+        raise RollbackException("obi split error, rollbacking view: "+str(e), o_view)
+
+    i_dms.record_command_line(command_line)
+    i_dms.close(force=True)
+
+    logger("info", "Done.")
+

python/obitools3/commands/stats.pyx

@@ -119,9 +119,12 @@ def mean(values, options):
 
 def variance(v):
     if len(v)==1: 
-        return 0 
+        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):
@@ -285,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("%d" %catcount[c]+"\t", end="")
-        print("%9d" %totcount[c])
+        print("%d" %totcount[c]+"\t")
 
     input[0].close(force=True)
     

python/obitools3/commands/taxonomy.pyx

@@ -0,0 +1,230 @@
+#cython: language_level=3
+
+from obitools3.apps.progress cimport ProgressBar  # @UnresolvedImport
+from obitools3.dms import DMS
+from obitools3.dms.view.view cimport View, Line_selection
+from obitools3.uri.decode import open_uri
+from obitools3.apps.optiongroups import addMinimalInputOption, addTaxonomyOption, addMinimalOutputOption, addNoProgressBarOption
+from obitools3.dms.view import RollbackException
+from obitools3.dms.column.column cimport Column
+from functools import reduce
+from obitools3.apps.config import logger
+from obitools3.utils cimport tobytes, str2bytes, tostr
+from io import BufferedWriter
+from obitools3.dms.capi.obiview cimport NUC_SEQUENCE_COLUMN, \
+                                        ID_COLUMN, \
+                                        DEFINITION_COLUMN, \
+                                        QUALITY_COLUMN, \
+                                        COUNT_COLUMN, \
+                                        TAXID_COLUMN
+from obitools3.dms.capi.obitypes cimport OBI_INT
+from obitools3.dms.capi.obitaxonomy cimport MIN_LOCAL_TAXID
+import time
+import math
+import sys
+
+from cpython.exc cimport PyErr_CheckSignals
+
+
+__title__="Add taxa with a new generated taxid to an NCBI taxonomy database"
+
+
+
+def addOptions(parser):
+
+    addMinimalInputOption(parser)
+    addTaxonomyOption(parser)
+    addMinimalOutputOption(parser)
+    addNoProgressBarOption(parser)
+
+    group=parser.add_argument_group('obi taxonomy specific options')
+
+    group.add_argument('-n', '--taxon-name-tag',
+                       action="store",
+                       dest="taxonomy:taxon_name_tag",
+                       metavar="<SCIENTIFIC_NAME_TAG>",
+                       default=b"SCIENTIFIC_NAME",
+                       help="Name of the tag giving the scientific name of the taxon "
+                            "(default: 'SCIENTIFIC_NAME').")
+
+#    group.add_argument('-g', '--try-genus-match',
+#                     action="store_true", dest="taxonomy:try_genus_match",
+#                     default=False,
+#                     help="Try matching the first word of <SCIENTIFIC_NAME_TAG> when can't find corresponding taxid for a taxon. "
+#                          "If there is a match it is added in the 'parent_taxid' tag. (Can be used by 'obi taxonomy' to add the taxon under that taxid).")
+
+    group.add_argument('-a', '--restricting-ancestor',
+                       action="store",
+                       dest="taxonomy:restricting_ancestor",
+                       metavar="<RESTRICTING_ANCESTOR>",
+                       default=None,
+                       help="Enables to restrict the addition of taxids under an ancestor specified by its taxid.")
+
+    group.add_argument('-t', '--taxid-tag',
+                       action="store",
+                       dest="taxonomy:taxid_tag",
+                       metavar="<TAXID_TAG>",
+                       default=b"TAXID",
+                       help="Name of the tag to store the new taxid "
+                            "(default: 'TAXID').")
+
+    group.add_argument('-l', '--log-file',
+                       action="store",
+                       dest="taxonomy:log_file",
+                       metavar="<LOG_FILE>",
+                       default='',
+                       help="Path to a log file to write informations about added taxids.")
+
+
+def run(config):
+
+    DMS.obi_atexit()
+
+    logger("info", "obi taxonomy")
+
+    # Open the input
+    input = open_uri(config['obi']['inputURI'])
+    if input is None:
+        raise Exception("Could not read input view")
+    i_dms = input[0]
+    i_view = input[1]
+    i_view_name = input[1].name
+
+    # Open the output: only the DMS, as the output view is going to be created by cloning the input view
+    # (could eventually be done via an open_uri() argument)
+    output = open_uri(config['obi']['outputURI'],
+                      input=False,
+                      dms_only=True)
+    if output is None:
+        raise Exception("Could not create output view")
+    o_dms = output[0]
+    output_0 = output[0]
+    o_view_name = output[1]
+
+    # stdout output: create temporary view
+    if type(output_0)==BufferedWriter:
+        o_dms = i_dms
+        i=0
+        o_view_name = b"temp"
+        while o_view_name in i_dms: # Making sure view name is unique in output DMS
+            o_view_name = o_view_name+b"_"+str2bytes(str(i))
+            i+=1
+        imported_view_name = o_view_name
+
+    # If the input and output DMS are not the same, import the input view in the output DMS before cloning it to modify it
+    # (could be the other way around: clone and modify in the input DMS then import the new view in the output DMS)
+    if i_dms != o_dms:
+        imported_view_name = i_view_name
+        i=0
+        while imported_view_name in o_dms:  # Making sure view name is unique in output DMS
+            imported_view_name = i_view_name+b"_"+str2bytes(str(i))
+            i+=1
+        View.import_view(i_dms.full_path[:-7], o_dms.full_path[:-7], i_view_name, imported_view_name)
+        i_view = o_dms[imported_view_name]
+
+    # Clone output view from input view
+    o_view = i_view.clone(o_view_name)
+    if o_view is None:
+        raise Exception("Couldn't create output view")
+    i_view.close()
+
+    # Open taxonomy
+    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]
+
+    # Initialize the progress bar
+    if config['obi']['noprogressbar'] == False:
+        pb = ProgressBar(len(o_view), config)
+    else:
+        pb = None
+
+    try:
+        if config['taxonomy']['log_file']:
+            logfile = open(config['taxonomy']['log_file'], 'w')
+        else:
+            logfile = sys.stdout
+        if 'restricting_ancestor' in config['taxonomy']:
+            res_anc = int(config['taxonomy']['restricting_ancestor'])
+        else:
+            res_anc = None
+        taxid_column_name = config['taxonomy']['taxid_tag']
+        parent_taxid_column_name = "PARENT_TAXID" # TODO macro
+        taxon_name_column_name = config['taxonomy']['taxon_name_tag']
+        taxid_column = Column.new_column(o_view, taxid_column_name, OBI_INT)
+        if parent_taxid_column_name in o_view:
+            parent_taxid_column = o_view[parent_taxid_column_name]
+        else:
+            parent_taxid_column = None
+        #parent_taxid_column = Column.new_column(o_view, parent_taxid_column_name, OBI_INT)
+        taxon_name_column = o_view[taxon_name_column_name]
+
+        for i in range(len(o_view)):
+            PyErr_CheckSignals()
+            #if pb is not None:
+            #    pb(i)
+            taxon_name = taxon_name_column[i]
+            taxon = taxo.get_taxon_by_name(taxon_name, res_anc)
+            if taxon is not None:
+                taxid_column[i] = taxon.taxid
+                if logfile:
+                    print(f"Found taxon '{tostr(taxon_name)}' already existing with taxid {taxid_column[i]}", file=logfile)
+            else: # try finding genus or other parent taxon from the first word
+                #print(i, o_view[i].id)
+                if parent_taxid_column is not None and parent_taxid_column[i] is not None:
+                    taxid_column[i] = taxo.add_taxon(taxon_name, 'species', parent_taxid_column[i])
+                    if logfile:
+                        print(f"Adding taxon '{tostr(taxon_name)}' under provided parent {parent_taxid_column[i]} with taxid {taxid_column[i]}", file=logfile)
+                else:
+                    taxon_name_sp = taxon_name.split(b" ")
+                    taxon = taxo.get_taxon_by_name(taxon_name_sp[0], res_anc)
+                    if taxon is not None:
+                        parent_taxid_column[i] = taxon.taxid
+                        taxid_column[i] = taxo.add_taxon(taxon_name, 'species', taxon.taxid)
+                        if logfile:
+                            print(f"Adding taxon '{tostr(taxon_name)}' under '{tostr(taxon.name)}' ({taxon.taxid}) with taxid {taxid_column[i]}", file=logfile)
+                    else:
+                        taxid_column[i] = taxo.add_taxon(taxon_name, 'species', res_anc)
+                        if logfile:
+                            print(f"Adding taxon '{tostr(taxon_name)}' under provided restricting ancestor {res_anc} with taxid {taxid_column[i]}", file=logfile)
+
+        taxo.write(taxo.name, update=True)
+        
+    except Exception, e:
+        raise RollbackException("obi taxonomy error, rollbacking view: "+str(e), o_view)
+
+    #if pb is not None:
+    #    pb(i, force=True)
+    #    print("", file=sys.stderr)
+
+    #logger("info", "\nTaxa already in the taxonomy: "+str(found_count)+"/"+str(len(o_view))+" ("+str(round(found_count*100.0/len(o_view), 2))+"%)")
+    #logger("info", "\nParent taxids found: "+str(parent_found_count)+"/"+str(len(o_view))+" ("+str(round(parent_found_count*100.0/len(o_view), 2))+"%)")
+    #logger("info", "\nTaxids not found: "+str(not_found_count)+"/"+str(len(o_view))+" ("+str(round(not_found_count*100.0/len(o_view), 2))+"%)")
+
+    # Save command config in View and DMS comments
+    command_line = " ".join(sys.argv[1:])
+    input_dms_name=[input[0].name]
+    input_view_name=[i_view_name]
+    if 'taxoURI' in config['obi'] and config['obi']['taxoURI'] is not None:
+        input_dms_name.append(config['obi']['taxoURI'].split("/")[-3])
+        input_view_name.append("taxonomy/"+config['obi']['taxoURI'].split("/")[-1])
+    o_view.write_config(config, "taxonomy", command_line, input_dms_name=input_dms_name, input_view_name=input_view_name)
+    o_dms.record_command_line(command_line)
+
+    #print("\n\nOutput view:\n````````````", file=sys.stderr)
+    #print(repr(o_view), file=sys.stderr)
+
+    # stdout output: write to buffer
+    if type(output_0)==BufferedWriter:
+        logger("info", "Printing to output...")
+        o_view.print_to_output(output_0, noprogressbar=config['obi']['noprogressbar'])
+        o_view.close()
+
+    # If the input and the output DMS are different or if stdout output, delete the temporary imported view used to create the final view
+    if i_dms != o_dms or type(output_0)==BufferedWriter:
+        View.delete_view(o_dms, imported_view_name)
+        o_dms.close(force=True)
+    i_dms.close(force=True)
+
+    logger("info", "Done.")

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/obitaxonomy.pxd

@@ -58,7 +58,7 @@ cdef extern from "obidms_taxonomy.h" nogil:
 
     OBIDMS_taxonomy_p obi_read_taxdump(const_char_p taxdump)
     
-    int obi_write_taxonomy(OBIDMS_p dms, OBIDMS_taxonomy_p tax, const_char_p tax_name)
+    int obi_write_taxonomy(OBIDMS_p dms, OBIDMS_taxonomy_p tax, const_char_p tax_name, bint update)
 
     int obi_close_taxonomy(OBIDMS_taxonomy_p taxonomy)
     

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.pxd

@@ -19,8 +19,8 @@ cdef class Taxonomy(OBIWrapper) :
     cpdef Taxon get_taxon_by_idx(self, int idx)
     cpdef Taxon get_taxon_by_taxid(self, int taxid)
     cpdef Taxon get_taxon_by_name(self, object taxon_name, object restricting_taxid=*)
-    cpdef write(self, object prefix)
+    cpdef write(self, object prefix, bint update=*)
-    cpdef int add_taxon(self, str name, str rank_name, int parent_taxid, int min_taxid=*)
+    cpdef int add_taxon(self, object name, object rank_name, int parent_taxid, int min_taxid=*)
     cpdef object get_species(self, int taxid)
     cpdef object get_genus(self, int taxid)
     cpdef object get_family(self, int taxid)

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)
@@ -146,7 +150,9 @@ cdef class Taxonomy(OBIWrapper) :
         taxo._ranks = []
         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
 
     
@@ -168,6 +174,7 @@ cdef class Taxonomy(OBIWrapper) :
 
 
     cpdef Taxon get_taxon_by_name(self, object taxon_name, object restricting_taxid=None):            
+        #print(taxon_name)
         taxon = self._name_dict.get(tobytes(taxon_name), None)
         if not taxon:
             return None
@@ -276,12 +283,12 @@ cdef class Taxonomy(OBIWrapper) :
             yield Taxon(taxon_capsule, self)
 
 
-    cpdef write(self, object prefix) :
+    cpdef write(self, object prefix, bint update=False) :
-        if obi_write_taxonomy(self._dms.pointer(), self.pointer(), tobytes(prefix)) < 0 :
+        if obi_write_taxonomy(self._dms.pointer(), self.pointer(), tobytes(prefix), update) < 0 :
             raise Exception("Error writing the taxonomy to binary files")
     
     
-    cpdef int add_taxon(self, str name, str rank_name, int parent_taxid, int min_taxid=10000000) :
+    cpdef int add_taxon(self, object name, object rank_name, int parent_taxid, int min_taxid=10000000) :
         cdef int taxid
         taxid = obi_taxo_add_local_taxon(self.pointer(), tobytes(name), tobytes(rank_name), parent_taxid, min_taxid)
         if taxid < 0 :
@@ -304,6 +311,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):
         """
@@ -313,15 +325,17 @@ cdef class Taxonomy(OBIWrapper) :
         cdef Taxon taxon
         try:
             taxon = self.get_taxon_by_taxid(taxid)
-        except:
+        except Exception as e:
-            raise StopIteration
+            print('\n'+e, file=sys.stderr)
+            return 
         if taxon is not None:
             while taxon.taxid != 1:
                 yield taxon
+                #print(taxon.taxid)
                 taxon = taxon.parent
             yield taxon
         else:
-            raise StopIteration
+            return
 
 
     def is_ancestor(self, int ancestor_taxid, int taxid):

python/obitools3/dms/view/view.pyx

@@ -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) :
@@ -799,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,17 +59,25 @@ 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:
-                        value = value.bytes()
+                        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:
                             value=list(value)
-                        value = str2bytes(str(bytes2str_object(value))) # genius programming
+                    value = str2bytes(str(bytes2str_object(value))) # genius programming
                 if value is not None:
                     lines.append(k + b"=" + value + b";")   
                 

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 bytes sep
+    cdef bint NAIntTo0
+    cdef bint metabaR
+    cdef bint ngsfilter

python/obitools3/format/tab.pyx

@@ -4,57 +4,81 @@ 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 object ktags
+        cdef list tags = [key for key in data]
+        
         line = []
-        
+        if self.tags != None and self.tags:
-        if self.first_line:
+            ktags = list(self.tags)
-            self.tags = [k for k in data.keys()]
+        else:
-        
+            ktags = list(set(tags))   
+                     
+        ktags.sort()
+                        
         if self.header and self.first_line:
-            for k in self.tags:
+            for k in ktags:
-                if isinstance(data.view[k], Column_multi_elts):
+                if k in tags:
-                    keys = data.view[k].keys()
+                    if self.metabaR:
-                    keys.sort()
+                        if k == b'NUC_SEQ':
-                    for k2 in keys:
+                            ktoprint = b'sequence'
-                        line.append(tobytes(k)+b':'+tobytes(k2))
+                        else:
-                else:
+                            ktoprint = k.lower()
-                    line.append(tobytes(k))
+                        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(ktoprint)+b':'+tobytes(k2))
+                    else:
+                        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:
-            value = data[k]
+            if k in tags:
-            if isinstance(data.view[k], Column_multi_elts):
+                value = data[k]
-                keys = data.view[k].keys()
+                if isinstance(data.view[k], Column_multi_elts):
-                keys.sort()
+                    keys = data.view[k].keys()
-                if value is None:  # all keys at None
+                    keys.sort()
-                    for k2 in keys: # TODO could be much more efficient
+                    if value is None:  # all keys at None
-                        line.append(self.NAString)
+                        for k2 in keys: # TODO could be much more efficient
-                else:
-                    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:
                             line.append(self.NAString)
-            else:
+                    else:
-                if value is not None:
+                        for k2 in keys: # TODO could be much more efficient
-                    line.append(str2bytes(str(bytes2str_object(value))))
+                            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:
-                    line.append(self.NAString)
+                    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)
                   	      	
         if self.header and self.first_line:
             r += self.sep.join(value for value in line)

python/obitools3/parsers/fastq.pyx

@@ -103,7 +103,11 @@ def fastqWithQualityIterator(lineiterator,
         yield seq
 
         read+=1
-        hline = next(i)
+        try:
+            hline = next(i)
+        except StopIteration:
+            return
+
     
 
 def fastqWithoutQualityIterator(lineiterator, 
@@ -174,5 +178,7 @@ def fastqWithoutQualityIterator(lineiterator,
         yield seq
     
         read+=1
-        hline = next(i)
+        try:
-        
+            hline = next(i)
+        except StopIteration:
+            return        

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)
+_seqMatcher    = re.compile(b'^ORIGIN.+(?=//\n)', re.DOTALL + re.M)
-_cleanSeq1     = re.compile(b'ORIGIN.+\n')
+_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,10 +155,10 @@ def genbankIterator_file(lineiterator,
         yield seq
         read+=1
     
-    # Last sequence
+    # Last sequence if not empty lines
-    seq = genbankParser(entry)
+    if entry.strip():
-    
+        seq = genbankParser(entry)
-    yield seq
+        yield seq
     
     free(entry)
 

python/obitools3/parsers/ngsfilter.pyx

@@ -48,24 +48,25 @@ 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)
+        if split_line:
-        tags = tags.split(b":")
+            tags = split_line.pop(2)
-        for t_idx in range(len(tags)):
+            tags = tags.split(b":")
-            if tags[t_idx]==b"-" or tags[t_idx]==b"None" or tags[t_idx]==b"":
+            for t_idx in range(len(tags)):
-                tags[t_idx] = nastring
+                if tags[t_idx]==b"-" or tags[t_idx]==b"None" or tags[t_idx]==b"":
-        if len(tags) == 1:          # Forward and reverse tags are the same
+                    tags[t_idx] = nastring
-            tags.append(tags[0])
+            if len(tags) == 1:          # Forward and reverse tags are the same
-        split_line.insert(2, tags[0])
+                tags.append(tags[0])
-        split_line.insert(3, tags[1])
+            split_line.insert(2, tags[0])
-        new_lines.append(out_sep.join(split_line[0:6]))
+            split_line.insert(3, tags[1])
-    
+            new_lines.append(out_sep.join(split_line[0:6]))
+        
     return tabIterator(iter(new_lines),
                        header = True,
                        sep = out_sep,

python/obitools3/parsers/tab.pyx

@@ -99,7 +99,10 @@ def tabIterator(lineiterator,
         
         read+=1
         
-        line = next(iterator)
+        try:
+            line = next(iterator)
+        except StopIteration:
+            return 
     
         
         

python/obitools3/uri/decode.pyx

@@ -173,7 +173,10 @@ def open_uri(uri,
              type newviewtype=View,
              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
@@ -277,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 'outputformat' in config['obi'] and 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')
@@ -299,11 +307,11 @@ def open_uri(uri,
                     format=config["obi"][formatkey]
                 except KeyError:
                     format=None
-                
+                       
         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:
@@ -427,7 +435,21 @@ def open_uri(uri,
                 nastring=tobytes(config["obi"][nakey])
             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])
@@ -462,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" or format==b"rdp":
+                if format==b"fasta" or format==b"silva" or format==b"rdp" or format == b"unite" or format == b"sintax":
                     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)
@@ -485,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)
@@ -521,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,
@@ -539,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,
@@ -557,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)

python/obitools3/utils.pyx

@@ -264,7 +264,7 @@ cdef obitype_t update_obitype(obitype_t obitype, object new_value) :
     if new_value == None or new_type==list or new_type==dict or new_type==tuple:
         return obitype
     
-    # TODO BOOL vers INT/FLOAT
+    # 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

python/obitools3/version.py

@@ -1,5 +1,5 @@
 major = 3
 minor = 0
-serial= '1b9'
+serial= '1b26'
 
 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/.gitignore

@@ -3,3 +3,4 @@
 /cmake_install.cmake
 /libcobitools3.dylib
 /Makefile
+/build/

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,
 								 "",

src/obi_ecopcr.c

@@ -365,8 +365,6 @@ static int print_seq(Obiview_p i_view, Obiview_p o_view,
 
 	int32_t i;
 
-	// TODO add check for primer longer than MAX_PAT_LEN (32)
-
 	// Get sequence id
 	seq_id = obi_get_str_with_elt_idx_and_col_p_in_view(i_view, i_id_column, i_idx, 0);
 
@@ -645,7 +643,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 +677,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;
@@ -749,6 +749,20 @@ int obi_ecopcr(const char* i_dms_name,
 	o1c = complementPattern(o1);
 	o2c = complementPattern(o2);
 
+	// check for primers equal or longer than MAX_PAT_LEN (32)
+	if (strlen(primer1) >= MAX_PAT_LEN)
+	{
+		obi_set_errno(OBI_ECOPCR_ERROR);
+		obidebug(1, "\nError: first primer is too long, needs to be < 32bp (%s)", primer1);
+		return -1;
+	}
+	if (strlen(primer2) >= MAX_PAT_LEN)
+	{
+		obi_set_errno(OBI_ECOPCR_ERROR);
+		obidebug(1, "\nError: second primer is too long, needs to be < 32bp (%s)", primer2);
+		return -1;
+	}
+
 	// Open input DMS
 	i_dms = obi_open_dms(i_dms_name, false);
 	if (i_dms == NULL)
@@ -965,8 +979,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/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);
+		strncpy(view_name, dp->d_name, last_dot_pos);
-		view_name[i] = '\0';
+		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;
 
@@ -1092,7 +1092,7 @@ static int write_ranks_idx(OBIDMS_p dms, OBIDMS_taxonomy_p tax, const char* taxo
 	free(taxonomy_path);
 
 	// Create file
-	file_descriptor = open(file_name, O_RDWR | O_CREAT | O_EXCL, 0777);
+	file_descriptor = open(file_name, O_RDWR | O_CREAT, 0777);
 	if (file_descriptor < 0)
 	{
 		obi_set_errno(OBI_TAXONOMY_ERROR);
@@ -1196,7 +1196,7 @@ static int write_taxonomy_idx(OBIDMS_p dms, OBIDMS_taxonomy_p tax, const char* t
 	free(taxonomy_path);
 
 	// Create file
-	file_descriptor = open(file_name, O_RDWR | O_CREAT | O_EXCL, 0777);
+	file_descriptor = open(file_name, O_RDWR | O_CREAT, 0777);
 	if (file_descriptor < 0)
 	{
 		obi_set_errno(OBI_TAXONOMY_ERROR);
@@ -1472,7 +1472,7 @@ static int write_names_idx(OBIDMS_p dms, OBIDMS_taxonomy_p tax, const char* taxo
 	free(taxonomy_path);
 
 	// Create file
-	file_descriptor = open(file_name, O_RDWR | O_CREAT | O_EXCL, 0777);
+	file_descriptor = open(file_name, O_RDWR | O_CREAT, 0777);
 	if (file_descriptor < 0)
 	{
 		obi_set_errno(OBI_TAXONOMY_ERROR);
@@ -1760,7 +1760,7 @@ static int write_merged_idx(OBIDMS_p dms, OBIDMS_taxonomy_p tax, const char* tax
 	free(taxonomy_path);
 
 	// Create file
-	file_descriptor = open(file_name, O_RDWR | O_CREAT | O_EXCL, 0777);
+	file_descriptor = open(file_name, O_RDWR | O_CREAT, 0777);
 	if (file_descriptor < 0)
 	{
 		obi_set_errno(OBI_TAXONOMY_ERROR);
@@ -3250,47 +3250,48 @@ OBIDMS_taxonomy_p obi_read_taxonomy(OBIDMS_p dms, const char* taxonomy_name, boo
 }
 
 
-int obi_write_taxonomy(OBIDMS_p dms, OBIDMS_taxonomy_p tax, const char* tax_name)
+int obi_write_taxonomy(OBIDMS_p dms, OBIDMS_taxonomy_p tax, const char* tax_name, bool update)
 {
 	char* taxonomy_path;
 
-	// Build the taxonomy directory path
+	if (!update) {
-	taxonomy_path = get_taxonomy_path(dms, tax_name);
+		// Build the taxonomy directory path
-	if (taxonomy_path == NULL)
+		taxonomy_path = get_taxonomy_path(dms, tax_name);
-		return -1;
+		if (taxonomy_path == NULL)
-
+			return -1;
-	// Try to create the directory
+		// Try to create the directory
-	if (mkdir(taxonomy_path, 00777) < 0)
+		if (mkdir(taxonomy_path, 00777) < 0)
-	{
+		{
-		if (errno == EEXIST)
+			if (errno == EEXIST)
-			obidebug(1, "\nA taxonomy already exists with this name.");
+				obidebug(1, "\nA taxonomy already exists with this name.");
-		obidebug(1, "\nProblem creating a new taxonomy directory");
+			obidebug(1, "\nProblem creating a new taxonomy directory");
+			free(taxonomy_path);
+			return -1;
+		}
 		free(taxonomy_path);
-		return -1;
 	}
 
-	free(taxonomy_path);
+	if (write_ranks_idx(dms, tax, tax_name) < 0)
-
+		return -1;
-    if (write_ranks_idx(dms, tax, tax_name) < 0)
+	if (write_taxonomy_idx(dms, tax, tax_name) < 0)
-        return -1;
+		return -1;
-    if (write_taxonomy_idx(dms, tax, tax_name) < 0)
+	if (write_names_idx(dms, tax, tax_name) < 0)
-        return -1;
+		return -1;
-    if (write_names_idx(dms, tax, tax_name) < 0)
+	if (write_merged_idx(dms, tax, tax_name) < 0)
-    	return -1;
+		return -1;
-    if (write_merged_idx(dms, tax, tax_name) < 0)
+	// Write preferred names if there are some
-    	return -1;
+	if (tax->preferred_names != NULL)
-    // Check if there are local taxa (if so last taxon is local)
+	{
+		if (write_preferred_names_idx(dms, tax, tax_name) < 0)
+			return -1;
+	}
+    // Write local taxa if there are some
     if ((tax->taxa)->local_count > 0)
     {
     	if (write_local_taxonomy_idx(dms, tax, tax_name) < 0)
     		return -1;
     }
-    // Write preferred names if there are some
+
-    if (tax->preferred_names != NULL)
-    {
-    	if (write_preferred_names_idx(dms, tax, tax_name) < 0)
-    		return -1;
-    }
 	return 0;
 }
 
@@ -3302,16 +3303,17 @@ int obi_close_taxonomy(OBIDMS_taxonomy_p taxonomy)
 	if (taxonomy)
 	{
 		// Update local informations (local taxa and preferred names) if there are any
-		if ((taxonomy->taxa)->local_count > 0)
+		// Done with write_taxo, edits all needed files. Only ldx file was edited in OBI1 but it led to issues. Discussable
-		{
+//		if ((taxonomy->taxa)->local_count > 0)
-			if (taxonomy->dms == NULL)
+//		{
-			{
+//			if (taxonomy->dms == NULL)
-				obi_set_errno(OBI_TAXONOMY_ERROR);
+//			{
-				obidebug(1, "\nError closing a taxonomy with local files but no DMS associated (probably read directly from taxdump)");		// TODO discuss
+//				obi_set_errno(OBI_TAXONOMY_ERROR);
-			}
+//				obidebug(1, "\nError closing a taxonomy with local files but no DMS associated (probably read directly from taxdump)");		// TODO discuss
-			if (write_local_taxonomy_idx(taxonomy->dms, taxonomy, taxonomy->tax_name) < 0)
+//			}
-				return -1;
+//			if (write_local_taxonomy_idx(taxonomy->dms, taxonomy, taxonomy->tax_name) < 0)
-		}
+//				return -1;
+//		}
 
 		// Write preferred names if there are some
 		if (taxonomy->preferred_names)
@@ -3377,9 +3379,10 @@ int obi_close_taxonomy(OBIDMS_taxonomy_p taxonomy)
 int obi_taxo_add_local_taxon(OBIDMS_taxonomy_p tax, const char* name, const char* rank_name, int32_t parent_taxid, int32_t min_taxid)
 {
 	int32_t    taxid;
+	int32_t    count;
 	ecotx_t*   taxon;
 	int        i;
-//	econame_t* name_struct;
+	econame_t* name_struct;
 
 	// Enlarge the structure memory for a new taxon
 	tax->taxa = (ecotxidx_t*) realloc(tax->taxa, sizeof(ecotxidx_t) + sizeof(ecotx_t) * (((tax->taxa)->count) + 1));
@@ -3441,42 +3444,65 @@ int obi_taxo_add_local_taxon(OBIDMS_taxonomy_p tax, const char* name, const char
 	((tax->taxa)->local_count)++;
 	(tax->taxa)->buffer_size = (tax->taxa)->count;
 
-//	// Add new name in names structure     // Commented because the new name was not added in the .ndx file in the OBITools1
+	// Add new name in names structure // On the OBI1, the new name was not added in the .ndx file but it could create issues
-//	// Allocate memory for new name
+	// Allocate memory for new name
-//	tax->names = (econameidx_t*) realloc(tax->names, sizeof(econameidx_t) + sizeof(econame_t) * ((tax->names)->count + 1));
+	tax->names = (econameidx_t*) realloc(tax->names, sizeof(econameidx_t) + sizeof(econame_t) * ((tax->names)->count + 1));
-//	if (tax->names == NULL)
+	if (tax->names == NULL)
-//	{
+	{
-//		obi_set_errno(OBI_MALLOC_ERROR);
+		obi_set_errno(OBI_MALLOC_ERROR);
-//		obidebug(1, "\nError reallocating memory for a taxonomy structure to add a new taxon");
+		obidebug(1, "\nError reallocating memory for a taxonomy structure to add a new taxon");
-//		return -1;
+		return -1;
-//	}
+	}
-//
+
-//	// Add new name
+	// Add new name
-//	name_struct = (tax->names)->names + ((tax->names)->count);
+	name_struct = (tax->names)->names + ((tax->names)->count);
-//	name_struct->name = (char*) malloc((strlen(name) + 1) * sizeof(char));
+	name_struct->name = (char*) malloc((strlen(name) + 1) * sizeof(char));
-//	if (name_struct->name == NULL)
+	if (name_struct->name == NULL)
-//	{
+	{
-//		obi_set_errno(OBI_MALLOC_ERROR);
+		obi_set_errno(OBI_MALLOC_ERROR);
-//		obidebug(1, "\nError allocating memory for a taxon name to add a new taxon");
+		obidebug(1, "\nError allocating memory for a taxon name to add a new taxon");
-//		return -1;
+		return -1;
-//	}
+	}
-//	strcpy(name_struct->name, name);
+	strcpy(name_struct->name, name);
-//	name_struct->class_name = (char*) malloc((strlen("scientific name") + 1) * sizeof(char));
+	name_struct->class_name = (char*) malloc((strlen("scientific name") + 1) * sizeof(char));
-//	if (name_struct->class_name == NULL)
+	if (name_struct->class_name == NULL)
-//	{
+	{
-//		obi_set_errno(OBI_MALLOC_ERROR);
+		obi_set_errno(OBI_MALLOC_ERROR);
-//		obidebug(1, "\nError allocating memory for a taxon class name to add a new taxon");
+		obidebug(1, "\nError allocating memory for a taxon class name to add a new taxon");
-//		return -1;
+		return -1;
-//	}
+	}
-//	strcpy(name_struct->class_name, "scientific name");
+	strcpy(name_struct->class_name, "scientific name");
-//	name_struct->is_scientific_name = true;
+	name_struct->is_scientific_name = true;
-//	name_struct->taxon = ((tax->taxa)->taxon) + ((tax->taxa)->count) - 1;
+	name_struct->taxon = ((tax->taxa)->taxon) + ((tax->taxa)->count) - 1;
-//
+
-//	// Sort names in alphabetical order
+	// Update name count
-//	qsort((tax->names)->names, (tax->names)->count, sizeof(econame_t), cmp_names);
+	((tax->names)->count)++;
-//
+
-//	// Update name count
+	// Sort names in alphabetical order
-//	((tax->names)->count)++;
+	qsort((tax->names)->names, (tax->names)->count, sizeof(econame_t), cmp_names);
+
+	// Add to merged index
+	tax->merged_idx = (ecomergedidx_t*) realloc(tax->merged_idx, sizeof(ecomergedidx_t) + sizeof(ecomerged_t) * ((tax->merged_idx)->count + 1));
+	if (tax->merged_idx == NULL)
+	{
+		obi_set_errno(OBI_MALLOC_ERROR);
+		obidebug(1, "\nError reallocating memory for a taxonomy structure");
+		return -1;
+	}
+
+	count = (tax->merged_idx)->count;
+	(tax->merged_idx)->count = count + 1;
+	(tax->merged_idx)->merged[count].taxid = taxid;
+	(tax->merged_idx)->merged[count].idx = taxon->idx;
+
+	//fprintf(stderr, "\nEntered in merged taxon.idx=%d", (tax->merged_idx)->merged[(tax->merged_idx)->count -1].idx);
+	//fprintf(stderr, "\nEntered in merged taxon.taxid=%d", (tax->merged_idx)->merged[(tax->merged_idx)->count -1].taxid);
+	//fprintf(stderr, "\nEntered in merged at %d", (tax->merged_idx)->count -1);
+	//taxon = obi_taxo_get_taxon_with_taxid(tax, taxid);
+	//fprintf(stderr, "\ntaxon=%x", taxon);
+	//fprintf(stderr, "\ntaxon.taxid=%d", taxon->taxid);
+	//fprintf(stderr, "\ntaxon.name=%s", taxon->name);
+	//fprintf(stderr, "\ntaxon.idx=%d\n\n", ((tax->merged_idx)->count));
 
 	return taxid;
 }
@@ -3547,11 +3573,12 @@ int obi_taxo_add_preferred_name_with_taxon(OBIDMS_taxonomy_p tax, ecotx_t* taxon
 	name_struct->is_scientific_name = false;
 	name_struct->taxon = taxon;
 
+	// Update preferred name count
+	((tax->preferred_names)->count)++;
+
 	// Sort preferred names in alphabetical order
 	qsort((tax->preferred_names)->names, (tax->preferred_names)->count, sizeof(econame_t), cmp_names);
 
-	// Update preferred name count
-	((tax->preferred_names)->count)++;
 
 	return 0;
 }
@@ -3669,8 +3696,10 @@ ecotx_t* obi_taxo_get_taxon_with_taxid(OBIDMS_taxonomy_p taxonomy, int32_t taxid
 	else if (indexed_taxon->idx == -1)
 		current_taxon = NULL;	// TODO discuss what to do when old deleted taxon
 	else
+	{
 		current_taxon = (taxonomy->taxa->taxon)+(indexed_taxon->idx);
-
+		//fprintf(stderr, "\n>>>idx %d, taxid %d<<<\n", indexed_taxon->idx, indexed_taxon->taxid);
+	}
 	return current_taxon;
 }
 

src/obidms_taxonomy.h

@@ -239,6 +239,7 @@ OBIDMS_taxonomy_p obi_read_taxonomy(OBIDMS_p dms, const char* taxonomy_name, boo
  * @param dms A pointer on the DMS to which the taxonomy belongs.
  * @param tax A pointer on the taxonomy structure.
  * @param tax_name The name (prefix) of the taxonomy.
+ * @param update Whether files should be rewritten or if it's a new taxonomy (set to true e.g. after adding local taxa).
  *
  * @returns An integer value indicating the success of the operation.
  * @retval 0 on success.
@@ -247,7 +248,7 @@ OBIDMS_taxonomy_p obi_read_taxonomy(OBIDMS_p dms, const char* taxonomy_name, boo
  * @since 2016
  * @author Celine Mercier (celine.mercier@metabarcoding.org)
  */
-int obi_write_taxonomy(OBIDMS_p dms, OBIDMS_taxonomy_p tax, const char* tax_name);
+int obi_write_taxonomy(OBIDMS_p dms, OBIDMS_taxonomy_p tax, const char* tax_name, bool update);
 
 
 /**

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);