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

ngsfilter: fixed sequence cutting when dealing with unaligned sequences.

Browse Source 
Could use optimization
This commit is contained in:
Celine Mercier
2019-12-12 17:01:31 +01:00
parent 81a179239c
commit 6bfd7441f3
1 changed files with 106 additions and 73 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@ -255,17 +255,12 @@ cdef tuple annotate(sequences, infos, no_tags, 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"):
@ -281,8 +276,6 @@ cdef tuple annotate(sequences, infos, no_tags, verbose=False):
# Try direct matching:
directmatch = []
first_matched_seq = None
second_matched_seq = None
for seq in sequences:
new_seq = True
pattern = 0
@ -301,43 +294,46 @@ cdef tuple annotate(sequences, infos, no_tags, 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]]
rev_prim = list(infos.keys())[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
@ -346,20 +342,48 @@ cdef tuple annotate(sequences, infos, no_tags, verbose=False):
(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 \
rev_prim != all_direct_matches[i][0]) :
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
@ -376,7 +400,7 @@ cdef tuple annotate(sequences, infos, no_tags, 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
@ -389,11 +413,11 @@ cdef tuple annotate(sequences, infos, no_tags, 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)
@ -403,42 +427,48 @@ cdef tuple annotate(sequences, infos, no_tags, 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] - directmatch[1][2])]
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 not no_tags:
@ -450,8 +480,8 @@ cdef tuple annotate(sequences, infos, no_tags, verbose=False):
if len(s)==1:
sample=s[0]
elif len(s)>1:
final_sequence[b'error']=b'Did not found reverse tag'
return False, final_sequence
sequences[0][b'error']=b'Did not found reverse tag'
return False, sequences[0]
else:
sample=None
else:
@ -460,21 +490,21 @@ cdef tuple annotate(sequences, infos, no_tags, verbose=False):
if len(s)==1:
sample=s[0]
elif len(s)>1:
final_sequence[b'error']=b'Did not found forward tag'
return False, final_sequence
sequences[0][b'error']=b'Did not found forward tag'
return False, sequences[0]
else:
sample=None
if sample is None:
final_sequence[b'error']=b"No tags found"
return False, final_sequence
sequences[0][b'error']=b"No tags found"
return False, sequences[0]
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):
@ -601,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)