obitools/obitools3
2
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 66 Wiki Activity

ngsfilter: fixed a bug (maybe 2) in the algo for the choice of the

Browse Source 
reverse primer when running on unaligned sequences
This commit is contained in:
Celine Mercier
2019-03-29 10:56:17 +01:00
parent 8e70bf1ee1
commit ceaafca427
1 changed files with 10 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
python/obitools3/commands/ngsfilter.pyx
View File

@ -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