API Reference

This section contains the complete API reference for SSIAMB.

CLI Module

ssiamb - SSI Ambiguous Site Detection Tool

Command-line interface for detecting ambiguous sites in bacterial genomes through mapping and variant calling approaches.

main(version=None, config=None, verbose=False, quiet=False, dry_run=False)

SSI Ambiguous Site Detection Tool.

Detect ambiguous sites in bacterial genomes through mapping and variant calling.

Source code in src/ssiamb/cli.py

@app.callback()
def main(
    version: Annotated[
        Optional[bool],
        typer.Option(
            "--version",
            "-V",
            help="Show version and exit",
            callback=version_callback,
            is_eager=True,
        ),
    ] = None,
    config: Annotated[
        Optional[Path],
        typer.Option(
            "--config",
            "-c",
            help="Path to custom configuration file",
            exists=True,
            file_okay=True,
            dir_okay=False,
            readable=True,
        ),
    ] = None,
    verbose: Annotated[
        bool,
        typer.Option(
            "--verbose", "-v", help="Enable verbose output"
        ),
    ] = False,
    quiet: Annotated[
        bool,
        typer.Option(
            "--quiet", "-q", help="Suppress non-error output"
        ),
    ] = False,
    dry_run: Annotated[
        bool,
        typer.Option(
            "--dry-run", help="Show what would be done without executing"
        ),
    ] = False,
) -> None:
    """
    SSI Ambiguous Site Detection Tool.

    Detect ambiguous sites in bacterial genomes through mapping and variant calling.
    """
    if verbose and quiet:
        console.print("[red]Error: --verbose and --quiet cannot be used together[/red]")
        raise typer.Exit(1)

    # Load configuration if specified
    if config:
        from .config import load_config
        try:
            load_config(config)
            if verbose:
                console.print(f"[green]Loaded configuration from: {config}[/green]")
        except Exception as e:
            console.print(f"[red]Error loading config file {config}: {e}[/red]")
            raise typer.Exit(1)

ref(ctx, r1, r2, reference=None, species=None, bracken=None, sample=None, output_dir=None, threads=4, mapper='minimap2', caller='bbtools', bbtools_mem=None, dp_min=10, maf_min=0.1, mapq=30, dp_cap=100, depth_tool=DepthTool.MOSDEPTH, require_pass=False, min_bracken_frac=0.7, min_bracken_reads=100000, ref_dir=None, on_fail='error', emit_vcf=False, emit_bed=False, emit_matrix=False, emit_per_contig=False, emit_multiqc=False, emit_provenance=False, tsv_mode=TSVMode.OVERWRITE, stdout=False)

Reference-mapping mode: map reads to a reference genome.

Maps paired-end reads to a reference genome to identify ambiguous sites. Reference can be provided directly, looked up by species, or selected from Bracken.

Source code in src/ssiamb/cli.py

@app.command()
def ref(
    ctx: typer.Context,
    r1: Annotated[
        Path,
        typer.Option(
            "--r1", help="Forward reads (FASTQ, can be gzipped)"
        ),
    ],
    r2: Annotated[
        Path,
        typer.Option(
            "--r2", help="Reverse reads (FASTQ, can be gzipped)"
        ),
    ],
    reference: Annotated[
        Optional[Path],
        typer.Option(
            "--reference", help="Reference genome FASTA file"
        ),
    ] = None,
    species: Annotated[
        Optional[str],
        typer.Option(
            "--species", help="Species name for reference lookup"
        ),
    ] = None,
    bracken: Annotated[
        Optional[Path],
        typer.Option(
            "--bracken", help="Bracken classification file"
        ),
    ] = None,
    sample: Annotated[
        Optional[str],
        typer.Option(
            "--sample", help="Sample name (inferred from filenames if not provided)"
        ),
    ] = None,
    output_dir: Annotated[
        Optional[Path],
        typer.Option(
            "--outdir", "-o", help="Output directory (default: current)"
        ),
    ] = None,
    threads: Annotated[
        int,
        typer.Option(
            "--threads", "-t", help="Number of threads"
        ),
    ] = 4,
    mapper: Annotated[
        str,
        typer.Option(
            "--mapper", help="Mapper to use"
        ),
    ] = "minimap2",
    caller: Annotated[
        str,
        typer.Option(
            "--caller", help="Variant caller to use"
        ),
    ] = "bbtools",
    bbtools_mem: Annotated[
        Optional[str],
        typer.Option(
            "--bbtools-mem", help="BBTools heap memory (e.g., '4g', '8g')"
        ),
    ] = None,
    dp_min: Annotated[
        int,
        typer.Option(
            "--dp-min", help="Minimum depth for ambiguous sites"
        ),
    ] = 10,
    maf_min: Annotated[
        float,
        typer.Option(
            "--maf-min", help="Minimum minor allele frequency"
        ),
    ] = 0.1,
    mapq: Annotated[
        int,
        typer.Option(
            "--mapq", help="Minimum mapping quality for depth analysis"
        ),
    ] = 30,
    dp_cap: Annotated[
        int,
        typer.Option(
            "--dp-cap", help="Maximum depth cap for variant analysis"
        ),
    ] = 100,
    depth_tool: Annotated[
        DepthTool,
        typer.Option(
            "--depth-tool", help="Tool for depth analysis"
        ),
    ] = DepthTool.MOSDEPTH,
    require_pass: Annotated[
        bool,
        typer.Option(
            "--require-pass", help="Only consider variants that pass caller filters"
        ),
    ] = False,
    min_bracken_frac: Annotated[
        float,
        typer.Option(
            "--min-bracken-frac", help="Minimum Bracken abundance fraction"
        ),
    ] = 0.70,
    min_bracken_reads: Annotated[
        int,
        typer.Option(
            "--min-bracken-reads", help="Minimum Bracken read count"
        ),
    ] = 100000,
    ref_dir: Annotated[
        Optional[Path],
        typer.Option(
            "--ref-dir", help="Admin reference directory"
        ),
    ] = None,
    on_fail: Annotated[
        str,
        typer.Option(
            "--on-fail", help="Action when reference resolution fails"
        ),
    ] = "error",
    emit_vcf: Annotated[
        bool,
        typer.Option(
            "--emit-vcf", help="Emit VCF file with ambiguous sites"
        ),
    ] = False,
    emit_bed: Annotated[
        bool,
        typer.Option(
            "--emit-bed", help="Emit BED file with ambiguous sites"
        ),
    ] = False,
    emit_matrix: Annotated[
        bool,
        typer.Option(
            "--emit-matrix", help="Emit variant matrix in TSV format"
        ),
    ] = False,
    emit_per_contig: Annotated[
        bool,
        typer.Option(
            "--emit-per-contig", help="Emit per-contig summary statistics"
        ),
    ] = False,
    emit_multiqc: Annotated[
        bool,
        typer.Option(
            "--emit-multiqc", help="Emit MultiQC-compatible metrics"
        ),
    ] = False,
    emit_provenance: Annotated[
        bool,
        typer.Option(
            "--emit-provenance", help="Emit JSON provenance file"
        ),
    ] = False,
    tsv_mode: Annotated[
        TSVMode,
        typer.Option(
            "--tsv-mode", help="TSV output mode"
        ),
    ] = TSVMode.OVERWRITE,
    stdout: Annotated[
        bool,
        typer.Option(
            "--stdout", help="Write summary to stdout instead of file"
        ),
    ] = False,
) -> None:
    """
    Reference-mapping mode: map reads to a reference genome.

    Maps paired-end reads to a reference genome to identify ambiguous sites.
    Reference can be provided directly, looked up by species, or selected from Bracken.
    """
    try:
        # Validate emit flags with stdout
        if stdout and any([emit_vcf, emit_bed, emit_matrix, emit_per_contig, emit_multiqc, emit_provenance]):
            console.print("[red]Error: --stdout cannot be used with --emit-* flags[/red]")
            raise typer.Exit(1)

        # Get global options from context
        dry_run = ctx.parent.params.get("dry_run", False) if ctx.parent else False

        # Create execution plan
        plan = create_run_plan(
            mode=Mode.REF,
            r1=r1,
            r2=r2,
            reference=reference,
            sample=sample,
            output_dir=output_dir,
            threads=threads,
            mapper=mapper,
            caller=caller,
            bbtools_mem=bbtools_mem,
            dp_min=dp_min,
            maf_min=maf_min,
            dp_cap=dp_cap,
            mapq=mapq,
            depth_tool=depth_tool.value,
            require_pass=require_pass,
            dry_run=dry_run,
            to_stdout=stdout,
            emit_vcf=emit_vcf,
            emit_bed=emit_bed,
            emit_matrix=emit_matrix,
            emit_per_contig=emit_per_contig,
            emit_multiqc=emit_multiqc,
            emit_provenance=emit_provenance,
            species=species,
            bracken=bracken,
            ref_dir=ref_dir,
            on_fail=on_fail,
            tsv_mode=tsv_mode,
            min_bracken_frac=min_bracken_frac,
            min_bracken_reads=min_bracken_reads,
        )

        # Execute plan
        result, provenance_record = execute_plan(plan, 
                            species=species, 
                            bracken=bracken, 
                            ref_dir=ref_dir, 
                            on_fail=on_fail,
                            min_bracken_frac=min_bracken_frac,
                            min_bracken_reads=min_bracken_reads)

        # If this was a dry run, we're done
        if dry_run:
            return

        # Handle provenance output
        if emit_provenance and provenance_record:
            from .provenance import write_provenance_json
            out_dir = output_dir or Path.cwd()
            provenance_path = out_dir / "run_provenance.json"
            write_provenance_json([provenance_record], provenance_path)
            console.print(f"Provenance written to {provenance_path}")

        if not stdout:
            console.print(f"[green]Reference-mapping completed for sample {result.sample}[/green]")
            total_sites = result.ambiguous_snv_count + result.ambiguous_indel_count + result.ambiguous_del_count
            console.print(f"Found {total_sites} ambiguous sites")

    except Exception as e:
        handle_exception_with_exit(e, "Reference-mapping mode failed")

self(ctx, r1, r2, assembly, sample=None, output_dir=None, threads=4, mapper='minimap2', caller='bbtools', bbtools_mem=None, dp_min=10, maf_min=0.1, mapq=30, dp_cap=100, depth_tool=DepthTool.MOSDEPTH, require_pass=False, emit_vcf=False, emit_bed=False, emit_matrix=False, emit_per_contig=False, emit_multiqc=False, emit_provenance=False, tsv_mode=TSVMode.OVERWRITE, stdout=False)

Self-mapping mode: map reads to their own assembly.

Maps paired-end reads to the provided assembly to identify ambiguous sites where the assembly may not represent the true sequence.

Source code in src/ssiamb/cli.py

@app.command()
def self(
    ctx: typer.Context,
    r1: Annotated[
        Path,
        typer.Option(
            "--r1", help="Forward reads (FASTQ, can be gzipped)"
        ),
    ],
    r2: Annotated[
        Path,
        typer.Option(
            "--r2", help="Reverse reads (FASTQ, can be gzipped)"
        ),
    ],
    assembly: Annotated[
        Path,
        typer.Option(
            "--assembly", help="Assembly FASTA file"
        ),
    ],
    sample: Annotated[
        Optional[str],
        typer.Option(
            "--sample", help="Sample name (inferred from filenames if not provided)"
        ),
    ] = None,
    output_dir: Annotated[
        Optional[Path],
        typer.Option(
            "--outdir", "-o", help="Output directory (default: current)"
        ),
    ] = None,
    threads: Annotated[
        int,
        typer.Option(
            "--threads", "-t", help="Number of threads"
        ),
    ] = 4,
    mapper: Annotated[
        str,
        typer.Option(
            "--mapper", help="Mapper to use"
        ),
    ] = "minimap2",
    caller: Annotated[
        str,
        typer.Option(
            "--caller", help="Variant caller to use"
        ),
    ] = "bbtools",
    bbtools_mem: Annotated[
        Optional[str],
        typer.Option(
            "--bbtools-mem", help="BBTools heap memory (e.g., '4g', '8g')"
        ),
    ] = None,
    dp_min: Annotated[
        int,
        typer.Option(
            "--dp-min", help="Minimum depth for ambiguous sites"
        ),
    ] = 10,
    maf_min: Annotated[
        float,
        typer.Option(
            "--maf-min", help="Minimum minor allele frequency"
        ),
    ] = 0.1,
    mapq: Annotated[
        int,
        typer.Option(
            "--mapq", help="Minimum mapping quality for depth analysis"
        ),
    ] = 30,
    dp_cap: Annotated[
        int,
        typer.Option(
            "--dp-cap", help="Maximum depth cap for variant analysis"
        ),
    ] = 100,
    depth_tool: Annotated[
        DepthTool,
        typer.Option(
            "--depth-tool", help="Tool for depth analysis"
        ),
    ] = DepthTool.MOSDEPTH,
    require_pass: Annotated[
        bool,
        typer.Option(
            "--require-pass", help="Only consider variants that pass caller filters"
        ),
    ] = False,
    emit_vcf: Annotated[
        bool,
        typer.Option(
            "--emit-vcf", help="Emit VCF file with ambiguous sites"
        ),
    ] = False,
    emit_bed: Annotated[
        bool,
        typer.Option(
            "--emit-bed", help="Emit BED file with ambiguous sites"
        ),
    ] = False,
    emit_matrix: Annotated[
        bool,
        typer.Option(
            "--emit-matrix", help="Emit variant matrix in TSV format"
        ),
    ] = False,
    emit_per_contig: Annotated[
        bool,
        typer.Option(
            "--emit-per-contig", help="Emit per-contig summary statistics"
        ),
    ] = False,
    emit_multiqc: Annotated[
        bool,
        typer.Option(
            "--emit-multiqc", help="Emit MultiQC-compatible metrics"
        ),
    ] = False,
    emit_provenance: Annotated[
        bool,
        typer.Option(
            "--emit-provenance", help="Emit JSON provenance file"
        ),
    ] = False,
    tsv_mode: Annotated[
        TSVMode,
        typer.Option(
            "--tsv-mode", help="TSV output mode"
        ),
    ] = TSVMode.OVERWRITE,
    stdout: Annotated[
        bool,
        typer.Option(
            "--stdout", help="Write summary to stdout instead of file"
        ),
    ] = False,
) -> None:
    """
    Self-mapping mode: map reads to their own assembly.

    Maps paired-end reads to the provided assembly to identify ambiguous sites
    where the assembly may not represent the true sequence.
    """
    try:
        # Validate emit flags with stdout
        if stdout and any([emit_vcf, emit_bed, emit_matrix, emit_per_contig, emit_multiqc, emit_provenance]):
            console.print("[red]Error: --stdout cannot be used with --emit-* flags[/red]")
            raise typer.Exit(1)

        # Get global options from context
        dry_run = ctx.parent.params.get("dry_run", False) if ctx.parent else False

        # Create execution plan
        plan = create_run_plan(
            mode=Mode.SELF,
            r1=r1,
            r2=r2,
            assembly=assembly,
            sample=sample,
            output_dir=output_dir,
            threads=threads,
            mapper=mapper,
            caller=caller,
            bbtools_mem=bbtools_mem,
            dp_min=dp_min,
            maf_min=maf_min,
            dp_cap=dp_cap,
            mapq=mapq,
            depth_tool=depth_tool.value,
            require_pass=require_pass,
            dry_run=dry_run,
            to_stdout=stdout,
            emit_vcf=emit_vcf,
            emit_bed=emit_bed,
            emit_matrix=emit_matrix,
            emit_per_contig=emit_per_contig,
            emit_multiqc=emit_multiqc,
            emit_provenance=emit_provenance,
            tsv_mode=tsv_mode,
        )

        # Execute plan
        result, provenance_record = execute_plan(plan)

        # If this was a dry run, we're done
        if dry_run:
            return

        # Handle provenance output
        if emit_provenance and provenance_record:
            from .provenance import write_provenance_json
            out_dir = output_dir or Path.cwd()
            provenance_path = out_dir / "run_provenance.json"
            write_provenance_json([provenance_record], provenance_path)
            console.print(f"Provenance written to {provenance_path}")

        if not stdout:
            console.print(f"[green]Self-mapping completed for sample {result.sample}[/green]")
            total_sites = result.ambiguous_snv_count + result.ambiguous_indel_count + result.ambiguous_del_count
            console.print(f"Found {total_sites} ambiguous sites")

    except Exception as e:
        handle_exception_with_exit(e, "Self-mapping mode failed")

summarize(ctx, vcf, bam, output=None, dp_min=10, maf_min=0.1, dp_cap=100, require_pass=False, emit_vcf=False, emit_bed=False, emit_matrix=False, emit_per_contig=False, emit_multiqc=False, emit_provenance=False, mode='combined', stdout=False)

Summarize VCF and BAM files to generate ambiguous site summary.

Analyzes a VCF file with BAM for denominator calculation to produce ambiguous site statistics.

Source code in src/ssiamb/cli.py

@app.command()
def summarize(
    ctx: typer.Context,
    vcf: Annotated[
        Path,
        typer.Option(
            "--vcf", help="VCF file to summarize"
        ),
    ],
    bam: Annotated[
        Path,
        typer.Option(
            "--bam", help="BAM file for denominator calculation"
        ),
    ],
    output: Annotated[
        Optional[Path],
        typer.Option(
            "--output", "-o", help="Output summary file"
        ),
    ] = None,
    dp_min: Annotated[
        int,
        typer.Option(
            "--dp-min", help="Minimum depth for ambiguous sites"
        ),
    ] = 10,
    maf_min: Annotated[
        float,
        typer.Option(
            "--maf-min", help="Minimum minor allele frequency"
        ),
    ] = 0.1,
    dp_cap: Annotated[
        int,
        typer.Option(
            "--dp-cap", help="Maximum depth cap for variant analysis"
        ),
    ] = 100,
    require_pass: Annotated[
        bool,
        typer.Option(
            "--require-pass", help="Only consider variants that pass caller filters"
        ),
    ] = False,
    emit_vcf: Annotated[
        bool,
        typer.Option(
            "--emit-vcf", help="Emit VCF file with ambiguous sites"
        ),
    ] = False,
    emit_bed: Annotated[
        bool,
        typer.Option(
            "--emit-bed", help="Emit BED file with ambiguous sites"
        ),
    ] = False,
    emit_matrix: Annotated[
        bool,
        typer.Option(
            "--emit-matrix", help="Emit variant matrix in TSV format"
        ),
    ] = False,
    emit_per_contig: Annotated[
        bool,
        typer.Option(
            "--emit-per-contig", help="Emit per-contig summary statistics"
        ),
    ] = False,
    emit_multiqc: Annotated[
        bool,
        typer.Option(
            "--emit-multiqc", help="Emit MultiQC-compatible metrics"
        ),
    ] = False,
    emit_provenance: Annotated[
        bool,
        typer.Option(
            "--emit-provenance", help="Emit JSON provenance file"
        ),
    ] = False,
    mode: Annotated[
        str,
        typer.Option(
            "--mode", help="Summary mode"
        ),
    ] = "combined",
    stdout: Annotated[
        bool,
        typer.Option(
            "--stdout", help="Write summary to stdout instead of file"
        ),
    ] = False,
) -> None:
    """
    Summarize VCF and BAM files to generate ambiguous site summary.

    Analyzes a VCF file with BAM for denominator calculation to produce
    ambiguous site statistics.
    """
    try:
        # Get global options from context
        dry_run = ctx.parent.params.get("dry_run", False) if ctx.parent else False

        if dry_run:
            console.print("[yellow]DRY RUN - Summarize mode plan:[/yellow]")
            console.print(f"  [bold]Input validation:[/bold]")
            console.print(f"    VCF file: {vcf} {'✓' if vcf.exists() else '✗ (missing)'}")
            console.print(f"    BAM file: {bam} {'✓' if bam.exists() else '✗ (missing)'}")
            console.print(f"  [bold]Sample name:[/bold] {vcf.stem.split('.')[0] if vcf.exists() else 'unknown'}")
            console.print(f"  [bold]Analysis plan:[/bold]")
            console.print(f"    1. Run mosdepth on BAM (MAPQ≥30, depth≥10, exclude duplicates)")
            console.print(f"    2. Parse VCF for ambiguous sites")
            console.print(f"    3. Count SNVs: dp_min={dp_min}, maf_min={maf_min}, dp_cap={dp_cap}")
            console.print(f"    4. Count indels and deletions for secondary metrics")
            console.print(f"    5. Calculate mapping rate from BAM")
            console.print(f"  [bold]Outputs planned:[/bold]")
            if stdout:
                console.print(f"    Summary: stdout")
            else:
                console.print(f"    Summary: {output or 'ambiguous_summary.tsv'}")
            if emit_vcf:
                console.print(f"    VCF: {vcf.stem}.ambiguous_sites.vcf.gz")
            if emit_bed:
                console.print(f"    BED: {vcf.stem}.ambiguous_sites.bed.gz")
            if emit_matrix:
                console.print(f"    Matrix: {vcf.stem}.variant_matrix.tsv.gz")
            if emit_per_contig:
                console.print(f"    Per-contig: {vcf.stem}.per_contig_summary.tsv")
            if emit_multiqc:
                console.print(f"    MultiQC: {vcf.stem}.multiqc.tsv")
            if emit_provenance:
                console.print(f"    Provenance: run_provenance.json")
            console.print(f"  [bold]Filters:[/bold] {'PASS-only' if require_pass else 'All variants'}")
            console.print("[green]Dry run completed - no files written[/green]")
            return

        console.print("[green]Running summarize mode...[/green]")
        console.print(f"VCF file: {vcf}")
        console.print(f"BAM file: {bam}")
        console.print(f"Output: {output or 'stdout' if stdout else 'auto'}")
        console.print(f"Mode: {mode}")
        console.print(f"Thresholds: dp_min={dp_min}, maf_min={maf_min}")

        # Call the summarize logic
        results = run_summarize(
            vcf=vcf,
            bam=bam,
            output=output,
            dp_min=dp_min,
            maf_min=maf_min,
            dp_cap=dp_cap,
            require_pass=require_pass,
            emit_vcf=emit_vcf,
            emit_bed=emit_bed,
            emit_matrix=emit_matrix,
            emit_per_contig=emit_per_contig,
            emit_multiqc=emit_multiqc,
            emit_provenance=emit_provenance,
            to_stdout=stdout,
        )

        if not stdout:
            console.print(f"[green]Summarize completed with {len(results)} results[/green]")

    except Exception as e:
        handle_exception_with_exit(e, "Summarize mode failed")

version_callback(value)

Print version and exit.

Source code in src/ssiamb/cli.py

def version_callback(value: bool) -> None:
    """Print version and exit."""
    if value:
        console.print(f"ssiamb version {__version__}")
        raise typer.Exit()

Configuration

Configuration management for ssiamb.

This module handles loading and merging configuration from: 1. Built-in defaults (config/defaults.yaml) 2. User-specified config files (--config) 3. Environment variables 4. Command-line overrides

SsiambConfig(thresholds, species_aliases, tools, output) dataclass

Complete ssiamb configuration.

This holds all configurable values that were previously hardcoded, allowing users to customize behavior via config files.

get_output_setting(key, default=None)

Get an output formatting setting.

Source code in src/ssiamb/config.py

def get_output_setting(self, key: str, default: Any = None) -> Any:
    """Get an output formatting setting."""
    return self.output.get(key, default)

get_species_alias(species)

Get species alias, returning original name if no alias exists.

Source code in src/ssiamb/config.py

def get_species_alias(self, species: str) -> str:
    """Get species alias, returning original name if no alias exists."""
    return self.species_aliases.get(species, species)

get_threshold(key, default=None)

Get a threshold value with fallback.

Source code in src/ssiamb/config.py

def get_threshold(self, key: str, default: Any = None) -> Any:
    """Get a threshold value with fallback."""
    return self.thresholds.get(key, default)

get_tool_setting(tool, key, default=None)

Get a tool-specific setting.

Source code in src/ssiamb/config.py

def get_tool_setting(self, tool: str, key: str, default: Any = None) -> Any:
    """Get a tool-specific setting."""
    return self.tools.get(tool, {}).get(key, default)

load(config_path=None) classmethod

Load configuration from files and environment.

Parameters:

Name	Type	Description	Default
config_path	Optional[Path]	Optional path to user config file	None

Returns:

Type	Description
SsiambConfig	Merged configuration object

Source code in src/ssiamb/config.py

@classmethod
def load(cls, config_path: Optional[Path] = None) -> SsiambConfig:
    """
    Load configuration from files and environment.

    Args:
        config_path: Optional path to user config file

    Returns:
        Merged configuration object
    """
    # Start with built-in defaults
    config = cls._load_defaults()

    # Overlay user config if provided
    if config_path and config_path.exists():
        user_config = cls._load_yaml(config_path)
        config = cls._merge_configs(config, user_config)

    # Apply environment variable overrides
    config = cls._apply_env_overrides(config)

    # Normalize species alias keys for consistent lookup
    if "species_aliases" in config:
        config["species_aliases"] = cls._normalize_species_aliases(config["species_aliases"])

    return cls(
        thresholds=config.get("thresholds", {}),
        species_aliases=config.get("species_aliases", {}),
        tools=config.get("tools", {}),
        output=config.get("output", {})
    )

get_config()

Get the global configuration instance.

Source code in src/ssiamb/config.py

def get_config() -> SsiambConfig:
    """Get the global configuration instance."""
    global _config
    if _config is None:
        _config = SsiambConfig.load()
    return _config

load_config(config_path=None)

Load and set configuration from file.

Source code in src/ssiamb/config.py

def load_config(config_path: Optional[Path] = None) -> SsiambConfig:
    """Load and set configuration from file."""
    config = SsiambConfig.load(config_path)
    set_config(config)
    return config

set_config(config)

Set the global configuration instance.

Source code in src/ssiamb/config.py

def set_config(config: SsiambConfig) -> None:
    """Set the global configuration instance."""
    global _config
    _config = config

Mapping

Mapping and index handling for minimap2 and bwa-mem2.

This module provides functionality for: - Building and managing indexes for reference sequences - Mapping paired-end FASTQ files to references - Generating sorted BAM files with proper read groups

ExternalToolError

Bases: Exception

Raised when external tools are missing or fail.

MappingError

Bases: Exception

Raised when mapping operations fail.

calculate_mapping_rate(bam_path)

Calculate mapping rate from BAM file using samtools stats.

Parameters:

Name	Type	Description	Default
bam_path	Path	Path to sorted BAM file	required

Returns:

Type	Description
float	Mapping rate as fraction (0.0-1.0)

Raises:

Type	Description
MappingError	If BAM file doesn't exist or samtools fails
ExternalToolError	If samtools is not available

Source code in src/ssiamb/mapping.py

def calculate_mapping_rate(bam_path: Path) -> float:
    """
    Calculate mapping rate from BAM file using samtools stats.

    Args:
        bam_path: Path to sorted BAM file

    Returns:
        Mapping rate as fraction (0.0-1.0)

    Raises:
        MappingError: If BAM file doesn't exist or samtools fails
        ExternalToolError: If samtools is not available
    """
    if not shutil.which("samtools"):
        raise ExternalToolError("samtools not found in PATH")

    if not bam_path.exists():
        raise MappingError(f"BAM file not found: {bam_path}")

    logger.debug(f"Calculating mapping rate for {bam_path}")

    try:
        # Run samtools stats to get read counts
        cmd = ["samtools", "stats", str(bam_path)]
        result = subprocess.run(
            cmd,
            capture_output=True,
            text=True,
            check=True
        )

        total_reads = 0
        mapped_reads = 0

        # Parse samtools stats output
        for line in result.stdout.splitlines():
            if line.startswith("SN\t"):
                parts = line.split("\t")
                if len(parts) >= 3:
                    metric = parts[1]
                    value_str = parts[2].strip()

                    # Extract numbers (some lines have additional text after the number)
                    try:
                        value = int(value_str.split()[0])
                    except (ValueError, IndexError):
                        continue

                    if "raw total sequences:" in metric:
                        total_reads = value
                    elif "reads mapped:" in metric:
                        mapped_reads = value

        if total_reads == 0:
            logger.warning(f"No reads found in BAM file {bam_path}")
            return 0.0

        mapping_rate = mapped_reads / total_reads
        logger.debug(f"Mapping rate: {mapped_reads}/{total_reads} = {mapping_rate:.4f}")

        return mapping_rate

    except subprocess.CalledProcessError as e:
        raise MappingError(f"samtools stats failed: {e.stderr}")
    except Exception as e:
        raise MappingError(f"Failed to calculate mapping rate: {e}")

check_external_tools()

Check availability and versions of external mapping tools.

Returns:

Type	Description
Dict[str, Dict[str, str]]	Dictionary mapping tool names to availability and version info.
Dict[str, Dict[str, str]]	Each tool entry contains:
Dict[str, Dict[str, str]]	'available': boolean status
Dict[str, Dict[str, str]]	'version': version string if available, or error message

Examples:

>>> tools = check_external_tools()
>>> if not tools['minimap2']['available']:
...     raise ExternalToolError("minimap2 not found")
>>> print(f"minimap2 version: {tools['minimap2']['version']}")

Source code in src/ssiamb/mapping.py

def check_external_tools() -> Dict[str, Dict[str, str]]:
    """
    Check availability and versions of external mapping tools.

    Returns:
        Dictionary mapping tool names to availability and version info.
        Each tool entry contains:
        - 'available': boolean status
        - 'version': version string if available, or error message

    Examples:
        >>> tools = check_external_tools()
        >>> if not tools['minimap2']['available']:
        ...     raise ExternalToolError("minimap2 not found")
        >>> print(f"minimap2 version: {tools['minimap2']['version']}")
    """
    tools = {}

    # Define tools and their version commands
    tool_commands = {
        'minimap2': ['minimap2', '--version'],
        'bwa-mem2': ['bwa-mem2', 'version'],
        'samtools': ['samtools', '--version']
    }

    for tool, version_cmd in tool_commands.items():
        tools[tool] = {'available': False, 'version': 'unknown'}

        # Check if tool is in PATH
        if shutil.which(tool) is None:
            tools[tool]['version'] = 'not found in PATH'
            logger.debug(f"Tool {tool}: not found")
            continue

        # Try to get version
        try:
            result = subprocess.run(
                version_cmd,
                capture_output=True,
                text=True,
                timeout=5  # Prevent hanging
            )
            if result.returncode == 0:
                tools[tool]['available'] = True
                # Extract version from output (first line typically)
                version_line = result.stdout.strip().split('\n')[0] if result.stdout else 'version check succeeded'
                tools[tool]['version'] = version_line
                logger.debug(f"Tool {tool}: available, version: {version_line}")
            else:
                tools[tool]['version'] = f'version check failed (exit {result.returncode})'
                logger.debug(f"Tool {tool}: found but version check failed")
        except subprocess.TimeoutExpired:
            tools[tool]['version'] = 'version check timed out'
            logger.debug(f"Tool {tool}: found but version check timed out")
        except Exception as e:
            tools[tool]['version'] = f'version check error: {e}'
            logger.debug(f"Tool {tool}: found but version check error: {e}")

    return tools

ensure_indexes_self(fasta_path, mapper)

Ensure index files exist for self-mode mapping, building them if missing.

This function builds index files next to the FASTA file if they don't exist. For minimap2, creates .mmi file. For bwa-mem2, creates .0123, .amb, .ann, .pac, .bwt.2bit.64 files.

Parameters:

Name	Type	Description	Default
fasta_path	Path	Path to reference FASTA file	required
mapper	Mapper	Mapper type (minimap2 or bwa-mem2)	required

Raises:

Type	Description
FileNotFoundError	If FASTA file doesn't exist
ExternalToolError	If required mapper tool is not available
MappingError	If index building fails

Examples:

>>> ensure_indexes_self(Path("ref.fasta"), Mapper.MINIMAP2)
# Creates ref.mmi if it doesn't exist

Source code in src/ssiamb/mapping.py

def ensure_indexes_self(fasta_path: Path, mapper: Mapper) -> None:
    """
    Ensure index files exist for self-mode mapping, building them if missing.

    This function builds index files next to the FASTA file if they don't exist.
    For minimap2, creates .mmi file. For bwa-mem2, creates .0123, .amb, .ann, 
    .pac, .bwt.2bit.64 files.

    Args:
        fasta_path: Path to reference FASTA file
        mapper: Mapper type (minimap2 or bwa-mem2)

    Raises:
        FileNotFoundError: If FASTA file doesn't exist
        ExternalToolError: If required mapper tool is not available
        MappingError: If index building fails

    Examples:
        >>> ensure_indexes_self(Path("ref.fasta"), Mapper.MINIMAP2)
        # Creates ref.mmi if it doesn't exist
    """
    if not fasta_path.exists():
        raise FileNotFoundError(f"Reference FASTA not found: {fasta_path}")

    # Check if indexes already exist
    if indexes_exist(fasta_path, mapper):
        logger.info(f"Index files already exist for {fasta_path} ({mapper.value})")
        return

    # Check tool availability
    tools = check_external_tools()
    tool_name = mapper.value  # Use the actual tool name without modification
    if not tools.get(tool_name, {}).get('available', False):
        version_info = tools.get(tool_name, {}).get('version', 'unknown')
        raise ExternalToolError(f"{tool_name} not available: {version_info}")

    logger.info(f"Building {mapper.value} index for {fasta_path}")

    try:
        if mapper == Mapper.MINIMAP2:
            _build_minimap2_index(fasta_path)
        elif mapper == Mapper.BWA_MEM2:
            _build_bwa_mem2_index(fasta_path)
        else:
            raise ValueError(f"Unsupported mapper: {mapper}")

    except subprocess.CalledProcessError as e:
        raise MappingError(f"Failed to build {mapper.value} index: {e}")

get_index_files(fasta_path, mapper)

Get expected index file paths for a given FASTA and mapper.

Parameters:

Name	Type	Description	Default
fasta_path	Path	Path to reference FASTA file	required
mapper	Mapper	Mapper type (minimap2 or bwa-mem2)	required

Returns:

Type	Description
List[Path]	List of expected index file paths

Raises:

Type	Description
ValueError	If mapper is not supported

Source code in src/ssiamb/mapping.py

def get_index_files(fasta_path: Path, mapper: Mapper) -> List[Path]:
    """
    Get expected index file paths for a given FASTA and mapper.

    Args:
        fasta_path: Path to reference FASTA file
        mapper: Mapper type (minimap2 or bwa-mem2)

    Returns:
        List of expected index file paths

    Raises:
        ValueError: If mapper is not supported
    """
    if mapper == Mapper.MINIMAP2:
        return [fasta_path.with_suffix('.mmi')]

    elif mapper == Mapper.BWA_MEM2:
        # BWA-MEM2 creates index files with the original filename plus extensions
        extensions = ['.0123', '.amb', '.ann', '.pac', '.bwt.2bit.64']
        return [Path(str(fasta_path) + ext) for ext in extensions]

    else:
        raise ValueError(f"Unsupported mapper: {mapper}")

index_bam(bam_path)

Index a BAM file using samtools.

Parameters:

Name	Type	Description	Default
bam_path	Path	Path to BAM file to index	required

Raises:

Type	Description
MappingError	If indexing fails

Source code in src/ssiamb/mapping.py

def index_bam(bam_path: Path) -> None:
    """
    Index a BAM file using samtools.

    Args:
        bam_path: Path to BAM file to index

    Raises:
        MappingError: If indexing fails
    """
    if not bam_path.exists():
        raise MappingError(f"BAM file not found: {bam_path}")

    index_path = bam_path.with_suffix('.bam.bai')

    logger.debug(f"Indexing BAM file: {bam_path}")

    try:
        cmd = ["samtools", "index", str(bam_path)]
        result = subprocess.run(
            cmd,
            capture_output=True,
            text=True,
            check=True
        )

        if not index_path.exists():
            raise MappingError(f"BAM index was not created: {index_path}")

        logger.debug(f"Created BAM index: {index_path}")

    except subprocess.CalledProcessError as e:
        raise MappingError(f"Failed to index BAM {bam_path}: {e.stderr}")

indexes_exist(fasta_path, mapper)

Check if all required index files exist for a given FASTA and mapper.

Parameters:

Name	Type	Description	Default
fasta_path	Path	Path to reference FASTA file	required
mapper	Mapper	Mapper type	required

Returns:

Type	Description
bool	True if all index files exist, False otherwise

Source code in src/ssiamb/mapping.py

def indexes_exist(fasta_path: Path, mapper: Mapper) -> bool:
    """
    Check if all required index files exist for a given FASTA and mapper.

    Args:
        fasta_path: Path to reference FASTA file
        mapper: Mapper type

    Returns:
        True if all index files exist, False otherwise
    """
    index_files = get_index_files(fasta_path, mapper)
    return all(idx_file.exists() for idx_file in index_files)

map_fastqs(mapper, fasta_path, r1_path, r2_path, sample_name, threads=4, output_path=None)

Map paired-end FASTQ files to reference and return sorted BAM.

Parameters:

Name	Type	Description	Default
mapper	Mapper	Mapper type (minimap2 or bwa-mem2)	required
fasta_path	Path	Path to reference FASTA file	required
r1_path	Path	Path to R1 FASTQ file	required
r2_path	Path	Path to R2 FASTQ file	required
sample_name	str	Sample name for read group	required
threads	int	Number of threads to use	4
output_path	Optional[Path]	Output BAM path (auto-generated if None)	None

Returns:

Type	Description
Path	Path to sorted BAM file

Raises:

Type	Description
FileNotFoundError	If input files don't exist
ExternalToolError	If required tools are not available
MappingError	If mapping fails

Examples:

>>> bam_path = map_fastqs(
...     Mapper.MINIMAP2,
...     Path("ref.fasta"),
...     Path("sample_R1.fastq.gz"),
...     Path("sample_R2.fastq.gz"),
...     "sample123"
... )

Source code in src/ssiamb/mapping.py

def map_fastqs(
    mapper: Mapper,
    fasta_path: Path,
    r1_path: Path,
    r2_path: Path,
    sample_name: str,
    threads: int = 4,
    output_path: Optional[Path] = None
) -> Path:
    """
    Map paired-end FASTQ files to reference and return sorted BAM.

    Args:
        mapper: Mapper type (minimap2 or bwa-mem2)
        fasta_path: Path to reference FASTA file
        r1_path: Path to R1 FASTQ file
        r2_path: Path to R2 FASTQ file  
        sample_name: Sample name for read group
        threads: Number of threads to use
        output_path: Output BAM path (auto-generated if None)

    Returns:
        Path to sorted BAM file

    Raises:
        FileNotFoundError: If input files don't exist
        ExternalToolError: If required tools are not available
        MappingError: If mapping fails

    Examples:
        >>> bam_path = map_fastqs(
        ...     Mapper.MINIMAP2,
        ...     Path("ref.fasta"),
        ...     Path("sample_R1.fastq.gz"),
        ...     Path("sample_R2.fastq.gz"),
        ...     "sample123"
        ... )
    """
    # Validate input files
    for file_path in [fasta_path, r1_path, r2_path]:
        if not file_path.exists():
            raise FileNotFoundError(f"Input file not found: {file_path}")

    # Check tool availability
    tools = check_external_tools()
    tool_name = mapper.value  # Use the actual tool name without modification
    if not tools.get(tool_name, {}).get('available', False) or not tools.get('samtools', {}).get('available', False):
        missing = [t for t, info in tools.items() 
                  if not info.get('available', False) and t in [tool_name, 'samtools']]
        raise ExternalToolError(f"Required tools not available: {missing}")

    # Ensure indexes exist
    ensure_indexes_self(fasta_path, mapper)

    # Generate output path if not provided
    if output_path is None:
        output_path = Path(f"{sample_name}.sorted.bam")

    logger.info(f"Mapping {r1_path.name} and {r2_path.name} to {fasta_path.name} using {mapper.value}")

    try:
        if mapper == Mapper.MINIMAP2:
            _map_with_minimap2(fasta_path, r1_path, r2_path, sample_name, threads, output_path)
        elif mapper == Mapper.BWA_MEM2:
            _map_with_bwa_mem2(fasta_path, r1_path, r2_path, sample_name, threads, output_path)
        else:
            raise ValueError(f"Unsupported mapper: {mapper}")

    except subprocess.CalledProcessError as e:
        raise MappingError(f"Mapping with {mapper.value} failed: {e}")

    if not output_path.exists():
        raise MappingError(f"Output BAM was not created: {output_path}")

    # Index the BAM file
    index_bam(output_path)

    logger.info(f"Created sorted BAM: {output_path}")
    return output_path

Variant Calling

Variant calling module for ssiamb.

This module implements variant calling using BBTools and bcftools pipelines. Supports both bacterial genome analysis with appropriate ploidy and quality settings.

VariantCallResult(vcf_path, caller, success, error_message=None, runtime_seconds=None) dataclass

Result of variant calling operation.

VariantCallingError

Bases: Exception

Raised when variant calling fails.

call_variants(bam_path, reference_path, output_vcf, caller, sample_name, threads=1, mapq_min=20, baseq_min=20, minallelefraction=0.0, bbtools_mem=None)

Run variant calling with the specified caller.

Parameters:

Name	Type	Description	Default
bam_path	Path	Input BAM file	required
reference_path	Path	Reference genome FASTA	required
output_vcf	Path	Output VCF file path	required
caller	Caller	Variant caller to use	required
sample_name	str	Sample name for VCF header	required
threads	int	Number of threads to use	1
mapq_min	int	Minimum mapping quality	20
baseq_min	int	Minimum base quality	20
minallelefraction	float	Minimum allele fraction (BBTools only)	0.0
bbtools_mem	Optional[str]	BBTools heap memory (e.g., '4g', '8g')	None

Returns:

Type	Description
VariantCallResult	VariantCallResult with execution details

Raises:

Type	Description
VariantCallingError	If caller tools are not available or calling fails

Source code in src/ssiamb/calling.py

def call_variants(
    bam_path: Path,
    reference_path: Path,
    output_vcf: Path,
    caller: Caller,
    sample_name: str,
    threads: int = 1,
    mapq_min: int = 20,
    baseq_min: int = 20,
    minallelefraction: float = 0.0,
    bbtools_mem: Optional[str] = None,
) -> VariantCallResult:
    """
    Run variant calling with the specified caller.

    Args:
        bam_path: Input BAM file
        reference_path: Reference genome FASTA
        output_vcf: Output VCF file path
        caller: Variant caller to use
        sample_name: Sample name for VCF header
        threads: Number of threads to use
        mapq_min: Minimum mapping quality
        baseq_min: Minimum base quality
        minallelefraction: Minimum allele fraction (BBTools only)
        bbtools_mem: BBTools heap memory (e.g., '4g', '8g')

    Returns:
        VariantCallResult with execution details

    Raises:
        VariantCallingError: If caller tools are not available or calling fails
    """
    # Check tool availability
    if not caller_tools_available(caller):
        raise VariantCallingError(f"Required tools for {caller.value} are not available")

    # Validate inputs
    if not bam_path.exists():
        raise VariantCallingError(f"BAM file not found: {bam_path}")
    if not reference_path.exists():
        raise VariantCallingError(f"Reference file not found: {reference_path}")

    # Run appropriate caller
    if caller == Caller.BBTOOLS:
        result = run_bbtools_calling(
            bam_path=bam_path,
            reference_path=reference_path,
            output_vcf=output_vcf,
            sample_name=sample_name,
            threads=threads,
            mapq_min=mapq_min,
            baseq_min=baseq_min,
            minallelefraction=minallelefraction,
            bbtools_mem=bbtools_mem,
        )
    elif caller == Caller.BCFTOOLS:
        result = run_bcftools_calling(
            bam_path=bam_path,
            reference_path=reference_path,
            output_vcf=output_vcf,
            sample_name=sample_name,
            threads=threads,
            mapq_min=mapq_min,
            baseq_min=baseq_min,
        )
    else:
        raise VariantCallingError(f"Unsupported caller: {caller}")

    # Raise exception if calling failed
    if not result.success:
        raise VariantCallingError(result.error_message or f"{caller.value} variant calling failed")

    return result

caller_tools_available(caller)

Check if all required tools for the specified caller are available.

Parameters:

Name	Type	Description	Default
caller	Caller	Variant caller to check	required

Returns:

Type	Description
bool	True if all required tools are available, False otherwise

Source code in src/ssiamb/calling.py

def caller_tools_available(caller: Caller) -> bool:
    """
    Check if all required tools for the specified caller are available.

    Args:
        caller: Variant caller to check

    Returns:
        True if all required tools are available, False otherwise
    """
    tools = check_caller_tools_detailed(caller)
    return all(tool_info.get('available', False) for tool_info in tools.values())

check_caller_tools(caller)

Check if required tools for the specified caller are available.

Parameters:

Name	Type	Description	Default
caller	Caller	Variant caller to check	required

Returns:

Type	Description
bool	True if all required tools are available, False otherwise

Source code in src/ssiamb/calling.py

def check_caller_tools(caller: Caller) -> bool:
    """
    Check if required tools for the specified caller are available.

    Args:
        caller: Variant caller to check

    Returns:
        True if all required tools are available, False otherwise
    """
    if caller == Caller.BBTOOLS:
        # Check for BBTools executables
        return (shutil.which("pileup.sh") is not None and 
                shutil.which("callvariants.sh") is not None)
    elif caller == Caller.BCFTOOLS:
        # Check for bcftools
        return shutil.which("bcftools") is not None
    else:
        raise ValueError(f"Unknown caller: {caller}")

check_caller_tools_detailed(caller)

Check availability and versions of required tools for the specified caller.

Parameters:

Name	Type	Description	Default
caller	Caller	Variant caller to check	required

Returns:

Type	Description
Dict[str, Dict[str, str]]	Dictionary mapping tool names to availability and version info.
Dict[str, Dict[str, str]]	Each tool entry contains:
Dict[str, Dict[str, str]]	'available': boolean status
Dict[str, Dict[str, str]]	'version': version string if available, or error message

Source code in src/ssiamb/calling.py

def check_caller_tools_detailed(caller: Caller) -> Dict[str, Dict[str, str]]:
    """
    Check availability and versions of required tools for the specified caller.

    Args:
        caller: Variant caller to check

    Returns:
        Dictionary mapping tool names to availability and version info.
        Each tool entry contains:
        - 'available': boolean status
        - 'version': version string if available, or error message
    """
    tools = {}

    if caller == Caller.BBTOOLS:
        # Check for BBTools executables
        for tool in ["pileup.sh", "callvariants.sh"]:
            tools[tool] = {'available': False, 'version': 'unknown'}

            if shutil.which(tool) is None:
                tools[tool]['version'] = 'not found in PATH'
                logger.debug(f"Tool {tool}: not found")
                continue

            # BBTools scripts typically don't have version commands, so just mark as available
            tools[tool]['available'] = True
            tools[tool]['version'] = 'BBTools (version check not supported)'
            logger.debug(f"Tool {tool}: available")

    elif caller == Caller.BCFTOOLS:
        # Check for bcftools
        tools['bcftools'] = {'available': False, 'version': 'unknown'}

        if shutil.which("bcftools") is None:
            tools['bcftools']['version'] = 'not found in PATH'
            logger.debug("Tool bcftools: not found")
        else:
            # Try to get bcftools version
            try:
                result = subprocess.run(
                    ['bcftools', '--version'],
                    capture_output=True,
                    text=True,
                    timeout=5
                )
                if result.returncode == 0:
                    tools['bcftools']['available'] = True
                    # Extract version from output (first line typically)
                    version_line = result.stdout.strip().split('\n')[0] if result.stdout else 'version check succeeded'
                    tools['bcftools']['version'] = version_line
                    logger.debug(f"Tool bcftools: available, version: {version_line}")
                else:
                    tools['bcftools']['version'] = f'version check failed (exit {result.returncode})'
                    logger.debug("Tool bcftools: found but version check failed")
            except subprocess.TimeoutExpired:
                tools['bcftools']['version'] = 'version check timed out'
                logger.debug("Tool bcftools: found but version check timed out")
            except Exception as e:
                tools['bcftools']['version'] = f'version check error: {e}'
                logger.debug(f"Tool bcftools: found but version check error: {e}")
    else:
        raise ValueError(f"Unknown caller: {caller}")

    return tools

get_available_callers()

Get list of available variant callers based on tool availability.

Returns:

Type	Description
List[Caller]	List of available Caller enum values

Source code in src/ssiamb/calling.py

def get_available_callers() -> List[Caller]:
    """
    Get list of available variant callers based on tool availability.

    Returns:
        List of available Caller enum values
    """
    available = []

    for caller in Caller:
        if caller_tools_available(caller):
            available.append(caller)

    return available

run_bbtools_calling(bam_path, reference_path, output_vcf, sample_name, threads=1, mapq_min=20, baseq_min=20, minallelefraction=0.0, bbtools_mem=None)

Run BBTools variant calling pipeline.

Executes: 1. callvariants.sh directly with BAM input (no pileup step needed)

Parameters:

Name	Type	Description	Default
bam_path	Path	Input BAM file	required
reference_path	Path	Reference genome FASTA	required
output_vcf	Path	Output VCF file path	required
sample_name	str	Sample name for VCF header	required
threads	int	Number of threads to use	1
mapq_min	int	Minimum mapping quality	20
baseq_min	int	Minimum base quality	20
minallelefraction	float	Minimum allele fraction for variant calling	0.0
bbtools_mem	Optional[str]	BBTools heap memory (e.g., '4g', '8g')	None

Returns:

Type	Description
VariantCallResult	VariantCallResult with execution details

Source code in src/ssiamb/calling.py

def run_bbtools_calling(
    bam_path: Path,
    reference_path: Path,
    output_vcf: Path,
    sample_name: str,
    threads: int = 1,
    mapq_min: int = 20,
    baseq_min: int = 20,
    minallelefraction: float = 0.0,
    bbtools_mem: Optional[str] = None,
) -> VariantCallResult:
    """
    Run BBTools variant calling pipeline.

    Executes:
    1. callvariants.sh directly with BAM input (no pileup step needed)

    Args:
        bam_path: Input BAM file
        reference_path: Reference genome FASTA
        output_vcf: Output VCF file path
        sample_name: Sample name for VCF header
        threads: Number of threads to use
        mapq_min: Minimum mapping quality
        baseq_min: Minimum base quality
        minallelefraction: Minimum allele fraction for variant calling
        bbtools_mem: BBTools heap memory (e.g., '4g', '8g')

    Returns:
        VariantCallResult with execution details
    """
    import time
    start_time = time.time()

    try:
        # Ensure output directory exists
        output_vcf.parent.mkdir(parents=True, exist_ok=True)

        # Run callvariants.sh directly with BAM input (simpler approach)
        callvariants_cmd = [
            "callvariants.sh",
            f"in={bam_path}",
            f"ref={reference_path}",
            f"vcf={output_vcf}",           # Use vcf= for VCF output
            "ploidy=1",                   # Haploid organism
            "clearfilters=t",             # Clear all filters to get raw variants
            f"minallelefraction={minallelefraction}",
            f"minavgmapq={mapq_min}",
            f"minquality={baseq_min}",
            f"threads={threads}",
        ]

        # Add memory setting if provided
        if bbtools_mem:
            callvariants_cmd.append(f"-Xmx{bbtools_mem}")

        logger.info(f"Running BBTools callvariants: {' '.join(map(str, callvariants_cmd))}")

        result = subprocess.run(
            callvariants_cmd,
            capture_output=True,
            text=True,
            check=False,
            timeout=3600  # 1 hour timeout
        )

        if result.returncode != 0:
            error_msg = f"BBTools callvariants failed (exit {result.returncode}): {result.stderr}"
            logger.error(error_msg)
            return VariantCallResult(
                vcf_path=output_vcf,
                caller=Caller.BBTOOLS,
                success=False,
                error_message=error_msg,
                runtime_seconds=time.time() - start_time
            )

        # Log stdout for debugging
        if result.stdout:
            logger.debug(f"BBTools callvariants stdout: {result.stdout.strip()}")

        # Verify output VCF was created
        if not output_vcf.exists() or output_vcf.stat().st_size == 0:
            error_msg = "BBTools callvariants completed but no VCF output found"
            logger.error(error_msg)
            return VariantCallResult(
                vcf_path=output_vcf,
                caller=Caller.BBTOOLS,
                success=False,
                error_message=error_msg,
                runtime_seconds=time.time() - start_time
            )

        logger.info(f"BBTools variant calling completed: {output_vcf}")
        return VariantCallResult(
            vcf_path=output_vcf,
            caller=Caller.BBTOOLS,
            success=True,
            runtime_seconds=time.time() - start_time
        )

    except subprocess.TimeoutExpired:
        error_msg = "BBTools variant calling timed out"
        logger.error(error_msg)
        return VariantCallResult(
            vcf_path=output_vcf,
            caller=Caller.BBTOOLS,
            success=False,
            error_message=error_msg,
            runtime_seconds=time.time() - start_time
        )
    except Exception as e:
        error_msg = f"BBTools variant calling failed with exception: {e}"
        logger.error(error_msg)
        return VariantCallResult(
            vcf_path=output_vcf,
            caller=Caller.BBTOOLS,
            success=False,
            error_message=error_msg,
            runtime_seconds=time.time() - start_time
        )

run_bcftools_calling(bam_path, reference_path, output_vcf, sample_name, threads=1, mapq_min=20, baseq_min=20)

Run bcftools variant calling pipeline.

Executes: 1. bcftools mpileup -q20 -Q20 -B -a AD,ADF,ADR,DP 2. bcftools call -m --ploidy 1

Parameters:

Name	Type	Description	Default
bam_path	Path	Input BAM file	required
reference_path	Path	Reference genome FASTA	required
output_vcf	Path	Output VCF file path	required
sample_name	str	Sample name for VCF header	required
threads	int	Number of threads to use	1
mapq_min	int	Minimum mapping quality	20
baseq_min	int	Minimum base quality	20

Returns:

Type	Description
VariantCallResult	VariantCallResult with execution details

Source code in src/ssiamb/calling.py

def run_bcftools_calling(
    bam_path: Path,
    reference_path: Path,
    output_vcf: Path,
    sample_name: str,
    threads: int = 1,
    mapq_min: int = 20,
    baseq_min: int = 20,
) -> VariantCallResult:
    """
    Run bcftools variant calling pipeline.

    Executes:
    1. bcftools mpileup -q20 -Q20 -B -a AD,ADF,ADR,DP
    2. bcftools call -m --ploidy 1

    Args:
        bam_path: Input BAM file
        reference_path: Reference genome FASTA
        output_vcf: Output VCF file path
        sample_name: Sample name for VCF header
        threads: Number of threads to use
        mapq_min: Minimum mapping quality
        baseq_min: Minimum base quality

    Returns:
        VariantCallResult with execution details
    """
    import time
    start_time = time.time()

    try:
        # Ensure output directory exists
        output_vcf.parent.mkdir(parents=True, exist_ok=True)

        # Combine mpileup and call in a pipeline
        mpileup_cmd = [
            "bcftools", "mpileup",
            "-Ou",                     # Uncompressed BCF output (required by spec)
            f"--threads", str(threads),
            f"-q", str(mapq_min),      # Minimum mapping quality
            f"-Q", str(baseq_min),     # Minimum base quality
            "-B",                      # Disable BAQ computation
            "--max-depth", "100000",   # Maximum depth (required by spec)
            "-a", "FORMAT/AD,ADF,ADR,DP",  # Annotations to include (fixed format)
            "-f", str(reference_path),  # Reference FASTA
            str(bam_path)              # Input BAM
        ]

        call_cmd = [
            "bcftools", "call",
            f"--threads", str(threads),
            "-m",                      # Multiallelic caller
            "--ploidy", "1",           # Haploid organism
            "--prior", "1.1e-3",       # Prior for novel mutation rate (required by spec)
            "-v",                      # Output only variants
            "-o", str(output_vcf)      # Output file
        ]

        logger.info(f"Running bcftools pipeline: {' '.join(mpileup_cmd)} | {' '.join(call_cmd)}")

        # Run mpileup | call pipeline
        mpileup_proc = subprocess.Popen(
            mpileup_cmd,
            stdout=subprocess.PIPE,
            stderr=subprocess.PIPE,
            text=True
        )

        call_proc = subprocess.Popen(
            call_cmd,
            stdin=mpileup_proc.stdout,
            stderr=subprocess.PIPE,
            text=True
        )

        # Close mpileup stdout to allow it to receive SIGPIPE
        if mpileup_proc.stdout:
            mpileup_proc.stdout.close()

        # Wait for both processes
        call_stderr = call_proc.communicate()[1]
        mpileup_stderr = mpileup_proc.communicate()[1]

        # Check return codes
        if mpileup_proc.returncode != 0:
            error_msg = f"bcftools mpileup failed (exit {mpileup_proc.returncode}): {mpileup_stderr}"
            logger.error(error_msg)
            return VariantCallResult(
                vcf_path=output_vcf,
                caller=Caller.BCFTOOLS,
                success=False,
                error_message=error_msg,
                runtime_seconds=time.time() - start_time
            )

        if call_proc.returncode != 0:
            error_msg = f"bcftools call failed (exit {call_proc.returncode}): {call_stderr}"
            logger.error(error_msg)
            return VariantCallResult(
                vcf_path=output_vcf,
                caller=Caller.BCFTOOLS,
                success=False,
                error_message=error_msg,
                runtime_seconds=time.time() - start_time
            )

        # Verify output VCF was created
        if not output_vcf.exists():
            error_msg = "bcftools call completed but no VCF output found"
            logger.error(error_msg)
            return VariantCallResult(
                vcf_path=output_vcf,
                caller=Caller.BCFTOOLS,
                success=False,
                error_message=error_msg,
                runtime_seconds=time.time() - start_time
            )

        logger.info(f"bcftools variant calling completed: {output_vcf}")
        return VariantCallResult(
            vcf_path=output_vcf,
            caller=Caller.BCFTOOLS,
            success=True,
            runtime_seconds=time.time() - start_time
        )

    except Exception as e:
        error_msg = f"bcftools variant calling failed with exception: {e}"
        logger.error(error_msg)
        return VariantCallResult(
            vcf_path=output_vcf,
            caller=Caller.BCFTOOLS,
            success=False,
            error_message=error_msg,
            runtime_seconds=time.time() - start_time
        )

Depth Analysis

Depth analysis using mosdepth for computing denominator metrics.

This module provides functionality for: - Running mosdepth on BAM files to compute depth statistics - Parsing mosdepth summary output files - Computing callable bases, genome length, breadth, and mean depth - Filtering contigs by size (≥500 bp) for consistent numerator/denominator

ContigDepthStats(name, length, bases_covered, mean_depth, breadth_10x) dataclass

Depth statistics for a single contig.

is_long_enough property

Check if contig meets minimum length threshold (≥500 bp).

DepthAnalysisError

Bases: Exception

Raised when depth analysis operations fail.

DepthSummary(callable_bases, genome_length, breadth_10x, mean_depth, total_contigs, included_contigs, contig_stats) dataclass

Summary of depth analysis results.

included_contig_names property

Set of contig names that meet length threshold.

analyze_depth(bam_path, output_dir, sample_name, mapq_threshold=30, depth_threshold=10, threads=4)

Complete depth analysis workflow: run mosdepth and parse results.

Parameters:

Name	Type	Description	Default
bam_path	Path	Path to input BAM file	required
output_dir	Path	Directory for output files	required
sample_name	str	Sample name for output prefix	required
mapq_threshold	int	Minimum mapping quality (default: 30)	30
depth_threshold	int	Depth threshold for breadth calculation (default: 10)	10
threads	int	Number of threads to use	4

Returns:

Type	Description
DepthSummary	DepthSummary object with computed statistics

Raises:

Type	Description
DepthAnalysisError	If analysis fails at any step

Source code in src/ssiamb/depth.py

def analyze_depth(
    bam_path: Path,
    output_dir: Path,
    sample_name: str,
    mapq_threshold: int = 30,
    depth_threshold: int = 10,
    threads: int = 4
) -> DepthSummary:
    """
    Complete depth analysis workflow: run mosdepth and parse results.

    Args:
        bam_path: Path to input BAM file
        output_dir: Directory for output files
        sample_name: Sample name for output prefix
        mapq_threshold: Minimum mapping quality (default: 30)
        depth_threshold: Depth threshold for breadth calculation (default: 10)
        threads: Number of threads to use

    Returns:
        DepthSummary object with computed statistics

    Raises:
        DepthAnalysisError: If analysis fails at any step
    """
    logger.info(f"Starting depth analysis for {sample_name}")

    # Create output directory if needed
    output_dir.mkdir(parents=True, exist_ok=True)

    # Run mosdepth
    output_prefix = output_dir / f"{sample_name}.depth"
    summary_file = run_mosdepth(
        bam_path=bam_path,
        output_prefix=output_prefix,
        mapq_threshold=mapq_threshold,
        threads=threads
    )

    # Parse results
    summary = parse_mosdepth_summary(summary_file)

    logger.info(f"Depth analysis complete for {sample_name}: "
               f"{summary.genome_length:,} bp genome, "
               f"{summary.breadth_10x:.2%} breadth ≥{depth_threshold}x")

    return summary

check_mosdepth_available()

Check if mosdepth is available in PATH.

Returns:

Type	Description
bool	True if mosdepth is available, False otherwise

Source code in src/ssiamb/depth.py

def check_mosdepth_available() -> bool:
    """
    Check if mosdepth is available in PATH.

    Returns:
        True if mosdepth is available, False otherwise
    """
    try:
        subprocess.run(['mosdepth', '--version'], 
                      capture_output=True, check=True)
        return True
    except (subprocess.CalledProcessError, FileNotFoundError):
        return False

get_depth_from_existing_summary(summary_file)

Parse existing mosdepth summary file without running mosdepth.

Useful for reusing existing depth analysis results.

Parameters:

Name	Type	Description	Default
summary_file	Path	Path to existing .mosdepth.summary.txt file	required

Returns:

Type	Description
DepthSummary	DepthSummary object with computed statistics

Source code in src/ssiamb/depth.py

def get_depth_from_existing_summary(summary_file: Path) -> DepthSummary:
    """
    Parse existing mosdepth summary file without running mosdepth.

    Useful for reusing existing depth analysis results.

    Args:
        summary_file: Path to existing .mosdepth.summary.txt file

    Returns:
        DepthSummary object with computed statistics
    """
    return parse_mosdepth_summary(summary_file)

list_included_contigs(summary_file, min_len=500)

Get set of contig names that meet minimum length threshold.

Parameters:

Name	Type	Description	Default
summary_file	Path	Path to .mosdepth.summary.txt file	required
min_len	int	Minimum contig length threshold	500

Returns:

Type	Description
Set[str]	Set of contig names that meet the length threshold

Source code in src/ssiamb/depth.py

def list_included_contigs(summary_file: Path, min_len: int = 500) -> Set[str]:
    """
    Get set of contig names that meet minimum length threshold.

    Args:
        summary_file: Path to .mosdepth.summary.txt file
        min_len: Minimum contig length threshold

    Returns:
        Set of contig names that meet the length threshold
    """
    if not summary_file.exists():
        logger.warning(f"Summary file not found: {summary_file}")
        return set()

    try:
        included_contigs = set()

        with open(summary_file, 'r') as f:
            for line in f:
                line = line.strip()
                if not line or line.startswith('#'):
                    continue

                fields = line.split('\t')
                if len(fields) < 4:
                    continue

                contig_name = fields[0]
                try:
                    contig_length = int(fields[1])
                    if contig_length >= min_len:
                        included_contigs.add(contig_name)
                except (ValueError, IndexError):
                    continue

        logger.debug(f"Found {len(included_contigs)} contigs >= {min_len} bp")
        return included_contigs

    except Exception as e:
        logger.error(f"Error reading mosdepth summary {summary_file}: {e}")
        return set()

parse_mosdepth_summary(summary_file)

Parse mosdepth summary file to extract depth statistics.

The mosdepth summary file has format: chrom length bases mean min max contig1 5000 4850 25.3 0 100 ... total 50000 48500 24.8 0 100

Parameters:

Name	Type	Description	Default
summary_file	Path	Path to .mosdepth.summary.txt file	required

Returns:

Type	Description
DepthSummary	DepthSummary object with computed statistics

Raises:

Type	Description
DepthAnalysisError	If file parsing fails
FileNotFoundError	If summary file doesn't exist

Source code in src/ssiamb/depth.py

def parse_mosdepth_summary(summary_file: Path) -> DepthSummary:
    """
    Parse mosdepth summary file to extract depth statistics.

    The mosdepth summary file has format:
    chrom	length	bases	mean	min	max
    contig1	5000	4850	25.3	0	100
    ...
    total	50000	48500	24.8	0	100

    Args:
        summary_file: Path to .mosdepth.summary.txt file

    Returns:
        DepthSummary object with computed statistics

    Raises:
        DepthAnalysisError: If file parsing fails
        FileNotFoundError: If summary file doesn't exist
    """
    if not summary_file.exists():
        raise FileNotFoundError(f"mosdepth summary file not found: {summary_file}")

    logger.info(f"Parsing mosdepth summary: {summary_file}")

    try:
        contig_stats = []
        total_line = None

        with open(summary_file, 'r') as f:
            # Skip header line
            header = f.readline().strip()
            if not header.startswith('chrom'):
                raise DepthAnalysisError(f"Unexpected header format: {header}")

            for line_num, line in enumerate(f, start=2):
                line = line.strip()
                if not line:
                    continue

                parts = line.split('\t')
                if len(parts) < 4:
                    raise DepthAnalysisError(f"Invalid line format at line {line_num}: {line}")

                chrom, length_str, bases_str, mean_str = parts[:4]

                try:
                    length = int(length_str)
                    bases_covered = int(bases_str)
                    mean_depth = float(mean_str)
                except ValueError as e:
                    raise DepthAnalysisError(f"Invalid numeric value at line {line_num}: {e}")

                if chrom == 'total':
                    total_line = (length, bases_covered, mean_depth)
                else:
                    # Calculate breadth for this contig
                    breadth_10x = bases_covered / length if length > 0 else 0.0

                    contig_stats.append(ContigDepthStats(
                        name=chrom,
                        length=length,
                        bases_covered=bases_covered,
                        mean_depth=mean_depth,
                        breadth_10x=breadth_10x
                    ))

        if not contig_stats:
            raise DepthAnalysisError("No contig data found in summary file")

        # Compute summary statistics for contigs ≥500bp
        long_contigs = [c for c in contig_stats if c.is_long_enough]

        if not long_contigs:
            logger.warning("No contigs ≥500bp found - all contigs are too short")
            genome_length = 0
            callable_bases = 0
            breadth_10x = 0.0
            mean_depth = 0.0
        else:
            genome_length = sum(c.length for c in long_contigs)
            callable_bases = sum(c.bases_covered for c in long_contigs)
            breadth_10x = callable_bases / genome_length if genome_length > 0 else 0.0

            # Weighted mean depth across long contigs
            total_depth = sum(c.mean_depth * c.length for c in long_contigs)
            mean_depth = total_depth / genome_length if genome_length > 0 else 0.0

        summary = DepthSummary(
            callable_bases=callable_bases,
            genome_length=genome_length,
            breadth_10x=breadth_10x,
            mean_depth=mean_depth,
            total_contigs=len(contig_stats),
            included_contigs=len(long_contigs),
            contig_stats=contig_stats
        )

        logger.info(f"Parsed depth summary: {summary.included_contigs}/{summary.total_contigs} contigs ≥500bp, "
                   f"{summary.callable_bases:,} callable bases, {summary.mean_depth:.1f}x mean depth")

        return summary

    except Exception as e:
        if isinstance(e, DepthAnalysisError):
            raise
        raise DepthAnalysisError(f"Failed to parse mosdepth summary: {e}")

run_mosdepth(bam_path, output_prefix, mapq_threshold=30, threads=4)

Run mosdepth on a BAM file to compute depth statistics.

Parameters:

Name	Type	Description	Default
bam_path	Path	Path to input BAM file	required
output_prefix	Path	Output prefix for mosdepth files	required
mapq_threshold	int	Minimum mapping quality (default: 30)	30
threads	int	Number of threads to use	4

Returns:

Type	Description
Path	Path to the generated summary file (.mosdepth.summary.txt)

Raises:

Type	Description
DepthAnalysisError	If mosdepth is not available or execution fails
FileNotFoundError	If BAM file doesn't exist

Source code in src/ssiamb/depth.py

def run_mosdepth(
    bam_path: Path,
    output_prefix: Path,
    mapq_threshold: int = 30,
    threads: int = 4
) -> Path:
    """
    Run mosdepth on a BAM file to compute depth statistics.

    Args:
        bam_path: Path to input BAM file
        output_prefix: Output prefix for mosdepth files
        mapq_threshold: Minimum mapping quality (default: 30)
        threads: Number of threads to use

    Returns:
        Path to the generated summary file (.mosdepth.summary.txt)

    Raises:
        DepthAnalysisError: If mosdepth is not available or execution fails
        FileNotFoundError: If BAM file doesn't exist
    """
    if not bam_path.exists():
        raise FileNotFoundError(f"BAM file not found: {bam_path}")

    if not check_mosdepth_available():
        raise DepthAnalysisError("mosdepth not found in PATH")

    # Construct mosdepth command
    cmd = [
        'mosdepth',
        '--mapq', str(mapq_threshold),
        '--no-per-base',  # Skip per-base output for efficiency
        '--threads', str(threads),
        str(output_prefix),
        str(bam_path)
    ]

    logger.info(f"Running mosdepth: {' '.join(cmd)}")

    try:
        result = subprocess.run(
            cmd, 
            capture_output=True, 
            text=True, 
            check=True
        )

        # Log any stderr output for debugging
        if result.stderr.strip():
            logger.debug(f"mosdepth stderr: {result.stderr.strip()}")

    except subprocess.CalledProcessError as e:
        error_msg = f"mosdepth failed with return code {e.returncode}"
        if e.stderr:
            error_msg += f": {e.stderr.strip()}"
        raise DepthAnalysisError(error_msg)

    # Check that expected output file was created
    summary_file = Path(str(output_prefix) + '.mosdepth.summary.txt')
    if not summary_file.exists():
        raise DepthAnalysisError(f"Expected mosdepth summary file not created: {summary_file}")

    logger.info(f"mosdepth completed successfully: {summary_file}")
    return summary_file

Quality Control

Quality control policies and warning thresholds for ssiamb.

This module implements QC policies that warn but do not fail runs, as specified in spec.md §8.

QCThresholds(min_breadth_10x=0.8, min_callable_bases=1000000, min_mapping_rate_ref=0.7) dataclass

QC warning thresholds as specified in spec.md §8.

QCWarning(metric, value, threshold, message) dataclass

A QC warning with message and metric values.

check_qc_metrics(breadth_10x, callable_bases, mapping_rate=None, mode=None, thresholds=None)

Check QC metrics against warning thresholds.

Parameters:

Name	Type	Description	Default
breadth_10x	float	Breadth of coverage at 10x depth (0.0-1.0)	required
callable_bases	int	Number of callable bases	required
mapping_rate	Optional[float]	Mapping rate (0.0-1.0), required for ref mode	None
mode	Optional[Mode]	Analysis mode (used to determine if mapping rate applies)	None
thresholds	Optional[QCThresholds]	QC thresholds (uses defaults if None)	None

Returns:

Type	Description
List[QCWarning]	List of QC warnings

Source code in src/ssiamb/qc.py

def check_qc_metrics(
    breadth_10x: float,
    callable_bases: int,
    mapping_rate: Optional[float] = None,
    mode: Optional[Mode] = None,
    thresholds: Optional[QCThresholds] = None
) -> List[QCWarning]:
    """
    Check QC metrics against warning thresholds.

    Args:
        breadth_10x: Breadth of coverage at 10x depth (0.0-1.0)
        callable_bases: Number of callable bases
        mapping_rate: Mapping rate (0.0-1.0), required for ref mode
        mode: Analysis mode (used to determine if mapping rate applies)
        thresholds: QC thresholds (uses defaults if None)

    Returns:
        List of QC warnings
    """
    if thresholds is None:
        thresholds = QCThresholds()

    warnings = []

    # Check breadth_10x threshold
    if breadth_10x < thresholds.min_breadth_10x:
        warnings.append(QCWarning(
            metric="breadth_10x",
            value=breadth_10x,
            threshold=thresholds.min_breadth_10x,
            message=f"Low breadth of coverage at 10x: {breadth_10x:.3f} < {thresholds.min_breadth_10x:.3f}"
        ))

    # Check callable bases threshold
    if callable_bases < thresholds.min_callable_bases:
        warnings.append(QCWarning(
            metric="callable_bases",
            value=callable_bases,
            threshold=thresholds.min_callable_bases,
            message=f"Low callable bases: {callable_bases:,} < {thresholds.min_callable_bases:,}"
        ))

    # Check mapping rate threshold (ref mode only)
    if mode == Mode.REF and mapping_rate is not None:
        if mapping_rate < thresholds.min_mapping_rate_ref:
            warnings.append(QCWarning(
                metric="mapping_rate",
                value=mapping_rate,
                threshold=thresholds.min_mapping_rate_ref,
                message=f"Low mapping rate in ref mode: {mapping_rate:.3f} < {thresholds.min_mapping_rate_ref:.3f}"
            ))

    return warnings

format_qc_warnings_for_summary(warnings)

Format QC warnings for inclusion in summary TSV qc_warnings field.

Parameters:

Name	Type	Description	Default
warnings	List[QCWarning]	List of QC warnings	required

Returns:

Type	Description
str	Semicolon-separated string of warning codes, or empty string if no warnings

Source code in src/ssiamb/qc.py

def format_qc_warnings_for_summary(warnings: List[QCWarning]) -> str:
    """
    Format QC warnings for inclusion in summary TSV qc_warnings field.

    Args:
        warnings: List of QC warnings

    Returns:
        Semicolon-separated string of warning codes, or empty string if no warnings
    """
    if not warnings:
        return ""

    # Create concise warning codes
    codes = []
    for warning in warnings:
        if warning.metric == "breadth_10x":
            codes.append(f"low_breadth_{warning.value:.2f}")
        elif warning.metric == "callable_bases":
            codes.append(f"low_callable_{warning.value:.0e}")
        elif warning.metric == "mapping_rate":
            codes.append(f"low_maprate_{warning.value:.2f}")

    return ";".join(codes)

log_qc_warnings(warnings)

Log QC warnings to the logger.

Parameters:

Name	Type	Description	Default
warnings	List[QCWarning]	List of QC warnings to log	required

Source code in src/ssiamb/qc.py

def log_qc_warnings(warnings: List[QCWarning]) -> None:
    """
    Log QC warnings to the logger.

    Args:
        warnings: List of QC warnings to log
    """
    if not warnings:
        logger.info("All QC metrics passed warning thresholds")
        return

    logger.warning(f"QC warnings detected ({len(warnings)} issues):")
    for warning in warnings:
        logger.warning(f"  - {warning.message}")

VCF Operations

VCF operations module for ssiamb.

This module implements VCF normalization, atomization, MAF extraction, variant classification, and grid-based counting for ambiguous sites.

AmbigGrid(dp_cap=100)

100×51 cumulative grid for ambiguous site counting.

Tracks sites by depth (0-100, with capping) and MAF bins (0-50, representing 0.00-0.50). MAF binning: bin = floor(100 * MAF), capped at 50.

Initialize grid.

Parameters:

Name	Type	Description	Default
dp_cap	int	Maximum depth value (higher values are capped)	100

Source code in src/ssiamb/vcf_ops.py

def __init__(self, dp_cap: int = 100):
    """
    Initialize grid.

    Args:
        dp_cap: Maximum depth value (higher values are capped)
    """
    self.dp_cap = dp_cap
    self.maf_bins = 51  # 0.00 to 0.50 in 0.01 increments

    # Grid: [depth][maf_bin] = count
    self.grid = np.zeros((dp_cap + 1, self.maf_bins), dtype=int)

add_site(depth, maf)

Add a site to the grid.

Parameters:

Name	Type	Description	Default
depth	int	Site depth (will be capped at dp_cap)	required
maf	float	Minor allele frequency (0.0 to 1.0)	required

Source code in src/ssiamb/vcf_ops.py

def add_site(self, depth: int, maf: float) -> None:
    """
    Add a site to the grid.

    Args:
        depth: Site depth (will be capped at dp_cap)
        maf: Minor allele frequency (0.0 to 1.0)
    """
    # Cap depth
    if depth > self.dp_cap:
        depth = self.dp_cap
        logger.debug(f"Depth capped at {self.dp_cap}")

    # Calculate MAF bin (floor of 100 * MAF)
    maf_bin = int(np.floor(100 * maf))

    # Cap MAF bin at 50 (representing 0.50)
    if maf_bin >= self.maf_bins:
        maf_bin = self.maf_bins - 1

    # Ensure non-negative
    depth = max(0, depth)
    maf_bin = max(0, maf_bin)

    self.grid[depth, maf_bin] += 1

build_cumulative()

Build cumulative counts matrix.

Returns:

Type	Description
ndarray	Cumulative grid where each cell contains count of sites
ndarray	with depth >= row_index and MAF >= col_index/100

Source code in src/ssiamb/vcf_ops.py

def build_cumulative(self) -> np.ndarray:
    """
    Build cumulative counts matrix.

    Returns:
        Cumulative grid where each cell contains count of sites 
        with depth >= row_index and MAF >= col_index/100
    """
    # Start from bottom-right and work backwards
    cumulative = np.zeros_like(self.grid)

    for dp in range(self.dp_cap, -1, -1):
        for maf_bin in range(self.maf_bins - 1, -1, -1):
            cumulative[dp, maf_bin] = self.grid[dp, maf_bin]

            # Add counts from higher depth (same MAF bin)
            if dp < self.dp_cap:
                cumulative[dp, maf_bin] += cumulative[dp + 1, maf_bin]

            # Add counts from higher MAF bin (same depth)
            if maf_bin < self.maf_bins - 1:
                cumulative[dp, maf_bin] += cumulative[dp, maf_bin + 1]

            # Subtract double-counted intersection
            if dp < self.dp_cap and maf_bin < self.maf_bins - 1:
                cumulative[dp, maf_bin] -= cumulative[dp + 1, maf_bin + 1]

    return cumulative

count_at(dp_min, maf_min)

Count sites meeting minimum thresholds.

Parameters:

Name	Type	Description	Default
dp_min	int	Minimum depth threshold	required
maf_min	float	Minimum MAF threshold	required

Returns:

Type	Description
int	Count of sites with depth >= dp_min and MAF >= maf_min

Source code in src/ssiamb/vcf_ops.py

def count_at(self, dp_min: int, maf_min: float) -> int:
    """
    Count sites meeting minimum thresholds.

    Args:
        dp_min: Minimum depth threshold
        maf_min: Minimum MAF threshold

    Returns:
        Count of sites with depth >= dp_min and MAF >= maf_min
    """
    cumulative = self.build_cumulative()

    # Convert MAF to bin (floor of 100 * MAF)
    maf_bin = int(np.floor(100 * maf_min))
    maf_bin = min(maf_bin, self.maf_bins - 1)

    # Cap depth
    dp_min = min(dp_min, self.dp_cap)

    return int(cumulative[dp_min, maf_bin])

to_wide_tsv(output_path)

Write cumulative grid to wide TSV format.

Per spec.md §4.4: "Wide table: rows depth=1..100; columns maf_0..maf_50"

Parameters:

Name	Type	Description	Default
output_path	Path	Output file path	required

Source code in src/ssiamb/vcf_ops.py

def to_wide_tsv(self, output_path: Path) -> None:
    """
    Write cumulative grid to wide TSV format.

    Per spec.md §4.4: "Wide table: rows depth=1..100; columns maf_0..maf_50"

    Args:
        output_path: Output file path
    """
    cumulative = self.build_cumulative()

    # Create header: depth, then MAF columns (0.00, 0.01, ..., 0.50)
    maf_headers = [f"{i/100:.2f}" for i in range(self.maf_bins)]
    header = ["depth"] + maf_headers

    with open(output_path, 'w') as f:
        # Write header
        f.write('\t'.join(header) + '\n')

        # Write data rows - start from depth=1, not depth=0 (per spec §4.4)
        for dp in range(1, self.dp_cap + 1):
            row = [str(dp)] + [str(cumulative[dp, maf_bin]) for maf_bin in range(self.maf_bins)]
            f.write('\t'.join(row) + '\n')

    logger.info(f"Cumulative grid written to {output_path}")

SiteRecord(chrom, pos, ref, alt, variant_class, depth, maf, original_filter) dataclass

Record for a genomic site with variant information.

VCFNormalizationResult(normalized_vcf_path, success, error_message=None, records_processed=0) dataclass

Result of VCF normalization operation.

VCFOperationError

Bases: Exception

Raised when VCF operations fail.

VariantClass

Bases: Enum

Variant classification types.

check_vcf_tools()

Check if required VCF processing tools are available.

Returns:

Type	Description
bool	True if all required tools are available, False otherwise

Source code in src/ssiamb/vcf_ops.py

def check_vcf_tools() -> bool:
    """
    Check if required VCF processing tools are available.

    Returns:
        True if all required tools are available, False otherwise
    """
    required_tools = ["bcftools", "bgzip", "tabix"]
    return all(shutil.which(tool) is not None for tool in required_tools)

classify_variant(ref, alt)

Classify variant type based on REF and ALT alleles.

Parameters:

Name	Type	Description	Default
ref	str	Reference allele	required
alt	str	Alternative allele	required

Returns:

Type	Description
VariantClass	VariantClass enum value

Source code in src/ssiamb/vcf_ops.py

def classify_variant(ref: str, alt: str) -> VariantClass:
    """
    Classify variant type based on REF and ALT alleles.

    Args:
        ref: Reference allele
        alt: Alternative allele

    Returns:
        VariantClass enum value
    """
    # Skip symbolic alleles and IUPAC codes
    if any(c in ref + alt for c in ["<", ">", "*", "N"]):
        return VariantClass.UNKNOWN

    # Skip non-ATCG characters
    valid_bases = set("ATCG")
    if not (set(ref.upper()) <= valid_bases and set(alt.upper()) <= valid_bases):
        return VariantClass.UNKNOWN

    ref_len = len(ref)
    alt_len = len(alt)

    if ref_len == alt_len == 1:
        return VariantClass.SNV
    elif alt_len > ref_len:
        return VariantClass.INS
    elif alt_len < ref_len:
        return VariantClass.DEL
    else:
        # Complex variants
        return VariantClass.UNKNOWN

count_ambiguous_sites(vcf_path, dp_min, maf_min, dp_cap=100, included_contigs=None, variant_classes=None)

Count ambiguous sites from normalized VCF.

Parameters:

Name	Type	Description	Default
vcf_path	Path	Path to normalized VCF file	required
dp_min	int	Minimum depth threshold	required
maf_min	float	Minimum MAF threshold	required
dp_cap	int	Maximum depth (higher values capped)	100
included_contigs	Optional[Set[str]]	Set of contigs to include	None
variant_classes	Optional[List[VariantClass]]	List of variant classes to count (default: [SNV])	None

Returns:

Type	Description
Tuple[int, AmbigGrid]	Tuple of (ambiguous_count, grid_object)

Source code in src/ssiamb/vcf_ops.py

def count_ambiguous_sites(vcf_path: Path, dp_min: int, maf_min: float, 
                         dp_cap: int = 100, 
                         included_contigs: Optional[Set[str]] = None,
                         variant_classes: Optional[List[VariantClass]] = None) -> Tuple[int, AmbigGrid]:
    """
    Count ambiguous sites from normalized VCF.

    Args:
        vcf_path: Path to normalized VCF file
        dp_min: Minimum depth threshold
        maf_min: Minimum MAF threshold  
        dp_cap: Maximum depth (higher values capped)
        included_contigs: Set of contigs to include
        variant_classes: List of variant classes to count (default: [SNV])

    Returns:
        Tuple of (ambiguous_count, grid_object)
    """
    if variant_classes is None:
        variant_classes = [VariantClass.SNV]

    grid = AmbigGrid(dp_cap=dp_cap)

    # Process all sites and build grid
    for site in parse_vcf_sites(vcf_path, included_contigs):
        if site.variant_class in variant_classes:
            grid.add_site(site.depth, site.maf)

    # Count ambiguous sites
    ambiguous_count = grid.count_at(dp_min, maf_min)

    logger.info(f"Found {ambiguous_count} ambiguous sites (dp>={dp_min}, maf>={maf_min})")

    return ambiguous_count, grid

emit_bed(normalized_vcf_path, output_path, dp_min, maf_min, sample_name, included_contigs=None)

Emit BED output with ambiguous sites in 0-based half-open coordinates.

BED format columns: chrom, start, end, name, score, strand, sample, variant_class, ref, alt, maf, dp, maf_bin, dp_cap

Parameters:

Name	Type	Description	Default
normalized_vcf_path	Path	Path to normalized input VCF	required
output_path	Path	Output BED path (will be bgzipped)	required
dp_min	int	Minimum depth threshold	required
maf_min	float	Minimum MAF threshold	required
sample_name	str	Sample name	required
included_contigs	Optional[Set[str]]	Set of contigs to include	None

Returns:

Type	Description
Path	Path to compressed, indexed BED file

Raises:

Type	Description
VCFOperationError	If emission fails

Source code in src/ssiamb/vcf_ops.py

def emit_bed(
    normalized_vcf_path: Path,
    output_path: Path,
    dp_min: int,
    maf_min: float,
    sample_name: str,
    included_contigs: Optional[Set[str]] = None
) -> Path:
    """
    Emit BED output with ambiguous sites in 0-based half-open coordinates.

    BED format columns:
    chrom, start, end, name, score, strand, sample, variant_class, ref, alt, maf, dp, maf_bin, dp_cap

    Args:
        normalized_vcf_path: Path to normalized input VCF
        output_path: Output BED path (will be bgzipped)
        dp_min: Minimum depth threshold
        maf_min: Minimum MAF threshold
        sample_name: Sample name
        included_contigs: Set of contigs to include

    Returns:
        Path to compressed, indexed BED file

    Raises:
        VCFOperationError: If emission fails
    """
    try:
        # Ensure output path has .bed.gz extension
        if not str(output_path).endswith('.bed.gz'):
            output_path = output_path.with_suffix('.bed.gz')

        output_path.parent.mkdir(parents=True, exist_ok=True)

        bed_records = []

        with pysam.VariantFile(str(normalized_vcf_path)) as input_vcf:
            for record in input_vcf:
                # Filter by included contigs
                if included_contigs is not None and record.chrom not in included_contigs:
                    continue

                # Skip if no ALT alleles
                if not record.alts:
                    continue

                # Process each ALT allele
                for alt_allele in record.alts:
                    # Skip if ref or alt is None
                    if record.ref is None or alt_allele is None:
                        continue

                    variant_class = classify_variant(record.ref, alt_allele)

                    # Skip unknown/symbolic variants
                    if variant_class == VariantClass.UNKNOWN:
                        continue

                    # Extract MAF and depth
                    maf = extract_maf_from_record(record)
                    if maf is None:
                        continue

                    # Get depth
                    depth = 0
                    if "DP" in record.info:
                        depth = record.info["DP"]
                    elif "AD" in record.format:
                        samples = list(record.samples)
                        if samples:
                            sample = record.samples[samples[0]]
                            ad = sample.get("AD")
                            if ad is not None:
                                depth = sum(ad)

                    # Check if passes thresholds
                    passes_thresholds = depth >= dp_min and maf >= maf_min

                    if passes_thresholds:
                        # Convert to 0-based coordinates
                        # VCF is 1-based, BED is 0-based half-open
                        start = record.start  # pysam already converts to 0-based

                        if variant_class == VariantClass.SNV:
                            # SNV: 1bp interval
                            end = start + 1
                        elif variant_class == VariantClass.INS:
                            # Insertion: 1bp anchor position
                            end = start + 1
                        elif variant_class == VariantClass.DEL:
                            # Deletion: span the deleted bases
                            if record.ref is not None:
                                end = start + len(record.ref)
                            else:
                                end = start + 1
                        else:
                            end = start + 1

                        # Create BED record
                        name = f"{record.chrom}:{record.pos}_{record.ref}>{alt_allele}"
                        score = min(1000, int(1000 * maf))  # Scale MAF to 0-1000
                        maf_bin = int(np.floor(100 * maf))

                        bed_record = [
                            record.chrom,           # chrom
                            str(start),             # start (0-based)
                            str(end),               # end (0-based half-open)
                            name,                   # name
                            str(score),             # score (0-1000)
                            ".",                    # strand (not applicable)
                            sample_name,            # sample
                            variant_class.value,    # variant_class
                            record.ref,             # ref
                            alt_allele,             # alt
                            f"{maf:.4f}",          # maf
                            str(depth),             # dp
                            str(maf_bin),          # maf_bin
                            "100"                   # dp_cap
                        ]

                        bed_records.append(bed_record)

        # Sort records by chromosome and position
        bed_records.sort(key=lambda x: (x[0], int(x[1])))

        # Write BED file
        temp_bed = output_path.with_suffix('.bed')
        with open(temp_bed, 'w') as f:
            # Write header comment
            f.write("# chrom\tstart\tend\tname\tscore\tstrand\tsample\tvariant_class\tref\talt\tmaf\tdp\tmaf_bin\tdp_cap\n")

            # Write records
            for record in bed_records:
                f.write('\t'.join(record) + '\n')

        # Compress with bgzip
        bgzip_cmd = ["bgzip", "-c", str(temp_bed)]
        with open(output_path, 'w') as f:
            subprocess.run(bgzip_cmd, stdout=f, check=True)

        # Remove temporary file
        temp_bed.unlink()

        # Index with tabix
        tabix_cmd = ["tabix", "-p", "bed", str(output_path)]
        subprocess.run(tabix_cmd, capture_output=True, text=True, check=True)

        logger.info(f"Emitted BED with {len(bed_records)} ambiguous sites: {output_path}")
        return output_path

    except Exception as e:
        error_msg = f"BED emission failed: {str(e)}"
        logger.error(error_msg)
        raise VCFOperationError(error_msg) from e

emit_matrix(grid, output_path, sample_name)

Emit variant matrix as compressed TSV.

Parameters:

Name	Type	Description	Default
grid	AmbigGrid	Populated AmbigGrid with variant counts	required
output_path	Path	Output path (will be ensured to end with .tsv.gz)	required
sample_name	str	Sample name for logging	required

Returns:

Type	Description
Path	Path to compressed matrix file

Raises:

Type	Description
VCFOperationError	If matrix emission fails

Source code in src/ssiamb/vcf_ops.py

def emit_matrix(
    grid: AmbigGrid,
    output_path: Path,
    sample_name: str
) -> Path:
    """
    Emit variant matrix as compressed TSV.

    Args:
        grid: Populated AmbigGrid with variant counts
        output_path: Output path (will be ensured to end with .tsv.gz)
        sample_name: Sample name for logging

    Returns:
        Path to compressed matrix file

    Raises:
        VCFOperationError: If matrix emission fails
    """
    try:
        # Ensure output path has .tsv.gz extension
        if not str(output_path).endswith('.tsv.gz'):
            output_path = output_path.with_suffix('.tsv.gz')

        output_path.parent.mkdir(parents=True, exist_ok=True)

        # Write to temporary uncompressed file first
        temp_tsv = output_path.with_suffix('.tsv')
        grid.to_wide_tsv(temp_tsv)

        # Compress with gzip
        import gzip
        with open(temp_tsv, 'rb') as f_in:
            with gzip.open(output_path, 'wb') as f_out:
                f_out.write(f_in.read())

        # Remove temporary file
        temp_tsv.unlink()

        logger.info(f"Emitted variant matrix for {sample_name}: {output_path}")
        return output_path

    except Exception as e:
        error_msg = f"Matrix emission failed: {str(e)}"
        logger.error(error_msg)
        raise VCFOperationError(error_msg) from e

emit_multiqc(sample_name, ambiguous_snv_count, breadth_10x, callable_bases, genome_length, dp_min, maf_min, mapper, caller, mode, output_path)

Emit MultiQC-compatible metrics TSV.

As per spec.md §4.6: sample, ambiguous_snv_count, ambiguous_snv_per_mb, breadth_10x, callable_bases, genome_length, dp_min, maf_min, mapper, caller, mode

Parameters:

Name	Type	Description	Default
sample_name	str	Sample identifier	required
ambiguous_snv_count	int	Count of ambiguous SNVs	required
breadth_10x	float	Fraction of genome with ≥10x coverage	required
callable_bases	int	Number of callable bases	required
genome_length	int	Total genome length	required
dp_min	int	Minimum depth threshold used	required
maf_min	float	Minimum MAF threshold used	required
mapper	str	Mapper tool name	required
caller	str	Caller tool name	required
mode	str	Analysis mode (self or ref)	required
output_path	Path	Output path (will be ensured to end with .tsv)	required

Returns:

Type	Description
Path	Path to emitted MultiQC TSV file

Raises:

Type	Description
VCFOperationError	If MultiQC emission fails

Source code in src/ssiamb/vcf_ops.py

def emit_multiqc(
    sample_name: str,
    ambiguous_snv_count: int,
    breadth_10x: float,
    callable_bases: int,
    genome_length: int,
    dp_min: int,
    maf_min: float,
    mapper: str,
    caller: str,
    mode: str,
    output_path: Path
) -> Path:
    """
    Emit MultiQC-compatible metrics TSV.

    As per spec.md §4.6: sample, ambiguous_snv_count, ambiguous_snv_per_mb, breadth_10x, 
    callable_bases, genome_length, dp_min, maf_min, mapper, caller, mode

    Args:
        sample_name: Sample identifier
        ambiguous_snv_count: Count of ambiguous SNVs 
        breadth_10x: Fraction of genome with ≥10x coverage
        callable_bases: Number of callable bases
        genome_length: Total genome length
        dp_min: Minimum depth threshold used
        maf_min: Minimum MAF threshold used
        mapper: Mapper tool name
        caller: Caller tool name
        mode: Analysis mode (self or ref)
        output_path: Output path (will be ensured to end with .tsv)

    Returns:
        Path to emitted MultiQC TSV file

    Raises:
        VCFOperationError: If MultiQC emission fails
    """
    try:
        # Ensure output path has .tsv extension
        if not str(output_path).endswith('.tsv'):
            output_path = output_path.with_suffix('.tsv')

        output_path.parent.mkdir(parents=True, exist_ok=True)

        # Calculate ambiguous SNVs per megabase
        # Per spec.md §4.1: ambiguous_snv_per_mb = ambiguous_snv_count * 1e6 / callable_bases
        ambiguous_snv_per_mb = (ambiguous_snv_count * 1_000_000) / callable_bases if callable_bases > 0 else 0.0

        with open(output_path, 'w') as f:
            # Write header
            header = [
                "sample", "ambiguous_snv_count", "ambiguous_snv_per_mb", "breadth_10x",
                "callable_bases", "genome_length", "dp_min", "maf_min", "mapper", "caller", "mode"
            ]
            f.write('\t'.join(header) + '\n')

            # Write data row
            row = [
                sample_name,
                str(ambiguous_snv_count),
                f"{ambiguous_snv_per_mb:.2f}",
                f"{breadth_10x:.4f}",
                str(callable_bases),
                str(genome_length),
                str(dp_min),
                f"{maf_min:.3f}",
                mapper,
                caller,
                mode
            ]
            f.write('\t'.join(row) + '\n')

        logger.info(f"Emitted MultiQC metrics for {sample_name}: {output_path}")
        return output_path

    except Exception as e:
        error_msg = f"MultiQC emission failed: {str(e)}"
        logger.error(error_msg)
        raise VCFOperationError(error_msg) from e

emit_per_contig(normalized_vcf_path, depth_summary_path, output_path, dp_min, maf_min, sample_name, included_contigs=None)

Emit per-contig summary statistics.

Parameters:

Name	Type	Description	Default
normalized_vcf_path	Path	Path to normalized VCF file	required
depth_summary_path	Path	Path to mosdepth summary file	required
output_path	Path	Output TSV path	required
dp_min	int	Minimum depth threshold	required
maf_min	float	Minimum MAF threshold	required
sample_name	str	Sample name	required
included_contigs	Optional[Set[str]]	Set of contigs to include	None

Returns:

Type	Description
Path	Path to per-contig summary file

Raises:

Type	Description
VCFOperationError	If emission fails

Source code in src/ssiamb/vcf_ops.py

def emit_per_contig(
    normalized_vcf_path: Path,
    depth_summary_path: Path,
    output_path: Path,
    dp_min: int,
    maf_min: float,
    sample_name: str,
    included_contigs: Optional[Set[str]] = None
) -> Path:
    """
    Emit per-contig summary statistics.

    Args:
        normalized_vcf_path: Path to normalized VCF file
        depth_summary_path: Path to mosdepth summary file
        output_path: Output TSV path
        dp_min: Minimum depth threshold
        maf_min: Minimum MAF threshold
        sample_name: Sample name
        included_contigs: Set of contigs to include

    Returns:
        Path to per-contig summary file

    Raises:
        VCFOperationError: If emission fails
    """
    try:
        from .depth import parse_mosdepth_summary

        output_path.parent.mkdir(parents=True, exist_ok=True)

        # Parse depth summary for contig-level stats
        depth_summary = parse_mosdepth_summary(depth_summary_path)

        # Count variants per contig
        contig_variant_counts = {}

        with pysam.VariantFile(str(normalized_vcf_path)) as vcf:
            for record in vcf:
                # Filter by included contigs
                if included_contigs is not None and record.chrom not in included_contigs:
                    continue

                # Skip if no ALT alleles
                if not record.alts:
                    continue

                # Initialize contig counts if needed
                if record.chrom not in contig_variant_counts:
                    contig_variant_counts[record.chrom] = {
                        'snv': 0,
                        'indel': 0,
                        'del': 0
                    }

                # Process each ALT allele
                for alt_allele in record.alts:
                    if record.ref is None or alt_allele is None:
                        continue

                    variant_class = classify_variant(record.ref, alt_allele)

                    if variant_class == VariantClass.UNKNOWN:
                        continue

                    # Extract MAF and depth
                    maf = extract_maf_from_record(record)
                    if maf is None:
                        continue

                    # Get depth
                    depth = 0
                    if "DP" in record.info:
                        depth = int(record.info["DP"])
                    elif "AD" in record.format:
                        samples = list(record.samples)
                        if samples:
                            sample = record.samples[samples[0]]
                            ad = sample.get("AD")
                            if ad is not None:
                                depth = sum(int(x) for x in ad)

                    # Check if passes thresholds
                    if depth >= dp_min and maf >= maf_min:
                        if variant_class == VariantClass.SNV:
                            contig_variant_counts[record.chrom]['snv'] += 1
                        elif variant_class == VariantClass.INS:
                            contig_variant_counts[record.chrom]['indel'] += 1
                        elif variant_class == VariantClass.DEL:
                            contig_variant_counts[record.chrom]['del'] += 1

        # Write per-contig summary
        with open(output_path, 'w') as f:
            # Write header
            header = [
                "sample", "contig", "length", "callable_bases_10x", "breadth_10x",
                "ambiguous_snv_count", "ambiguous_indel_count", "ambiguous_del_count",
                "ambiguous_snv_per_mb", "mean_depth"
            ]
            f.write('\t'.join(header) + '\n')

            # Write data for each contig in depth summary
            for contig_stat in depth_summary.contig_stats:
                contig_name = contig_stat.name
                length = contig_stat.length
                callable_bases = contig_stat.bases_covered
                mean_depth = contig_stat.mean_depth
                breadth_10x = contig_stat.breadth_10x

                # Get variant counts for this contig
                counts = contig_variant_counts.get(contig_name, {'snv': 0, 'indel': 0, 'del': 0})
                snv_count = counts['snv']
                indel_count = counts['indel']
                del_count = counts['del']

                # Calculate SNV per Mb - use callable_bases, not total length
                # Per spec.md §4.1: ambiguous_snv_per_mb = ambiguous_snv_count * 1e6 / callable_bases
                snv_per_mb = (snv_count * 1_000_000) / callable_bases if callable_bases > 0 else 0.0

                row = [
                    sample_name,
                    contig_name,
                    str(length),
                    str(callable_bases),
                    f"{breadth_10x:.4f}",
                    str(snv_count),
                    str(indel_count),
                    str(del_count),
                    f"{snv_per_mb:.2f}",
                    f"{mean_depth:.2f}"
                ]
                f.write('\t'.join(row) + '\n')

        logger.info(f"Emitted per-contig summary for {sample_name}: {output_path}")
        return output_path

    except Exception as e:
        error_msg = f"Per-contig summary emission failed: {str(e)}"
        logger.error(error_msg)
        raise VCFOperationError(error_msg) from e

emit_vcf(normalized_vcf_path, output_path, dp_min, maf_min, sample_name, require_pass=False, included_contigs=None)

Emit VCF output with only records passing ambiguous site thresholds.

Parameters:

Name	Type	Description	Default
normalized_vcf_path	Path	Path to normalized input VCF	required
output_path	Path	Output VCF path (will be bgzipped)	required
dp_min	int	Minimum depth threshold	required
maf_min	float	Minimum MAF threshold	required
sample_name	str	Sample name for output	required
require_pass	bool	If True, include ORIG_FILTER info field	False
included_contigs	Optional[Set[str]]	Set of contigs to include	None

Returns:

Type	Description
Path	Path to compressed, indexed VCF file

Raises:

Type	Description
VCFOperationError	If emission fails

Source code in src/ssiamb/vcf_ops.py

def emit_vcf(
    normalized_vcf_path: Path,
    output_path: Path,
    dp_min: int,
    maf_min: float,
    sample_name: str,
    require_pass: bool = False,
    included_contigs: Optional[Set[str]] = None
) -> Path:
    """
    Emit VCF output with only records passing ambiguous site thresholds.

    Args:
        normalized_vcf_path: Path to normalized input VCF
        output_path: Output VCF path (will be bgzipped)
        dp_min: Minimum depth threshold
        maf_min: Minimum MAF threshold
        sample_name: Sample name for output
        require_pass: If True, include ORIG_FILTER info field
        included_contigs: Set of contigs to include

    Returns:
        Path to compressed, indexed VCF file

    Raises:
        VCFOperationError: If emission fails
    """
    try:
        # Ensure output path has .vcf.gz extension
        if not str(output_path).endswith('.vcf.gz'):
            output_path = output_path.with_suffix('.vcf.gz')

        output_path.parent.mkdir(parents=True, exist_ok=True)

        with pysam.VariantFile(str(normalized_vcf_path)) as input_vcf:
            # Create output VCF with enhanced header
            header = input_vcf.header.copy()

            # Add custom INFO fields
            header.add_line('##INFO=<ID=AMBIG,Number=0,Type=Flag,Description="Site passes ambiguous thresholds">')
            header.add_line('##INFO=<ID=MAF,Number=1,Type=Float,Description="Minor allele frequency">')
            header.add_line('##INFO=<ID=MAF_BIN,Number=1,Type=Integer,Description="MAF bin (floor(100*MAF))">')
            header.add_line('##INFO=<ID=DP_CAP,Number=1,Type=Integer,Description="Depth cap used (100)">')
            header.add_line('##INFO=<ID=VARIANT_CLASS,Number=1,Type=String,Description="Variant class (SNV/INS/DEL)">')

            if not require_pass:
                header.add_line('##INFO=<ID=ORIG_FILTER,Number=1,Type=String,Description="Original FILTER value">')

            # Add FILTER definitions
            header.add_line('##FILTER=<ID=PASS,Description="All filters passed">')

            with pysam.VariantFile(str(output_path), 'w', header=header) as output_vcf:
                records_written = 0

                for record in input_vcf:
                    # Filter by included contigs
                    if included_contigs is not None and record.chrom not in included_contigs:
                        continue

                    # Skip if no ALT alleles
                    if not record.alts:
                        continue

                    # Process each ALT allele (though normalization should make these single)
                    for alt_allele in record.alts:
                        # Skip if ref or alt is None
                        if record.ref is None or alt_allele is None:
                            continue

                        variant_class = classify_variant(record.ref, alt_allele)

                        # Skip unknown/symbolic variants
                        if variant_class == VariantClass.UNKNOWN:
                            continue

                        # Extract MAF and depth
                        maf = extract_maf_from_record(record)
                        if maf is None:
                            continue

                        # Get depth
                        depth = 0
                        if "DP" in record.info:
                            depth = record.info["DP"]
                        elif "AD" in record.format:
                            samples = list(record.samples)
                            if samples:
                                sample = record.samples[samples[0]]
                                ad = sample.get("AD")
                                if ad is not None:
                                    depth = sum(ad)

                        # Check if passes thresholds
                        passes_thresholds = depth >= dp_min and maf >= maf_min

                        if passes_thresholds:
                            # Create output record
                            new_record = output_vcf.new_record(
                                contig=record.chrom,
                                start=record.start,
                                stop=record.stop,
                                alleles=(record.ref, alt_allele),
                                id=record.id,
                                qual=record.qual
                            )

                            # Copy existing INFO fields
                            for key, value in record.info.items():
                                new_record.info[key] = value

                            # Add our custom INFO fields
                            new_record.info["AMBIG"] = True
                            new_record.info["MAF"] = round(maf, 4)
                            new_record.info["MAF_BIN"] = int(np.floor(100 * maf))
                            new_record.info["DP_CAP"] = 100
                            new_record.info["VARIANT_CLASS"] = variant_class.value

                            if not require_pass:
                                filter_keys = list(record.filter.keys())
                                orig_filter = ";".join(str(k) for k in filter_keys) if filter_keys else "PASS"
                                new_record.info["ORIG_FILTER"] = orig_filter

                            # Set FILTER to PASS
                            new_record.filter.clear()
                            new_record.filter.add("PASS")

                            # Copy FORMAT fields
                            for sample in record.samples:
                                for key in record.format:
                                    new_record.samples[sample][key] = record.samples[sample][key]

                            output_vcf.write(new_record)
                            records_written += 1

        # Index the compressed VCF
        tabix_cmd = ["tabix", "-p", "vcf", str(output_path)]
        subprocess.run(tabix_cmd, capture_output=True, text=True, check=True)

        logger.info(f"Emitted VCF with {records_written} ambiguous sites: {output_path}")
        return output_path

    except Exception as e:
        error_msg = f"VCF emission failed: {str(e)}"
        logger.error(error_msg)
        raise VCFOperationError(error_msg) from e

extract_maf_from_record(record)

Extract minor allele frequency from VCF record using precedence: AD → DP4 → AF.

Parameters:

Name	Type	Description	Default
record		pysam VariantRecord	required

Returns:

Type	Description
Optional[float]	MAF as float between 0 and 1, or None if cannot be determined

Source code in src/ssiamb/vcf_ops.py

def extract_maf_from_record(record) -> Optional[float]:
    """
    Extract minor allele frequency from VCF record using precedence: AD → DP4 → AF.

    Args:
        record: pysam VariantRecord

    Returns:
        MAF as float between 0 and 1, or None if cannot be determined
    """
    try:
        # Precedence 1: FORMAT/AD (Allelic Depth)
        if "AD" in record.format:
            samples = list(record.samples)
            if samples:
                sample = record.samples[samples[0]]  # Single sample
                ad = sample.get("AD")
                dp = sample.get("DP")  # Total depth

                if ad is not None and dp is not None:
                    # Handle different AD formats
                    if isinstance(ad, (list, tuple)) and len(ad) >= 2:
                        # Standard format: [ref_depth, alt_depth, ...]
                        ref_depth = int(ad[0])
                        alt_depth = sum(int(x) for x in ad[1:])  # Sum all alt alleles
                        total_depth = ref_depth + alt_depth
                    elif isinstance(ad, (int, str)):
                        # BBTools format: AD is just alt_depth, use DP for total
                        alt_depth = int(ad)
                        total_depth = int(dp)
                        ref_depth = total_depth - alt_depth
                    else:
                        # Skip if AD format is unexpected
                        alt_depth = ref_depth = total_depth = 0

                    if total_depth > 0 and ref_depth >= 0 and alt_depth >= 0:
                        ref_freq = ref_depth / total_depth
                        alt_freq = alt_depth / total_depth
                        return min(ref_freq, alt_freq)  # MAF is the minor frequency

        # Precedence 2: INFO/DP4 (legacy format from bcftools/samtools)
        if "DP4" in record.info:
            dp4 = record.info["DP4"]
            if len(dp4) == 4:
                ref_depth = int(dp4[0]) + int(dp4[1])  # Forward + reverse ref
                alt_depth = int(dp4[2]) + int(dp4[3])  # Forward + reverse alt
                total_depth = ref_depth + alt_depth
                if total_depth > 0:
                    ref_freq = ref_depth / total_depth
                    alt_freq = alt_depth / total_depth
                    return min(ref_freq, alt_freq)

        # Precedence 3: INFO/AF (Allele Frequency)
        if "AF" in record.info:
            af = record.info["AF"]
            if isinstance(af, (list, tuple)):
                # Multi-allelic: compute site-level MAF
                if len(af) > 0:
                    max_af = max(af)
                    ref_freq = 1.0 - sum(af)
                    return min(ref_freq, max_af)
            else:
                # Single alt allele
                af_val = float(af)
                ref_freq = 1.0 - af_val
                return min(ref_freq, af_val)

        # No MAF data available
        warnings.warn(f"No MAF data found for variant at {record.chrom}:{record.pos}")
        return None

    except (ValueError, TypeError, KeyError) as e:
        warnings.warn(f"Error extracting MAF from {record.chrom}:{record.pos}: {e}")
        return None

normalize_and_split(vcf_in, reference, output_dir=None)

Normalize VCF to reference and decompose into primitive records.

Uses bcftools norm to: 1. Normalize variants to reference (-f REF) 2. Split multiallelic sites (-m -both)
3. Atomize variants (--atomize) 4. Compress and index output

Parameters:

Name	Type	Description	Default
vcf_in	Path	Input VCF file path	required
reference	Path	Reference FASTA file path	required
output_dir	Optional[Path]	Output directory (defaults to input directory)	None

Returns:

Type	Description
Path	Path to normalized, compressed VCF file

Raises:

Type	Description
VCFOperationError	If normalization fails

Source code in src/ssiamb/vcf_ops.py

def normalize_and_split(vcf_in: Path, reference: Path, 
                       output_dir: Optional[Path] = None) -> Path:
    """
    Normalize VCF to reference and decompose into primitive records.

    Uses bcftools norm to:
    1. Normalize variants to reference (-f REF)
    2. Split multiallelic sites (-m -both)  
    3. Atomize variants (--atomize)
    4. Compress and index output

    Args:
        vcf_in: Input VCF file path
        reference: Reference FASTA file path
        output_dir: Output directory (defaults to input directory)

    Returns:
        Path to normalized, compressed VCF file

    Raises:
        VCFOperationError: If normalization fails
    """
    if not check_vcf_tools():
        raise VCFOperationError(
            "Required VCF tools not found. Need: bcftools, bgzip, tabix"
        )

    if not vcf_in.exists():
        raise VCFOperationError(f"Input VCF file not found: {vcf_in}")

    if not reference.exists():
        raise VCFOperationError(f"Reference file not found: {reference}")

    # Determine output path
    if output_dir is None:
        output_dir = vcf_in.parent
    output_dir.mkdir(parents=True, exist_ok=True)

    normalized_vcf = output_dir / f"{vcf_in.stem}.normalized.vcf.gz"

    try:
        # Build bcftools norm command
        cmd = [
            "bcftools", "norm",
            "-f", str(reference),
            "-m", "-both",
            "--atomize",
            "-cw",  # Check REF alleles and warn about issues (instead of error)
            "-d", "exact",  # Remove exact duplicates (replaces deprecated -D)
            "-O", "z",  # Compress output
            "-o", str(normalized_vcf),
            str(vcf_in)
        ]

        logger.debug(f"Running bcftools norm: {' '.join(cmd)}")

        result = subprocess.run(
            cmd,
            capture_output=True,
            text=True,
            check=True
        )

        if result.stderr:
            logger.debug(f"bcftools norm stderr: {result.stderr}")

        # Index the compressed VCF
        tabix_cmd = ["tabix", "-p", "vcf", str(normalized_vcf)]
        logger.debug(f"Indexing VCF: {' '.join(tabix_cmd)}")

        subprocess.run(
            tabix_cmd,
            capture_output=True,
            text=True,
            check=True
        )

        logger.info(f"VCF normalized and indexed: {normalized_vcf}")
        return normalized_vcf

    except subprocess.CalledProcessError as e:
        error_msg = f"VCF normalization failed: {e.stderr if e.stderr else str(e)}"
        logger.error(error_msg)
        raise VCFOperationError(error_msg) from e
    except Exception as e:
        error_msg = f"Unexpected error during VCF normalization: {str(e)}"
        logger.error(error_msg)
        raise VCFOperationError(error_msg) from e

parse_vcf_sites(vcf_path, included_contigs=None)

Parse normalized VCF and yield site records with MAF and variant classification.

Parameters:

Name	Type	Description	Default
vcf_path	Path	Path to VCF file (can be compressed)	required
included_contigs	Optional[Set[str]]	Set of contig names to include (None = include all)	None

Yields:

Type	Description
SiteRecord	SiteRecord objects for each variant site

Source code in src/ssiamb/vcf_ops.py

def parse_vcf_sites(vcf_path: Path, included_contigs: Optional[Set[str]] = None) -> Iterator[SiteRecord]:
    """
    Parse normalized VCF and yield site records with MAF and variant classification.

    Args:
        vcf_path: Path to VCF file (can be compressed)
        included_contigs: Set of contig names to include (None = include all)

    Yields:
        SiteRecord objects for each variant site
    """
    try:
        with pysam.VariantFile(str(vcf_path)) as vcf:
            for record in vcf:
                # Filter by included contigs
                if included_contigs is not None and record.chrom not in included_contigs:
                    continue

                # Skip if no ALT alleles
                if not record.alts:
                    continue

                # For multi-allelic sites (after normalization should be rare),
                # process each ALT allele separately
                for alt_allele in record.alts:
                    variant_class = classify_variant(record.ref, alt_allele)

                    # Skip unknown/symbolic variants
                    if variant_class == VariantClass.UNKNOWN:
                        continue

                    # Extract MAF
                    maf = extract_maf_from_record(record)
                    if maf is None:
                        continue

                    # Get depth from DP field or calculate from AD
                    depth = 0
                    if "DP" in record.info:
                        depth = record.info["DP"]
                    elif "AD" in record.format:
                        samples = list(record.samples)
                        if samples:
                            sample = record.samples[samples[0]]
                            ad = sample.get("AD")
                            if ad is not None:
                                depth = sum(ad)

                    # Get original filter
                    orig_filter = ";".join(record.filter.keys()) if record.filter.keys() else "PASS"

                    yield SiteRecord(
                        chrom=record.chrom,
                        pos=record.pos,
                        ref=record.ref,
                        alt=alt_allele,
                        variant_class=variant_class,
                        depth=depth,
                        maf=maf,
                        original_filter=orig_filter
                    )

    except Exception as e:
        error_msg = f"Error parsing VCF file {vcf_path}: {str(e)}"
        logger.error(error_msg)
        raise VCFOperationError(error_msg) from e

I/O Utilities

I/O utilities for ssiamb.

This module provides file handling utilities including sample name validation, TSV writing with atomic operations, and file hashing.

SampleNameError

Bases: ValueError

Raised when sample name is invalid.

TSVWriteError

Bases: Exception

Raised when TSV writing fails.

compute_md5(file_path)

Compute MD5 hash of a file.

Parameters:

Name	Type	Description	Default
file_path	Path	Path to file	required

Returns:

Type	Description
str	MD5 hash as hex string

Source code in src/ssiamb/io_utils.py

def compute_md5(file_path: Path) -> str:
    """
    Compute MD5 hash of a file.

    Args:
        file_path: Path to file

    Returns:
        MD5 hash as hex string
    """
    hash_md5 = hashlib.md5()
    with open(file_path, "rb") as f:
        for chunk in iter(lambda: f.read(4096), b""):
            hash_md5.update(chunk)
    return hash_md5.hexdigest()

file_lock(file_path)

Context manager for file locking (best-effort cross-platform).

Parameters:

Name	Type	Description	Default
file_path	Path	Path to lock	required

Source code in src/ssiamb/io_utils.py

@contextmanager
def file_lock(file_path: Path):
    """
    Context manager for file locking (best-effort cross-platform).

    Args:
        file_path: Path to lock
    """
    lock_file = file_path.with_suffix(file_path.suffix + ".lock")

    try:
        with open(lock_file, "w") as f:
            try:
                fcntl.flock(f.fileno(), fcntl.LOCK_EX | fcntl.LOCK_NB)
                yield
            except (OSError, IOError):
                # Lock failed - either file is locked or system doesn't support locks
                # For now, just proceed (best effort)
                yield
    except (OSError, IOError):
        # Can't create lock file - just proceed
        yield
    finally:
        try:
            lock_file.unlink(missing_ok=True)
        except (OSError, IOError):
            pass

infer_sample_name(r1, r2=None, vcf=None, bam=None)

Infer sample name from input filenames.

Tries to extract a common prefix from R1/R2 filenames, removing common suffixes like _R1, _R2, .fastq, .gz, etc.

Parameters:

Name	Type	Description	Default
r1	Path	Forward reads file	required
r2	Optional[Path]	Reverse reads file (optional)	None
vcf	Optional[Path]	VCF file (optional)	None
bam	Optional[Path]	BAM file (optional)	None

Returns:

Type	Description
str	Inferred sample name

Raises:

Type	Description
SampleNameError	If sample name cannot be inferred

Source code in src/ssiamb/io_utils.py

def infer_sample_name(
    r1: Path,
    r2: Optional[Path] = None,
    vcf: Optional[Path] = None,
    bam: Optional[Path] = None,
) -> str:
    """
    Infer sample name from input filenames.

    Tries to extract a common prefix from R1/R2 filenames, removing
    common suffixes like _R1, _R2, .fastq, .gz, etc.

    Args:
        r1: Forward reads file
        r2: Reverse reads file (optional)
        vcf: VCF file (optional)
        bam: BAM file (optional)

    Returns:
        Inferred sample name

    Raises:
        SampleNameError: If sample name cannot be inferred
    """
    # Start with R1 filename (remove .gz first if present)
    name = r1.name
    if name.endswith('.gz'):
        name = name[:-3]

    # Remove .fastq extension
    if name.endswith('.fastq'):
        name = name[:-6]
    elif name.endswith('.fq'):
        name = name[:-3]

    # Remove R1 suffix
    for suffix in ["_R1", "_1", ".R1", ".1"]:
        if name.endswith(suffix):
            name = name[:-len(suffix)]
            break

    # If R2 provided, ensure it has a compatible name
    if r2:
        r2_name = r2.name
        if r2_name.endswith('.gz'):
            r2_name = r2_name[:-3]

        # Remove .fastq extension
        if r2_name.endswith('.fastq'):
            r2_name = r2_name[:-6]
        elif r2_name.endswith('.fq'):
            r2_name = r2_name[:-3]

        # Remove R2 suffix
        for suffix in ["_R2", "_2", ".R2", ".2"]:
            if r2_name.endswith(suffix):
                r2_name = r2_name[:-len(suffix)]
                break

        # Check if R1 and R2 have the same base name
        if name != r2_name:
            raise SampleNameError(
                f"Cannot infer sample name: R1 stem '{name}' and "
                f"R2 stem '{r2_name}' do not match"
            )

    # Validate the inferred name
    if not name:
        raise SampleNameError("Cannot infer sample name: filename too short")

    try:
        validate_sample_name(name)
        return name
    except SampleNameError as e:
        raise SampleNameError(f"Cannot infer valid sample name from '{r1}': {e}")

validate_sample_name(sample)

Validate sample name according to ssiamb rules.

Parameters:

Name	Type	Description	Default
sample	str	Sample name to validate	required

Returns:

Type	Description
str	Validated sample name

Raises:

Type	Description
SampleNameError	If sample name is invalid

Source code in src/ssiamb/io_utils.py

def validate_sample_name(sample: str) -> str:
    """
    Validate sample name according to ssiamb rules.

    Args:
        sample: Sample name to validate

    Returns:
        Validated sample name

    Raises:
        SampleNameError: If sample name is invalid
    """
    if not sample:
        raise SampleNameError("Sample name cannot be empty")

    if not SAMPLE_NAME_PATTERN.match(sample):
        raise SampleNameError(
            f"Sample name '{sample}' is invalid. "
            "Must be 1-64 characters containing only letters, numbers, "
            "periods, underscores, and hyphens."
        )

    return sample

write_tsv_summary(output_path, rows, mode=TSVMode.OVERWRITE)

Write summary rows to TSV file with atomic operations.

Parameters:

Name	Type	Description	Default
output_path	Path	Output TSV file path	required
rows	List[SummaryRow]	List of SummaryRow objects to write	required
mode	TSVMode	Write mode (overwrite, append, or fail if exists)	OVERWRITE

Raises:

Type	Description
TSVWriteError	If writing fails
FileExistsError	If file exists and mode is FAIL

Source code in src/ssiamb/io_utils.py

def write_tsv_summary(
    output_path: Path,
    rows: List[SummaryRow],
    mode: TSVMode = TSVMode.OVERWRITE,
) -> None:
    """
    Write summary rows to TSV file with atomic operations.

    Args:
        output_path: Output TSV file path
        rows: List of SummaryRow objects to write
        mode: Write mode (overwrite, append, or fail if exists)

    Raises:
        TSVWriteError: If writing fails
        FileExistsError: If file exists and mode is FAIL
    """
    if not rows:
        raise TSVWriteError("No rows to write")

    # Check file existence and mode
    file_exists = output_path.exists()

    if mode == TSVMode.FAIL and file_exists:
        raise FileExistsError(f"Output file already exists: {output_path}")

    # Determine if we need to write header
    write_header = (mode == TSVMode.OVERWRITE) or (mode == TSVMode.APPEND and not file_exists)

    # For atomic writes, use temporary file
    with file_lock(output_path):
        if mode == TSVMode.APPEND and file_exists:
            # Direct append
            with open(output_path, "a", newline="") as f:
                writer = csv.DictWriter(f, fieldnames=rows[0].to_dict().keys(), delimiter="\t")
                for row in rows:
                    writer.writerow(row.to_dict())
        else:
            # Atomic write using temporary file
            temp_file = None
            try:
                # Create temporary file in same directory
                temp_fd, temp_path = tempfile.mkstemp(
                    suffix=".tmp",
                    prefix=output_path.name + ".",
                    dir=output_path.parent,
                )
                temp_file = Path(temp_path)

                with os.fdopen(temp_fd, "w", newline="") as f:
                    writer = csv.DictWriter(f, fieldnames=rows[0].to_dict().keys(), delimiter="\t")

                    if write_header:
                        writer.writeheader()

                    for row in rows:
                        writer.writerow(row.to_dict())

                # Atomic rename
                temp_file.replace(output_path)
                temp_file = None  # Successfully moved

            except Exception as e:
                if temp_file and temp_file.exists():
                    temp_file.unlink()
                raise TSVWriteError(f"Failed to write TSV: {e}") from e

write_tsv_to_stdout(rows)

Write summary rows to stdout in TSV format.

Parameters:

Name	Type	Description	Default
rows	List[SummaryRow]	List of SummaryRow objects to write	required

Source code in src/ssiamb/io_utils.py

def write_tsv_to_stdout(rows: List[SummaryRow]) -> None:
    """
    Write summary rows to stdout in TSV format.

    Args:
        rows: List of SummaryRow objects to write
    """
    import sys

    if not rows:
        return

    writer = csv.DictWriter(sys.stdout, fieldnames=rows[0].to_dict().keys(), delimiter="\t")
    writer.writeheader()

    for row in rows:
        writer.writerow(row.to_dict())

Models

Core data models for ssiamb.

This module defines the primary data structures used throughout the application, including configuration objects, output records, and enums for various options.

Caller

Bases: Enum

Variant calling tool.

DepthTool

Bases: Enum

Depth analysis tool.

Mapper

Bases: Enum

Read mapping tool.

Mode

Bases: Enum

Analysis mode.

OnFail

Bases: Enum

Action when reference resolution fails.

Paths(r1, r2, assembly=None, reference=None, output_dir=Path('.'), sample=None, bam=None, vcf=None, depth_summary=None) dataclass

File paths for input and output.

__post_init__()

Validate paths and resolve output directory.

Source code in src/ssiamb/models.py

def __post_init__(self) -> None:
    """Validate paths and resolve output directory."""
    if not self.r1.exists():
        raise FileNotFoundError(f"R1 file not found: {self.r1}")
    if not self.r2.exists():
        raise FileNotFoundError(f"R2 file not found: {self.r2}")

    # Ensure output directory exists
    self.output_dir.mkdir(parents=True, exist_ok=True)

Provenance(ssiamb_version, command_line, timestamp, input_files=dict(), tool_versions=dict(), hostname='', working_dir='') dataclass

Provenance information for reproducibility.

Tracks tool versions, command lines, input file hashes, etc.

to_dict()

Convert to dictionary for JSON output.

Source code in src/ssiamb/models.py

def to_dict(self) -> Dict[str, Any]:
    """Convert to dictionary for JSON output."""
    return {
        "ssiamb_version": self.ssiamb_version,
        "command_line": self.command_line,
        "timestamp": self.timestamp,
        "input_files": self.input_files,
        "tool_versions": self.tool_versions,
        "hostname": self.hostname,
        "working_dir": self.working_dir,
    }

RunPlan(mode, sample, paths, thresholds, mapper=Mapper.MINIMAP2, caller=Caller.BBTOOLS, depth_tool=DepthTool.MOSDEPTH, bbtools_mem=None, require_pass=False, ref_source='unknown', ref_label='unknown', emit_vcf=False, emit_bed=False, emit_matrix=False, emit_per_contig=False, emit_provenance=False, emit_multiqc=False, to_stdout=False, tsv_mode=TSVMode.OVERWRITE, threads=1, dry_run=False) dataclass

Complete execution plan for a ssiamb run.

This captures all resolved inputs, tools, and output specifications before execution begins.

get_sample_prefix()

Get prefix for output files.

Source code in src/ssiamb/models.py

def get_sample_prefix(self) -> str:
    """Get prefix for output files."""
    return f"{self.sample}_{self.mode.value}"

SummaryRow(sample, mode, mapper, caller, dp_min, maf_min, dp_cap, denom_policy, callable_bases, genome_length, breadth_10x, ambiguous_snv_count, ambiguous_snv_per_mb, ambiguous_indel_count, ambiguous_del_count, ref_label, ref_accession, bracken_species, bracken_frac, bracken_reads, alias_used, reused_bam, reused_vcf, runtime_sec, tool_version) dataclass

Single row of the ambiguous_summary.tsv output.

Based on spec.md §4.1 - Primary Output Format. Schema: sample, mode, mapper, caller, dp_min, maf_min, dp_cap, denom_policy, callable_bases, genome_length, breadth_10x, ambiguous_snv_count, ambiguous_snv_per_mb, ambiguous_indel_count, ambiguous_del_count, ref_label, ref_accession, bracken_species, bracken_frac, bracken_reads, alias_used, reused_bam, reused_vcf, runtime_sec, tool_version

to_dict()

Convert to dictionary for TSV output.

Source code in src/ssiamb/models.py

def to_dict(self) -> Dict[str, Any]:
    """Convert to dictionary for TSV output."""
    return {
        "sample": self.sample,
        "mode": self.mode,
        "mapper": self.mapper,
        "caller": self.caller,
        "dp_min": self.dp_min,
        "maf_min": self.maf_min,
        "dp_cap": self.dp_cap,
        "denom_policy": self.denom_policy,
        "callable_bases": self.callable_bases,
        "genome_length": self.genome_length,
        "breadth_10x": f"{self.breadth_10x:.4f}",
        "ambiguous_snv_count": self.ambiguous_snv_count,
        "ambiguous_snv_per_mb": f"{self.ambiguous_snv_per_mb:.2f}",
        "ambiguous_indel_count": self.ambiguous_indel_count,
        "ambiguous_del_count": self.ambiguous_del_count,
        "ref_label": self.ref_label,
        "ref_accession": self.ref_accession,
        "bracken_species": self.bracken_species,
        "bracken_frac": self.bracken_frac,
        "bracken_reads": self.bracken_reads,
        "alias_used": self.alias_used,
        "reused_bam": self.reused_bam,
        "reused_vcf": self.reused_vcf,
        "runtime_sec": self.runtime_sec,
        "tool_version": self.tool_version,
    }

TSVMode

Bases: Enum

TSV output mode.

Thresholds(dp_min=None, maf_min=None, dp_cap=None, mapq_min=None, baseq_min=None) dataclass

Thresholds for ambiguous site detection.

__post_init__()

Load defaults from configuration and validate threshold values.

Source code in src/ssiamb/models.py

def __post_init__(self) -> None:
    """Load defaults from configuration and validate threshold values."""
    # Import here to avoid circular imports
    from .config import get_config

    # Load defaults from config if not explicitly set
    config = get_config()
    if self.dp_min is None:
        self.dp_min = config.get_threshold("dp_min", 10)
    if self.maf_min is None:
        self.maf_min = config.get_threshold("maf_min", 0.1)
    if self.dp_cap is None:
        self.dp_cap = config.get_threshold("dp_cap", 100)
    if self.mapq_min is None:
        self.mapq_min = config.get_threshold("mapq_min", 20)
    if self.baseq_min is None:
        self.baseq_min = config.get_threshold("baseq_min", 20)

    # Validate threshold values (now they should all be set)
    assert self.dp_min is not None
    assert self.maf_min is not None
    assert self.dp_cap is not None
    assert self.mapq_min is not None
    assert self.baseq_min is not None

    if self.dp_min < 1:
        raise ValueError("dp_min must be >= 1")
    if not 0.0 <= self.maf_min <= 1.0:
        raise ValueError("maf_min must be between 0.0 and 1.0")
    if self.dp_cap < self.dp_min:
        raise ValueError("dp_cap must be >= dp_min")
    if self.mapq_min < 0:
        raise ValueError("mapq_min must be >= 0")
    if self.baseq_min < 0:
        raise ValueError("baseq_min must be >= 0")

from_config(config=None, **overrides) classmethod

Create Thresholds from configuration with optional overrides.

Parameters:

Name	Type	Description	Default
config	Optional['SsiambConfig']	Configuration object (uses global if None)	None
**overrides		Explicit threshold values to override	{}

Returns:

Type	Description
'Thresholds'	Thresholds instance with config defaults and overrides applied

Source code in src/ssiamb/models.py

@classmethod
def from_config(cls, config: Optional["SsiambConfig"] = None, **overrides) -> "Thresholds":
    """
    Create Thresholds from configuration with optional overrides.

    Args:
        config: Configuration object (uses global if None)
        **overrides: Explicit threshold values to override

    Returns:
        Thresholds instance with config defaults and overrides applied
    """
    if config is None:
        from .config import get_config
        config = get_config()

    # Start with config defaults
    kwargs = {
        "dp_min": config.get_threshold("dp_min", 10),
        "maf_min": config.get_threshold("maf_min", 0.1),
        "dp_cap": config.get_threshold("dp_cap", 100),
        "mapq_min": config.get_threshold("mapq_min", 20),
        "baseq_min": config.get_threshold("baseq_min", 20),
    }

    # Apply overrides
    kwargs.update(overrides)

    return cls(**kwargs)

Provenance

Provenance tracking for ssiamb runs.

This module handles collection and output of detailed per-sample provenance information when --emit-provenance is used.

ProvenanceCallerParams(exact_cmdlines) dataclass

Variant caller parameters.

ProvenanceContigFilters(min_contig_len=500) dataclass

Contig filtering settings.

ProvenanceCounts(ambiguous_snv_count=None, ambiguous_indel_count=None, ambiguous_del_count=None, callable_bases=None, genome_length=None) dataclass

Analysis counts.

ProvenanceDuplicatePolicy(denominator_excludes_dups=True, bam_had_dups_flag=None) dataclass

Duplicate handling policy.

ProvenanceExtrasEmitted(vcf=False, bed=False, matrix=False, per_contig=False) dataclass

Optional outputs emitted.

ProvenanceGridCell(depth, maf_bin) dataclass

Grid cell used for analysis.

ProvenanceInput(path, md5=None) dataclass

Input file information with MD5 hash.

ProvenanceMappingStats(total_reads=None, mapped_reads=None, map_rate=None, mean_depth=None, breadth_1x=None, breadth_10x=None) dataclass

Mapping statistics.

ProvenanceRecord(tool_version, python_version, conda_env, started_at, finished_at, runtime_sec, threads, sample, mode, mapper, caller, thresholds, denom_policy, depth_tool, inputs, reference_info, species_selection=None, mapping_stats=None, duplicate_policy=None, contig_filters=None, caller_params=None, counts=None, grid_cell_used=None, extras_emitted=None, warnings=None) dataclass

Complete provenance record for a sample.

ProvenanceReferenceInfo(species_requested=None, alias_applied=None, species_final=None, fasta_path=None, fasta_md5=None, source=None) dataclass

Reference genome information.

ProvenanceSpeciesSelection(bracken_species=None, bracken_frac=None, bracken_reads=None, thresholds=None, on_fail=None) dataclass

Species selection from Bracken.

calculate_md5(file_path)

Calculate MD5 hash of a file.

Source code in src/ssiamb/provenance.py

def calculate_md5(file_path: Path) -> Optional[str]:
    """Calculate MD5 hash of a file."""
    hash_md5 = hashlib.md5()
    try:
        with open(file_path, "rb") as f:
            for chunk in iter(lambda: f.read(4096), b""):
                hash_md5.update(chunk)
        return hash_md5.hexdigest()
    except (FileNotFoundError, PermissionError):
        return None

create_provenance_record(sample, mode, started_at, finished_at, threads, mapper, caller, dp_min, maf_min, dp_cap, denom_policy, depth_tool, emit_vcf, emit_bed, emit_matrix, emit_per_contig, r1=None, r2=None, assembly=None, reference=None, bam=None, vcf=None, species=None, mapping_stats=None, species_selection=None, counts=None, warnings=None, reference_path=None, reference_species=None)

Create a complete provenance record for a sample.

Source code in src/ssiamb/provenance.py

def create_provenance_record(
    sample: str,
    mode: Mode,
    started_at: datetime,
    finished_at: datetime,
    threads: int,
    mapper: str,
    caller: str,
    dp_min: int,
    maf_min: float,
    dp_cap: int,
    denom_policy: str,
    depth_tool: str,
    emit_vcf: bool,
    emit_bed: bool,
    emit_matrix: bool,
    emit_per_contig: bool,
    r1: Optional[Path] = None,
    r2: Optional[Path] = None,
    assembly: Optional[Path] = None,
    reference: Optional[Path] = None,
    bam: Optional[Path] = None,
    vcf: Optional[Path] = None,
    species: Optional[str] = None,
    mapping_stats: Optional[ProvenanceMappingStats] = None,
    species_selection: Optional[ProvenanceSpeciesSelection] = None,
    counts: Optional[ProvenanceCounts] = None,
    warnings: Optional[List[QCWarning]] = None,
    reference_path: Optional[Path] = None,
    reference_species: Optional[str] = None,
) -> ProvenanceRecord:
    """Create a complete provenance record for a sample."""

    # Calculate runtime
    runtime_sec = (finished_at - started_at).total_seconds()

    # Prepare inputs
    inputs = {}
    if r1:
        inputs["r1"] = {"path": str(r1), "md5": calculate_md5(r1)}
    if r2:
        inputs["r2"] = {"path": str(r2), "md5": calculate_md5(r2)}
    if assembly:
        inputs["assembly"] = {"path": str(assembly), "md5": calculate_md5(assembly)}
    if reference:
        inputs["reference"] = {"path": str(reference), "md5": calculate_md5(reference)}
    if bam:
        inputs["bam"] = {"path": str(bam), "md5": calculate_md5(bam)}
    if vcf:
        inputs["vcf"] = {"path": str(vcf), "md5": calculate_md5(vcf)}

    # Reference info
    reference_info = {}
    if mode == Mode.REF:
        reference_info["species_requested"] = species
        reference_info["species_final"] = reference_species
        if reference_path:
            reference_info["fasta_path"] = str(reference_path)
            reference_info["fasta_md5"] = calculate_md5(reference_path)

    # Grid cell used
    grid_cell = {
        "depth": dp_min,
        "maf_bin": int(maf_min * 100)  # Floor of 100 * maf_min
    }

    # Extras emitted
    extras = {
        "vcf": emit_vcf,
        "bed": emit_bed,
        "matrix": emit_matrix,
        "per_contig": emit_per_contig
    }

    # Convert warnings to string list
    warning_messages = []
    if warnings:
        warning_messages = [f"{w.metric}: {w.message}" for w in warnings]

    return ProvenanceRecord(
        tool_version=get_tool_version(),
        python_version=f"{sys.version_info.major}.{sys.version_info.minor}.{sys.version_info.micro}",
        conda_env=get_conda_env(),
        started_at=started_at.isoformat(),
        finished_at=finished_at.isoformat(),
        runtime_sec=runtime_sec,
        threads=threads,
        sample=sample,
        mode=mode.value,
        mapper=mapper,
        caller=caller,
        thresholds={
            "dp_min": dp_min,
            "maf_min": maf_min,
            "dp_cap": dp_cap
        },
        denom_policy=denom_policy,
        depth_tool=depth_tool,
        inputs=inputs,
        reference_info=reference_info,
        species_selection=asdict(species_selection) if species_selection else None,
        mapping_stats=asdict(mapping_stats) if mapping_stats else None,
        duplicate_policy=asdict(ProvenanceDuplicatePolicy()) if mode == Mode.SELF else None,
        contig_filters=asdict(ProvenanceContigFilters(min_contig_len=500)),
        caller_params={"exact_cmdlines": []},  # TODO: Collect actual command lines
        counts=asdict(counts) if counts else {},
        grid_cell_used=grid_cell,
        extras_emitted=extras,
        warnings=warning_messages
    )

get_conda_env()

Get the name of the current conda environment.

Source code in src/ssiamb/provenance.py

def get_conda_env() -> Optional[str]:
    """Get the name of the current conda environment."""
    conda_env = os.environ.get('CONDA_DEFAULT_ENV')
    if conda_env and conda_env != 'base':
        return conda_env
    return None

get_tool_version()

Get the current tool version.

Source code in src/ssiamb/provenance.py

def get_tool_version() -> str:
    """Get the current tool version."""
    try:
        from importlib.metadata import version
        return version('ssiamb')
    except Exception:
        return "unknown"

write_provenance_json(records, output_path)

Write provenance records to JSON file.

Source code in src/ssiamb/provenance.py

def write_provenance_json(records: List[ProvenanceRecord], output_path: Path) -> None:
    """Write provenance records to JSON file."""
    with open(output_path, 'w') as f:
        json.dump([asdict(record) for record in records], f, indent=2)