fixed alignment bug, tweaked config, fixed plot warnings

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"

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

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
 

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()

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,18 +112,18 @@ 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
             # following amplicon.
             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)