From ceaafca4271b62ad31b3084474cb3897afcc931a Mon Sep 17 00:00:00 2001
From: Celine Mercier <celine.mercier@metabarcoding.org>
Date: Fri, 29 Mar 2019 10:56:17 +0100
Subject: [PATCH] ngsfilter: fixed a bug (maybe 2) in the algo for the choice
 of the reverse primer when running on unaligned sequences

---
 python/obitools3/commands/ngsfilter.pyx | 16 ++++++++++------
 1 file changed, 10 insertions(+), 6 deletions(-)

diff --git a/python/obitools3/commands/ngsfilter.pyx b/python/obitools3/commands/ngsfilter.pyx
index 9375d9a..073855d 100755
--- a/python/obitools3/commands/ngsfilter.pyx
+++ b/python/obitools3/commands/ngsfilter.pyx
@@ -283,7 +283,8 @@ cdef tuple annotate(sequences, infos, verbose=False):
             if pattern == MAX_PATTERN:
                 new_seq = True
                 pattern = 0
-            directmatch.append((p, p(seq, same_sequence=not new_seq, pattern=pattern), seq))
+            # Saving original primer as 4th member of the tuple to serve as correct key in infos dict even if it might be reversed complemented (not here)
+            directmatch.append((p, p(seq, same_sequence=not new_seq, pattern=pattern), seq, p))
             new_seq = False
             pattern+=1
     
@@ -329,10 +330,11 @@ cdef tuple annotate(sequences, infos, verbose=False):
     # Keep only paired reverse primer
     infos = infos[directmatch[0]]    
     
-    # If not aligned, look for other match in already computed match (choose the one that makes the biggest amplicon)
+    # If not aligned, look for other match in already computed matches (choose the one that makes the biggest amplicon)
     if not_aligned:
         i=1
-        while all_direct_matches[i][1] is None and all_direct_matches[i][0].forward and i<len(all_direct_matches):
+        # 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]):
             i+=1
         if i < len(all_direct_matches):
             reversematch = all_direct_matches[i]
@@ -341,7 +343,7 @@ cdef tuple annotate(sequences, infos, verbose=False):
     
     # 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[1] is None):
+    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
         else:
@@ -360,7 +362,9 @@ cdef tuple annotate(sequences, infos, verbose=False):
                 primer=p.revcomp
             else:
                 primer=p
-            reversematch.append((primer, primer(sequence_to_match, same_sequence=not new_seq, pattern=pattern, begin=begin)))
+            # 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))
             new_seq = False
             pattern+=1
         # Choose match closer to the end of the sequence
@@ -416,7 +420,7 @@ cdef tuple annotate(sequences, infos, verbose=False):
         if tags[1] is not None:
             final_sequence[b'reverse_tag'] = tags[1]
 
-        samples = infos[reversematch[0]]
+        samples = infos[reversematch[3]]
         
     if not directmatch[0].forward and not not_aligned:   # don't reverse complement if not_aligned
         final_sequence = final_sequence.reverse_complement