bedpull 0.3.0

mod cli;

use anyhow::{Context, Result};
use bedpull::paf::PafIndex;
use bedpull::reads::{BamConfig, StitchConfig, get_bam_reads, get_cram_reads, get_paf_reads};
use bedpull::utils::{read_bed, write_fasta_record, write_fastq_record};
use clap::Parser;
use noodles::bam;
use noodles::cram;
use noodles::fasta;
use noodles::fasta::repository::adapters::IndexedReader as FastaIndexedReader;
use std::collections::HashSet;
use std::fs::File;
use std::fs::OpenOptions;
use std::io::{BufReader, BufWriter, Write};
use std::path::{Path, PathBuf};

/// Replace empty `VN:` fields in @PG and @RG header lines.
/// Some samtools versions write `VN:\t` (empty value), which noodles rejects.
fn sanitize_bam_header_text(text: &str) -> String {
    text.lines()
        .map(|line| {
            if !line.starts_with("@PG") && !line.starts_with("@RG") {
                return line.to_string();
            }
            line.split('\t')
                .map(|f| if f == "VN:" { "VN:unknown" } else { f })
                .collect::<Vec<_>>()
                .join("\t")
        })
        .collect::<Vec<_>>()
        .join("\n")
}

/// Open the BAM file with a raw BGZF reader, extract and sanitize the SAM header
/// text, then parse it into a `sam::Header`.  Used as a fallback when noodles'
/// strict parser rejects the header (e.g. empty `VN:` in `@PG` records).
fn read_bam_header_lenient(path: &Path) -> Result<noodles::sam::Header> {
    use std::io::Read;

    let file = File::open(path).context("failed to open BAM file")?;
    let mut reader = noodles::bgzf::io::Reader::new(file);

    let mut magic = [0u8; 4];
    reader
        .read_exact(&mut magic)
        .context("failed to read BAM magic bytes")?;
    if &magic != b"BAM\x01" {
        anyhow::bail!("not a valid BAM file (bad magic)");
    }

    let mut len_buf = [0u8; 4];
    reader
        .read_exact(&mut len_buf)
        .context("failed to read BAM header length")?;
    let l_text = u32::from_le_bytes(len_buf) as usize;

    let mut raw_text = vec![0u8; l_text];
    reader
        .read_exact(&mut raw_text)
        .context("failed to read BAM header text")?;

    let text = std::str::from_utf8(&raw_text)
        .context("BAM header text is not valid UTF-8")?
        .trim_end_matches('\0');

    let sanitized = sanitize_bam_header_text(text);

    sanitized
        .parse::<noodles::sam::Header>()
        .map_err(|e| anyhow::anyhow!("failed to parse sanitized BAM header: {e}"))
}

use cli::Opts;

fn effective_flanks(opts: &Opts) -> (usize, usize) {
    cli::resolve_flanks(opts.flanks, opts.lflank, opts.rflank)
}

/// Recover a region's original 0-based `(start, end)` bounds.
///
/// `read_bed` stores both bounds shifted by `+1` to work around `Position`
/// being `NonZeroUsize` (see its docs) — this is the one place that shift
/// gets undone, via `usize::from(position) - 1`, for all three extraction
/// modes.
fn region_bounds(region: &noodles::core::Region, region_name: &str) -> Result<(usize, usize)> {
    let region_start = region
        .interval()
        .start()
        .map(usize::from)
        .map(|v| v - 1)
        .ok_or_else(|| anyhow::anyhow!("BED region '{}' has unbounded start", region_name))?;
    let region_end = region
        .interval()
        .end()
        .map(usize::from)
        .map(|v| v - 1)
        .ok_or_else(|| anyhow::anyhow!("BED region '{}' has unbounded end", region_name))?;
    Ok((region_start, region_end))
}

/// Build a header suffix describing bases missing from a `--partial` read whose
/// alignment didn't fully span the requested (region ± flank) window, e.g.
/// `|missing_left=12bp|missing_right=8bp`. Returns an empty string when the read
/// covers the full window (the normal, non-partial case).
fn missing_bases_suffix(
    desired_start: usize,
    desired_end: usize,
    ref_start: usize,
    ref_end: usize,
) -> String {
    let mut suffix = String::new();
    if ref_start > desired_start {
        suffix.push_str(&format!("|missing_left={}bp", ref_start - desired_start));
    }
    if ref_end < desired_end {
        suffix.push_str(&format!("|missing_right={}bp", desired_end - ref_end));
    }
    suffix
}

fn bam_config(opts: &Opts) -> BamConfig {
    BamConfig {
        min_mapq: opts.min_mapq,
        include_secondary: opts.include_secondary,
        include_supplementary: opts.include_supplementary,
        partial: opts.partial,
        min_partial_coverage: opts.min_partial_coverage,
        min_region_quality: opts.min_region_quality,
    }
}

fn hap_output_path(base: &Path, hap: u8) -> PathBuf {
    let mut name = base
        .file_stem()
        .unwrap_or_default()
        .to_string_lossy()
        .into_owned();
    name.push_str(&format!(".h{hap}"));
    if let Some(ext) = base.extension() {
        name.push('.');
        name.push_str(&ext.to_string_lossy());
    }
    base.parent().unwrap_or(Path::new(".")).join(name)
}

fn open_writer(path: &Path) -> Result<BufWriter<File>> {
    let f = OpenOptions::new()
        .write(true)
        .create(true)
        .truncate(true)
        .open(path)
        .with_context(|| format!("failed to open output file: {}", path.display()))?;
    Ok(BufWriter::new(f))
}

/// Open a writer for an optional `--<flag>`-style path option: `None` if the option
/// wasn't set (still `"None"`), stdout if set to `-`, otherwise a truncated file.
fn open_optional_writer(path: &Path) -> Result<Option<Box<dyn std::io::Write>>> {
    if path.to_str() == Some("None") {
        return Ok(None);
    }
    if cli::is_stdout(path) {
        return Ok(Some(Box::new(BufWriter::new(std::io::stdout()))));
    }
    let f = OpenOptions::new()
        .write(true)
        .create(true)
        .truncate(true)
        .open(path)
        .with_context(|| format!("failed to open file: {}", path.display()))?;
    Ok(Some(Box::new(BufWriter::new(f))))
}

/// Write one unmapped-region entry to `--unmapped` output: a `#reason` comment line
/// followed by the input BED record (chrom, start, end, name), mirroring liftOver's
/// own unmapped-file convention.
fn write_unmapped_region(
    writer: &mut dyn std::io::Write,
    chr: &str,
    region_start: usize,
    region_end: usize,
    region_name: &str,
    reason: &str,
) -> Result<()> {
    writeln!(writer, "#{reason}").context("failed to write unmapped reason comment")?;
    writeln!(writer, "{chr}\t{region_start}\t{region_end}\t{region_name}")
        .context("failed to write unmapped BED record")?;
    Ok(())
}

fn main() -> Result<()> {
    let opts: Opts = Opts::parse();
    if opts.debug {
        eprintln!("{:#?}", opts);
    }
    crate::cli::check_option_values(&opts)?;
    crate::cli::check_inputs_exist(&opts)?;

    if opts.debug {
        eprintln!("Reading bed file");
    }
    let regions = read_bed(&opts.bed, opts.debug)?;

    let mut read_writer: Box<dyn std::io::Write> = if cli::is_stdout(&opts.output) {
        Box::new(BufWriter::new(std::io::stdout()))
    } else {
        let output_file = OpenOptions::new()
            .write(true)
            .create(true)
            .truncate(true)
            .open(&opts.output)
            .with_context(|| format!("failed to open output file: {}", opts.output.display()))?;
        Box::new(BufWriter::new(output_file))
    };

    if opts.bam.to_str() != Some("None") {
        eprintln!("BAM mode");
        eprintln!("Extracting sequences");
        extract_from_bam(&opts, regions, read_writer.as_mut())?;
    } else if opts.cram.to_str() != Some("None") {
        eprintln!("CRAM mode");
        eprintln!("Extracting sequences");
        extract_from_cram(&opts, regions, read_writer.as_mut())?;
    } else if opts.paf.to_str() != Some("None") && opts.query_ref.to_str() != Some("None") {
        eprintln!("PAF mode");
        eprintln!("Extracting sequences");
        extract_from_paf(&opts, regions, read_writer.as_mut())?;
    }

    eprintln!("Done");
    Ok(())
}

pub fn extract_from_bam(
    opts: &Opts,
    regions: Vec<(noodles::core::Region, String, String)>,
    read_writer: &mut dyn std::io::Write,
) -> Result<()> {
    let mut seen: HashSet<String> = HashSet::new();

    // Open per-haplotype writers when --hap_split is set.
    let mut hap_writers: Option<[BufWriter<File>; 3]> = None;
    if opts.hap_split {
        hap_writers = Some([
            open_writer(&hap_output_path(&opts.output, 0))?,
            open_writer(&hap_output_path(&opts.output, 1))?,
            open_writer(&hap_output_path(&opts.output, 2))?,
        ]);
    }

    // Open the BAM once and reuse it for every region — reopening + re-parsing the
    // header per region dominates runtime once a BED file has more than a handful
    // of regions.
    let mut reader = bam::io::indexed_reader::Builder::default()
        .build_from_path(&opts.bam)
        .context("failed to open BAM file")?;
    let header = reader
        .read_header()
        .or_else(|e| {
            if e.kind() == std::io::ErrorKind::InvalidData {
                eprintln!(
                    "Warning: BAM header has non-standard fields (e.g. empty VN: in @PG \
                     records — produced by some samtools versions). Retrying with lenient parser."
                );
                read_bam_header_lenient(&opts.bam)
                    .map_err(|ae| std::io::Error::new(std::io::ErrorKind::InvalidData, ae))
            } else {
                Err(e)
            }
        })
        .context("failed to read BAM header")?;

    let mut unmapped_writer = open_optional_writer(&opts.unmapped)?;

    for (region, region_name, chr) in regions.iter() {
        if opts.debug {
            eprintln!("===============================");
            eprintln!("Analysing region: {}, {}", region, region_name);
            eprintln!("===============================");
        }

        let (region_start, region_end) = region_bounds(region, region_name)?;

        if region.name().contains(&b'#') {
            let reason = "region skipped (chromosome name contains '#')";
            eprintln!("Region {} has a #, skipping", region_name);
            if let Some(w) = unmapped_writer.as_mut() {
                write_unmapped_region(
                    w.as_mut(),
                    chr,
                    region_start,
                    region_end,
                    region_name,
                    reason,
                )?;
            }
            continue;
        }

        let query = reader
            .query(&header, region)
            .context("BAM region query failed")?;

        // find all reads that map to region
        // apply filters (full length, quality, etc)
        // cut out sequence (optionally qstring too and do quality calculation)
        let (lflank, rflank) = effective_flanks(opts);
        let (overlapping_reads, candidates_seen) =
            get_bam_reads(&bam_config(opts), query, region, lflank, rflank)?;
        if overlapping_reads.is_empty() {
            let reason = if candidates_seen == 0 {
                "no overlapping reads found".to_string()
            } else {
                format!(
                    "{candidates_seen} candidate read(s) found but all were filtered out (--min_mapq/--include_secondary/--include_supplementary/--partial/--min_partial_coverage/--min_region_quality)"
                )
            };
            eprintln!(
                "No reads found for region in bam file. Skipping region: {}",
                region_name
            );
            if let Some(w) = unmapped_writer.as_mut() {
                write_unmapped_region(
                    w.as_mut(),
                    chr,
                    region_start,
                    region_end,
                    region_name,
                    &reason,
                )?;
            }
            continue;
        }
        // write to fasta or fastq
        let desired_start = region_start.saturating_sub(lflank);
        let desired_end = region_end + rflank;
        let matched_count = overlapping_reads.len();
        let mut written_count = 0usize;
        for (name, subseq, subqual, ref_start, ref_end, hap) in overlapping_reads {
            if opts.dedup && !seen.insert(name.clone()) {
                continue;
            }
            written_count += 1;
            let hap_suffix = if hap > 0 {
                format!("|h{}", hap)
            } else {
                String::new()
            };
            let missing_suffix =
                missing_bases_suffix(desired_start, desired_end, ref_start, ref_end);
            let head = format!(
                "{}|{}:{}-{}|{}{}{}",
                name, chr, region_start, region_end, region_name, hap_suffix, missing_suffix
            );
            let seq_str =
                std::str::from_utf8(&subseq).context("BAM sequence contains invalid UTF-8")?;
            let writer: &mut dyn std::io::Write = match hap_writers.as_mut() {
                Some(writers) => match hap {
                    1 => &mut writers[1],
                    2 => &mut writers[2],
                    _ => {
                        if hap > 2 {
                            eprintln!("Warning: unexpected HP tag value {hap}, routing to h0");
                        }
                        &mut writers[0]
                    }
                },
                None => read_writer,
            };
            if opts.fastq {
                write_fastq_record(writer, &head, seq_str, &subqual)
                    .context("failed to write FASTQ record")?;
            } else {
                write_fasta_record(writer, &head, seq_str)
                    .context("failed to write FASTA record")?;
            }
        }
        if written_count == 0 {
            let reason = format!(
                "{matched_count} matching read(s) found but all were already emitted for another region (--dedup)"
            );
            if let Some(w) = unmapped_writer.as_mut() {
                write_unmapped_region(
                    w.as_mut(),
                    chr,
                    region_start,
                    region_end,
                    region_name,
                    &reason,
                )?;
            }
        }
    }
    Ok(())
}

pub fn extract_from_cram(
    opts: &Opts,
    regions: Vec<(noodles::core::Region, String, String)>,
    read_writer: &mut dyn std::io::Write,
) -> Result<()> {
    let mut seen: HashSet<String> = HashSet::new();

    // Open per-haplotype writers when --hap_split is set.
    let mut hap_writers: Option<[BufWriter<File>; 3]> = None;
    if opts.hap_split {
        hap_writers = Some([
            open_writer(&hap_output_path(&opts.output, 0))?,
            open_writer(&hap_output_path(&opts.output, 1))?,
            open_writer(&hap_output_path(&opts.output, 2))?,
        ]);
    }

    // Build a reference sequence repository from --reference if provided.
    let reference_repo = if opts.reference.to_str() != Some("None") {
        let indexed = fasta::io::indexed_reader::Builder::default()
            .build_from_path(&opts.reference)
            .with_context(|| {
                format!(
                    "failed to open reference FASTA: {}",
                    opts.reference.display()
                )
            })?;
        fasta::Repository::new(FastaIndexedReader::new(indexed))
    } else {
        // Whether this CRAM actually needs an external reference isn't
        // reliably detectable up front without deeper container/codec
        // inspection than noodles exposes here, so rather than guess, warn:
        // a reference-compressed CRAM read without --reference won't error
        // loudly, it'll just decode wrong (empty/garbled) sequence.
        eprintln!(
            "Warning: --cram given without --reference. If this CRAM was compressed \
             against an external reference (the common case), sequences will decode \
             incorrectly rather than error — pass --reference <fasta> if extracted \
             sequences look empty or wrong. Only CRAMs with embedded sequences (see \
             docs/src/cram-mode.md) can safely omit --reference."
        );
        fasta::Repository::default()
    };

    // Open the CRAM once and reuse it for every region — reopening + re-parsing the
    // header per region dominates runtime once a BED file has more than a handful
    // of regions.
    let mut reader = cram::io::indexed_reader::Builder::default()
        .set_reference_sequence_repository(reference_repo.clone())
        .build_from_path(&opts.cram)
        .context("failed to open CRAM file")?;
    let header = reader.read_header().context("failed to read CRAM header")?;

    let mut unmapped_writer = open_optional_writer(&opts.unmapped)?;

    for (region, region_name, chr) in regions.iter() {
        if opts.debug {
            eprintln!("===============================");
            eprintln!("Analysing region: {}, {}", region, region_name);
            eprintln!("===============================");
        }

        let (region_start, region_end) = region_bounds(region, region_name)?;

        if region.name().contains(&b'#') {
            let reason = "region skipped (chromosome name contains '#')";
            eprintln!("Region {} has a #, skipping", region_name);
            if let Some(w) = unmapped_writer.as_mut() {
                write_unmapped_region(
                    w.as_mut(),
                    chr,
                    region_start,
                    region_end,
                    region_name,
                    reason,
                )?;
            }
            continue;
        }

        let query = reader
            .query(&header, region)
            .context("CRAM region query failed")?;

        let (lflank, rflank) = effective_flanks(opts);
        let (overlapping_reads, candidates_seen) =
            get_cram_reads(&bam_config(opts), query, region, lflank, rflank)?;
        if overlapping_reads.is_empty() {
            let reason = if candidates_seen == 0 {
                "no overlapping reads found".to_string()
            } else {
                format!(
                    "{candidates_seen} candidate read(s) found but all were filtered out (--min_mapq/--include_secondary/--include_supplementary/--partial/--min_partial_coverage/--min_region_quality)"
                )
            };
            eprintln!(
                "No reads found for region in CRAM file. Skipping region: {}",
                region_name
            );
            if let Some(w) = unmapped_writer.as_mut() {
                write_unmapped_region(
                    w.as_mut(),
                    chr,
                    region_start,
                    region_end,
                    region_name,
                    &reason,
                )?;
            }
            continue;
        }

        let desired_start = region_start.saturating_sub(lflank);
        let desired_end = region_end + rflank;
        let matched_count = overlapping_reads.len();
        let mut written_count = 0usize;
        for (name, subseq, subqual, ref_start, ref_end, hap) in overlapping_reads {
            if opts.dedup && !seen.insert(name.clone()) {
                continue;
            }
            written_count += 1;
            let hap_suffix = if hap > 0 {
                format!("|h{}", hap)
            } else {
                String::new()
            };
            let missing_suffix =
                missing_bases_suffix(desired_start, desired_end, ref_start, ref_end);
            let head = format!(
                "{}|{}:{}-{}|{}{}{}",
                name, chr, region_start, region_end, region_name, hap_suffix, missing_suffix
            );
            let seq_str =
                std::str::from_utf8(&subseq).context("CRAM sequence contains invalid UTF-8")?;
            let writer: &mut dyn std::io::Write = match hap_writers.as_mut() {
                Some(writers) => match hap {
                    1 => &mut writers[1],
                    2 => &mut writers[2],
                    _ => {
                        if hap > 2 {
                            eprintln!("Warning: unexpected HP tag value {hap}, routing to h0");
                        }
                        &mut writers[0]
                    }
                },
                None => read_writer,
            };
            if opts.fastq {
                write_fastq_record(writer, &head, seq_str, &subqual)
                    .context("failed to write FASTQ record")?;
            } else {
                write_fasta_record(writer, &head, seq_str)
                    .context("failed to write FASTA record")?;
            }
        }
        if written_count == 0 {
            let reason = format!(
                "{matched_count} matching read(s) found but all were already emitted for another region (--dedup)"
            );
            if let Some(w) = unmapped_writer.as_mut() {
                write_unmapped_region(
                    w.as_mut(),
                    chr,
                    region_start,
                    region_end,
                    region_name,
                    &reason,
                )?;
            }
        }
    }
    Ok(())
}

pub fn extract_from_paf(
    opts: &Opts,
    regions: Vec<(noodles::core::Region, String, String)>,
    read_writer: &mut dyn std::io::Write,
) -> Result<()> {
    let mut seen: HashSet<String> = HashSet::new();

    // Open per-haplotype writers when --hap_split is set.
    let mut hap_writers: Option<[BufWriter<File>; 3]> = None;
    if opts.hap_split {
        hap_writers = Some([
            open_writer(&hap_output_path(&opts.output, 0))?,
            open_writer(&hap_output_path(&opts.output, 1))?,
            open_writer(&hap_output_path(&opts.output, 2))?,
        ]);
    }

    // Open --bed_out writer if requested.
    let mut bed_out_writer: Option<Box<dyn std::io::Write>> =
        if opts.bed_out.to_str() != Some("None") {
            if cli::is_stdout(&opts.bed_out) {
                Some(Box::new(BufWriter::new(std::io::stdout())))
            } else {
                let f = OpenOptions::new()
                    .write(true)
                    .create(true)
                    .truncate(true)
                    .open(&opts.bed_out)
                    .with_context(|| {
                        format!("failed to open bed_out file: {}", opts.bed_out.display())
                    })?;
                Some(Box::new(BufWriter::new(f)))
            }
        } else {
            None
        };

    // Build or load index
    let paf_path = opts
        .paf
        .to_str()
        .context("PAF path contains invalid UTF-8")?;
    let query_ref = opts
        .query_ref
        .to_str()
        .context("query_ref path contains invalid UTF-8")?;
    let index_path = format!("{}.idx", paf_path);
    let index = if opts.use_paf_index {
        if std::path::Path::new(&index_path).exists() {
            eprintln!("Loading PAF index from {}", index_path);
            PafIndex::load(&index_path)
                .with_context(|| format!("failed to load PAF index from {}", index_path))?
        } else {
            eprintln!("Building PAF index...");
            let index = PafIndex::build(paf_path).context("failed to build PAF index")?;
            index
                .save(&index_path)
                .with_context(|| format!("failed to save PAF index to {}", index_path))?;
            eprintln!("Index saved to {}", index_path);
            index
        }
    } else {
        PafIndex::build(paf_path).context("failed to build PAF index")?
    };

    // Open the PAF file and query FASTA once and reuse them for every region/record —
    // reopening per record dominates runtime once a BED file has more than a handful
    // of regions.
    let mut paf_reader = BufReader::new(
        File::open(paf_path).with_context(|| format!("failed to open PAF file: {paf_path}"))?,
    );
    let mut fasta_reader = fasta::io::indexed_reader::Builder::default()
        .build_from_path(query_ref)
        .with_context(|| format!("failed to open query_ref FASTA: {query_ref}"))?;

    let mut unmapped_writer = open_optional_writer(&opts.unmapped)?;

    let stitch = StitchConfig {
        enabled: opts.stitch_records,
        max_gap: opts.max_stitch_gap,
    };

    // for each region, get paf regions and extract sequences
    for (region, region_name, chr) in regions.iter() {
        if opts.debug {
            eprintln!("===============================");
            eprintln!("Analysing region: {}, {}", region, region_name);
            eprintln!("===============================");
        }

        let (region_start, region_end) = region_bounds(region, region_name)?;

        if region.name().contains(&b'#') {
            let reason = "region skipped (chromosome name contains '#')";
            eprintln!("Region {} has a #, skipping", region_name);
            if let Some(w) = unmapped_writer.as_mut() {
                write_unmapped_region(
                    w.as_mut(),
                    chr,
                    region_start,
                    region_end,
                    region_name,
                    reason,
                )?;
            }
            continue;
        }

        let (lflank, rflank) = effective_flanks(opts);

        // Query the flank-expanded window, not just the raw region — otherwise a
        // record that only overlaps within the flank zone (relevant to both plain
        // flank extraction and --stitch_records chain members) never reaches
        // get_paf_reads at all, since it's filtered out here first.
        let query_start = region_start.saturating_sub(lflank);
        let query_end = region_end + rflank;
        let overlapping_entries = index.query(chr, query_start, query_end);
        if opts.debug {
            eprintln!("Found {} overlapping alignments", overlapping_entries.len());
        }

        let reads = get_paf_reads(
            &mut paf_reader,
            &mut fasta_reader,
            &overlapping_entries,
            region_start,
            region_end,
            lflank,
            rflank,
            stitch,
            opts.debug,
        )?;
        if reads.is_empty() {
            let reason = if overlapping_entries.is_empty() {
                "no overlapping alignments found".to_string()
            } else {
                format!(
                    "{} overlapping alignment(s) found but none produced output (missing CIGAR, no valid overlap, or invalid coordinates)",
                    overlapping_entries.len()
                )
            };
            eprintln!(
                "No overlapping alignments produced output for region in PAF file. Skipping region: {}",
                region_name
            );
            if let Some(w) = unmapped_writer.as_mut() {
                write_unmapped_region(
                    w.as_mut(),
                    chr,
                    region_start,
                    region_end,
                    region_name,
                    &reason,
                )?;
            }
            continue;
        }

        let matched_count = reads.len();
        let mut written_count = 0usize;
        for (sequence, query_name, query_start, query_end, strand, hap) in reads {
            if opts.dedup && !seen.insert(query_name.clone()) {
                continue;
            }
            written_count += 1;
            if let Some(bed_writer) = bed_out_writer.as_mut() {
                writeln!(
                    bed_writer,
                    "{}\t{}\t{}\t{}\t0\t{}",
                    query_name, query_start, query_end, region_name, strand
                )
                .context("failed to write BED record")?;
            }

            let hap_suffix = if hap > 0 {
                format!("|h{}", hap)
            } else {
                String::new()
            };
            let header = format!(
                "{}|{}:{}-{}|{}|{}:{}-{}|{}{}",
                query_name,
                chr,
                region_start,
                region_end,
                region_name,
                query_name,
                query_start,
                query_end,
                strand,
                hap_suffix
            );
            let writer: &mut dyn std::io::Write = match hap_writers.as_mut() {
                Some(writers) => match hap {
                    1 => &mut writers[1],
                    2 => &mut writers[2],
                    _ => {
                        if hap > 2 {
                            eprintln!("Warning: unexpected HP tag value {hap}, routing to h0");
                        }
                        &mut writers[0]
                    }
                },
                None => read_writer,
            };
            write_fasta_record(writer, &header, &sequence)
                .context("failed to write FASTA record")?;
        }
        if written_count == 0 {
            let reason = format!(
                "{matched_count} matching alignment(s) found but all were already emitted for another region (--dedup)"
            );
            if let Some(w) = unmapped_writer.as_mut() {
                write_unmapped_region(
                    w.as_mut(),
                    chr,
                    region_start,
                    region_end,
                    region_name,
                    &reason,
                )?;
            }
        }
    }
    Ok(())
}

#[cfg(test)]
mod tests {
    use crate::cli::resolve_flanks;
    use bedpull::ToCigarOps;
    use bedpull::paf::{PafIndex, read_paf_record_at_offset};
    use bedpull::utils::get_read_cuts;

    const PAF_PATH: &str = "examples/hg002pat_to_hs1.rfc1_only.paf";

    // RFC1 BED region: chr4 39318077-39318136 (59 bp reference span)
    // This alignment has a 520 bp insertion inside the region, so bedpull
    // should extract 579 bp from the query — see README for context.
    const RFC1_TARGET_START: usize = 31058861; // PAF field 8  (align_start)
    const RFC1_REGION_START: usize = 39318077; // BED start
    const RFC1_REGION_END: usize = 39318136; // BED end
    const RFC1_EXPECTED_BP: usize = 579;

    #[test]
    fn region_bounds_recovers_zero_start_without_erroring() {
        // Regression test: BED start == 0 (the first base of a chromosome, an
        // extremely common, valid coordinate) used to be forwarded straight
        // into noodles' NonZeroUsize-backed Position, which can't represent
        // 0 — read_bed errored, aborting the whole run. read_bed now stores
        // both bounds shifted by +1; region_bounds must recover the exact
        // original 0-based values, including 0 itself.
        let dir = tempfile::tempdir().unwrap();
        let bed_path = dir.path().join("zero_start.bed");
        std::fs::write(&bed_path, "chr1\t0\t1000\tZERO_START\n").unwrap();

        let regions = bedpull::utils::read_bed(&bed_path, false)
            .expect("read_bed should not error on start=0");
        assert_eq!(regions.len(), 1);
        let (region, name, _chr) = &regions[0];
        assert_eq!(name, "ZERO_START");

        let (start, end) =
            super::region_bounds(region, name).expect("region_bounds should recover bounds");
        assert_eq!((start, end), (0, 1000));
    }

    #[test]
    fn paf_index_build_finds_rfc1_region() {
        let idx = PafIndex::build(PAF_PATH).expect("failed to build index");
        let hits = idx.query("chr4", RFC1_REGION_START, RFC1_REGION_END);
        assert_eq!(hits.len(), 1, "expected exactly one alignment over RFC1");
    }

    #[test]
    fn paf_record_at_offset_zero_parses_correctly() {
        let r = read_paf_record_at_offset(PAF_PATH, 0).expect("failed to read record");
        assert_eq!(r.query_name, "chr4_PATERNAL");
        assert_eq!(r.target_name, "chr4");
        assert_eq!(r.target_start, RFC1_TARGET_START);
        assert!(r.cigar.is_some(), "expected cg:Z: tag");
    }

    #[test]
    fn get_read_cuts_rfc1_captures_insertion() {
        use noodles::sam::alignment::record::cigar::op::Kind;
        let r = read_paf_record_at_offset(PAF_PATH, 0).unwrap();
        let cigar_str = r.cigar.as_deref().expect("no CIGAR");
        let ops = cigar_str.to_cigar_ops().expect("valid CIGAR from PAF file");
        // align_end = target_start + reference bases consumed (M/=/X/D/N).
        let ref_len: usize = ops
            .iter()
            .filter(|op| {
                matches!(
                    op.kind,
                    Kind::Match
                        | Kind::SequenceMatch
                        | Kind::SequenceMismatch
                        | Kind::Deletion
                        | Kind::Skip
                )
            })
            .map(|op| op.len)
            .sum();
        let align_end = RFC1_TARGET_START + ref_len;
        let cuts = get_read_cuts(
            &ops,
            RFC1_TARGET_START,
            align_end,
            RFC1_REGION_START,
            RFC1_REGION_END,
        );
        let extracted_len = cuts.read_end - cuts.read_start;
        assert_eq!(
            extracted_len, RFC1_EXPECTED_BP,
            "expected {RFC1_EXPECTED_BP} bp (59 bp ref span + 520 bp insertion); got {extracted_len}"
        );
    }

    #[test]
    fn resolve_flanks_zero_is_identity() {
        assert_eq!(resolve_flanks(0, 0, 0), (0, 0));
    }
}