feat(utils): Ensure reference first in alignment

README.md

@@ -1,10 +1,62 @@
 # BarcodeForge
+
 [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.15654220.svg)](https://doi.org/10.5281/zenodo.15654220)
 
+A CLI tool for generating pathogen-specific barcodes for Freyja.
+
+## Installation
+
+BarcodeForge requires **Python >= 3.10** (tested on 3.10–3.13).
+
+### Conda (recommended)
 
-https://anaconda.org/bioconda/barcodeforge
+BarcodeForge is available on [Bioconda](https://anaconda.org/bioconda/barcodeforge):
 
-To install run:
 ```
 conda install bioconda::barcodeforge
 ```
+
+### From source (GitHub)
+
+Install the latest version directly from the GitHub repository:
+
+```
+pip install git+https://github.com/gp201/BarcodeForge.git
+```
+
+Or clone the repository and install it locally:
+
+```
+git clone https://github.com/gp201/BarcodeForge.git
+cd BarcodeForge
+pip install .
+```
+
+For development (includes test and formatting tools):
+
+```
+pip install -e ".[dev-requirements,test-requirements]"
+```
+
+### Additional dependencies
+
+The following tools are required and must be installed separately via Conda
+(from the `bioconda` and `conda-forge` channels):
+
+```
+conda install -c conda-forge -c bioconda usher ucsc-fatovcf
+```
+
+## Usage
+
+Once installed, run the tool with:
+
+```
+barcodeforge --help
+```
+
+For example, to generate a barcode:
+
+```
+barcodeforge barcode reference.fasta aligned.fasta tree.nwk lineages.tsv --threads 4 --prefix RSVa
+```

barcodeforge/cli.py

@@ -4,6 +4,7 @@
 from .utils import (
     resolve_tree_format,
     run_subprocess_command,
+    ensure_reference_is_first_in_alignment,
     print_error,
     print_success,
     print_info,
@@ -109,6 +110,16 @@ def barcode(
     intermediate_dir = "barcodeforge_workdir"
     os.makedirs(intermediate_dir, exist_ok=True)
 
+    # faToVcf uses the first sequence in the alignment as the reference for calling
+    # variants. If the reference genome is not already first, prepend it and use the
+    # corrected alignment (written to the workdir) for the rest of the pipeline.
+    alignment = ensure_reference_is_first_in_alignment(
+        reference_genome,
+        alignment,
+        os.path.join(intermediate_dir, "alignment_with_reference.fasta"),
+        is_debug,
+    )
+
     # Run faToVcf command using utility function
     fatovcf_output_vcf = os.path.join(intermediate_dir, "aligned.vcf")
     fatovcf_cmd = ["faToVcf", alignment, fatovcf_output_vcf]

barcodeforge/utils.py

@@ -2,6 +2,7 @@
 import subprocess
 import rich_click as click
 from rich.console import Console
+from Bio import SeqIO
 
 console = Console()
 
@@ -71,6 +72,75 @@ def resolve_tree_format(
     return resolved_format
 
 
+def ensure_reference_is_first_in_alignment(
+    reference_genome_path: str,
+    alignment_path: str,
+    output_alignment_path: str,
+    debug: bool,
+) -> str:
+    """
+    Ensures the reference genome is the first sequence in the alignment.
+
+    faToVcf treats the first sequence in the alignment as the reference. If the
+    reference genome is not already the first sequence, a corrected alignment is
+    written to ``output_alignment_path`` with the reference genome prepended (and
+    any existing record sharing the reference id removed to avoid duplicates), and
+    that path is returned. If the reference is already first, the original
+    alignment path is returned unchanged.
+
+    Args:
+        reference_genome_path (str): Path to the reference genome FASTA file.
+        alignment_path (str): Path to the alignment FASTA file.
+        output_alignment_path (str): Path to write the corrected alignment to if
+            the reference genome needs to be added.
+        debug (bool): If True, prints debug information.
+    Returns:
+        str: Path to the alignment to use downstream (the original path if no
+            change was needed, otherwise ``output_alignment_path``).
+    Raises:
+        click.Abort: If the alignment is empty, or the reference genome length
+            does not match the aligned sequence length.
+    """
+    ref = SeqIO.read(reference_genome_path, "fasta")
+    records = list(SeqIO.parse(alignment_path, "fasta"))
+
+    if not records:
+        print_error(f"Alignment file '{alignment_path}' contains no sequences.")
+        raise click.Abort()
+
+    if records[0].id == ref.id:
+        if debug:
+            print_debug(
+                f"Reference genome '{ref.id}' is already the first sequence in the alignment."
+            )
+        return alignment_path
+
+    print_warning(
+        f"The first sequence in the alignment ('{records[0].id}') is not the "
+        f"reference genome ('{ref.id}'). Prepending the reference genome to the alignment."
+    )
+
+    alignment_width = len(records[0].seq)
+    if len(ref.seq) != alignment_width:
+        print_error(
+            f"Cannot add the reference genome to the alignment: the reference length "
+            f"({len(ref.seq)}) does not match the aligned sequence length "
+            f"({alignment_width}). faToVcf requires all sequences to be the same length. "
+            f"Please provide a reference that is aligned to the same coordinates as the alignment."
+        )
+        raise click.Abort()
+
+    # Drop any existing record that shares the reference id to avoid duplicate
+    # sample names, then place the reference genome first.
+    corrected = [ref] + [rec for rec in records if rec.id != ref.id]
+    SeqIO.write(corrected, output_alignment_path, "fasta")
+    print_success(
+        f"Wrote corrected alignment with reference genome '{ref.id}' as the first "
+        f"sequence to {output_alignment_path}"
+    )
+    return output_alignment_path
+
+
 def run_subprocess_command(
     cmd: list[str],
     debug: bool,

tests/test_cli.py

@@ -32,9 +32,11 @@ def temp_files(tmp_path):
     lineages = tmp_path / "lineages.tsv"
 
     create_dummy_file(ref_genome, ">ref_genome_id\nAGCTAGCTAGCTAGCT")
+    # The reference genome must be the first sequence in the alignment, since
+    # faToVcf treats the first sequence as the reference.
     create_dummy_file(
         alignment,
-        ">seq1\nAGCTAGCTAGCTAGCT\n>seq2\nAGCTAGCTAGCTCGCT\n>seq3\nAGCTAGCTAGCTAGGT",
+        ">ref_genome_id\nAGCTAGCTAGCTAGCT\n>seq1\nAGCTAGCTAGCTAGCT\n>seq2\nAGCTAGCTAGCTCGCT\n>seq3\nAGCTAGCTAGCTAGGT",
     )
     create_dummy_file(tree, "((seq1:0.1,seq2:0.1):0.05,seq3:0.15);")
     create_dummy_file(lineages, "clade\tsequences\nlineageA\tseq1,seq2\nlineageB\tseq3")
@@ -537,6 +539,60 @@ def test_barcode_command_missing_file(runner, temp_files):
     assert "Invalid value for 'ALIGNMENT'" in result.output
 
 
+def test_barcode_command_reference_not_first_in_alignment(runner, temp_files, mocker):
+    # Alignment whose first sequence is NOT the reference genome: the reference
+    # should be prepended and the corrected alignment used for the rest of the run.
+    create_dummy_file(
+        temp_files["alignment"],
+        ">seq1\nAGCTAGCTAGCTAGCT\n>seq2\nAGCTAGCTAGCTCGCT\n>seq3\nAGCTAGCTAGCTAGGT",
+    )
+    mock_run_subp = mocker.patch(
+        "barcodeforge.cli.run_subprocess_command", return_value=True
+    )
+    mocker.patch("barcodeforge.cli.resolve_tree_format", return_value="newick")
+    mocker.patch("barcodeforge.cli.convert_nexus_to_newick")
+    mock_process_reroot = mocker.patch("barcodeforge.cli.process_and_reroot_lineages")
+    mocker.patch("barcodeforge.cli.create_barcodes_from_lineage_paths")
+    mocker.patch("barcodeforge.cli.create_barcode_plot")
+    mocker.patch.object(barcodeforge.utils, "console", MagicMock(spec=Console))
+
+    intermediate_dir = "barcodeforge_workdir"
+    corrected_alignment_fn = f"{intermediate_dir}/alignment_with_reference.fasta"
+    aligned_vcf_fn = f"{intermediate_dir}/aligned.vcf"
+
+    args = [
+        "barcode",
+        temp_files["ref_genome"],
+        temp_files["alignment"],
+        temp_files["tree"],
+        temp_files["lineages"],
+    ]
+    try:
+        result = runner.invoke(cli, args, catch_exceptions=False)
+        assert result.exit_code == 0, f"CLI failed: {result.output}"
+
+        # faToVcf and rerooting should use the corrected alignment, not the input.
+        mock_run_subp.assert_any_call(
+            ["faToVcf", corrected_alignment_fn, aligned_vcf_fn],
+            False,
+            success_message=ANY,
+            error_message_prefix=ANY,
+        )
+        assert (
+            mock_process_reroot.call_args.kwargs["sequences_fasta_path"]
+            == corrected_alignment_fn
+        )
+
+        # The corrected alignment should have the reference genome as its first record.
+        from Bio import SeqIO
+
+        records = list(SeqIO.parse(corrected_alignment_fn, "fasta"))
+        assert records[0].id == "ref_genome_id"
+        assert [r.id for r in records] == ["ref_genome_id", "seq1", "seq2", "seq3"]
+    finally:
+        shutil.rmtree(intermediate_dir, ignore_errors=True)
+
+
 def test_extract_auspice_data_command(runner, tmp_path):
     json_path = tmp_path / "tree.json"
     json_data = {

tests/test_utils.py

@@ -1,15 +1,91 @@
 import pytest
 import click
 import subprocess
+from pathlib import Path
 from unittest.mock import MagicMock, patch, call
 from rich.console import Console
 import barcodeforge.utils
 from barcodeforge.utils import (
     sortFun,
     resolve_tree_format,
     run_subprocess_command,
+    ensure_reference_is_first_in_alignment,
     STYLES,
 )
+from Bio import SeqIO
+
+
+def _write_fasta(path, content):
+    with open(path, "w") as f:
+        f.write(content)
+    return str(path)
+
+
+def test_ensure_reference_first_already_first(tmp_path):
+    ref = _write_fasta(tmp_path / "ref.fasta", ">ref\nACGT")
+    aln = _write_fasta(tmp_path / "aln.fasta", ">ref\nACGT\n>seq1\nACGA")
+    out = str(tmp_path / "out.fasta")
+    # Reference already first: returns the original path and writes nothing.
+    result = ensure_reference_is_first_in_alignment(ref, aln, out, debug=False)
+    assert result == aln
+    assert not Path(out).exists()
+
+
+def test_ensure_reference_first_prepends_when_absent(tmp_path):
+    ref = _write_fasta(tmp_path / "ref.fasta", ">ref\nACGT")
+    aln = _write_fasta(tmp_path / "aln.fasta", ">seq1\nACGA\n>seq2\nACGC")
+    out = str(tmp_path / "out.fasta")
+    result = ensure_reference_is_first_in_alignment(ref, aln, out, debug=False)
+    assert result == out
+    records = list(SeqIO.parse(out, "fasta"))
+    assert [r.id for r in records] == ["ref", "seq1", "seq2"]
+    assert str(records[0].seq) == "ACGT"
+
+
+def test_ensure_reference_first_dedupes_existing_reference(tmp_path):
+    ref = _write_fasta(tmp_path / "ref.fasta", ">ref\nACGT")
+    # Reference present but not first; should be moved to the front, not duplicated.
+    aln = _write_fasta(tmp_path / "aln.fasta", ">seq1\nACGA\n>ref\nACGT")
+    out = str(tmp_path / "out.fasta")
+    result = ensure_reference_is_first_in_alignment(ref, aln, out, debug=False)
+    assert result == out
+    records = list(SeqIO.parse(out, "fasta"))
+    assert [r.id for r in records] == ["ref", "seq1"]
+
+
+def test_ensure_reference_first_length_mismatch_aborts(tmp_path):
+    ref = _write_fasta(tmp_path / "ref.fasta", ">ref\nACG")
+    aln = _write_fasta(tmp_path / "aln.fasta", ">seq1\nACGA\n>seq2\nACGC")
+    out = str(tmp_path / "out.fasta")
+    mock_console = MagicMock(spec=Console)
+    with patch.object(barcodeforge.utils, "console", mock_console):
+        with pytest.raises(click.Abort):
+            ensure_reference_is_first_in_alignment(ref, aln, out, debug=False)
+    assert not Path(out).exists()
+    printed = " ".join(c.args[0] for c in mock_console.print.call_args_list)
+    assert "does not match" in printed
+
+
+def test_ensure_reference_first_empty_alignment_aborts(tmp_path):
+    ref = _write_fasta(tmp_path / "ref.fasta", ">ref\nACGT")
+    aln = _write_fasta(tmp_path / "aln.fasta", "")
+    out = str(tmp_path / "out.fasta")
+    with pytest.raises(click.Abort):
+        ensure_reference_is_first_in_alignment(ref, aln, out, debug=False)
+
+
+def test_ensure_reference_first_debug_output(tmp_path):
+    ref = _write_fasta(tmp_path / "ref.fasta", ">ref\nACGT")
+    aln = _write_fasta(tmp_path / "aln.fasta", ">ref\nACGT\n>seq1\nACGA")
+    out = str(tmp_path / "out.fasta")
+    mock_console = MagicMock(spec=Console)
+    with patch.object(barcodeforge.utils, "console", mock_console):
+        ensure_reference_is_first_in_alignment(ref, aln, out, debug=True)
+    mock_console.print.assert_any_call(
+        "[DEBUG] Reference genome 'ref' is already the first sequence in the alignment.",
+        style=STYLES["debug"],
+        markup=False,
+    )
 
 
 def test_sortFun():