ecopcr: fixed and improved the options to keep nuclotides around the

amplicon

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/ecopcr.pyx

@@ -107,14 +107,20 @@ def addOptions(parser):
                        help="Defines the method used for estimating the Tm (melting temperature) between the primers and their corresponding "
                             "target sequences. SANTALUCIA: 1, or OWCZARZY: 2. Default: 1.")
 
+    group.add_argument('--keep-primers', '-p',
+                       action="store_true", 
+                       dest="ecopcr:keep-primers",
+                       default=False,
+                       help="Whether to keep the primers attached to the output sequences (default: the primers are cut out).")
+
     group.add_argument('--keep-nucs', '-D',
                        action="store", 
                        dest="ecopcr:keep-nucs",
-                       metavar="<INTEGER>",
+                       metavar="<N>",
                        type=int,
                        default=0,
-                       help="Keeps the specified number of nucleotides on each side of the in silico amplified sequences, "
-                            "(already including the amplified DNA fragment plus the two target sequences of the primers).")
+                       help="Keeps N nucleotides on each side of the in silico amplified sequences, "
+                            "not including the primers (implying that primers are automatically kept if N > 0).")
 
     group.add_argument('--kingdom-mode', '-k',
                        action="store_true", 
@@ -185,7 +191,7 @@ def run(config):
                   config['ecopcr']['min-length'], config['ecopcr']['max-length'], \
                   restrict_to_taxids_p, ignore_taxids_p, \
                   config['ecopcr']['circular'], config['ecopcr']['salt-concentration'], config['ecopcr']['salt-correction-method'], \
-                  config['ecopcr']['keep-nucs'], config['ecopcr']['kingdom-mode']) < 0:
+                  config['ecopcr']['keep-nucs'], config['ecopcr']['keep-primers'], config['ecopcr']['kingdom-mode']) < 0:
         raise Exception("Error running ecopcr")
 
     # Save command config in DMS comments

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/dms/capi/obiecopcr.pxd

@@ -23,6 +23,7 @@ cdef extern from "obi_ecopcr.h" nogil:
                    double salt_concentration,
                    int salt_correction_method,
                    int keep_nucleotides,
+                   bint keep_primers,
                    bint kingdom_mode)
 
     

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-beta4'
 
 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"

src/obi_ecopcr.c

@@ -77,7 +77,8 @@ static int create_output_columns(Obiview_p o_view, bool kingdom_mode);
  * @param err2 The number of errors in the second primer.
  * @param strand The DNA strand direction of the amplicon (R(everse) or D(irect)).
  * @param kingdom_mode Whether the kingdom or the superkingdom informations should be printed to the output.
- * @param keep_nucleotides Number of nucleotides kept on each side of the amplicon.
+ * @param keep_nucleotides Number of nucleotides kept on each side of the amplicon (not including the primers if they are kept).
+ * @param keep_primers Whether to keep the primers.
  * @param i_id_column A pointer on the input sequence identifier column.
  * @param o_id_column A pointer on the output sequence identifier column.
  * @param o_ori_seq_len_column A pointer on the original sequence length column.
@@ -124,6 +125,7 @@ static int print_seq(Obiview_p i_view, Obiview_p o_view,
 					 int32_t err1, int32_t err2,
 					 char strand, bool kingdom_mode,
 					 int keep_nucleotides,
+					 bool keep_primers,
 					 OBIDMS_column_p i_id_column, OBIDMS_column_p o_id_column, OBIDMS_column_p o_ori_seq_len_column,
 					 OBIDMS_column_p o_amplicon_column, OBIDMS_column_p o_amplicon_length_column,
 					 OBIDMS_column_p o_taxid_column, OBIDMS_column_p o_rank_column, OBIDMS_column_p o_name_column,
@@ -328,6 +330,7 @@ static int print_seq(Obiview_p i_view, Obiview_p o_view,
 					 int32_t err1, int32_t err2,
 					 char strand, bool kingdom_mode,
 					 int keep_nucleotides,
+					 bool keep_primers,
 					 OBIDMS_column_p i_id_column, OBIDMS_column_p o_id_column, OBIDMS_column_p o_ori_seq_len_column,
 					 OBIDMS_column_p o_amplicon_column, OBIDMS_column_p o_amplicon_length_column,
 					 OBIDMS_column_p o_taxid_column, OBIDMS_column_p o_rank_column, OBIDMS_column_p o_name_column,
@@ -382,7 +385,7 @@ static int print_seq(Obiview_p i_view, Obiview_p o_view,
 		oligo2[o1->patlen] = 0;
 		error2 = err1;
 
-		if (keep_nucleotides == 0)
+		if (!keep_primers)
 			amplicon+=o2->patlen;
 		else
 		{
@@ -401,7 +404,7 @@ static int print_seq(Obiview_p i_view, Obiview_p o_view,
 		oligo2[o2->patlen] = 0;
 		error2 = err2;
 
-		if (keep_nucleotides==0)
+		if (!keep_primers)
 			amplicon+=o1->patlen;
 		else
 		{
@@ -411,16 +414,11 @@ static int print_seq(Obiview_p i_view, Obiview_p o_view,
 	}
 
 	ecoComplementSequence(oligo2);
-	if (keep_nucleotides == 0)
+	if (!keep_primers)
 		amplicon[amplicon_len]=0;
 	else
 	{
 		amplicon_len = ldelta+rdelta+amplicon_len;
-		for (i=0; i<ldelta; i++)
-			amplicon[i]|=32;
-		for (i=1; i<=rdelta; i++)
-			amplicon[amplicon_len-i]|=32;
-
 		amplicon[amplicon_len] = 0;
 	}
 
@@ -659,6 +657,7 @@ int obi_ecopcr(const char* i_dms_name,
 			   double salt,
 			   int saltmethod,
 			   int keep_nucleotides,
+			   bool keep_primers,
 			   bool kingdom_mode)
 {
 
@@ -717,6 +716,9 @@ int obi_ecopcr(const char* i_dms_name,
 
 	signal(SIGINT, sig_handler);
 
+	if (keep_nucleotides > 0)
+		keep_primers = true;
+
 	if (circular)
 	{
 		circular = strlen(primer1);
@@ -1076,6 +1078,7 @@ int obi_ecopcr(const char* i_dms_name,
 														  	  erri, errj,
 															  'D', kingdom_mode,
 															  keep_nucleotides,
+															  keep_primers,
 															  i_id_column, o_id_column, o_ori_seq_len_column,
 															  o_amplicon_column, o_amplicon_length_column,
 															  o_taxid_column, o_rank_column, o_name_column,
@@ -1163,6 +1166,7 @@ int obi_ecopcr(const char* i_dms_name,
 														  	  erri, errj,
 															  'R', kingdom_mode,
 															  keep_nucleotides,
+															  keep_primers,
 															  i_id_column, o_id_column, o_ori_seq_len_column,
 															  o_amplicon_column, o_amplicon_length_column,
 															  o_taxid_column, o_rank_column, o_name_column,

src/obi_ecopcr.h

@@ -93,8 +93,8 @@
  * @param salt_concentration The salt concentration used for estimating the Tm.
  * @param salt_correction_method The method used for estimating the Tm (melting temperature) between the primers and their corresponding
  *                              target sequences. SANTALUCIA: 1, or OWCZARZY: 2.
- * @param keep_nucleotides The number of nucleotides to keep on each side of the in silico amplified sequences
- *                         (already including the amplified DNA fragment plus the two target sequences of the primers).
+ * @param keep_nucleotides The number of nucleotides to keep on each side of the in silico amplified sequences, not including primers (primers automatically entirely kept if > 0).
+ * @param keep_primers Whether primers are kept attached to the output sequences.
  * @param kingdom_mode Whether the kingdom or the superkingdom informations should be printed to the output.
  *
  * @returns A value indicating the success of the operation.
@@ -121,6 +121,7 @@ int obi_ecopcr(const char* i_dms_name,
 			   double salt_concentration,
 			   int salt_correction_method,
 			   int keep_nucleotides,
+			   bool keep_primers,
 			   bool kingdom_mode);
 
 #endif /* OBI_ECOPCR_H_ */