Gracefully handle low-coverage runs with explicit failure output (instead of crashing) [fix #109]

src/optitype/cli.py

@@ -173,7 +173,7 @@ def run(
       # With custom settings
       optitype run -i reads.fq --dna -o results/ --solver cbc --threads 8
     """
-    from optitype.pipeline import PipelineConfig, load_config, load_mapper_config, run_pipeline
+    from optitype.pipeline import LowCoverageError, PipelineConfig, load_config, load_mapper_config, run_pipeline
 
     # Resolve razers3 path (--razers3 option or search PATH)
     if razers3 is None:
@@ -237,6 +237,9 @@ def run(
         click.echo(f"Results written to: {result.output_csv}")
         click.echo(f"Coverage plot: {result.output_plot}")
 
+    except LowCoverageError as e:
+        click.secho(f"WARNING: {e}", fg="yellow", err=True)
+        return
     except Exception as e:
         raise click.ClickException(str(e)) from None
 

src/optitype/io/readers.py

@@ -190,11 +190,13 @@ def pysam_to_dataframe(samfile: str) -> tuple[pd.DataFrame, pd.DataFrame]:
 
     logger.debug("\n %s %d reads loaded. Creating dataframe...", elapsed(), len(hits))
 
-    pos_df = pd.DataFrame.from_dict(hits, orient="index")
-    pos_df.columns = sam.references[:]
+    pos_df = pd.DataFrame.from_dict(hits, orient="index", columns=sam.references[:])
 
-    details_df = pd.DataFrame.from_dict(read_details, orient="index")
-    details_df.columns = ["mismatches", "read_length"]
+    details_df = pd.DataFrame.from_dict(
+        read_details,
+        orient="index",
+        columns=["mismatches", "read_length"],
+    )
 
     if hit_counter > 0:
         logger.debug(

src/optitype/pipeline.py

@@ -68,6 +68,29 @@ class PipelineResult:
     output_plot: str
 
 
+class LowCoverageError(RuntimeError):
+    """Raised when HLA typing cannot proceed due to insufficient mapped HLA reads."""
+
+
+def _write_low_coverage_result(out_csv: str, message: str) -> None:
+    """Write a structured failure TSV so callers get an explicit output artifact."""
+    failed = pd.DataFrame([
+        {
+            "A1": None,
+            "A2": None,
+            "B1": None,
+            "B2": None,
+            "C1": None,
+            "C2": None,
+            "Reads": 0,
+            "Objective": None,
+            "Status": "FAILED",
+            "Message": message,
+        }
+    ])
+    failed.to_csv(out_csv, sep="\t", index=False)
+
+
 def get_num_threads(configured_threads: int) -> int:
     """Get the number of threads to use, capped by available CPUs."""
     try:
@@ -324,6 +347,15 @@ def run_pipeline(
         pos2, read_details2 = pysam_to_dataframe(bam_paths[1])
         binary2 = np.sign(pos2)
 
+        if pos.empty and pos2.empty:
+            msg = (
+                "HLA typing was not possible due to low coverage: no mapped HLA reads "
+                "were found in either input read file."
+            )
+            logger.error(msg)
+            _write_low_coverage_result(out_csv, msg)
+            raise LowCoverageError(f"{msg} See {out_csv}.")
+
         if not bam_input and config.delete_bam:
             os.remove(bam_paths[0])
             os.remove(bam_paths[1])
@@ -361,15 +393,42 @@ def cut_last_char(x): return x[:-1]
     else:
         pos, read_details = pysam_to_dataframe(bam_paths[0])
 
+        if pos.empty:
+            msg = (
+                "HLA typing was not possible due to low coverage: no mapped HLA reads "
+                "were found in the input file."
+            )
+            logger.error(msg)
+            _write_low_coverage_result(out_csv, msg)
+            raise LowCoverageError(f"{msg} See {out_csv}.")
+
         if not bam_input and config.delete_bam:
             os.remove(bam_paths[0])
 
         binary = np.sign(pos)
 
+    if binary.shape[0] == 0:
+        msg = (
+            "HLA typing was not possible due to low coverage: no usable mapped HLA reads "
+            "remain after read pairing/selection."
+        )
+        logger.error(msg)
+        _write_low_coverage_result(out_csv, msg)
+        raise LowCoverageError(f"{msg} See {out_csv}.")
+
     # Filter to frequent alleles
     alleles_to_keep = [col for col in binary.columns if _is_frequent(col, table)]
     binary = binary[alleles_to_keep]
 
+    if binary.shape[1] == 0:
+        msg = (
+            "HLA typing was not possible due to low coverage: no informative HLA allele "
+            "hits remained after filtering."
+        )
+        logger.error(msg)
+        _write_low_coverage_result(out_csv, msg)
+        raise LowCoverageError(f"{msg} See {out_csv}.")
+
     logger.debug("%s Temporary pruning of identical rows and columns", elapsed())
 
     unique_col, _representing = ht.prune_identical_alleles(binary, report_groups=True)
@@ -381,6 +440,15 @@ def cut_last_char(x): return x[:-1]
 
     minimal_alleles = ht.prune_overshadowed_alleles(temp_pruned)
 
+    if len(minimal_alleles) == 0:
+        msg = (
+            "HLA typing was not possible due to low coverage: no minimal informative "
+            "alleles could be retained."
+        )
+        logger.error(msg)
+        _write_low_coverage_result(out_csv, msg)
+        raise LowCoverageError(f"{msg} See {out_csv}.")
+
     logger.debug("%s Keeping only the minimal number of required alleles %s", elapsed(), minimal_alleles.shape)
 
     binary = binary[minimal_alleles]

tests/test_cli.py

@@ -6,6 +6,7 @@
 from click.testing import CliRunner
 
 from optitype.cli import main
+from optitype.pipeline import LowCoverageError
 
 
 def test_cli_help():
@@ -75,3 +76,29 @@ def test_cli_run_yara_missing_binary(tmp_path):
         )
     assert result.exit_code != 0
     assert "yara_mapper" in result.output
+
+
+def test_cli_run_low_coverage_reports_reason_and_result_path(tmp_path):
+    """CLI should surface low-coverage reason and result TSV path to users."""
+    fq = tmp_path / "reads.fq"
+    fq.write_text("@read1\nACGT\n+\nIIII\n")
+
+    outdir = tmp_path / "out"
+    prefix = "lowcov"
+    out_csv = outdir / f"{prefix}_result.tsv"
+    out_csv.parent.mkdir(parents=True, exist_ok=True)
+    out_csv.write_text("Status\tMessage\nFAILED\tlow coverage\n")
+
+    msg = f"HLA typing was not possible due to low coverage. See {out_csv}."
+
+    runner = CliRunner()
+    with patch("optitype.pipeline.run_pipeline", side_effect=LowCoverageError(msg)):
+        result = runner.invoke(
+            main,
+            ["run", "-i", str(fq), "--dna", "-o", str(outdir), "--prefix", prefix],
+        )
+
+    assert result.exit_code == 0
+    assert "warning" in result.output.lower()
+    assert "low coverage" in result.output.lower()
+    assert str(out_csv) in result.output

tests/test_io.py

@@ -4,6 +4,7 @@
 import pandas as pd
 
 from optitype.io.data import get_data_path, load_reference_data, get_reference_fasta
+from optitype.io import readers
 
 
 def test_get_data_path():
@@ -48,3 +49,30 @@ def test_get_reference_fasta_invalid():
     """Test that invalid seq_type raises ValueError."""
     with pytest.raises(ValueError):
         get_reference_fasta("invalid")
+
+
+def test_pysam_to_dataframe_empty_alignments(monkeypatch):
+    """Empty BAM/SAM input should return empty DataFrames with valid columns."""
+
+    class DummySam:
+        header = {"PG": [{"ID": "yara", "CL": ""}]}
+        nreferences = 3
+        references = ["HLA:A*01:01", "HLA:B*07:02", "HLA:C*07:02"]
+
+        def __iter__(self):
+            return iter(())
+
+    class DummyPysam:
+        @staticmethod
+        def AlignmentFile(_samfile, _mode):
+            return DummySam()
+
+    monkeypatch.setattr(readers, "PYSAM_AVAILABLE", True)
+    monkeypatch.setattr(readers, "pysam", DummyPysam(), raising=False)
+
+    pos_df, details_df = readers.pysam_to_dataframe("empty.bam")
+
+    assert pos_df.empty
+    assert list(pos_df.columns) == ["HLA:A*01:01", "HLA:B*07:02", "HLA:C*07:02"]
+    assert details_df.empty
+    assert list(details_df.columns) == ["mismatches", "read_length"]

tests/test_pipeline.py

@@ -12,6 +12,7 @@
 from optitype.io.data import load_reference_data
 from optitype.pipeline import (
     PipelineConfig,
+    LowCoverageError,
     run_pipeline,
     load_config,
     load_mapper_config,
@@ -174,6 +175,33 @@ def test_invalid_enumerate(self, tmp_path):
         with pytest.raises(ValueError, match="enumerate"):
             run_pipeline([str(f)], "dna", str(tmp_path), enumerate_count=0)
 
+    def test_low_coverage_writes_failure_tsv(self, tmp_path):
+        out_prefix = "lowcov"
+        out_csv = tmp_path / f"{out_prefix}_result.tsv"
+
+        with patch(
+            "optitype.pipeline.pysam_to_dataframe",
+            return_value=(
+                pd.DataFrame(columns=["HLA:A*01:01"]),
+                pd.DataFrame(columns=["mismatches", "read_length"]),
+            ),
+        ):
+            with pytest.raises(LowCoverageError, match="low coverage"):
+                run_pipeline(
+                    input_files=["sample.bam"],
+                    seq_type="dna",
+                    output_dir=str(tmp_path),
+                    prefix=out_prefix,
+                    config=PipelineConfig(),
+                    verbose=False,
+                )
+
+        assert out_csv.exists()
+        failed = pd.read_csv(out_csv, sep="\t")
+        assert failed.loc[0, "Status"] == "FAILED"
+        assert "low coverage" in failed.loc[0, "Message"]
+        assert int(failed.loc[0, "Reads"]) == 0
+
 
 class TestMapperConfig:
     """Tests for mapper config loading."""