v2.2.1: fixed minor bugs in longitudinal tracking

CITATION.cff

@@ -1,8 +1,8 @@
 cff-version: 1.2.0
 message: "If you use this software, please cite it as below."
 title: "vartracker"
-version: "2.2.0"
-date-released: "2026-05-05"
+version: "2.2.1"
+date-released: "2026-05-06"
 license: MIT
 repository-code: "https://github.com/charlesfoster/vartracker"
 url: "https://github.com/charlesfoster/vartracker"
@@ -22,7 +22,7 @@ preferred-citation:
     - family-names: Foster
       given-names: Charles
   title: "vartracker"
-  version: "2.2.0"
+  version: "2.2.1"
   doi: "10.5281/zenodo.18452274"
   url: "https://github.com/charlesfoster/vartracker"
-  date-released: "2026-05-05"
+  date-released: "2026-05-06"

README.md

@@ -167,7 +167,7 @@ Docker is a self-contained reproducible option. If you publish the image, record
 set it when running to include it in the run manifest:
 
 ```bash
-export VARTRACKER_CONTAINER_IMAGE=ghcr.io/your-org/vartracker:2.0.0
+export VARTRACKER_CONTAINER_IMAGE=ghcr.io/your-org/vartracker:2.2.1
 export VARTRACKER_CONTAINER_DIGEST=sha256:...
 ```
 
@@ -253,8 +253,10 @@ LoFreq primer-overlap rescue:
 - `bam` and `end-to-end` therefore run LoFreq with `--no-default-filter`, then apply the normal `lofreq filter` step so standard LoFreq PASS calls are unchanged.
 - With the default `--lofreq-primer-rescue auto`, the rescue step runs only when `--primer-bed` is supplied. In other words, `auto` means "use primer rescue when an amplicon primer scheme has been explicitly provided."
 - In `end-to-end` mode, the same `--primer-bed` is used for `samtools ampliconclip` and for rescue. In `bam` mode, vartracker does not clip the input BAMs; the primer BED is used only to identify primer-overlap sites for rescue.
-- Rescue candidates must be single-ALT SNPs that overlap a primer interval, fail the default LoFreq filter, and pass conservative near-fixed thresholds (`AF>=0.95`, `DP>=100`, `DP4 alt count>=95`, `QUAL>=100`, `DP4 ref count<=20`) with one-sided alternate-strand support. Indels, multi-ALT records, lower-frequency variants, and non-primer-overlap variants are not rescued by this rule.
+- Rescue candidates must be single-ALT SNPs that overlap a primer interval, fail LoFreq's default strand-bias filtering, and pass conservative near-fixed thresholds (`AF>=0.95`, `DP>=100`, `DP4 alt count>=95`, `QUAL>=100`, `DP4 ref count<=20`, minor ALT strand fraction `<=0.05`). Indels, multi-ALT records, lower-frequency variants, non-primer-overlap variants, and variants filtered for non-strand-bias reasons are not rescued by this rule.
+- The raw LoFreq calls are retained as `<sample>_variants.raw.vcf.gz` and listed in the updated spreadsheet as `raw_vcf`.
 - Rescued variants are marked with `FILTER=RESCUED_PRIMER_OVERLAP`, `INFO/PRIMER_OVERLAP`, and `INFO/RESCUED_BY=overlap_primer_interval`; per-sample details are written to `<sample>_variants.rescued.tsv` and listed in the updated spreadsheet as `lofreq_rescued_tsv`.
+- Variants called by raw LoFreq but filtered out of the final VCF are written to `<sample>_variants.filtered_out.tsv` with the LoFreq filter reason and core metrics. This is useful for auditing high-frequency calls that fail strand-bias or other LoFreq filters.
 - Use `--lofreq-primer-rescue off` to disable rescue even when a primer BED is supplied, or `--lofreq-primer-rescue on` to require rescue and fail if `--primer-bed` is missing. The rescue thresholds can be adjusted with the `--lofreq-rescue-*` options.
 
 Example amplicon run with primer rescue:
@@ -284,10 +286,11 @@ Mode-specific expectations:
 - **End-to-end mode** requires `reads1` (and optionally `reads2`); remaining fields are generated.
 
 The `bam` and `end-to-end` workflows also write two consensus FASTA columns to
-the updated Snakemake spreadsheet, plus the LoFreq rescue audit column:
+the updated Snakemake spreadsheet, plus LoFreq audit columns:
 `consensus` for a simple consensus, `iupac_consensus` for an IUPAC-aware
-consensus, and `lofreq_rescued_tsv` for the per-sample primer-overlap rescue
-table. SNPs below
+consensus, `raw_vcf` for raw LoFreq calls, `lofreq_rescued_tsv` for the
+per-sample primer-overlap rescue table, and `lofreq_filtered_out_tsv` for raw
+LoFreq records excluded from the final VCF. SNPs below
 `--consensus-snp-min-af` are ignored, SNPs from `--consensus-snp-min-af` up to
 `--consensus-snp-thresh` stay as reference bases in the simple consensus
 and become REF+ALT ambiguity codes in the IUPAC consensus, and SNPs at or above
@@ -509,7 +512,9 @@ vartracker produces several output files:
 
 - **results.csv**: Comprehensive variant analysis with all metrics
 - **results_metadata.json**: Output schema version and results metadata
+- **`<sample>_variants.raw.vcf.gz`** (`bam`/`end-to-end`): Raw LoFreq calls before default filtering and primer-overlap rescue
 - **`<sample>_variants.rescued.tsv`** (`bam`/`end-to-end`): LoFreq primer-overlap rescue audit table, empty when rescue is disabled or no variants are rescued
+- **`<sample>_variants.filtered_out.tsv`** (`bam`/`end-to-end`): Raw LoFreq calls excluded from the final VCF, including filter reason and call metrics
 - **new_mutations.csv**: Mutations not present in the first sample
 - **persistent_new_mutations.csv**: New mutations that persist to the final sample
 - **cumulative_mutations.pdf**: Plot showing mutation accumulation over time
@@ -566,7 +571,7 @@ The pipeline performs the following analysis:
 When using vartracker, please cite the software release you used. Citation metadata is provided
 in `CITATION.cff`, and GitHub releases are archived on Zenodo.
 
-- Foster, C. (2026). *vartracker* (Version 2.2.0). Zenodo. https://doi.org/10.5281/zenodo.18452274
+- Foster, C. (2026). *vartracker* (Version 2.2.1). Zenodo. https://doi.org/10.5281/zenodo.18452274
 
 Note: the DOI above is the Zenodo concept DOI for all versions; a version-specific DOI is minted by Zenodo after each GitHub release.
 

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "hatchling.build"
 
 [project]
 name = "vartracker"
-version = "2.2.0"
+version = "2.2.1"
 authors = [
     {name = "Dr Charles Foster"},
 ]

tests/test_analysis.py

@@ -133,18 +133,75 @@ def test_prepare_variant_heatmap_matrix_orders_variants_by_genome():
     expected_long_label = f"nsp2:{head}+{middle}{tail}\n(G1946GT)"
 
     expected_index = [
-        "nsp2:T809=\n(C1059T)",
+        "nsp2:T629=\n(C1059T)",
         expected_long_label,
         "S:D215G\n(A22206G)",
     ]
 
     assert list(matrix.index) == expected_index
-    assert matrix.loc["nsp2:T809=\n(C1059T)", "P0"] == 0.0
-    assert matrix.loc["nsp2:T809=\n(C1059T)", "P1"] == 0.5
+    assert matrix.loc["nsp2:T629=\n(C1059T)", "P0"] == 0.0
+    assert matrix.loc["nsp2:T629=\n(C1059T)", "P1"] == 0.5
     assert matrix.loc[expected_long_label, "P1"] == 0.8
     assert matrix.loc["S:D215G\n(A22206G)", "P1"] == 1.0
 
 
+def test_prepare_variant_heatmap_matrix_normalises_starred_synonymous_label():
+    table = pd.DataFrame(
+        [
+            {
+                "gene": "ORF1ab",
+                "amino_acid_consequence": "924F",
+                "nsp_aa_change": "",
+                "type_of_change": "*synonymous",
+                "type_of_variant": "snp",
+                "alt_freq": "0.5 / 0.0",
+                "samples": "P0 / P1",
+                "variant": "C3037T",
+                "start": 3037,
+            },
+            {
+                "gene": "ORF1ab",
+                "amino_acid_consequence": "924F",
+                "nsp_aa_change": "",
+                "type_of_change": "synonymous",
+                "type_of_variant": "snp",
+                "alt_freq": "0.0 / 0.6",
+                "samples": "P0 / P1",
+                "variant": "C3037T",
+                "start": 3037,
+            },
+        ]
+    )
+
+    matrix = _prepare_variant_heatmap_matrix(table, ["P0", "P1"], 0.2, 0.3)
+
+    assert list(matrix.index) == ["nsp3_PLpro:F106=\n(C3037T)"]
+    assert matrix.loc["nsp3_PLpro:F106=\n(C3037T)", "P0"] == 0.5
+    assert matrix.loc["nsp3_PLpro:F106=\n(C3037T)", "P1"] == 0.6
+
+
+def test_prepare_variant_heatmap_matrix_repairs_stale_stop_gained_nsp_label():
+    table = pd.DataFrame(
+        [
+            {
+                "gene": "ORF1ab",
+                "amino_acid_consequence": "L889*",
+                "nsp_aa_change": "nsp3_PLpro:71L",
+                "type_of_change": "stop_gained",
+                "type_of_variant": "snp",
+                "alt_freq": "0.052",
+                "samples": "P0",
+                "variant": "T2931A",
+                "start": 2931,
+            },
+        ]
+    )
+
+    matrix = _prepare_variant_heatmap_matrix(table, ["P0"], 0.0, 0.0)
+
+    assert list(matrix.index) == ["nsp3_PLpro:L71*\n(T2931A)"]
+
+
 def test_prepare_variant_heatmap_matrix_excludes_selected_consequence_types():
     table = pd.DataFrame(
         [
@@ -329,7 +386,7 @@ def test_process_joint_variants_only_adds_single_joint_prefix(tmp_path):
                 "aa1_weight": "",
                 "aa2_weight": "",
                 "weight_difference": "",
-                "type_of_change": "joint_joint_frameshift",
+                "type_of_change": "joint_*frameshift",
             },
             {
                 "start": 101,

tests/test_analysis_launcher.py

@@ -50,6 +50,7 @@ def test_validate_lofreq_primer_rescue_on_requires_primer_bed():
             min_alt_count=95,
             min_qual=100,
             max_ref_count=20,
+            max_minor_alt_fraction=0.05,
         )
 
 
@@ -63,4 +64,17 @@ def test_validate_lofreq_primer_rescue_rejects_invalid_threshold():
             min_alt_count=95,
             min_qual=100,
             max_ref_count=20,
+            max_minor_alt_fraction=0.05,
+        )
+
+    with pytest.raises(ValueError, match="between 0 and 1"):
+        _validate_lofreq_primer_rescue(
+            "auto",
+            None,
+            min_af=0.95,
+            min_dp=100,
+            min_alt_count=95,
+            min_qual=100,
+            max_ref_count=20,
+            max_minor_alt_fraction=1.5,
         )

tests/test_constants.py

@@ -0,0 +1,12 @@
+"""Tests for reference-coordinate formatting helpers."""
+
+from __future__ import annotations
+
+from vartracker.constants import reformat_csq_notation
+
+
+def test_reformat_csq_notation_preserves_orf1ab_stop_gained_nsp_change():
+    reformatted, nsp_change = reformat_csq_notation("ORF1ab", "889L>889*")
+
+    assert reformatted == "L889*"
+    assert nsp_change == "nsp3_PLpro:L71*"

tests/test_lofreq_primer_rescue.py

@@ -1,7 +1,10 @@
 import subprocess
 from pathlib import Path
 
-from vartracker.lofreq_primer_rescue import rescue_lofreq_primer_variants
+from vartracker.lofreq_primer_rescue import (
+    lofreq_filter_with_audit,
+    rescue_lofreq_primer_variants,
+)
 
 
 def test_rescue_lofreq_primer_variants_keeps_pass_and_rescues_overlap(
@@ -15,7 +18,10 @@ def test_rescue_lofreq_primer_variants_keeps_pass_and_rescues_overlap(
         '##INFO=<ID=DP4,Number=4,Type=Integer,Description="Strand depths">\n'
         "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n"
         "chr1\t10\t.\tA\tG\t200\t.\tAF=0.98;DP=120;DP4=1,0,118,0\n"
-        "chr1\t30\t.\tC\tT\t150\t.\tAF=0.40;DP=120;DP4=30,30,30,30\n",
+        "chr1\t15\t.\tT\tC\t210\t.\tAF=0.98;DP=125;DP4=0,5,114,5\n"
+        "chr1\t18\t.\tG\tT\t220\t.\tAF=0.98;DP=125;DP4=0,4,115,5\n"
+        "chr1\t30\t.\tC\tT\t150\t.\tAF=0.40;DP=120;DP4=30,30,30,30\n"
+        "chr1\t50\t.\tT\tC\t160\t.\tAF=0.96;DP=130;DP4=2,2,60,65\n",
         encoding="utf-8",
     )
     primers = tmp_path / "primers.bed"
@@ -27,11 +33,14 @@ def test_rescue_lofreq_primer_variants_keeps_pass_and_rescues_overlap(
         '##INFO=<ID=DP4,Number=4,Type=Integer,Description="Strand depths">\n'
         "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n"
         "chr1\t10\t.\tA\tG\t200\tstrandbias\tAF=0.98;DP=120;DP4=1,0,118,0\n"
+        "chr1\t15\t.\tT\tC\t210\tsb_fdr\tAF=0.98;DP=125;DP4=0,5,114,5\n"
+        "chr1\t18\t.\tG\tT\t220\tlowqual\tAF=0.98;DP=125;DP4=0,4,115,5\n"
         "chr1\t30\t.\tC\tT\t150\tPASS\tAF=0.40;DP=120;DP4=30,30,30,30\n"
+        "chr1\t50\t.\tT\tC\t160\tsb_fdr\tAF=0.96;DP=130;DP4=2,2,60,65\n"
     )
 
     def fake_run(cmd, check):
-        assert cmd[:4] == ["lofreq", "filter", "-i", str(raw_vcf)]
+        assert cmd[:5] == ["lofreq", "filter", "--print-all", "-i", str(raw_vcf)]
         assert check is True
         output_path = Path(cmd[cmd.index("-o") + 1])
         output_path.write_text(default_filtered_vcf, encoding="utf-8")
@@ -41,21 +50,26 @@ def fake_run(cmd, check):
 
     output_vcf = tmp_path / "final.vcf"
     rescued_tsv = tmp_path / "rescued.tsv"
+    filtered_out_tsv = tmp_path / "filtered_out.tsv"
 
     result = rescue_lofreq_primer_variants(
         raw_vcf=raw_vcf,
         primers_bed=primers,
         output_vcf=output_vcf,
         rescued_tsv=rescued_tsv,
+        filtered_out_tsv=filtered_out_tsv,
     )
 
     assert result.normal_passed == 1
-    assert result.rescued == 1
-    assert result.discarded == 0
+    assert result.rescued == 2
+    assert result.discarded == 2
+    assert result.filtered_out == 2
 
     output_text = output_vcf.read_text(encoding="utf-8")
     assert "##FILTER=<ID=RESCUED_PRIMER_OVERLAP," in output_text
     assert "chr1\t10\t.\tA\tG\t200\tRESCUED_PRIMER_OVERLAP" in output_text
+    assert "chr1\t15\t.\tT\tC\t210\tRESCUED_PRIMER_OVERLAP" in output_text
+    assert "chr1\t18\t.\tG\tT" not in output_text
     assert "PRIMER_OVERLAP;RESCUED_BY=overlap_primer_interval" in output_text
     assert "chr1\t30\t.\tC\tT\t150\tPASS" in output_text
 
@@ -64,3 +78,61 @@ def fake_run(cmd, check):
     assert rescued_lines[1].startswith(
         "A10G\tstrandbias\toverlap_primer_interval\tAF=0.98"
     )
+    assert rescued_lines[2].startswith("T15C\tsb_fdr\toverlap_primer_interval\tAF=0.98")
+    assert "minor_alt_fraction=0.0420168" in rescued_lines[2]
+
+    filtered_out_lines = filtered_out_tsv.read_text(encoding="utf-8").splitlines()
+    assert filtered_out_lines[0] == "variant\treason_filtered\tmetrics"
+    assert filtered_out_lines[1].startswith("G18T\tlowqual\tAF=0.98")
+    assert filtered_out_lines[2].startswith("T50C\tsb_fdr\tAF=0.96")
+
+
+def test_lofreq_filter_with_audit_writes_pass_vcf_and_filtered_table(
+    monkeypatch, tmp_path
+):
+    raw_vcf = tmp_path / "raw.vcf"
+    raw_vcf.write_text(
+        "##fileformat=VCFv4.2\n"
+        '##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">\n'
+        '##INFO=<ID=DP,Number=1,Type=Integer,Description="Depth">\n'
+        '##INFO=<ID=DP4,Number=4,Type=Integer,Description="Strand depths">\n'
+        "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n"
+        "chr1\t10\t.\tA\tG\t200\t.\tAF=0.98;DP=120;DP4=1,0,118,0\n"
+        "chr1\t50\t.\tT\tC\t160\t.\tAF=0.96;DP=130;DP4=2,2,60,65\n",
+        encoding="utf-8",
+    )
+    default_filtered_vcf = (
+        "##fileformat=VCFv4.2\n"
+        '##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">\n'
+        '##INFO=<ID=DP,Number=1,Type=Integer,Description="Depth">\n'
+        '##INFO=<ID=DP4,Number=4,Type=Integer,Description="Strand depths">\n'
+        "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n"
+        "chr1\t10\t.\tA\tG\t200\tPASS\tAF=0.98;DP=120;DP4=1,0,118,0\n"
+        "chr1\t50\t.\tT\tC\t160\tsb_fdr\tAF=0.96;DP=130;DP4=2,2,60,65\n"
+    )
+
+    def fake_run(cmd, check):
+        assert cmd[:5] == ["lofreq", "filter", "--print-all", "-i", str(raw_vcf)]
+        assert check is True
+        output_path = Path(cmd[cmd.index("-o") + 1])
+        output_path.write_text(default_filtered_vcf, encoding="utf-8")
+        return subprocess.CompletedProcess(cmd, 0)
+
+    monkeypatch.setattr("vartracker.lofreq_primer_rescue.subprocess.run", fake_run)
+
+    output_vcf = tmp_path / "final.vcf"
+    filtered_out_tsv = tmp_path / "filtered_out.tsv"
+    result = lofreq_filter_with_audit(
+        raw_vcf=raw_vcf,
+        output_vcf=output_vcf,
+        filtered_out_tsv=filtered_out_tsv,
+    )
+
+    assert result.normal_passed == 1
+    assert result.filtered_out == 1
+    output_text = output_vcf.read_text(encoding="utf-8")
+    assert "chr1\t10\t.\tA\tG\t200\tPASS" in output_text
+    assert "chr1\t50\t.\tT\tC" not in output_text
+    filtered_out_lines = filtered_out_tsv.read_text(encoding="utf-8").splitlines()
+    assert filtered_out_lines[0] == "variant\treason_filtered\tmetrics"
+    assert filtered_out_lines[1].startswith("T50C\tsb_fdr\tAF=0.96")

tests/test_main.py

@@ -54,12 +54,15 @@ def test_bam_parser_accepts_lofreq_primer_rescue_options():
             "off",
             "--lofreq-rescue-min-af",
             "0.9",
+            "--lofreq-rescue-max-minor-alt-fraction",
+            "0.1",
         ]
     )
 
     assert args.primer_bed == "primers.bed"
     assert args.lofreq_primer_rescue == "off"
     assert args.lofreq_rescue_min_af == 0.9
+    assert args.lofreq_rescue_max_minor_alt_fraction == 0.1
 
 
 def test_drop_exact_duplicate_result_rows_removes_only_exact_duplicates(capsys):
@@ -1072,6 +1075,7 @@ def fake_vcf(args):
     assert recorded["workflow_kwargs"]["ampliconclip_tolerance"] == 2
     assert recorded["workflow_kwargs"]["lofreq_primer_rescue"] == "auto"
     assert recorded["workflow_kwargs"]["lofreq_rescue_min_af"] == 0.95
+    assert recorded["workflow_kwargs"]["lofreq_rescue_max_minor_alt_fraction"] == 0.05
     assert recorded["vcf_input"] == str(updated_csv)
     assert modes_checked == ["e2e"]
 

tests/test_snakemake_workflow.py

@@ -33,12 +33,30 @@ def test_snakemake_rules_write_logs_under_outdir():
     assert "--tolerance {params.tolerance}" in snakefile
     assert "--no-default-filter" in snakefile
     assert "_variants.raw.vcf.gz" in snakefile
+    assert "_variants.raw.vcf.gz.tbi" in snakefile
+    assert (
+        'vcf_raw = temp(f"{OUTDIR}/{{sample}}/{{sample}}_variants.raw.vcf.gz")'
+        not in snakefile
+    )
+    assert (
+        'vcf_raw = f"{OUTDIR}/{{sample}}/{{sample}}_variants.raw.vcf.gz"' in snakefile
+    )
+    assert (
+        'vcf_raw_tbi = f"{OUTDIR}/{{sample}}/{{sample}}_variants.raw.vcf.gz.tbi"'
+        in snakefile
+    )
+    assert "_variants.filtered_out.tsv" in snakefile
     assert "LOFREQ_PRIMER_RESCUE_ENABLED" in snakefile
+    assert "LOFREQ_RESCUE_MAX_MINOR_ALT_FRACTION" in snakefile
+    assert "max_minor_alt_fraction=LOFREQ_RESCUE_MAX_MINOR_ALT_FRACTION" in snakefile
+    assert "lofreq_filter_with_audit" in snakefile
     assert "rescue_lofreq_primer_variants" in snakefile
     assert "_variants.rescued.tsv" in snakefile
     assert "_validate_primer_bed_reference(PRIMER_BED, REF)" in snakefile
     assert "_consensus.fasta" in snakefile
     assert "_iupac_consensus.fasta" in snakefile
     assert "df['consensus']" in snakefile
     assert "df['iupac_consensus']" in snakefile
+    assert "df['raw_vcf']" in snakefile
     assert "df['lofreq_rescued_tsv']" in snakefile
+    assert "df['lofreq_filtered_out_tsv']" in snakefile