Compare commits

...

10 Commits

7 changed files with 99 additions and 82 deletions

View File

@ -205,19 +205,25 @@ def run(config):
if type(entries) == list: if type(entries) == list:
forward = entries[0] forward = entries[0]
reverse = entries[1] reverse = entries[1]
aligner = Kmer_similarity(forward, \ if len(forward) == 0 or len(reverse) == 0:
view2=reverse, \ aligner = None
kmer_size=config['alignpairedend']['kmersize'], \ else:
reversed_column=None) aligner = Kmer_similarity(forward, \
view2=reverse, \
kmer_size=config['alignpairedend']['kmersize'], \
reversed_column=None)
else: else:
aligner = Kmer_similarity(entries, \ if len(entries) == 0:
column2=entries[REVERSE_SEQUENCE_COLUMN], \ aligner = None
qual_column2=entries[REVERSE_QUALITY_COLUMN], \ else:
kmer_size=config['alignpairedend']['kmersize'], \ aligner = Kmer_similarity(entries, \
reversed_column=entries[b'reversed']) # column created by the ngsfilter tool column2=entries[REVERSE_SEQUENCE_COLUMN], \
qual_column2=entries[REVERSE_QUALITY_COLUMN], \
kmer_size=config['alignpairedend']['kmersize'], \
reversed_column=entries[b'reversed']) # column created by the ngsfilter tool
ba = alignmentIterator(entries, aligner) ba = alignmentIterator(entries, aligner)
i = 0 i = 0
for ali in ba: for ali in ba:
@ -251,7 +257,7 @@ def run(config):
pb(i, force=True) pb(i, force=True)
print("", file=sys.stderr) print("", file=sys.stderr)
if kmer_ali : if kmer_ali and aligner is not None:
aligner.free() aligner.free()
# Save command config in View and DMS comments # Save command config in View and DMS comments

View File

@ -322,7 +322,7 @@ cdef tuple annotate(sequences, infos, no_tags, verbose=False):
sequences[0] = sequences[0][directmatch[1][2]:] sequences[0] = sequences[0][directmatch[1][2]:]
else: else:
sequences[1] = sequences[1][directmatch[1][2]:] sequences[1] = sequences[1][directmatch[1][2]:]
sequences[0][REVERSE_SEQUENCE_COLUMN] = sequences[1].seq # used by alignpairedend tool sequences[0][REVERSE_SEQUENCE_COLUMN] = sequences[1].seq # used by alignpairedend tool
sequences[0][REVERSE_QUALITY_COLUMN] = sequences[1].quality # used by alignpairedend tool sequences[0][REVERSE_QUALITY_COLUMN] = sequences[1].quality # used by alignpairedend tool
if directmatch[0].forward: if directmatch[0].forward:
@ -369,7 +369,7 @@ cdef tuple annotate(sequences, infos, no_tags, verbose=False):
sequences[0] = sequences[0][:r[1]] sequences[0] = sequences[0][:r[1]]
else: else:
sequences[1] = sequences[1][:r[1]] sequences[1] = sequences[1][:r[1]]
sequences[0][REVERSE_SEQUENCE_COLUMN] = sequences[1].seq # used by alignpairedend tool sequences[0][REVERSE_SEQUENCE_COLUMN] = sequences[1].seq # used by alignpairedend tool
sequences[0][REVERSE_QUALITY_COLUMN] = sequences[1].quality # used by alignpairedend tool sequences[0][REVERSE_QUALITY_COLUMN] = sequences[1].quality # used by alignpairedend tool
# do the same on the other seq # do the same on the other seq
if first_match_first_seq: if first_match_first_seq:
@ -394,7 +394,7 @@ cdef tuple annotate(sequences, infos, no_tags, verbose=False):
seq_to_match = sequences[0] seq_to_match = sequences[0]
reversematch = [] reversematch = []
# Compute begin # Compute begin
begin=directmatch[1][2]+1 # end of match + 1 on the same sequence #begin=directmatch[1][2]+1 # end of match + 1 on the same sequence -- No, already cut out forward primer
# Try reverse matching on the other sequence: # Try reverse matching on the other sequence:
new_seq = True new_seq = True
pattern = 0 pattern = 0
@ -408,7 +408,7 @@ cdef tuple annotate(sequences, infos, no_tags, verbose=False):
primer=p 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 # 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) # (3rd member already used by directmatch)
reversematch.append((primer, primer(seq_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=0), None, p))
new_seq = False new_seq = False
pattern+=1 pattern+=1
# Choose match closer to the end of the sequence # Choose match closer to the end of the sequence
@ -645,6 +645,7 @@ def run(config):
g = 0 g = 0
u = 0 u = 0
i = 0
no_tags = config['ngsfilter']['notags'] no_tags = config['ngsfilter']['notags']
try: try:
for i in range(entries_len): for i in range(entries_len):

View File

@ -354,8 +354,8 @@ cdef uniq_sequences(View_NUC_SEQS view, View_NUC_SEQS o_view, ProgressBar pb, di
key = mergedKeys[k] key = mergedKeys[k]
merged_col_name = mergedKeys_m[k] merged_col_name = mergedKeys_m[k]
if merged_infos[merged_col_name]['nb_elts'] == 1: # if merged_infos[merged_col_name]['nb_elts'] == 1:
raise Exception("Can't merge information from a tag with only one element (e.g. one sample ; don't use -m option)") # raise Exception("Can't merge information from a tag with only one element (e.g. one sample ; don't use -m option)")
if merged_col_name in view: if merged_col_name in view:
i_col = view[merged_col_name] i_col = view[merged_col_name]

View File

@ -8,7 +8,7 @@ Created on feb 20th 2018
import types import types
from obitools3.utils cimport __etag__ from obitools3.utils cimport __etag__
from obitools3.utils cimport str2bytes
def tabIterator(lineiterator, def tabIterator(lineiterator,
bint header = False, bint header = False,
@ -75,7 +75,7 @@ def tabIterator(lineiterator,
continue continue
else: else:
# TODO ??? default column names? like R? # TODO ??? default column names? like R?
keys = [i for i in range(len(line.split(sep)))] keys = [str2bytes(str(i)) for i in range(len(line.split(sep)))]
while skipped < skip : while skipped < skip :
line = next(iterator) line = next(iterator)

View File

@ -53,7 +53,11 @@ def entryIteratorFactory(lineiterator,
i = iterator i = iterator
first=next(i) try:
first=next(i)
except StopIteration:
first=""
pass
format=b"tabular" format=b"tabular"

View File

@ -1,5 +1,5 @@
major = 3 major = 3
minor = 0 minor = 0
serial= '0b43' serial= '1b4'
version ="%d.%d.%s" % (major,minor,serial) version ="%d.%d.%s" % (major,minor,serial)

View File

@ -229,6 +229,8 @@ int obi_clean(const char* dms_name,
return -1; return -1;
} }
seq_count = (i_view->infos)->line_count;
// Open the sequence column // Open the sequence column
if (strcmp((i_view->infos)->view_type, VIEW_TYPE_NUC_SEQS) == 0) if (strcmp((i_view->infos)->view_type, VIEW_TYPE_NUC_SEQS) == 0)
iseq_column = obi_view_get_column(i_view, NUC_SEQUENCE_COLUMN); iseq_column = obi_view_get_column(i_view, NUC_SEQUENCE_COLUMN);
@ -245,7 +247,7 @@ int obi_clean(const char* dms_name,
} }
// Open the sample column if there is one // Open the sample column if there is one
if ((strcmp(sample_column_name, "") == 0) || (sample_column_name == NULL)) if ((strcmp(sample_column_name, "") == 0) || (sample_column_name == NULL) || (seq_count == 0))
{ {
fprintf(stderr, "Info: No sample information provided, assuming one sample.\n"); fprintf(stderr, "Info: No sample information provided, assuming one sample.\n");
sample_column = obi_view_get_column(i_view, COUNT_COLUMN); sample_column = obi_view_get_column(i_view, COUNT_COLUMN);
@ -340,66 +342,67 @@ int obi_clean(const char* dms_name,
return -1; return -1;
} }
// Build kmer tables if (seq_count > 0)
ktable = hash_seq_column(i_view, iseq_column, 0);
if (ktable == NULL)
{ {
obi_set_errno(OBI_CLEAN_ERROR); // Build kmer tables
obidebug(1, "\nError building kmer tables before aligning"); ktable = hash_seq_column(i_view, iseq_column, 0);
return -1; if (ktable == NULL)
} {
obi_set_errno(OBI_CLEAN_ERROR);
obidebug(1, "\nError building kmer tables before aligning");
return -1;
}
seq_count = (i_view->infos)->line_count; // Allocate arrays for sample counts otherwise reading in mapped files takes longer
complete_sample_count_array = (int*) malloc(seq_count * sample_count * sizeof(int));
// Allocate arrays for sample counts otherwise reading in mapped files takes longer if (complete_sample_count_array == NULL)
complete_sample_count_array = (int*) malloc(seq_count * sample_count * sizeof(int)); {
if (complete_sample_count_array == NULL) obi_set_errno(OBI_MALLOC_ERROR);
{ obidebug(1, "\nError allocating memory for the array of sample counts, size: %lld", seq_count * sample_count * sizeof(int));
obi_set_errno(OBI_MALLOC_ERROR); return -1;
obidebug(1, "\nError allocating memory for the array of sample counts, size: %lld", seq_count * sample_count * sizeof(int)); }
return -1;
}
for (samp=0; samp < sample_count; samp++)
{
for (k=0; k<seq_count; k++)
complete_sample_count_array[k+(samp*seq_count)] = obi_get_int_with_elt_idx_and_col_p_in_view(i_view, sample_column, k, samp);
}
// Allocate arrays for blobs otherwise reading in mapped files takes longer
blob_array = (Obi_blob_p*) malloc(seq_count * sizeof(Obi_blob_p));
if (blob_array == NULL)
{
obi_set_errno(OBI_MALLOC_ERROR);
obidebug(1, "\nError allocating memory for the array of blobs");
return -1;
}
for (k=0; k<seq_count; k++)
{
blob_array[k] = obi_get_blob_with_elt_idx_and_col_p_in_view(i_view, iseq_column, k, 0);
}
// Allocate alignment result array (byte at 0 if not aligned yet,
// 1 if sequence at index has a similarity above the threshold with the current sequence,
// 2 if sequence at index has a similarity below the threshold with the current sequence)
alignment_result_array = (byte_t*) calloc(seq_count, sizeof(byte_t));
if (alignment_result_array == NULL)
{
obi_set_errno(OBI_MALLOC_ERROR);
obidebug(1, "\nError allocating memory for alignment result array");
return -1;
}
// Initialize all sequences to singletons or NA if no sequences in that sample
for (k=0; k<seq_count; k++)
{
for (samp=0; samp < sample_count; samp++) for (samp=0; samp < sample_count; samp++)
{ {
if (obi_get_int_with_elt_idx_and_col_p_in_view(i_view, sample_column, k, samp) != OBIInt_NA) // Only initialize samples where there are some sequences for (k=0; k<seq_count; k++)
complete_sample_count_array[k+(samp*seq_count)] = obi_get_int_with_elt_idx_and_col_p_in_view(i_view, sample_column, k, samp);
}
// Allocate arrays for blobs otherwise reading in mapped files takes longer
blob_array = (Obi_blob_p*) malloc(seq_count * sizeof(Obi_blob_p));
if (blob_array == NULL)
{
obi_set_errno(OBI_MALLOC_ERROR);
obidebug(1, "\nError allocating memory for the array of blobs");
return -1;
}
for (k=0; k<seq_count; k++)
{
blob_array[k] = obi_get_blob_with_elt_idx_and_col_p_in_view(i_view, iseq_column, k, 0);
}
// Allocate alignment result array (byte at 0 if not aligned yet,
// 1 if sequence at index has a similarity above the threshold with the current sequence,
// 2 if sequence at index has a similarity below the threshold with the current sequence)
alignment_result_array = (byte_t*) calloc(seq_count, sizeof(byte_t));
if (alignment_result_array == NULL)
{
obi_set_errno(OBI_MALLOC_ERROR);
obidebug(1, "\nError allocating memory for alignment result array");
return -1;
}
// Initialize all sequences to singletons or NA if no sequences in that sample
for (k=0; k<seq_count; k++)
{
for (samp=0; samp < sample_count; samp++)
{ {
if (obi_set_char_with_elt_idx_and_col_p_in_view(o_view, status_column, k, samp, 's') < 0) if (obi_get_int_with_elt_idx_and_col_p_in_view(i_view, sample_column, k, samp) != OBIInt_NA) // Only initialize samples where there are some sequences
{ {
obidebug(1, "\nError initializing all sequences to singletons"); if (obi_set_char_with_elt_idx_and_col_p_in_view(o_view, status_column, k, samp, 's') < 0)
return -1; {
obidebug(1, "\nError initializing all sequences to singletons");
return -1;
}
} }
} }
} }
@ -551,17 +554,20 @@ int obi_clean(const char* dms_name,
} }
} }
free_kmer_tables(ktable, seq_count); if (seq_count > 0)
free(complete_sample_count_array); {
free(blob_array); free_kmer_tables(ktable, seq_count);
free(alignment_result_array); free(complete_sample_count_array);
free(blob_array);
free(alignment_result_array);
}
fprintf(stderr, "\n"); fprintf(stderr, "\n");
if (stop) if (stop)
return -1; return -1;
if (heads_only) if (heads_only && (seq_count > 0))
{ {
line_selection = malloc((((o_view->infos)->line_count) + 1) * sizeof(index_t)); line_selection = malloc((((o_view->infos)->line_count) + 1) * sizeof(index_t));
if (line_selection == NULL) if (line_selection == NULL)
@ -635,7 +641,7 @@ int obi_clean(const char* dms_name,
} }
// Flag the end of the line selection // Flag the end of the line selection
if (heads_only) if (heads_only && (seq_count > 0))
line_selection[l] = -1; line_selection[l] = -1;
// Create new view with line selection if heads only // Create new view with line selection if heads only