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

ngsfilter: made more robust and practical to use with empty tags

Browse Source
This commit is contained in:
Celine Mercier
2019-11-29 15:21:08 +01:00
parent 8bc249b2f4
commit 53f18316b0
3 changed files with 52 additions and 40 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
python/obitools3/apps/optiongroups/__init__.py
View File

@ -55,7 +55,7 @@ def __addImportInputOption(optionManager):
action="store_const", dest="obi:inputformat", action="store_const", dest="obi:inputformat",
default=None, default=None,
const=b'ngsfilter', 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', group.add_argument('--ecopcr-result-input',
action="store_const", dest="obi:inputformat", action="store_const", dest="obi:inputformat",

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

@ -56,6 +56,11 @@ def addOptions(parser):
type=str, type=str,
default=None, default=None,
help="URI to the view used to store the sequences unassigned to any sample") 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', group.add_argument('-e','--error',
action="store", dest="ngsfilter: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 i=0
for p in info_view: for p in info_view:
forward=Primer(p[b'forward_primer'], 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, True,
max_errors=max_errors, max_errors=max_errors,
verbose=verbose, verbose=verbose,
@ -178,7 +183,7 @@ cdef read_info_view(info_view, max_errors=2, verbose=False, not_aligned=False):
infos[forward]=fp infos[forward]=fp
reverse=Primer(p[b'reverse_primer'], 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, False,
max_errors=max_errors, max_errors=max_errors,
verbose=verbose, verbose=verbose,
@ -213,10 +218,11 @@ cdef read_info_view(info_view, max_errors=2, verbose=False, not_aligned=False):
rpp=rp.get(cf,{}) rpp=rp.get(cf,{})
rp[cf]=rpp rp[cf]=rpp
tags = (p[b'forward_tag'] if p[b'forward_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 p[b'reverse_tag']!=b'-' 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) "Tag pair %s is already used with primer pairs: (%s,%s)" % (str(tags),forward,reverse)
# Save additional data # 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 return infos, primer_list
cdef tuple annotate(sequences, infos, verbose=False): cdef tuple annotate(sequences, infos, no_tags, verbose=False):
def sortMatch(match): def sortMatch(match):
if match[1] is None: if match[1] is None:
@ -430,35 +436,35 @@ cdef tuple annotate(sequences, infos, verbose=False):
final_sequence[b'reversed'] = True # used by the alignpairedend tool (in kmer_similarity.c) final_sequence[b'reversed'] = True # used by the alignpairedend tool (in kmer_similarity.c)
sample=None sample=None
if not no_tags:
if tags[0] is not None: # Direct tag known if tags[0] is not None: # Direct tag known
if tags[1] is not None: # Reverse tag known if tags[1] is not None: # Reverse tag known
sample = samples.get(tags, None) sample = samples.get(tags, None)
else: # Only direct tag known else: # Only direct tag known
s=[samples[x] for x in samples if x[0]==tags[0]] s=[samples[x] for x in samples if x[0]==tags[0]]
if len(s)==1: if len(s)==1:
sample=s[0] sample=s[0]
elif len(s)>1: elif len(s)>1:
final_sequence[b'error']=b'multiple samples match tags' final_sequence[b'error']=b'multiple samples match tags'
return False, final_sequence return False, final_sequence
else: else:
sample=None sample=None
else: else:
if tags[1] is not None: # Only reverse tag known if tags[1] is not None: # Only reverse tag known
s=[samples[x] for x in samples if x[1]==tags[1]] s=[samples[x] for x in samples if x[1]==tags[1]]
if len(s)==1: if len(s)==1:
sample=s[0] sample=s[0]
elif len(s)>1: elif len(s)>1:
final_sequence[b'error']=b'multiple samples match tags' final_sequence[b'error']=b'multiple samples match tags'
return False, final_sequence return False, final_sequence
else: else:
sample=None sample=None
if sample is None:
final_sequence[b'error']=b"Cannot assign sequence to a sample"
return False, final_sequence
if sample is None: final_sequence.update(sample)
final_sequence[b'error']=b"Cannot assign sequence to a sample"
return False, final_sequence
final_sequence.update(sample)
if not not_aligned: if not not_aligned:
final_sequence[b'seq_length']=len(final_sequence) final_sequence[b'seq_length']=len(final_sequence)
@ -572,6 +578,7 @@ def run(config):
g = 0 g = 0
u = 0 u = 0
no_tags = config['ngsfilter']['notags']
try: try:
for i in range(entries_len): for i in range(entries_len):
PyErr_CheckSignals() PyErr_CheckSignals()
@ -580,7 +587,7 @@ def run(config):
modseq = [Nuc_Seq.new_from_stored(forward[i]), Nuc_Seq.new_from_stored(reverse[i])] modseq = [Nuc_Seq.new_from_stored(forward[i]), Nuc_Seq.new_from_stored(reverse[i])]
else: else:
modseq = [Nuc_Seq.new_from_stored(entries[i])] modseq = [Nuc_Seq.new_from_stored(entries[i])]
good, oseq = annotate(modseq, infos) good, oseq = annotate(modseq, infos, no_tags)
if good: if good:
o_view[g].set(oseq.id, oseq.seq, definition=oseq.definition, quality=oseq.quality, tags=oseq) o_view[g].set(oseq.id, oseq.seq, definition=oseq.definition, quality=oseq.quality, tags=oseq)
g+=1 g+=1
@ -596,9 +603,10 @@ def run(config):
# Save command config in View and DMS comments # Save command config in View and DMS comments
command_line = " ".join(sys.argv[1:]) 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) 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) if unidentified is not None:
# Add comment about unidentified seqs unidentified.write_config(config, "ngsfilter", command_line, input_dms_name=input_dms_name, input_view_name=input_view_name)
unidentified.comments["info"] = "View containing sequences categorized as unidentified by the ngsfilter command" # 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) output[0].record_command_line(command_line)
#print("\n\nOutput view:\n````````````", file=sys.stderr) #print("\n\nOutput view:\n````````````", file=sys.stderr)
@ -607,7 +615,8 @@ def run(config):
input[0].close() input[0].close()
output[0].close() output[0].close()
info_input[0].close() info_input[0].close()
unidentified_input[0].close() if unidentified is not None:
unidentified_input[0].close()
aligner.free() aligner.free()
logger("info", "Done.") logger("info", "Done.")

3
python/obitools3/parsers/ngsfilter.pyx Executable file → Normal file
View File

@ -57,6 +57,9 @@ def ngsfilterIterator(lineiterator,
split_line = line.split() split_line = line.split()
tags = split_line.pop(2) tags = split_line.pop(2)
tags = tags.split(b":") 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 if len(tags) == 1: # Forward and reverse tags are the same
tags.append(tags[0]) tags.append(tags[0])
split_line.insert(2, tags[0]) split_line.insert(2, tags[0])