added fasta parsing for compatible primers and revisited primer dimer evaluation

docs/how_varvamp_works.md

@@ -56,6 +56,24 @@ To search for the best amplicon, varVAMP uses a graph based approach.
 7. Test amplicons for their [minimal free energy](https://en.wikipedia.org/wiki/Gibbs_free_energy) at their lowest primer temperature with [`seqfold`](https://github.com/Lattice-Automation/seqfold) and filter to avoid secondary structures. Amplicons with large potential deletions (>QAMPLICON_DEL_CUTOFF) will be ignored. Smaller deletions will be accepted.
 8. Take the best qPCR schemes of overlapping schemes.
 
+
+#### Multi-processing
+varVAMP can use multiple cores at some steps in the workflow. If you have performance issues
+at the following steps it might be worth increasing the number of cores.
+
+1. All workflows:
+- BLAST search: Each amplicon is searched for off-targets in the BLAST database.
+- Primer dimer search against a user-provided list of primers.
+- Primer search: Each kmer is evaluated as a potential primer.
+
+2. Tiled workflow:
+- Final primer dimer solve.
+
+3. qPCR workflow:
+- Probe search: Each kmer is evaluated as a potential probe.
+- Amplicon search: Each potential amplicon is checked for its parameters.
+- Amplicon deltaG calculation: Each amplicon is checked for potential secondary structures.
+
 #### Penalty calculation
 
 ```python

varvamp/command.py

@@ -73,7 +73,7 @@ def get_args(sysargs):
         par.add_argument(
             "-th",
             "--threads",
-            help="number of threads for BLAST search and deltaG calculations",
+            help="number of threads",
             metavar="1",
             type=int,
             default=1
@@ -85,6 +85,14 @@ def get_args(sysargs):
             type=str,
             default="varVAMP"
         )
+        par.add_argument(
+            "--compatible-primers",
+            metavar="None",
+            type=str,
+            default=None,
+            help="FASTA primer file with which new primers should not form dimers. Sequences >40 nt are ignored. Can significantly increase runtime."
+        )
+
     for par in (SINGLE_parser, TILED_parser):
         par.add_argument(
             "-t",
@@ -137,7 +145,6 @@ def get_args(sysargs):
     QPCR_parser.add_argument(
         "-t",
         "--threshold",
-        required=True,
         metavar="0.9",
         type=float,
         default=0.9,
@@ -184,9 +191,13 @@ def shared_workflow(args, log_file):
     """
     # start varvamp
     logging.varvamp_progress(log_file, mode=args.mode)
-
     # read in alignment and preprocess
     preprocessed_alignment = alignment.preprocess(args.input[0])
+    # read in external primer sequences with which new primers should not form dimers
+    if args.compatible_primers is not None:
+        compatible_primers = primers.parse_primer_fasta(args.compatible_primers)
+    else:
+        compatible_primers = None
     # check alignment length and number of gaps and report if its significantly more/less than expected
     logging.check_alignment_length(preprocessed_alignment, log_file)
     logging.check_gaped_sequences(preprocessed_alignment, log_file)
@@ -276,7 +287,8 @@ def shared_workflow(args, log_file):
     left_primer_candidates, right_primer_candidates = primers.find_primers(
         kmers,
         ambiguous_consensus,
-        alignment_cleaned
+        alignment_cleaned,
+        args.threads
     )
     for primer_type, primer_candidates in [("+", left_primer_candidates), ("-", right_primer_candidates)]:
         if not primer_candidates:
@@ -292,7 +304,22 @@ def shared_workflow(args, log_file):
         progress_text=f"{len(left_primer_candidates)} fw and {len(right_primer_candidates)} rv potential primers"
     )
 
-    return alignment_cleaned, majority_consensus, ambiguous_consensus, primer_regions, left_primer_candidates, right_primer_candidates, potential_primer_regions
+    # filter primers against user-provided list of compatible primers, can use multi-processing
+    if compatible_primers is not None:
+        left_primer_candidates = primers.filter_non_dimer_candidates(
+            left_primer_candidates, compatible_primers, args.threads
+        )
+        right_primer_candidates = primers.filter_non_dimer_candidates(
+            right_primer_candidates, compatible_primers, args.threads
+        )
+        logging.varvamp_progress(
+            log_file,
+            progress=0.65,
+            job="Filtering primers against provided primers.",
+            progress_text=f"{len(left_primer_candidates)} fw and {len(right_primer_candidates)} rv primers after filtering"
+        )
+
+    return alignment_cleaned, majority_consensus, ambiguous_consensus, primer_regions, left_primer_candidates, right_primer_candidates, potential_primer_regions, compatible_primers
 
 
 def single_and_tiled_shared_workflow(args, left_primer_candidates, right_primer_candidates, potential_primer_regions, data_dir, log_file):
@@ -323,8 +350,7 @@ def single_and_tiled_shared_workflow(args, left_primer_candidates, right_primer_
     )
     if not amplicons:
         logging.raise_error(
-            "no amplicons found. Increase the max amplicon length or \
-            number of ambiguous bases or lower threshold!\n",
+            "no amplicons found. Increase the max amplicon length or number of ambiguous bases or lower threshold!\n",
             log_file,
             exit=True
         )
@@ -383,21 +409,9 @@ def tiled_workflow(args, amplicons, left_primer_candidates, right_primer_candida
         amplicon_graph
     )
 
-    # check for dimers
-    dimers_not_solved = scheme.check_and_solve_heterodimers(
-        amplicon_scheme,
-        left_primer_candidates,
-        right_primer_candidates,
-        all_primers)
-    if dimers_not_solved:
-        logging.raise_error(
-            f"varVAMP found {len(dimers_not_solved)} primer dimers without replacements. Check the dimer file and perform the PCR for incomaptible amplicons in a sperate reaction.",
-            log_file
-        )
-        reporting.write_dimers(results_dir, dimers_not_solved)
-
     # evaluate coverage
-    # ATTENTION: Genome coverage of the scheme might still change slightly through resolution of primer dimers, but this potential, minor inaccuracy is currently accepted.
+    # ATTENTION: Genome coverage of the scheme might still change slightly through resolution of primer dimers,
+    # but this potential, minor inaccuracy is currently accepted.
     percent_coverage = round(coverage/len(ambiguous_consensus)*100, 2)
     logging.varvamp_progress(
         log_file,
@@ -414,10 +428,37 @@ def tiled_workflow(args, amplicons, left_primer_candidates, right_primer_candida
             "\t - relax primer settings (not recommended)\n",
             log_file
         )
+
+    # check for dimers
+    dimers_not_solved, n_initial_dimers = scheme.check_and_solve_heterodimers(
+        amplicon_scheme,
+        left_primer_candidates,
+        right_primer_candidates,
+        all_primers,
+        args.threads
+    )
+
+    # report dimers solve
+    if n_initial_dimers > 0 and not dimers_not_solved:
+        logging.varvamp_progress(
+            log_file,
+            progress=0.95,
+            job="Trying to solve primer dimers.",
+            progress_text=f"all dimers (n={n_initial_dimers}) could be resolved"
+        )
+    elif dimers_not_solved:
+        logging.varvamp_progress(
+            log_file,
+            progress=0.95,
+            job="Trying to solve primer dimers.",
+            progress_text=f"{len(dimers_not_solved)}/{n_initial_dimers} dimers could not be resolved"
+        )
+        reporting.write_dimers(results_dir, dimers_not_solved)
+
     return amplicon_scheme
 
 
-def qpcr_workflow(args, data_dir, alignment_cleaned, ambiguous_consensus, majority_consensus, left_primer_candidates, right_primer_candidates, log_file):
+def qpcr_workflow(args, data_dir, alignment_cleaned, ambiguous_consensus, majority_consensus, left_primer_candidates, right_primer_candidates, compatible_primers, log_file):
     """
     part of the workflow specific for the tiled mode
     """
@@ -428,7 +469,7 @@ def qpcr_workflow(args, data_dir, alignment_cleaned, ambiguous_consensus, majori
     )
     if not probe_regions:
         logging.raise_error(
-            "no regions that fullfill probe criteria! lower threshold or increase number of ambiguous chars in probe\n",
+            "no regions that fulfill probe criteria! lower threshold or increase number of ambiguous chars in probe\n",
             log_file,
             exit=True
         )
@@ -440,7 +481,7 @@ def qpcr_workflow(args, data_dir, alignment_cleaned, ambiguous_consensus, majori
         config.QPROBE_SIZES
     )
     # find potential probes
-    qpcr_probes = qpcr.get_qpcr_probes(probe_kmers, ambiguous_consensus, alignment_cleaned)
+    qpcr_probes = qpcr.get_qpcr_probes(probe_kmers, ambiguous_consensus, alignment_cleaned, args.threads)
     if not qpcr_probes:
         logging.raise_error(
             "no qpcr probes found\n",
@@ -454,8 +495,21 @@ def qpcr_workflow(args, data_dir, alignment_cleaned, ambiguous_consensus, majori
         progress_text=f"{len(qpcr_probes)} potential qPCR probes"
     )
 
+    # filter primers against non-dimer sequences if provided
+    if compatible_primers is not None:
+        qpcr_probes = primers.filter_non_dimer_candidates(
+            qpcr_probes, compatible_primers, args.threads)
+        logging.varvamp_progress(
+            log_file,
+            progress=0.75,
+            job="Filtering probes against provided primer sequences.",
+            progress_text=f"{len(qpcr_probes)} potential qPCR probes after filtering"
+        )
+
     # find unique amplicons with a low penalty and an internal probe
-    qpcr_scheme_candidates = qpcr.find_qcr_schemes(qpcr_probes, left_primer_candidates, right_primer_candidates, majority_consensus, ambiguous_consensus)
+    qpcr_scheme_candidates = qpcr.find_qcr_schemes(
+        qpcr_probes, left_primer_candidates, right_primer_candidates, majority_consensus, ambiguous_consensus, args.threads
+    )
     if not qpcr_scheme_candidates:
         logging.raise_error(
             "no qPCR scheme candidates found. lower threshold or increase number of ambiguous chars in primer and/or probe\n",
@@ -516,7 +570,7 @@ def main():
         blast.check_BLAST_installation(log_file)
 
     # mode unspecific part of the workflow
-    alignment_cleaned, majority_consensus, ambiguous_consensus, primer_regions, left_primer_candidates, right_primer_candidates, potential_primer_regions = shared_workflow(args, log_file)
+    alignment_cleaned, majority_consensus, ambiguous_consensus, primer_regions, left_primer_candidates, right_primer_candidates, potential_primer_regions, compatible_primers = shared_workflow(args, log_file)
 
     # write files that are shared in all modes
     reporting.write_regions_to_bed(primer_regions, args.name, data_dir)
@@ -600,6 +654,7 @@ def main():
             majority_consensus,
             left_primer_candidates,
             right_primer_candidates,
+            compatible_primers,
             log_file
         )
 

varvamp/scripts/default_config.py

@@ -9,7 +9,7 @@
     'PCR_DNA_CONC', 'PCR_DNTP_CONC', 'PCR_DV_CONC', 'PCR_MV_CONC',
     'PRIMER_3_PENALTY', 'PRIMER_GC_END', 'PRIMER_GC_PENALTY',
     'PRIMER_GC_RANGE', 'PRIMER_HAIRPIN', 'PRIMER_MAX_BASE_PENALTY',
-    'PRIMER_MAX_DIMER_TMP', 'PRIMER_MAX_DINUC_REPEATS', 'PRIMER_MAX_POLYX',
+    'PRIMER_MAX_DIMER_TMP', 'PRIMER_MAX_DIMER_DELTAG', 'PRIMER_MAX_DINUC_REPEATS', 'PRIMER_MAX_POLYX',
     'PRIMER_MIN_3_WITHOUT_AMB', 'PRIMER_PERMUTATION_PENALTY',
     'PRIMER_SIZES', 'PRIMER_SIZE_PENALTY',
     'PRIMER_TMP', 'PRIMER_TM_PENALTY',
@@ -30,7 +30,10 @@
 PRIMER_HAIRPIN = 47  # max melting temp for secondary structure
 PRIMER_GC_END = (1, 3)  # min/max GCs in the last 5 bases of the 3' end
 PRIMER_MIN_3_WITHOUT_AMB = 3  # min len of 3' without ambiguous charaters
-PRIMER_MAX_DIMER_TMP = 47  # max melting temp for dimers (homo- or heterodimers)
+# primer dimer constraints
+PRIMER_MAX_DIMER_TMP = 35  # max melting temp for dimers, lower temperature means more stringent filtering
+PRIMER_MAX_DIMER_DELTAG = -9000  # max allowed gibbs free energy for dimer formation, higher values mean more stringent filtering
+END_OVERLAP = 5  # maximum allowed nt overlap between ends of primers to be considered a dimer
 
 # QPCR parameters
 # basic probe parameters
@@ -42,7 +45,6 @@
 QPRIMER_DIFF = 2  # maximal temperature diff of qPCR primers
 QPROBE_TEMP_DIFF = (5, 10)  # min/max temp diff between probe and primers
 QPROBE_DISTANCE = (4, 15)  # min/max distance to the primer on the same strand
-END_OVERLAP = 5  # maximum allowed nt overlap between the ends of probe and primer
 QAMPLICON_LENGTH = (70, 200)  # min/max length of the qPCR amplicon
 QAMPLICON_GC = (40, 60)  # GC min/max of the qPCR amplicon
 QAMPLICON_DEL_CUTOFF = 4  # consider regions of the alignment for deltaG calculation if they have smaller deletions than cutoff

varvamp/scripts/logging.py

@@ -107,13 +107,12 @@ def raise_arg_errors(args, log_file):
             "degeneracy. Consider reducing.",
             log_file
         )
-    if args.database is not None:
-        if args.threads < 1:
-            raise_error(
-                "number of threads cannot be smaller than 1.",
-                log_file,
-                exit=True
-            )
+    if args.threads < 1:
+        raise_error(
+            "number of threads cannot be smaller than 1.",
+            log_file,
+            exit=True
+        )
     if args.mode in ("tiled", "single"):
         if args.opt_length > args.max_length:
             raise_error(
@@ -281,7 +280,7 @@ def confirm_config(args, log_file):
 
     # check if all variables exists
     all_vars = [
-        # arg independent TILED, SINGLE mode
+        # arg independent all modes
         (
             "PRIMER_TMP",
             "PRIMER_GC_RANGE",
@@ -291,6 +290,7 @@ def confirm_config(args, log_file):
             "PRIMER_MAX_DINUC_REPEATS",
             "PRIMER_GC_END",
             "PRIMER_MAX_DIMER_TMP",
+            "PRIMER_MAX_DIMER_DELTAG",
             "PRIMER_MIN_3_WITHOUT_AMB",
             "PCR_MV_CONC",
             "PCR_DV_CONC",
@@ -301,7 +301,8 @@ def confirm_config(args, log_file):
             "PRIMER_SIZE_PENALTY",
             "PRIMER_MAX_BASE_PENALTY",
             "PRIMER_3_PENALTY",
-            "PRIMER_PERMUTATION_PENALTY"
+            "PRIMER_PERMUTATION_PENALTY",
+            "END_OVERLAP"
         ),
         # arg independent QPCR mode
         (
@@ -312,7 +313,6 @@ def confirm_config(args, log_file):
             "QPRIMER_DIFF",
             "QPROBE_TEMP_DIFF",
             "QPROBE_DISTANCE",
-            "END_OVERLAP",
             "QAMPLICON_LENGTH",
             "QAMPLICON_GC",
             "QAMPLICON_DEL_CUTOFF"
@@ -408,7 +408,7 @@ def confirm_config(args, log_file):
         ("max base penalty", config.PRIMER_MAX_BASE_PENALTY),
         ("primer permutation penalty", config.PRIMER_PERMUTATION_PENALTY),
         ("qpcr flanking primer difference", config.QPRIMER_DIFF),
-        ("probe and primer end overlap", config.END_OVERLAP),
+        ("end overlap", config.END_OVERLAP),
         ("qpcr deletion size still considered for deltaG calculation", config.QAMPLICON_DEL_CUTOFF),
         ("maximum difference between primer and blast db", config.BLAST_MAX_DIFF),
         ("multiplier of the maximum length for non-specific amplicons", config.BLAST_SIZE_MULTI),
@@ -446,15 +446,20 @@ def confirm_config(args, log_file):
         )
     if config.PRIMER_HAIRPIN < 0:
         raise_error(
-            "decreasing hairpin melting temp to negative values "
-            "will influence successful primer search!",
+            "decreasing hairpin melting temp to negative values will influence successful primer search!",
+            log_file
+        )
+    if config.PRIMER_MAX_DIMER_TMP < 21:
+        raise_error(
+            "there is no need to set max dimer melting temp below room temperature.",
             log_file
         )
-    if config.PRIMER_MAX_DIMER_TMP < 0:
+    if config.PRIMER_MAX_DIMER_DELTAG > -6000:
         raise_error(
-            "there is no need to set max dimer melting temp below 0.",
+            "primer interactions with deltaG values higher than -6000 are generally considered unproblematic.",
             log_file
         )
+
     if config.PRIMER_MAX_BASE_PENALTY < 8:
         raise_error(
             "decreasing the base penalty will filter out more primers.",
@@ -595,7 +600,6 @@ def goodbye_message():
     messages = [
         "Thank you. Come again.",
         ">Placeholder for your advertisement<",
-        "Make primers great again!",
         "Ciao cacao!",
         "And now lets pray to the PCR gods.",
         "**bibobibobop** task finished",
@@ -611,6 +615,19 @@ def goodbye_message():
         "Task failed successfully.",
         "Never gonna give you up, never gonna let you down.",
         "Have you tried turning it off and on again?",
+        "Just try it. PCR is magic!",
+        "One goat was sacrificed for this primer design to work.",
+        "You seem trustworthy. Here's a cookie!",
+        "Owwww yeah, primers done!",
+        "These primers were designed without AI assistance.",
+        "Research is fun (if you ignore the pipetting).",
+        "Balance your primers, balance your life.",
+        "Keep calm and PCR on.",
+        "In silico we trust.",
+        "May the primers be with you.",
+        "Designing primers like a boss!",
+        "Primer design completed. Time for a break!",
+        "Eureka! Your primers are ready.",
         "Look, I am your primer scheme.",
         "Quod erat demonstrandum.",
         "Miau?",