-
Notifications
You must be signed in to change notification settings - Fork 69
Description
Description
When using PRIMER_TASK=generic with sequences where all primers should be rejected for high self-complementarity, the EXPLAIN output shows an incorrect "ok 1" count instead of "ok 0". This occurs due to a lazy evaluation optimization where one primer's self-complementarity is never calculated if no valid primer pairs can be formed.
This can make reporting stats screwy for repetitive regions and throw off unit testing. Otherwise, I can't imagine a functional impact. There may be a slight performance impact as the error causes all primer pairs to be evaluated, but they fail fast since the pairs don't exist.
Steps to Reproduce
# Minimal test case
echo 'SEQUENCE_ID=gatc_bug
SEQUENCE_TEMPLATE=GATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATC
PRIMER_TASK=generic
PRIMER_MIN_SIZE=20
PRIMER_MAX_SIZE=20
PRIMER_EXPLAIN_FLAG=1
PRIMER_THERMODYNAMIC_OLIGO_ALIGNMENT=0
PRIMER_PRODUCT_SIZE_RANGE=30-52
=' | primer3_coreNote: Only affects PRIMER_TASK=generic (not check_primers or pick_primer_list)
Note: The "ok 1" primer is never returned (PRIMER_LEFT_NUM_RETURNED=0)
Expected Behavior
PRIMER_LEFT_EXPLAIN=considered 23, low tm 6, high any compl 17, ok 0
PRIMER_RIGHT_EXPLAIN=considered 23, low tm 6, ok 17
PRIMER_PAIR_EXPLAIN=considered 0, ok 0
Actual Behavior
PRIMER_LEFT_EXPLAIN=considered 23, low tm 6, high any compl 16, ok 1
PRIMER_RIGHT_EXPLAIN=considered 23, low tm 6, ok 17
PRIMER_PAIR_EXPLAIN=considered 23, unacceptable product size 23, ok 0
Note:
- LEFT shows "ok 1" when all 17 non-low-tm primers have self_any scores of 18.0-20.0 (well above the 8.0 threshold)
- PAIR shows "considered 23" even though no valid left primers exist (the "ok 1" is a ghost primer that was never fully evaluated)
Root Cause
The issue is caused by lazy evaluation in primer3_core when self-complementarity weights are 0:
In libprimer3.cc lines 3550-3587:
if ((must_use
|| pa->file_flag // P3_FILE_FLAG check here
|| retval->output_type == primer_list
|| po_args->weights.compl_any // weight = 0 by default
|| po_args->weights.compl_end) // weight = 0 by default
&& pa->thermodynamic_oligo_alignment==0) {
// Calculate self-complementarity
oligo_compl(h, po_args, stats, dpal_arg_to_use, ...);
} else {
// LAZY EVALUATION BUG: Skip calculation entirely
h->self_any = h->self_end = ALIGN_SCORE_UNDEF; // line 3584
}When P3_FILE_FLAG=0 (default):
- Condition evaluates to FALSE (all weights are 0, file_flag is 0)
- Takes else branch:
self_any = ALIGN_SCORE_UNDEF(never calculated) - Explain counters not updated:
okstays at 1,compl_anynot incremented
When P3_FILE_FLAG=1:
- Condition evaluates to TRUE (pa->file_flag is true)
- Calls
oligo_compl()which calculates self_any - In
oligo_compl()at line 5318:ostats->ok--properly decrements ok count - Result: Correct "ok 0" in explain output
Demonstration
1. Melting Temperature Analysis for ALL 23 Primers:
# Using oligotm to calculate Tm for each primer position:
oligotm GATCGATCGATCGATCGATC # Position 0: Output: 56.580423 (REJECTED < 57°C)
oligotm ATCGATCGATCGATCGATCG # Position 1: Output: 57.363116 (OK)
oligotm TCGATCGATCGATCGATCGA # Position 2: Output: 59.312195 (OK)
oligotm CGATCGATCGATCGATCGAT # Position 3: Output: 57.363116 (OK)
oligotm GATCGATCGATCGATCGATC # Position 4: Output: 56.580423 (REJECTED < 57°C)
oligotm ATCGATCGATCGATCGATCG # Position 5: Output: 57.363116 (OK)
oligotm TCGATCGATCGATCGATCGA # Position 6: Output: 59.312195 (OK)
oligotm CGATCGATCGATCGATCGAT # Position 7: Output: 57.363116 (OK)
oligotm GATCGATCGATCGATCGATC # Position 8: Output: 56.580423 (REJECTED < 57°C)
oligotm ATCGATCGATCGATCGATCG # Position 9: Output: 57.363116 (OK)
oligotm TCGATCGATCGATCGATCGA # Position 10: Output: 59.312195 (OK)
oligotm CGATCGATCGATCGATCGAT # Position 11: Output: 57.363116 (OK)
oligotm GATCGATCGATCGATCGATC # Position 12: Output: 56.580423 (REJECTED < 57°C)
oligotm ATCGATCGATCGATCGATCG # Position 13: Output: 57.363116 (OK)
oligotm TCGATCGATCGATCGATCGA # Position 14: Output: 59.312195 (OK)
oligotm CGATCGATCGATCGATCGAT # Position 15: Output: 57.363116 (OK)
oligotm GATCGATCGATCGATCGATC # Position 16: Output: 56.580423 (REJECTED < 57°C)
oligotm ATCGATCGATCGATCGATCG # Position 17: Output: 57.363116 (OK)
oligotm TCGATCGATCGATCGATCGA # Position 18: Output: 59.312195 (OK)
oligotm CGATCGATCGATCGATCGAT # Position 19: Output: 57.363116 (OK)
oligotm GATCGATCGATCGATCGATC # Position 20: Output: 56.580423 (REJECTED < 57°C)
oligotm ATCGATCGATCGATCGATCG # Position 21: Output: 57.363116 (OK)
oligotm TCGATCGATCGATCGATCGA # Position 22: Output: 59.312195 (OK)
# Summary: 6 primers rejected for low Tm, 17 primers have acceptable Tm2. Remaining 17 all have High Self-Complementarity with primer3_core:
# For loop to test all 17 primers with acceptable Tm
SEQUENCE="GATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATC"
OK_TM_POSITIONS=(1 2 3 5 6 7 9 10 11 13 14 15 17 18 19 21 22)
for pos in "${OK_TM_POSITIONS[@]}"; do
primer=${SEQUENCE:$pos:20}
echo "Position $pos: $primer"
echo "SEQUENCE_TEMPLATE=$SEQUENCE
PRIMER_TASK=check_primers
PRIMER_LEFT_INPUT=$primer
PRIMER_THERMODYNAMIC_OLIGO_ALIGNMENT=0
PRIMER_MAX_SELF_ANY=8.0
=" | primer3_core | grep -E "(PRIMER_LEFT_NUM_RETURNED|EXPLAIN)"
echo
done
# Output for all 17 primers:
# PRIMER_LEFT_NUM_RETURNED=0Workaround
Force complete evaluation by adding P3_FILE_FLAG=1
# With workaround flags - shows CORRECT output
echo 'SEQUENCE_ID=gatc_bug_workaround
SEQUENCE_TEMPLATE=GATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATC
PRIMER_TASK=generic
PRIMER_MIN_SIZE=20
PRIMER_MAX_SIZE=20
PRIMER_EXPLAIN_FLAG=1
PRIMER_PRODUCT_SIZE_RANGE=30-52
P3_FILE_FLAG=1
=' | primer3_core
# Output with workaround:
# PRIMER_LEFT_EXPLAIN=considered 23, low tm 6, high any compl 17, ok 0
# This is the CORRECT output - all 17 primers rejected for high self_any
# PRIMER_PAIR_EXPLAIN=considered 0, ok 0
# No ghost primer pairs consideredP3_FILE_FLAG=1 forces complete evaluation by making the condition at line 3550 in libprimer3.cc evaluate to TRUE:
if ((must_use || pa->file_flag || ...)) // pa->file_flag is now TRUE
oligo_compl(h, ...); // Self-complementarity IS calculatedWithout P3_FILE_FLAG (when it's 0), the condition is FALSE and the code skips to:
else
h->self_any = ALIGN_SCORE_UNDEF; // Never calculatedThis bypasses the lazy evaluation optimization that causes incorrect EXPLAIN statistics.
Note: P3_FILE_FLAG=1 requires SEQUENCE_ID to be set.
Possible fixes
Option 1: Force evaluation when EXPLAIN_FLAG is set
if ((must_use
|| pa->file_flag
|| pa->explain_flag // Add this condition
|| retval->output_type == primer_list
|| po_args->weights.compl_any
|| po_args->weights.compl_end)
&& pa->thermodynamic_oligo_alignment==0) {Option 2: Update statistics during lazy evaluation
} else {
h->self_any = h->self_end = ALIGN_SCORE_UNDEF;
// Add: Check if primer would fail and update stats
if (pa->explain_flag) {
// Calculate self_any just for statistics
double temp_self_any = align(oligo_seq, revc_oligo_seq, ...);
if (temp_self_any > po_args->max_self_any) {
ostats->compl_any++;
ostats->ok--;
}
}
}System Information
- primer3 Version: 2.6.1
- Operating System: Linux
- Architecture: x86_64