diff --git a/varvamp/__init__.py b/varvamp/__init__.py index 19c42d0..1ae9eb8 100644 --- a/varvamp/__init__.py +++ b/varvamp/__init__.py @@ -1,3 +1,3 @@ """Tool to design amplicons for highly variable virusgenomes""" _program = "varvamp" -__version__ = "0.9.3" +__version__ = "0.9.4" diff --git a/varvamp/scripts/alignment.py b/varvamp/scripts/alignment.py index 5b95eb9..1bbda28 100644 --- a/varvamp/scripts/alignment.py +++ b/varvamp/scripts/alignment.py @@ -171,7 +171,7 @@ def clean_gaps(alignment, gaps_to_mask): stop = region[0] masked_seq_temp = sequence[1][start:stop] # check if the deletion is at the start - if len(masked_seq_temp) == 0: + if start == 0 and len(masked_seq_temp) == 0: masked_seq = mask # check if deletion is not at start elif start == 0 and len(masked_seq_temp) != 0: @@ -181,7 +181,7 @@ def clean_gaps(alignment, gaps_to_mask): masked_seq = masked_seq + mask + masked_seq_temp start = region[1]+1 if max(gaps_to_mask)[1] < len(sequence[1])-1: - # append the last gaps if it is not + # append the last seq if no gap is at # the end of the sequence start = max(gaps_to_mask)[1] stop = len(sequence[1])-1 diff --git a/varvamp/scripts/config.py b/varvamp/scripts/config.py index dc91200..3c7b9ec 100644 --- a/varvamp/scripts/config.py +++ b/varvamp/scripts/config.py @@ -5,7 +5,7 @@ # CAN BE CHANGED, DO NOT DELETE # basic primer parameters -PRIMER_TMP = (57, 63, 60) # melting temperatur (min, max, opt) +PRIMER_TMP = (56, 63, 60) # melting temperatur (min, max, opt) PRIMER_GC_RANGE = (35, 65, 50) # gc (min, max, opt) PRIMER_SIZES = (18, 24, 21) # size (min, max, opt) PRIMER_MAX_POLYX = 3 # max number of polyx repeats @@ -40,7 +40,7 @@ PRIMER_TM_PENALTY = 2 # temperature penalty PRIMER_GC_PENALTY = 0.2 # gc penalty PRIMER_SIZE_PENALTY = 0.5 # size penalty -PRIMER_MAX_BASE_PENALTY = 8 # max base penalty for a primer +PRIMER_MAX_BASE_PENALTY = 10 # max base penalty for a primer PRIMER_3_PENALTY = (32, 16, 8, 4, 2) # penalties for 3' mismatches PRIMER_PERMUTATION_PENALTY = 0.1 # penalty for the number of permutations diff --git a/varvamp/scripts/reporting.py b/varvamp/scripts/reporting.py index e149d14..81e0ef7 100644 --- a/varvamp/scripts/reporting.py +++ b/varvamp/scripts/reporting.py @@ -557,6 +557,11 @@ def per_base_mismatch_plot(dir, amplicon_scheme, threshold, mode="SANGER/TILED") ax[idx].xaxis.set_ticklabels(x, rotation=45) ax[idx].set_ylabel(ylabel="freq of sequences") ax[idx].set_xlabel("position") - ax[idx].set_ylim(0, 1-threshold) + # set yaxis lims reasonably + if threshold <= 0.9: + ax[idx].set_ylim(0, 1-threshold) + else: + ax[idx].set_ylim(0, 0.1) # - to pdf pdf.savefig(fig, bbox_inches='tight') + plt.close() diff --git a/varvamp/scripts/scheme.py b/varvamp/scripts/scheme.py index 43da059..52f91cc 100644 --- a/varvamp/scripts/scheme.py +++ b/varvamp/scripts/scheme.py @@ -4,6 +4,7 @@ # BUILT-INS import heapq +import math # varVAMP from varvamp.scripts import config, primers @@ -78,7 +79,7 @@ def find_amplicons(all_primers, opt_len, max_len): continue # calculate length dependend amplicon costs as the cumulative primer # score multiplied by the fold length of the optimal length. - amplicon_costs = (right_primer[3] + left_primer[3])*(amplicon_length/opt_len) + amplicon_costs = (right_primer[3] + left_primer[3])*math.exp(amplicon_length/opt_len) amplicon_name = "amplicon_"+str(amplicon_number) amplicon_dict[amplicon_name] = [ left_primer[1], # start @@ -111,8 +112,8 @@ def create_amplicon_graph(amplicons, min_overlap): current_amplicon = amplicons[current] start = current_amplicon[0] + current_amplicon[4]/2 stop = current_amplicon[1] - min_overlap - for next in amplicons: - next_amplicon = amplicons[next] + for next_amp in amplicons: + next_amplicon = amplicons[next_amp] # check if the next amplicon lies within the start/stop range of # the current amplicon and if its non-overlapping part is large # enough to ensure space for a primer and the min overlap of the @@ -120,9 +121,9 @@ def create_amplicon_graph(amplicons, min_overlap): if not all((start <= next_amplicon[0] <= stop, next_amplicon[1] > current_amplicon[1] + next_amplicon[4]/2)): continue if current not in amplicon_graph: - amplicon_graph[current] = {next: next_amplicon[5]} + amplicon_graph[current] = {next_amp: next_amplicon[5]} else: - amplicon_graph[current][next] = next_amplicon[5] + amplicon_graph[current][next_amp] = next_amplicon[5] # return a graph object return Graph(nodes, amplicon_graph)