.

Author: vsc20958

README.md

@@ -9,7 +9,7 @@ RecMpox is a command-line tool that **flags potential recombination events** in
 1. **References are required**: RecMpox compares your genomes against two reference sequences (for example, Clade Ia vs. Ib, or Ib vs. IIb), because recombination can only occur between two distinct lineages.
 2. **Alignment and diagnostic SNPs**: The two reference genomes are aligned using [Squirrel](https://github.com/aineniamh/squirrel), so that the same genomic positions correspond across all sequences. RecMpox then identifies positions where the two references differ at the same coordinates. These positions are defined as diagnostic SNPs, because they distinguish between the reference lineages. Positions where the references are identical are ignored, as they do not provide information for detecting recombination.
 3. **Consensus genome classification**: Your consensus genomes are aligned to the same references. At each diagnostic SNP, the base is classified as matching reference 1, reference 2, or other (e.g., gaps or ambiguous bases).
-4. **Flagging potential recombinants**: If both references contribute at least 5% of the diagnostic positions in a genome, RecMpox flags it as a potential recombinant, since no single lineage clearly dominates.
+4. **Flagging potential recombinants**: If both references contribute at least 10% of the diagnostic positions in a genome, RecMpox flags it as a potential recombinant, since no single lineage clearly dominates.
 5. **Recombination tracts and breakpoints**: By examining the pattern of reference matches along the genome, RecMpox infers recombination tracts and identifies their breakpoints (start and end positions). To reduce false positives, runs of fewer than 2 consecutive diagnostic SNPs are ignored.
 6 **Outputs**:
    - TSV file: or each genome, reports the number and proportion of diagnostic SNPs matching each reference, the resulting recombinant flag, and summary statistics used for tract inference.
@@ -98,6 +98,7 @@ recmpox -i accessions.txt -o output -ref Ia,Ib   # one accession per line or com
 - `-ref1_g`, `-ref2_g`: Genotype labels for TSV/HTML (default from `-ref` or accession)
 - `-include-indels`: Include diagnostic indels (default: SNPs only)
 - `-min-indel-size`: Min indel length (bp) when using `-include-indels` (default: 100)
+- `-m, --minor-ref-pct`: Minor reference % threshold for calling "potential recombinant" (default: 10). Increase to be more conservative (e.g. 15, 20).
 - `-t, --threads`: Number of threads
 - `-q, --quiet`: Log to file only
 
@@ -126,8 +127,8 @@ Intermediate files (e.g. diagnostic_snps.txt, Squirrel outputs) are written unde
 
 ## Interpretation
 
-- **No recombinant**: One ref dominates (minor ref < 5% of diagnostic sites).
-- **Potential recombinant**: Both refs contribute ≥5% (minor ref % ≥ 5%). The HTML report shows recombination tracts (beginning/end of each tract) and breakpoints between tracts. A single tract means the genome is entirely one clade (no recombination).
+- **No recombinant**: One ref dominates (minor ref < 10% of diagnostic sites).
+- **Potential recombinant**: Both refs contribute ≥10% (minor ref % ≥ 10%). The HTML report shows recombination tracts (beginning/end of each tract) and breakpoints between tracts. A single tract means the genome is entirely one clade (no recombination).
 - **High pct_other**: Many Ns, gaps, or non-ref bases at diagnostic sites (poor coverage or alignment).
 
 ## HTML Output example

conda-recipes/meta.yaml

@@ -51,7 +51,7 @@ about:
     It does not confirm recombination; it flags genomes that may be recombinant for further investigation.
     RecMpox takes two references (e.g. Clade Ia vs Ib, or IIa vs IIb), aligns them with Squirrel,
     finds diagnostic SNP positions where the references differ, then classifies each consensus genome
-    at those positions. Genomes where both refs contribute at least 5% are flagged as potential
+    at those positions. Genomes where both refs contribute at least 10% are flagged as potential
     recombinants. Recombination tracts and breakpoints are inferred along the genome.
   doc_url: https://github.com/DaanJansen94/RecMpox/blob/main/README.md
   dev_url: https://github.com/DaanJansen94/RecMpox

recmpox/recmpox.py

@@ -93,8 +93,8 @@ def resolve_ref(spec: str, work_dir: Path, label: str) -> Optional[Path]:
 
 
 # Minor ref % = smaller of pct_ref1, pct_ref2 over ALL diagnostic sites. When >= threshold, flag as potential recombinant.
-# Recombinant threshold: 5% for all (intra- and inter-clade)
-MINOR_REF_PCT_THRESHOLD = 5.0
+# Default recombinant threshold: 10% for all (intra- and inter-clade)
+MINOR_REF_PCT_THRESHOLD = 10.0
 
 # Default NCBI accessions per clade/subclade (used when -ref X,Y is given)
 REF_DEFAULTS = {
@@ -1220,6 +1220,15 @@ def main() -> None:
     optional = parser.add_argument_group("optional arguments")
     parser.add_argument("-h", "-help", "--help", action="help", help="show this help message and exit")
     optional.add_argument("--version", action="version", version=f"RecMpox v{__version__}")
+    optional.add_argument(
+        "-m",
+        "--minor-ref-pct",
+        dest="minor_ref_pct",
+        type=float,
+        default=MINOR_REF_PCT_THRESHOLD,
+        metavar="",
+        help=f"Minor reference %% threshold for calling 'potential recombinant' (default: {MINOR_REF_PCT_THRESHOLD:g}).",
+    )
     required.add_argument("-i", "-input", dest="input", type=Path, default=None, metavar="", help="FASTA file, directory of .fa/.fasta/.fna, .txt file of accessions (one per line or comma-separated), NCBI accession, or comma-separated accessions (e.g. -i ACC1,ACC2 or -i accessions.txt)")
     required.add_argument("-ref", dest="ref", type=str, default=None, metavar="", help="Reference pair: two comma-separated labels among Ia, Ib, IIa, IIb (e.g. Ia,Ib or Ib,IIb). Uses built-in defaults. Either -ref or both -ref1 and -ref2 are required.")
     required.add_argument("-ref1", type=str, default=None, metavar="", help="First reference: FASTA path or NCBI accession; overrides ref1 when using -ref. Required if -ref is not used.")
@@ -1239,6 +1248,11 @@ def main() -> None:
 
     args = parser.parse_args()
 
+    if getattr(args, "minor_ref_pct", None) is None:
+        args.minor_ref_pct = MINOR_REF_PCT_THRESHOLD
+    if args.minor_ref_pct < 0 or args.minor_ref_pct > 100:
+        parser.error("--minor-ref-pct must be between 0 and 100")
+
     if args.input is None:
         parser.error("-i/-input is required")
 
@@ -1294,7 +1308,7 @@ def _normalize_ref_label(s: str) -> str:
     if squirrel_clade is None:
         squirrel_clade = _infer_squirrel_clade(ref1_label, ref2_label)
     is_intra_clade = squirrel_clade is not None
-    minor_threshold = MINOR_REF_PCT_THRESHOLD
+    minor_threshold = float(getattr(args, "minor_ref_pct", MINOR_REF_PCT_THRESHOLD))
     if squirrel_clade == "cladei":
         logger.info("Inferred Squirrel --clade cladei from ref1_g/ref2_g (Clade I)")
     elif squirrel_clade == "cladeii":
@@ -1526,8 +1540,8 @@ def row(r: Dict[str, Any]) -> str:
         logger.info("Wrote %s (%d potential recombinant samples)", out_sites_tsv, len(rec_samples))
 
     recombinant_threshold_note = (
-        "A 5% threshold is used for all recombinant calls: "
-        "when minor ref % ≥ 5%, the sample is flagged as potential recombinant."
+        f"A {minor_threshold:g}% threshold is used for all recombinant calls: "
+        f"when minor ref % ≥ {minor_threshold:g}%, the sample is flagged as potential recombinant."
     )
     other_explanation = (
         "%% other = diagnostic sites where the query neither matched %s nor %s (different base; at SNPs, gap/N count as other)."