Optimize amplicon overlap calculations and fix off-by-one errors

... in all interval clculations

varvamp/scripts/primers.py

@@ -386,13 +386,13 @@ def find_best_primers(left_primer_candidates, right_primer_candidates):
         primer_candidates.sort(key=lambda x: (x[3], x[1]))
         # ini everything with the primer with the lowest penalty
         to_retain = [primer_candidates[0]]
-        primer_ranges = list(range(primer_candidates[0][1], primer_candidates[0][2]+1))
+        primer_ranges = list(range(primer_candidates[0][1], primer_candidates[0][2]))
         primer_set = set(primer_ranges)
 
         for primer in primer_candidates:
             # get the thirds of the primer, only consider the middle
             thirds_len = int((primer[2] - primer[1])/3)
-            primer_positions = list(range(primer[1] + thirds_len, primer[2] - thirds_len + 1))
+            primer_positions = list(range(primer[1] + thirds_len, primer[2] - thirds_len))
             # check if none of the nucleotides of the next primer
             # are already covered by a better primer
             if not any(x in primer_positions for x in primer_set):

varvamp/scripts/qpcr.py

@@ -211,13 +211,13 @@ def assess_amplicons(left_subset, right_subset, qpcr_probes, probe, majority_con
             if "LEFT" in probe:
                 if not qpcr_probes[probe][1] in range(
                         left_primer[2] + config.QPROBE_DISTANCE[0],
-                        left_primer[2] + config.QPROBE_DISTANCE[1] + 1
+                        left_primer[2] + config.QPROBE_DISTANCE[1]
                 ):
                     continue
             elif "RIGHT" in probe:
                 if not right_primer[1] in range(
                         qpcr_probes[probe][2] + config.QPROBE_DISTANCE[0],
-                        qpcr_probes[probe][2] + config.QPROBE_DISTANCE[1] + 1
+                        qpcr_probes[probe][2] + config.QPROBE_DISTANCE[1]
 
                 ):
                     continue
@@ -297,20 +297,17 @@ def process_single_amplicon_deltaG(amplicon, majority_consensus):
     Process a single amplicon to test its deltaG and apply filtering.
     This function will be called concurrently by multiple threads.
     """
-    start = amplicon["LEFT"][1]
-    stop = amplicon["RIGHT"][2]
-    seq = majority_consensus[start:stop]
+    seq = majority_consensus[amplicon["LEFT"][1]:amplicon["RIGHT"][2]]
     seq = seq.replace("N", "")
     seq = seq.replace("n", "")
-    amp_positions = list(range(start, stop + 1))
     # check if the amplicon overlaps with an amplicon that was previously
     # found and had a high enough deltaG
     min_temp = min((primers.calc_temp(amplicon["LEFT"][0]),
                     primers.calc_temp(amplicon["RIGHT"][0])))
     # calculate deltaG at the minimal primer temp
     amplicon["deltaG"] = seqfold.dg(seq, min_temp)
 
-    return amp_positions, amplicon
+    return amplicon
 
 
 def test_amplicon_deltaG_parallel(qpcr_schemes_candidates, majority_consensus, n_to_test, deltaG_cutoff, n_threads):
@@ -320,7 +317,6 @@ def test_amplicon_deltaG_parallel(qpcr_schemes_candidates, majority_consensus, n
     for processing amplicons in parallel.
     """
     final_amplicons = []
-    amplicon_set = set()
 
     # Create a pool of processes to handle the concurrent processing
     with multiprocessing.Pool(processes=n_threads) as pool:
@@ -334,15 +330,20 @@ def test_amplicon_deltaG_parallel(qpcr_schemes_candidates, majority_consensus, n
         # process amplicons concurrently
         results = pool.starmap(process_single_amplicon_deltaG, [(amp, majority_consensus) for amp in amplicons])
         # Process the results
-        for amp_positions, amp in results:
+        retained_ranges = []
+        for amp in results:
             # check if the amplicon overlaps with an amplicon that was previously
             # found and had a high enough deltaG
-            if any(x in amp_positions for x in amplicon_set):
+            if amp["deltaG"] <= deltaG_cutoff:
                 continue
-            # and if this passes cutoff make a dict entry and do not allow further
-            # amplicons in that region (they will have a lower penalty)
-            if amp["deltaG"] > deltaG_cutoff:
+            amp_range = range(amp["LEFT"][1], amp["RIGHT"][2])
+            overlaps_retained = False
+            for r in retained_ranges:
+                if amp_range.start < r.stop and r.start < amp_range.stop:
+                    overlaps_retained = True
+                    break
+            if not overlaps_retained:
                 final_amplicons.append(amp)
-                amplicon_set.update(amp_positions)
+                retained_ranges.append(amp_range)
 
     return final_amplicons

varvamp/scripts/scheme.py

@@ -319,15 +319,15 @@ def get_overlapping_primers(dimer, left_primer_candidates, right_primer_candidat
         overlapping_primers_temp = []
         thirds_len = int((primer[3][2] - primer[3][1]) / 3)
         # get the middle third of the primer (here are the previously excluded primers)
-        overlap_set = set(range(primer[3][1] + thirds_len, primer[3][2] - thirds_len + 1))
+        overlap_set = set(range(primer[3][1] + thirds_len, primer[3][2] - thirds_len))
         # check in which list to look for them
         if "RIGHT" in primer[2]:
             primers_to_test = right_primer_candidates
         else:
             primers_to_test = left_primer_candidates
         # and check this list for all primers that overlap
         for potential_new in primers_to_test:
-            primer_positions = list(range(potential_new[1], potential_new[2]+1))
+            primer_positions = list(range(potential_new[1], potential_new[2]))
             if not any(x in primer_positions for x in overlap_set):
                 continue
             overlapping_primers_temp.append((primer[0], primer[1], primer[2], potential_new))
@@ -399,15 +399,13 @@ def find_single_amplicons(amplicons, all_primers, n):
     # find lowest non-overlapping
     for amp in sorted_amplicons:
         overlaps_retained = False
-        amp_range = range(amp["LEFT"][1], amp["RIGHT"][2]+1)
+        amp_range = range(amp["LEFT"][1], amp["RIGHT"][2])
         for r in retained_ranges:
-            if amp_range.start in r or amp_range.stop in r or r.start in amp_range or r.stop in amp_range:
+            if amp_range.start < r.stop and r.start < amp_range.stop:
                 overlaps_retained = True
                 break
         if not overlaps_retained:
-            retained_ranges.append(
-                range(amp["LEFT"][1], amp["RIGHT"][2]+1)
-            )
+            retained_ranges.append(amp_range)
             to_retain.append(amp)
             if len(to_retain) == n:
                 break