Tagging version beta3

Could use optimization

python/obitools3/apps/optiongroups/__init__.py

@@ -55,7 +55,7 @@ def __addImportInputOption(optionManager):
                      action="store_const", dest="obi:inputformat",
                      default=None,
                      const=b'ngsfilter',
-                     help="Input file is an ngsfilter file")
+                     help="Input file is an ngsfilter file. If not using tags, use ':' or 'None:None' or '-:-' or any combination")
 
     group.add_argument('--ecopcr-result-input',
                      action="store_const", dest="obi:inputformat",

python/obitools3/commands/clean_dms.pyx

@@ -21,7 +21,10 @@ def run(config):
     
     logger("info", "obi clean_dms")
     
-    if obi_clean_dms(tobytes(config['obi']['inputURI'])) < 0 :
+    dms_path = tobytes(config['obi']['inputURI'])
+    if b'.obidms' in dms_path:
+        dms_path = dms_path.split(b'.obidms')[0]
+    if obi_clean_dms(dms_path) < 0 :
         raise Exception("Error cleaning DMS", config['obi']['inputURI'])
 
     logger("info", "Done.")

python/obitools3/commands/ngsfilter.pyx

@@ -56,6 +56,11 @@ def addOptions(parser):
                      type=str,
                      default=None,
                      help="URI to the view used to store the sequences unassigned to any sample")
+
+    group.add_argument('--no-tags',
+                     action="store_true", dest="ngsfilter:notags",
+                     default=False,
+                     help="Use this option if your experiment does not use tags to identify samples")
     
     group.add_argument('-e','--error',
                      action="store", dest="ngsfilter:error",
@@ -167,7 +172,7 @@ cdef read_info_view(info_view, max_errors=2, verbose=False, not_aligned=False):
     i=0
     for p in info_view:
         forward=Primer(p[b'forward_primer'],
-                       len(p[b'forward_tag']) if p[b'forward_tag']!=b'-' else None,
+                       len(p[b'forward_tag']) if (b'forward_tag' in p and p[b'forward_tag']!=None) else None,
                        True,
                        max_errors=max_errors,
                        verbose=verbose,
@@ -178,7 +183,7 @@ cdef read_info_view(info_view, max_errors=2, verbose=False, not_aligned=False):
         infos[forward]=fp
         
         reverse=Primer(p[b'reverse_primer'],
-                       len(p[b'reverse_tag']) if p[b'reverse_tag']!=b'-' else None,
+                       len(p[b'reverse_tag']) if (b'reverse_tag' in p and p[b'reverse_tag']!=None) else None,
                        False,
                        max_errors=max_errors,
                        verbose=verbose,
@@ -213,10 +218,11 @@ cdef read_info_view(info_view, max_errors=2, verbose=False, not_aligned=False):
             rpp=rp.get(cf,{})
             rp[cf]=rpp
 
-        tags = (p[b'forward_tag'] if p[b'forward_tag']!=b'-' else None,
-                p[b'reverse_tag'] if p[b'reverse_tag']!=b'-' else None)
+        tags = (p[b'forward_tag'] if (b'forward_tag' in p and p[b'forward_tag']!=None) else None,
+                p[b'reverse_tag'] if (b'reverse_tag' in p and p[b'reverse_tag']!=None) else None)
         
-        assert tags not in dpp, \
+        if tags != (None, None):
+            assert tags not in dpp, \
                "Tag pair %s is already used with primer pairs: (%s,%s)" % (str(tags),forward,reverse)
         
         # Save additional data
@@ -234,7 +240,7 @@ cdef read_info_view(info_view, max_errors=2, verbose=False, not_aligned=False):
     return infos, primer_list
     
 
-cdef tuple annotate(sequences, infos, verbose=False):
+cdef tuple annotate(sequences, infos, no_tags, verbose=False):
         
     def sortMatch(match):
         if match[1] is None:
@@ -249,17 +255,12 @@ cdef tuple annotate(sequences, infos, verbose=False):
             return match[1][1]
     
     not_aligned = len(sequences) > 1
-    sequenceF = sequences[0]
-    sequenceR = None
-    if not not_aligned: 
-        final_sequence = sequenceF
-    else:
-        final_sequence = sequenceF.clone()   # TODO maybe not cloning and then deleting quality tags is more efficient
+    sequences[0] = sequences[0].clone()
 
     if not_aligned:
-        sequenceR = sequences[1]
-        final_sequence[REVERSE_SEQ_COLUMN_NAME] = sequenceR.seq             # used by alignpairedend tool
-        final_sequence[REVERSE_QUALITY_COLUMN_NAME] = sequenceR.quality     # used by alignpairedend tool
+        sequences[1] = sequences[1].clone()
+        sequences[0][REVERSE_SEQ_COLUMN_NAME] = sequences[1].seq             # used by alignpairedend tool
+        sequences[0][REVERSE_QUALITY_COLUMN_NAME] = sequences[1].quality     # used by alignpairedend tool
 
     for seq in sequences:
         if hasattr(seq, "quality_array"): 
@@ -275,8 +276,6 @@ cdef tuple annotate(sequences, infos, verbose=False):
     
     # Try direct matching:
     directmatch = []
-    first_matched_seq = None
-    second_matched_seq = None
     for seq in sequences:
         new_seq = True
         pattern = 0
@@ -295,60 +294,96 @@ cdef tuple annotate(sequences, infos, verbose=False):
     directmatch = directmatch[0] if directmatch[0][1] is not None else None
 
     if directmatch is None:
-        final_sequence[b'error']=b'No primer match'
-        return False, final_sequence
+        if not_aligned:
+            sequences[0][REVERSE_SEQ_COLUMN_NAME] = sequences[1].seq             # used by alignpairedend tool
+            sequences[0][REVERSE_QUALITY_COLUMN_NAME] = sequences[1].quality     # used by alignpairedend tool
+        sequences[0][b'error']=b'No primer match'
+        return False, sequences[0]
 
-    first_matched_seq = directmatch[2]
-    if id(first_matched_seq) == id(sequenceF) and not_aligned:
-        second_matched_seq = sequenceR
+    if id(directmatch[2]) == id(sequences[0]):
+        first_match_first_seq = True
     else:
-        second_matched_seq = sequenceF
-   
-    match = first_matched_seq[directmatch[1][1]:directmatch[1][2]]
+        first_match_first_seq = False
+       
+    match = directmatch[2][directmatch[1][1]:directmatch[1][2]]
     
     if not not_aligned:
-        final_sequence[b'seq_length_ori']=len(final_sequence)
+        sequences[0][b'seq_length_ori']=len(sequences[0])
     
-    if not not_aligned or id(first_matched_seq) == id(sequenceF):
-        final_sequence = final_sequence[directmatch[1][2]:]
+    if not not_aligned or first_match_first_seq:
+        sequences[0] = sequences[0][directmatch[1][2]:]
     else:
-        cut_seq = sequenceR[directmatch[1][2]:]
-        final_sequence[REVERSE_SEQ_COLUMN_NAME] = cut_seq.seq           # used by alignpairedend tool
-        final_sequence[REVERSE_QUALITY_COLUMN_NAME] = cut_seq.quality   # used by alignpairedend tool
+        sequences[1] = sequences[1][directmatch[1][2]:]
+        sequences[0][REVERSE_SEQ_COLUMN_NAME] = sequences[1].seq           # used by alignpairedend tool
+        sequences[0][REVERSE_QUALITY_COLUMN_NAME] = sequences[1].quality   # used by alignpairedend tool
     
     if directmatch[0].forward:
-        final_sequence[b'direction']=b'forward'
-        final_sequence[b'forward_errors']=directmatch[1][0]
-        final_sequence[b'forward_primer']=directmatch[0].raw
-        final_sequence[b'forward_match']=match.seq
+        sequences[0][b'direction']=b'forward'
+        sequences[0][b'forward_errors']=directmatch[1][0]
+        sequences[0][b'forward_primer']=directmatch[0].raw
+        sequences[0][b'forward_match']=match.seq
         
     else:
-        final_sequence[b'direction']=b'reverse'
-        final_sequence[b'reverse_errors']=directmatch[1][0]
-        final_sequence[b'reverse_primer']=directmatch[0].raw
-        final_sequence[b'reverse_match']=match.seq
+        sequences[0][b'direction']=b'reverse'
+        sequences[0][b'reverse_errors']=directmatch[1][0]
+        sequences[0][b'reverse_primer']=directmatch[0].raw
+        sequences[0][b'reverse_match']=match.seq
 
     # Keep only paired reverse primer
-    infos = infos[directmatch[0]]    
-    
+    infos = infos[directmatch[0]]
+    reverse_primer = list(infos.keys())[0]
+    direct_primer = directmatch[0]
+        
     # If not aligned, look for other match in already computed matches (choose the one that makes the biggest amplicon)
     if not_aligned:
         i=1
         # TODO comment
-        while i<len(all_direct_matches) and (all_direct_matches[i][1] is None or all_direct_matches[i][0].forward == directmatch[0].forward or all_direct_matches[i][0] == directmatch[0]):
+        while i<len(all_direct_matches) and \
+            (all_direct_matches[i][1] is None or \
+             all_direct_matches[i][0].forward == directmatch[0].forward or \
+             all_direct_matches[i][0] == directmatch[0] or \
+             reverse_primer != all_direct_matches[i][0]) :
             i+=1
         if i < len(all_direct_matches):
             reversematch = all_direct_matches[i]
         else:
             reversematch = None
+        
+    # Cut reverse primer out of 1st matched seq if it contains it, because if it's also in the other sequence, the next step will "choose" only the one on the other sequence
+    if not_aligned:
+        # do it on same seq
+        if first_match_first_seq:
+            r = reverse_primer.revcomp(sequences[0])
+        else:
+            r = reverse_primer.revcomp(sequences[1]) 
+        if r is not None: # found
+            if first_match_first_seq :
+                sequences[0] = sequences[0][:r[1]]
+            else:
+                sequences[1] = sequences[1][:r[1]]
+                sequences[0][REVERSE_SEQ_COLUMN_NAME] = sequences[1].seq           # used by alignpairedend tool
+                sequences[0][REVERSE_QUALITY_COLUMN_NAME] = sequences[1].quality   # used by alignpairedend tool
+        # do the same on the other seq
+        if first_match_first_seq: 
+            r = direct_primer.revcomp(sequences[1])
+        else:
+            r = direct_primer.revcomp(sequences[0])
+        if r is not None: # found
+            if first_match_first_seq:
+                sequences[1] = sequences[1][:r[1]]
+            else:
+                sequences[0] = sequences[0][:r[1]] 
+                sequences[0][REVERSE_SEQ_COLUMN_NAME] = sequences[1].seq
+                sequences[0][REVERSE_QUALITY_COLUMN_NAME] = sequences[1].quality
+    
     
     # Look for other primer in the other direction on the sequence, or
     # If sequences are not already aligned and reverse primer not found in most likely sequence (the one without the forward primer), try matching on the same sequence than the first match (primer in the other direction)
     if not not_aligned or (not_aligned and (reversematch is None or reversematch[1] is None)):
-        if not not_aligned:
-            sequence_to_match = second_matched_seq
+        if not_aligned and first_match_first_seq:
+            seq_to_match = sequences[1]
         else:
-            sequence_to_match = first_matched_seq
+            seq_to_match = sequences[0]
         reversematch = []
         # Compute begin
         begin=directmatch[1][2]+1  # end of match + 1 on the same sequence
@@ -365,7 +400,7 @@ cdef tuple annotate(sequences, infos, verbose=False):
                 primer=p
             # Saving original primer as 4th member of the tuple to serve as correct key in infos dict even if it might have been reversed complemented
             # (3rd member already used by directmatch)
-            reversematch.append((primer, primer(sequence_to_match, same_sequence=not new_seq, pattern=pattern, begin=begin), None, p))
+            reversematch.append((primer, primer(seq_to_match, same_sequence=not new_seq, pattern=pattern, begin=begin), None, p))
             new_seq = False
             pattern+=1
         # Choose match closer to the end of the sequence
@@ -378,11 +413,11 @@ cdef tuple annotate(sequences, infos, verbose=False):
             message = b'No reverse primer match'
         else:
             message = b'No direct primer match'
-        final_sequence[b'error']=message
-        return False, final_sequence
+        sequences[0][b'error']=message
+        return False, sequences[0]
     
     if reversematch is None:
-        final_sequence[b'status']=b'partial'
+        sequences[0][b'status']=b'partial'
         
         if directmatch[0].forward:
             tags=(directmatch[1][3],None)
@@ -392,78 +427,84 @@ cdef tuple annotate(sequences, infos, verbose=False):
         samples = infos[None]
         
     else:
-        final_sequence[b'status']=b'full'
+        sequences[0][b'status']=b'full'
+        
+        if not not_aligned or first_match_first_seq:
+            match = sequences[0][reversematch[1][1]:reversematch[1][2]]
+        else:
+            match = sequences[1][reversematch[1][1]:reversematch[1][2]]
         
-        match = second_matched_seq[reversematch[1][1]:reversematch[1][2]]
         match = match.reverse_complement
         
-        if not not_aligned or id(second_matched_seq) == id(sequenceF):
-            final_sequence = final_sequence[0:reversematch[1][1]]
-        else:
-            cut_seq = sequenceR[reversematch[1][2]:]
+        if not not_aligned:
+            sequences[0] = sequences[0][0:reversematch[1][1]]
+        elif first_match_first_seq:
+            sequences[1] = sequences[1][reversematch[1][2]:]
             if not directmatch[0].forward:
-                cut_seq = cut_seq.reverse_complement
-            final_sequence[REVERSE_SEQ_COLUMN_NAME] = cut_seq.seq           # used by alignpairedend tool
-            final_sequence[REVERSE_QUALITY_COLUMN_NAME] = cut_seq.quality   # used by alignpairedend tool
-        
+                sequences[1] = sequences[1].reverse_complement
+            sequences[0][REVERSE_SEQ_COLUMN_NAME] = sequences[1].seq           # used by alignpairedend tool
+            sequences[0][REVERSE_QUALITY_COLUMN_NAME] = sequences[1].quality   # used by alignpairedend tool
+        else:
+            sequences[0] = sequences[0][reversematch[1][2]:]
+            
         if directmatch[0].forward:
             tags=(directmatch[1][3], reversematch[1][3])
-            final_sequence[b'reverse_errors'] = reversematch[1][0]
-            final_sequence[b'reverse_primer'] = reversematch[0].raw
-            final_sequence[b'reverse_match'] = match.seq
+            sequences[0][b'reverse_errors'] = reversematch[1][0]
+            sequences[0][b'reverse_primer'] = reversematch[0].raw
+            sequences[0][b'reverse_match'] = match.seq
         
         else:
             tags=(reversematch[1][3], directmatch[1][3])
-            final_sequence[b'forward_errors'] = reversematch[1][0]
-            final_sequence[b'forward_primer'] = reversematch[0].raw
-            final_sequence[b'forward_match'] = match.seq
+            sequences[0][b'forward_errors'] = reversematch[1][0]
+            sequences[0][b'forward_primer'] = reversematch[0].raw
+            sequences[0][b'forward_match'] = match.seq
                 
         if tags[0] is not None:
-            final_sequence[b'forward_tag'] = tags[0]
+            sequences[0][b'forward_tag'] = tags[0]
         if tags[1] is not None:
-            final_sequence[b'reverse_tag'] = tags[1]
+            sequences[0][b'reverse_tag'] = tags[1]
 
         samples = infos[reversematch[3]]
         
     if not directmatch[0].forward:
-        final_sequence = final_sequence.reverse_complement
-        final_sequence[b'reversed'] = True   # used by the alignpairedend tool (in kmer_similarity.c)
+        sequences[0] = sequences[0].reverse_complement
+        sequences[0][b'reversed'] = True   # used by the alignpairedend tool (in kmer_similarity.c)
 
     sample=None
-
-    if tags[0] is not None:                                    # Direct  tag known
-        if tags[1] is not None:                                # Reverse tag known
-            sample = samples.get(tags, None)             
-        else:                                                   # Only direct tag known
-            s=[samples[x] for x in samples if x[0]==tags[0]]
-            if len(s)==1:
-                sample=s[0]
-            elif len(s)>1:
-                final_sequence[b'error']=b'multiple samples match tags'
-                return False, final_sequence
-            else:
-                sample=None
-    else: 
-        if tags[1] is not None:                                 # Only reverse tag known
-            s=[samples[x] for x in samples if x[1]==tags[1]]
-            if len(s)==1:
-                sample=s[0]
-            elif len(s)>1:
-                final_sequence[b'error']=b'multiple samples match tags'
-                return False, final_sequence
-            else:
-                sample=None
+    if not no_tags:
+        if tags[0] is not None:                                    # Direct  tag known
+            if tags[1] is not None:                                # Reverse tag known
+                sample = samples.get(tags, None)             
+            else:                                                   # Only direct tag known
+                s=[samples[x] for x in samples if x[0]==tags[0]]
+                if len(s)==1:
+                    sample=s[0]
+                elif len(s)>1:
+                    sequences[0][b'error']=b'Did not found reverse tag'
+                    return False, sequences[0]
+                else:
+                    sample=None
+        else: 
+            if tags[1] is not None:                                 # Only reverse tag known
+                s=[samples[x] for x in samples if x[1]==tags[1]]
+                if len(s)==1:
+                    sample=s[0]
+                elif len(s)>1:
+                    sequences[0][b'error']=b'Did not found forward tag'
+                    return False, sequences[0]
+                else:
+                    sample=None
+        
+        if sample is None:
+            sequences[0][b'error']=b"No tags found"
+            return False, sequences[0]
     
-    if sample is None:
-        final_sequence[b'error']=b"Cannot assign sequence to a sample"
-        return False, final_sequence
-    
-    final_sequence.update(sample)
+        sequences[0].update(sample)
     
     if not not_aligned:
-        final_sequence[b'seq_length']=len(final_sequence)
+        sequences[0][b'seq_length']=len(sequences[0])
     
-    return True, final_sequence
+    return True, sequences[0]
 
 
 def run(config):
@@ -567,11 +608,13 @@ def run(config):
         Column.new_column(o_view, REVERSE_SEQ_COLUMN_NAME, OBI_SEQ)
         Column.new_column(o_view, REVERSE_QUALITY_COLUMN_NAME, OBI_QUAL, associated_column_name=REVERSE_SEQ_COLUMN_NAME, associated_column_version=o_view[REVERSE_SEQ_COLUMN_NAME].version)
         
-        Column.new_column(unidentified, REVERSE_SEQ_COLUMN_NAME, OBI_SEQ)
-        Column.new_column(unidentified, REVERSE_QUALITY_COLUMN_NAME, OBI_QUAL, associated_column_name=REVERSE_SEQ_COLUMN_NAME, associated_column_version=unidentified[REVERSE_SEQ_COLUMN_NAME].version)
+        if unidentified is not None:
+            Column.new_column(unidentified, REVERSE_SEQ_COLUMN_NAME, OBI_SEQ)
+            Column.new_column(unidentified, REVERSE_QUALITY_COLUMN_NAME, OBI_QUAL, associated_column_name=REVERSE_SEQ_COLUMN_NAME, associated_column_version=unidentified[REVERSE_SEQ_COLUMN_NAME].version)
     
     g = 0
     u = 0
+    no_tags = config['ngsfilter']['notags']
     try:
         for i in range(entries_len):
             PyErr_CheckSignals()
@@ -580,7 +623,7 @@ def run(config):
                 modseq = [Nuc_Seq.new_from_stored(forward[i]), Nuc_Seq.new_from_stored(reverse[i])]
             else:
                 modseq = [Nuc_Seq.new_from_stored(entries[i])]
-            good, oseq = annotate(modseq, infos)
+            good, oseq = annotate(modseq, infos, no_tags)
             if good:
                 o_view[g].set(oseq.id, oseq.seq, definition=oseq.definition, quality=oseq.quality, tags=oseq)
                 g+=1
@@ -588,7 +631,10 @@ def run(config):
                 unidentified[u].set(oseq.id, oseq.seq, definition=oseq.definition, quality=oseq.quality, tags=oseq)
                 u+=1
     except Exception, e:
-        raise RollbackException("obi ngsfilter error, rollbacking views: "+str(e), o_view, unidentified)
+        if unidentified is not None:
+            raise RollbackException("obi ngsfilter error, rollbacking views: "+str(e), o_view, unidentified)
+        else:
+            raise RollbackException("obi ngsfilter error, rollbacking view: "+str(e), o_view)
     
     pb(i, force=True)
     print("", file=sys.stderr)
@@ -596,9 +642,10 @@ def run(config):
     # Save command config in View and DMS comments
     command_line = " ".join(sys.argv[1:])
     o_view.write_config(config, "ngsfilter", command_line, input_dms_name=input_dms_name, input_view_name=input_view_name)
-    unidentified.write_config(config, "ngsfilter", command_line, input_dms_name=input_dms_name, input_view_name=input_view_name)
-    # Add comment about unidentified seqs
-    unidentified.comments["info"] = "View containing sequences categorized as unidentified by the ngsfilter command"
+    if unidentified is not None:
+        unidentified.write_config(config, "ngsfilter", command_line, input_dms_name=input_dms_name, input_view_name=input_view_name)
+        # Add comment about unidentified seqs
+        unidentified.comments["info"] = "View containing sequences categorized as unidentified by the ngsfilter command"
     output[0].record_command_line(command_line)
     
     #print("\n\nOutput view:\n````````````", file=sys.stderr)
@@ -607,7 +654,8 @@ def run(config):
     input[0].close()
     output[0].close()
     info_input[0].close()
-    unidentified_input[0].close()
+    if unidentified is not None:
+        unidentified_input[0].close()
     aligner.free()
 
     logger("info", "Done.")

python/obitools3/commands/uniq.pyx

@@ -23,7 +23,9 @@ from cpython.exc cimport PyErr_CheckSignals
 
 
 __title__="Group sequence records together"
- 
+
+
+REVERSE_QUALITY_COLUMN_NAME = b"REVERSE_QUALITY"   # TODO from ngsfilter, move to C
 
  
 def addOptions(parser):
@@ -491,9 +493,11 @@ cdef uniq_sequences(View_NUC_SEQS view, View_NUC_SEQS o_view, ProgressBar pb, li
 
         o_idx += 1
     
-    # Deletes quality column if there is one because the matching between sequence and quality will be broken (quality set to NA when sequence not)
+    # Deletes quality columns if there is one because the matching between sequence and quality will be broken (quality set to NA when sequence not)
     if QUALITY_COLUMN in view:
         o_view.delete_column(QUALITY_COLUMN)
+    if REVERSE_QUALITY_COLUMN_NAME in view:
+        o_view.delete_column(REVERSE_QUALITY_COLUMN_NAME)
     
     if taxonomy is not None:
         print("")  # TODO because in the middle of progress bar. Better solution?

python/obitools3/parsers/ngsfilter.pyx

@@ -57,6 +57,9 @@ def ngsfilterIterator(lineiterator,
         split_line = line.split()
         tags = split_line.pop(2)
         tags = tags.split(b":")
+        for t_idx in range(2):
+            if tags[t_idx]==b"-" or tags[t_idx]==b"None" or tags[t_idx]==b"":
+                tags[t_idx] = nastring
         if len(tags) == 1:          # Forward and reverse tags are the same
             tags.append(tags[0])
         split_line.insert(2, tags[0])

python/obitools3/version.py

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

setup.py

@@ -16,6 +16,8 @@ from distutils.extension import Extension
 
 from distutils.dist import Distribution as ori_Distribution
 
+from python.obitools3.version import version
+
 
 class Distribution(ori_Distribution):
     
@@ -83,7 +85,7 @@ def findPackage(root,base=None):
   
 
 PACKAGE     = "OBITools3"
-VERSION     = "3.0.0-beta1"
+VERSION     = version
 AUTHOR      = 'Celine Mercier'
 EMAIL       = 'celine.mercier@metabarcoding.org'
 URL         = "http://metabarcoding.org/obitools3"