gwseq_io/bbi/

convert.rs

//! bedGraph / WIG → bigWig, BED → bigBed.
//!
//! Streaming throughout: nothing is sorted or spooled,
//! so a conversion of any size holds a megabyte of input and one open section.
//! Input that is not already pooled by chromosome and ordered raises, naming
//! the line — which is what `sort -k1,1 -k2,2n` is for.

use std::io::{BufRead, Read as _};
use std::path::Path;

use indexmap::IndexMap;

use crate::bbi::header::{BbiKind, BED_FIELD_NAMES};
use crate::bbi::writer::{BbiWriter, BbiWriterOptions};
use crate::error::{Error, Result};
use crate::genomic::ChrMap;
use crate::progress::{CancelFlag, ProgressFn, ProgressTracker};

/// Values gathered before they are handed to the writer in one call.
const VALUE_BATCH: usize = 65536;
/// Lines between two progress reports.
const PROGRESS_INTERVAL: u64 = 65536;
/// A line longer than this is a binary file being read as text, not a record.
const MAX_LINE_SIZE: usize = 16 << 20;

/// What an ordering complaint is asking for. Appended to anything the *writer*
/// refuses, which is where an out-of-order or overlapping input surfaces — the
/// converter itself cannot tell that from a malformed line, and the caller's
/// next move is the same either way.
const ORDER_HINT: &str = "input must be pooled by chromosome and sorted by start, \
                          eg with `sort -k1,1 -k2,2n`";

#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum TextFormat {
    BedGraph,
    Wig,
    Bed,
}

impl TextFormat {
    pub fn as_str(self) -> &'static str {
        match self {
            TextFormat::BedGraph => "bedgraph",
            TextFormat::Wig => "wig",
            TextFormat::Bed => "bed",
        }
    }
}

#[derive(Debug, Clone)]
pub struct ConvertResult {
    pub format: TextFormat,
    pub line_count: u64,
    pub item_count: u64,
    pub skipped_count: u64,
    pub clipped_count: u64,
    pub chr_sizes: Vec<(String, i64)>,
}

/// Reads lines and says which one it is on, so a failure names it.
struct LineReader {
    input: crate::source::TextInput,
    path: String,
    size: u64,
    consumed: u64,
    line_number: u64,
    buffer: Vec<u8>,
}

impl LineReader {
    fn open(path: &Path) -> Result<Self> {
        let (input, size) = crate::source::open_text(path)?;
        Ok(Self {
            input,
            path: path.to_string_lossy().into_owned(),
            size,
            consumed: 0,
            line_number: 0,
            buffer: Vec::with_capacity(4096),
        })
    }

    /// The next line, with its newline and any carriage return before it
    /// stripped. `None` at the end of the file; a final line with no newline of
    /// its own is still returned.
    /// The next line into `out`, `false` at the end of the file.
    ///
    /// Into the caller's buffer rather than back as a `&str`, so the loop can
    /// hold the line and still call `fail`/`guard` — which borrow the reader
    /// for the path and the line number. Returning a borrow meant a
    /// `to_string()` per line to break that conflict, which on a multi-GB
    /// bedGraph is an allocation per record; one buffer reused for the file is
    /// a `memcpy` instead.
    fn read_line_into(&mut self, out: &mut String) -> Result<bool> {
        let Some(line) = self.next_line()? else {
            return Ok(false);
        };
        // Two moves of the same bytes would be one too many, but the borrow
        // checker cannot see that `line` comes out of `self` — and this is the
        // copy the `to_string` was making anyway, minus the allocation.
        out.clear();
        out.push_str(line);
        Ok(true)
    }

    fn next_line(&mut self) -> Result<Option<&str>> {
        self.buffer.clear();
        // Capped as it is read, not after: `read_until` on a binary file with
        // no newline in it grows the buffer to the whole file before anything
        // gets to complain about the length. One byte over the limit is enough
        // to tell that it was exceeded.
        let read = self
            .input
            .by_ref()
            .take(MAX_LINE_SIZE as u64 + 1)
            .read_until(b'\n', &mut self.buffer)
            .map_err(|e| Error::io(&self.path, e))?;
        if read == 0 {
            return Ok(None);
        }
        // Counted before the trim, since this is progress through the file.
        self.consumed += read as u64;
        self.line_number += 1;
        if self.buffer.len() > MAX_LINE_SIZE {
            return Err(Error::format(
                &self.path,
                format!(
                    "line {} is longer than {MAX_LINE_SIZE} bytes",
                    self.line_number
                ),
            ));
        }
        while matches!(self.buffer.last(), Some(b'\n' | b'\r')) {
            self.buffer.pop();
        }
        // Lossy rather than strict, and lossy *per byte*: a bed's name column
        // carries whatever the annotation had in it, and refusing a file for
        // one bad byte in a comment would be worse than passing the replacement
        // character on. Replacing the whole line with the one replacement
        // character — which is what `unwrap_or` did here — is worse still: a
        // `# café` in Latin-1 became a line that is neither a record nor a
        // declaration, and a bed with one such byte in a name column became a
        // line with one column.
        //
        // The fix-up rewrites the buffer rather than staging the string
        // elsewhere, so this still hands back a borrow of one place. It costs a
        // second validation pass on every line, which beside the allocation and
        // the split this line is about to go through is nothing.
        if std::str::from_utf8(&self.buffer).is_err() {
            self.buffer = String::from_utf8_lossy(&self.buffer)
                .into_owned()
                .into_bytes();
        }
        Ok(Some(
            std::str::from_utf8(&self.buffer).expect("the buffer above is valid utf-8 or replaced"),
        ))
    }

    /// Wrap an error with the line it happened on, which is what makes a
    /// conversion failure actionable.
    ///
    /// `Format`, not `InvalidArgument`: everything that reaches here is a
    /// complaint about the *input file* — a column that will not parse, a
    /// record out of order, a coordinate past the end of its chromosome — and
    /// the hierarchy's own definition puts that under `InvalidFile`, which is
    /// what a caller catching "the file is bad" reaches for. `InvalidArgument`
    /// is for what the caller asked for, and Python surfaces it as a
    /// `ValueError`, which a malformed bedGraph is not. The rendered message is
    /// unchanged: `Format` prints as `{path}: {what}`.
    ///
    /// Every helper that parses a field is wrapped through here, so their own
    /// variants never reach a caller.
    fn fail(&self, message: impl std::fmt::Display) -> Error {
        Error::format(&self.path, format!("line {}: {message}", self.line_number))
    }

    /// As `fail`, and says what an ordering complaint is asking for. Used
    /// around the calls into the writer, which is where one surfaces.
    fn guard(&self, message: impl std::fmt::Display) -> Error {
        Error::format(
            &self.path,
            format!("line {}: {message}\n{ORDER_HINT}", self.line_number),
        )
    }
}

// -- line shapes -----------------------------------------------------------

fn is_blank(c: char) -> bool {
    c == ' ' || c == '\t'
}

/// Split on runs of spaces and tabs, leading and trailing ones dropped.
///
/// What bedGraph and WIG are read with. UCSC's readers split them this way, and
/// neither format has a column that can hold a space, so a run of whitespace is
/// always a separator and never data.
fn split_blanks(line: &str) -> Vec<&str> {
    line.split(is_blank).filter(|f| !f.is_empty()).collect()
}

/// Split on single tabs, so an empty column stays an empty field.
///
/// What BED is read with. A BED is tab-delimited and its name column is allowed
/// to hold spaces — UCSC's own tables are full of them — so splitting a BED on
/// whitespace would cut a name in half.
fn split_tabs(line: &str) -> Vec<&str> {
    line.split('\t').collect()
}

/// Drop a trailing empty column that is a formatting artefact rather than a
/// field.
///
/// A line ending on a tab and a line whose last column is genuinely empty look
/// identical on their own, so this decides by the shape the file has already
/// shown. `expected` is the column count the first record fixed — or the one
/// `fields=` declared, when the caller said — and a line carrying exactly one
/// more than that, empty, ended on a tab. A line already the right width keeps
/// its empty last column, which is what it is.
///
/// Dropping unconditionally, as this once did, made a BED whose `name` column
/// is empty on some records and not on others fail with "has N columns and the
/// first one had N+1" — the pop applied to some lines and not to others.
fn trim_trailing_tab(fields: &mut Vec<&str>, expected: Option<usize>) {
    if fields.len() < 2 || !fields.last().is_some_and(|f| f.is_empty()) {
        return;
    }
    match expected {
        // Nothing to compare against yet: the first record of a file with no
        // declared fields, where a trailing tab is the likelier reading.
        None => {
            fields.pop();
        }
        Some(width) if fields.len() == width + 1 => {
            fields.pop();
        }
        Some(_) => {}
    }
}

/// True when the line begins with `token`, case-insensitively, followed by a
/// space, a tab or the end.
fn starts_with_token(line: &str, token: &str) -> bool {
    let bytes = line.as_bytes();
    let token = token.as_bytes();
    if bytes.len() < token.len() {
        return false;
    }
    if !bytes[..token.len()].eq_ignore_ascii_case(token) {
        return false;
    }
    bytes.len() == token.len() || bytes[token.len()] == b' ' || bytes[token.len()] == b'\t'
}

/// True when a line carries no data: blank, a comment, or one of the `track`
/// and `browser` declarations a text track is wrapped in.
fn is_skipped_line(line: &str) -> bool {
    let trimmed = line.trim_start_matches(is_blank);
    trimmed.is_empty()
        || trimmed.starts_with('#')
        || starts_with_token(trimmed, "track")
        || starts_with_token(trimmed, "browser")
}

fn is_declaration(field: &str) -> bool {
    starts_with_token(field, "fixedstep") || starts_with_token(field, "variablestep")
}

/// True when the four fields have the shape of a bedGraph record: two whole
/// numbers and a number after the chromosome.
///
/// Only ever asked of the first line carrying data, to settle the format. Once
/// that is settled the fields are parsed for real, and a later record that does
/// not parse is an error rather than a change of mind.
fn is_bedgraph_record(fields: &[&str]) -> bool {
    fields[1].parse::<i64>().is_ok()
        && fields[2].parse::<i64>().is_ok()
        && fields[3].parse::<f64>().is_ok()
}

/// True when the fields are the body of a WIG section: one number for a
/// fixedStep, or a position and a number for a variableStep.
///
/// Only ever asked of a line that has settled neither format, to tell a WIG
/// missing its declaration from a file that is not a WIG at all.
fn is_orphan_wig_data(fields: &[&str]) -> bool {
    !fields.is_empty() && fields.len() <= 2 && fields.iter().all(|f| f.parse::<f64>().is_ok())
}

/// Which format the input is comes from its **content**, not its name: the
/// first line that is neither blank, a comment, nor a `track` or `browser`
/// declaration decides. A `fixedStep` or `variableStep` line makes it a WIG,
/// four columns a bedGraph, anything else is refused. The decision is made
/// once, so a file holding both is refused too.
pub fn sniff_format(first_data_line: &str) -> Result<TextFormat> {
    let fields = split_blanks(first_data_line);
    if fields.first().is_some_and(|f| is_declaration(f)) {
        return Ok(TextFormat::Wig);
    }
    if fields.len() == 4 && is_bedgraph_record(&fields) {
        return Ok(TextFormat::BedGraph);
    }
    if is_orphan_wig_data(&fields) {
        return Err(Error::invalid(
            "wig data before any fixedStep or variableStep declaration",
        ));
    }
    Err(Error::invalid(format!(
        "\"{first_data_line}\" is neither a bedgraph record (chr, start, end, value) \
         nor a wig declaration (fixedStep or variableStep), so the format of the input \
         cannot be told"
    )))
}

// -- the wig declaration ---------------------------------------------------

/// The declaration a WIG section is written under (the UCSC wiggle format).
///
/// WIG coordinates are 1-based, bedGraph 0-based half-open. `step` and `span`
/// both default to 1; a declaration with no `chrom`, or a `fixedStep` with no
/// `start`, is an error rather than a guess — a WIG with no start has no
/// coordinates at all.
#[derive(Debug, Clone, Default)]
pub struct WigDeclaration {
    pub fixed_step: bool,
    pub chr: String,
    /// 0-based, converted from the 1-based coordinate the file carries.
    /// Meaningless for a variableStep.
    pub start: i64,
    pub step: i64,
    pub span: i64,
}

pub fn parse_wig_declaration(line: &str) -> Result<WigDeclaration> {
    let fields = split_blanks(line);
    let mut declaration = WigDeclaration {
        fixed_step: fields
            .first()
            .is_some_and(|f| starts_with_token(f, "fixedstep")),
        step: 1,
        span: 1,
        ..Default::default()
    };
    let (mut has_chr, mut has_start) = (false, false);
    for field in &fields[1..] {
        let Some((key, value)) = field.split_once('=') else {
            return Err(Error::invalid(format!(
                "\"{field}\" is not a key=value of a wig declaration"
            )));
        };
        let number = |what: &str| -> Result<i64> {
            value
                .parse::<i64>()
                .map_err(|_| Error::invalid(format!("could not read \"{value}\" as a {what}")))
        };
        match key.to_ascii_lowercase().as_str() {
            "chrom" => {
                declaration.chr = value.to_string();
                has_chr = true;
            }
            "start" => {
                let start = number("start")?;
                if start < 1 {
                    return Err(Error::invalid(format!(
                        "start {start} is not a 1-based coordinate"
                    )));
                }
                declaration.start = start - 1;
                has_start = true;
            }
            "step" => declaration.step = number("step")?,
            "span" => declaration.span = number("span")?,
            other => {
                return Err(Error::invalid(format!(
                    "{other} is not a wig declaration key (chrom, start, step, span)"
                )))
            }
        }
    }
    if !has_chr {
        return Err(Error::invalid("wig declaration has no chrom"));
    }
    if declaration.fixed_step && !has_start {
        return Err(Error::invalid("fixedStep declaration has no start"));
    }
    if declaration.step <= 0 {
        return Err(Error::invalid(format!(
            "step {} must be positive",
            declaration.step
        )));
    }
    if declaration.span <= 0 {
        return Err(Error::invalid(format!(
            "span {} must be positive",
            declaration.span
        )));
    }
    Ok(declaration)
}

// -- the binning sink ------------------------------------------------------

/// Takes values from the reader and hands them to the writer, binning them
/// first when a bin size was asked for.
///
/// With a bin size the writer is handed bins rather than values, so this
/// carries the ordering and clipping contracts itself. Otherwise overlapping
/// input averages into a bin and comes out as a number nothing in the file
/// says, and a value hanging over a chromosome opens a bin past it and fails
/// only at the flush, naming a position rather than its line. The wording is
/// the writer's, so the same input reports the same way either way.
struct ValueSink {
    bin_size: i64,
    declared: Option<ChrMap>,
    chr: String,
    /// The open bin, and the run of closed ones waiting to go out.
    bin: Option<i64>,
    bin_sum: f64,
    bin_covered: i64,
    run_start_bin: Option<i64>,
    run_values: Vec<f32>,
    item_count: u64,
    clipped_count: u64,
    /// End of the last value taken on this chromosome, and the size the
    /// declaration gives it. Both reset by `set_chr`.
    last_end: i64,
    chr_size: Option<i64>,
}

impl ValueSink {
    fn new(bin_size: i64, declared: Option<ChrMap>) -> Self {
        Self {
            bin_size,
            declared,
            chr: String::new(),
            bin: None,
            bin_sum: 0.0,
            bin_covered: 0,
            run_start_bin: None,
            run_values: Vec::new(),
            item_count: 0,
            clipped_count: 0,
            last_end: 0,
            chr_size: None,
        }
    }

    fn binning(&self) -> bool {
        self.bin_size > 0
    }

    /// Point the sink at a chromosome, closing what the last one left open.
    ///
    /// The declared size is resolved here for its *error*, which the writer
    /// would raise the same but later: with a bin size it sees no value until a
    /// bin closes and a run flushes, by which point the line carrying the name
    /// is thousands behind.
    fn set_chr(&mut self, writer: &mut BbiWriter, chr: &str) -> Result<()> {
        if self.chr == chr {
            return Ok(());
        }
        self.finish(writer)?;
        self.chr.clear();
        self.chr.push_str(chr);
        self.last_end = 0;
        self.chr_size = None;
        if let Some(declared) = &self.declared {
            self.chr_size = Some(declared.resolve(chr)?.size);
        }
        Ok(())
    }

    fn add(
        &mut self,
        writer: &mut BbiWriter,
        chr: &str,
        start: i64,
        end: i64,
        value: f32,
    ) -> Result<()> {
        self.set_chr(writer, chr)?;
        self.item_count += 1;
        if !self.binning() {
            let chr = std::mem::take(&mut self.chr);
            let result = writer.write_value(&chr, start, end, value);
            self.chr = chr;
            return result;
        }
        self.accumulate(writer, start, end, value)
    }

    fn add_run(
        &mut self,
        writer: &mut BbiWriter,
        chr: &str,
        start: i64,
        span: i64,
        values: &[f32],
    ) -> Result<()> {
        if values.is_empty() {
            return Ok(());
        }
        self.set_chr(writer, chr)?;
        self.item_count += values.len() as u64;
        if !self.binning() {
            let chr = std::mem::take(&mut self.chr);
            let result = writer.write_values(&chr, start, span, values);
            self.chr = chr;
            return result;
        }
        for (i, value) in values.iter().enumerate() {
            let i = i as i64;
            self.accumulate(writer, start + span * i, start + span * (i + 1), *value)?;
        }
        Ok(())
    }

    /// Spread one value over the bins it overlaps.
    fn accumulate(
        &mut self,
        writer: &mut BbiWriter,
        start: i64,
        mut end: i64,
        value: f32,
    ) -> Result<()> {
        if start < 0 {
            return Err(Error::invalid(format!("start {start} is negative")));
        }
        if end <= start {
            return Err(Error::invalid(format!(
                "end {end} is not past the start {start}"
            )));
        }
        if start < self.last_end {
            return Err(Error::invalid(format!(
                "{}:{start}-{end} starts before the end {} of the previous value, values \
                 must be added in order and without overlap",
                self.chr, self.last_end
            )));
        }
        if let Some(size) = self.chr_size {
            if end > size {
                if start >= size {
                    return Err(Error::invalid(format!(
                        "{}:{start}-{end} starts past the end of {}, which is {size} bases long",
                        self.chr, self.chr
                    )));
                }
                self.clipped_count += 1;
                end = size;
            }
        }
        self.last_end = end;
        let mut index = start / self.bin_size;
        while index * self.bin_size < end {
            if Some(index) != self.bin {
                self.close_bin(writer)?;
                self.bin = Some(index);
                self.bin_sum = 0.0;
                self.bin_covered = 0;
            }
            let overlap = end.min((index + 1) * self.bin_size) - start.max(index * self.bin_size);
            self.bin_sum += value as f64 * overlap as f64;
            self.bin_covered += overlap;
            index += 1;
        }
        Ok(())
    }

    /// Emit the open bin, into the run when it can extend it.
    fn close_bin(&mut self, writer: &mut BbiWriter) -> Result<()> {
        let Some(index) = self.bin.take() else {
            return Ok(());
        };
        let (covered, sum) = (self.bin_covered, self.bin_sum);
        self.bin_sum = 0.0;
        self.bin_covered = 0;
        if covered <= 0 {
            return Ok(());
        }
        let value = (sum / covered as f64) as f32;
        if self
            .run_start_bin
            .is_some_and(|first| index != first + self.run_values.len() as i64)
        {
            self.flush_run(writer)?;
        }
        if self.run_start_bin.is_none() {
            self.run_start_bin = Some(index);
        }
        self.run_values.push(value);
        if self.run_values.len() >= VALUE_BATCH {
            self.flush_run(writer)?;
        }
        Ok(())
    }

    fn flush_run(&mut self, writer: &mut BbiWriter) -> Result<()> {
        let Some(first) = self.run_start_bin.take() else {
            return Ok(());
        };
        if self.run_values.is_empty() {
            return Ok(());
        }
        let start = first * self.bin_size;
        let values = std::mem::take(&mut self.run_values);
        let chr = std::mem::take(&mut self.chr);
        let result = writer.write_values(&chr, start, self.bin_size, &values);
        self.chr = chr;
        self.run_values = values;
        self.run_values.clear();
        result
    }

    /// Close the open bin and the open run, at the end of the input.
    fn finish(&mut self, writer: &mut BbiWriter) -> Result<()> {
        self.close_bin(writer)?;
        self.flush_run(writer)
    }
}

// -- convert_to_bigwig -----------------------------------------------------

/// Convert a bedGraph or WIG file into a bigWig.
///
/// Which of the two the input is comes from its content, not its name. A
/// `bin_size` of `None` writes the values as they stand; a positive one
/// averages them into bins of that width, weighting each value by the bases it
/// covers.
pub fn convert_to_bigwig(
    input: &Path,
    output: &Path,
    bin_size: Option<i64>,
    mut options: BbiWriterOptions,
    progress: Option<ProgressFn>,
    cancel: Option<CancelFlag>,
) -> Result<ConvertResult> {
    let bin_size = bin_size.unwrap_or(0);
    if bin_size < 0 {
        return Err(Error::invalid(format!(
            "bin_size {bin_size} must not be negative"
        )));
    }
    options.kind = BbiKind::BigWig;
    let declared = options.chr_sizes.clone();

    let mut reader = LineReader::open(input)?;
    let mut writer = BbiWriter::create(&output.to_string_lossy(), options)?;
    let mut sink = ValueSink::new(bin_size, declared);
    let tracker = ProgressTracker::with_callback(reader.size, progress);

    let mut format: Option<TextFormat> = None;
    let mut line_count = 0u64;
    let mut declaration = WigDeclaration::default();
    let mut declared_yet = false;
    // The run of contiguous fixedStep values being gathered, and where it
    // started. Only a fixedStep whose step is its span can build one: anything
    // else leaves gaps or overlaps, which a run has neither of.
    let mut run: Vec<f32> = Vec::new();
    let mut run_start = 0i64;
    let mut reported = 0u64;

    // The body is a closure so the writer can be abandoned on any failure: a
    // conversion promises its caller either a whole file or none, and a
    // `BbiWriter` left to drop would close itself and leave a valid bigWig
    // holding however much had been read.
    let mut line = String::new();
    let result = (|| -> Result<()> {
        while reader.read_line_into(&mut line)? {
            line_count += 1;
            if line_count % PROGRESS_INTERVAL == 0 {
                tracker.add(reader.consumed - reported);
                reported = reader.consumed;
                // Checked here rather than per line: the callback is the only
                // place a cancellation can come from, so nothing can have
                // changed in between.
                if cancel.as_ref().is_some_and(CancelFlag::is_cancelled) {
                    return Err(Error::invalid("conversion cancelled"));
                }
            }
            if is_skipped_line(&line) {
                continue;
            }
            let fields = split_blanks(&line);
            if fields.is_empty() {
                continue;
            }
            let declaration_line = is_declaration(fields[0]);

            if format.is_none() {
                format = Some(sniff_format(&line).map_err(|e| reader.fail(e))?);
            }

            if format == Some(TextFormat::Wig) {
                if declaration_line {
                    flush_run(&mut sink, &mut writer, &declaration, run_start, &mut run)
                        .map_err(|e| reader.guard(e))?;
                    declaration = parse_wig_declaration(&line).map_err(|e| reader.fail(e))?;
                    declared_yet = true;
                    continue;
                }
                if !declared_yet {
                    return Err(
                        reader.fail("wig data before any fixedStep or variableStep declaration")
                    );
                }
                if declaration.fixed_step {
                    if fields.len() != 1 {
                        return Err(reader.fail(format!(
                            "fixedStep data has {} columns, not 1",
                            fields.len()
                        )));
                    }
                    let value = parse_f32(fields[0]).map_err(|e| reader.fail(e))?;
                    // A contiguous run is what write_values takes in one
                    // extend; anything else is handed over a value at a time,
                    // which the writer encodes just as tightly but has to walk.
                    if declaration.step == declaration.span {
                        if run.is_empty() {
                            run_start = declaration.start;
                        }
                        run.push(value);
                        if run.len() >= VALUE_BATCH {
                            flush_run(&mut sink, &mut writer, &declaration, run_start, &mut run)
                                .map_err(|e| reader.guard(e))?;
                        }
                    } else {
                        sink.add(
                            &mut writer,
                            &declaration.chr,
                            declaration.start,
                            declaration.start + declaration.span,
                            value,
                        )
                        .map_err(|e| reader.guard(e))?;
                    }
                    declaration.start += declaration.step;
                } else {
                    if fields.len() != 2 {
                        return Err(reader.fail(format!(
                            "variableStep data has {} columns, not 2",
                            fields.len()
                        )));
                    }
                    let start = parse_i64(fields[0], "position").map_err(|e| reader.fail(e))?;
                    let value = parse_f32(fields[1]).map_err(|e| reader.fail(e))?;
                    if start < 1 {
                        return Err(reader.fail(format!("position {start} is not 1-based")));
                    }
                    sink.add(
                        &mut writer,
                        &declaration.chr,
                        start - 1,
                        start - 1 + declaration.span,
                        value,
                    )
                    .map_err(|e| reader.guard(e))?;
                }
                continue;
            }

            if declaration_line {
                return Err(
                    reader.fail("wig declaration in what has been read as a bedgraph so far")
                );
            }
            if fields.len() != 4 {
                return Err(reader.fail(format!(
                    "bedgraph record has {} columns, not 4",
                    fields.len()
                )));
            }
            let start = parse_i64(fields[1], "start").map_err(|e| reader.fail(e))?;
            let end = parse_i64(fields[2], "end").map_err(|e| reader.fail(e))?;
            let value = parse_f32(fields[3]).map_err(|e| reader.fail(e))?;
            sink.add(&mut writer, fields[0], start, end, value)
                .map_err(|e| reader.guard(e))?;
        }

        flush_run(&mut sink, &mut writer, &declaration, run_start, &mut run)
            .map_err(|e| reader.guard(e))?;
        sink.finish(&mut writer).map_err(|e| reader.guard(e))?;
        writer.close()
    })();

    if let Err(error) = result {
        writer.abandon();
        return Err(error);
    }
    tracker.done_report();

    Ok(ConvertResult {
        // An input holding no data at all still produces a valid, empty bigWig;
        // it simply never settled on a format.
        format: format.unwrap_or(TextFormat::BedGraph),
        line_count,
        item_count: sink.item_count,
        skipped_count: writer.skipped_count(),
        // Only one of the two sees any given value — with a bin size the writer
        // is handed bins and never the values they came from — so the report
        // adds them.
        clipped_count: writer.clipped_count() + sink.clipped_count,
        chr_sizes: writer.chr_sizes(),
    })
}

fn flush_run(
    sink: &mut ValueSink,
    writer: &mut BbiWriter,
    declaration: &WigDeclaration,
    run_start: i64,
    run: &mut Vec<f32>,
) -> Result<()> {
    if run.is_empty() {
        return Ok(());
    }
    let result = sink.add_run(writer, &declaration.chr, run_start, declaration.span, run);
    run.clear();
    result
}

fn parse_i64(text: &str, what: &str) -> Result<i64> {
    text.parse()
        .map_err(|_| Error::invalid(format!("could not read \"{text}\" as a {what}")))
}

fn parse_f32(text: &str) -> Result<f32> {
    text.parse::<f64>()
        .map(|v| v as f32)
        .map_err(|_| Error::invalid(format!("could not read \"{text}\" as a number")))
}

// -- convert_to_bigbed -----------------------------------------------------

/// The types the standard BED columns carry, parallel to `BED_FIELD_NAMES`.
///
/// itemRgb, blockSizes and blockStarts hold comma-separated lists — autoSql
/// calls them `uint[3]` and `int[]` — and this library writes only the four
/// scalar types, so they go out as strings. A reader gets the text either way.
const BED_FIELD_STANDARD_TYPES: &[&str] = &[
    "string", "uint", "uint", "string", "uint", "string", "uint", "uint", "string", "uint",
    "string", "string",
];

/// Names and types for a BED of `col_count` columns that says nothing about
/// itself.
///
/// The standard BED columns as far as the file has them, then `field13` and up
/// — the same naming a bigBed carrying no autoSql is read with, so a BED
/// converted here and a bigBed read there describe their columns alike.
pub fn default_bed_fields(col_count: usize) -> IndexMap<String, String> {
    (0..col_count)
        .map(|index| {
            if index < BED_FIELD_NAMES.len() {
                (
                    BED_FIELD_NAMES[index].to_string(),
                    BED_FIELD_STANDARD_TYPES[index].to_string(),
                )
            } else {
                (format!("field{}", index + 1), "string".to_string())
            }
        })
        .collect()
}

/// Convert a BED file into a bigBed.
///
/// Entries may overlap and nest, which is the ordinary shape of a BED, so only
/// their starts have to be in order. Chromosomes still have to be pooled.
pub fn convert_to_bigbed(
    input: &Path,
    output: &Path,
    mut options: BbiWriterOptions,
    progress: Option<ProgressFn>,
    cancel: Option<CancelFlag>,
) -> Result<ConvertResult> {
    options.kind = BbiKind::BigBed;
    let mut reader = LineReader::open(input)?;
    let tracker = ProgressTracker::with_callback(reader.size, progress);

    // The writer needs the columns at construction, since a bigBed's autoSql is
    // written before any record is. When they were not given they come from the
    // first record, so the writer is opened lazily — which also means a BED that
    // turns out to be malformed leaves no file behind at all.
    let mut writer: Option<BbiWriter> = None;
    let mut declared_fields = std::mem::take(&mut options.fields);
    let mut col_count = 0usize;
    let mut line_count = 0u64;
    let mut reported = 0u64;
    let mut values: IndexMap<String, String> = IndexMap::new();

    let mut line = String::new();
    let result = (|| -> Result<()> {
        while reader.read_line_into(&mut line)? {
            line_count += 1;
            if line_count % PROGRESS_INTERVAL == 0 {
                tracker.add(reader.consumed - reported);
                reported = reader.consumed;
                // Checked here rather than per line: the callback is the only
                // place a cancellation can come from, so nothing can have
                // changed in between.
                if cancel.as_ref().is_some_and(CancelFlag::is_cancelled) {
                    return Err(Error::invalid("conversion cancelled"));
                }
            }
            if is_skipped_line(&line) {
                continue;
            }
            let mut fields = split_tabs(&line);
            trim_trailing_tab(
                &mut fields,
                if writer.is_none() {
                    (!declared_fields.is_empty()).then(|| declared_fields.len())
                } else {
                    Some(col_count)
                },
            );
            if fields.len() < 3 {
                return Err(reader.fail(format!(
                    "bed record has {} tab-separated columns, and needs at least 3 \
                     (chrom, chromStart, chromEnd)",
                    fields.len()
                )));
            }

            if writer.is_none() {
                col_count = fields.len();
                if declared_fields.is_empty() {
                    declared_fields = default_bed_fields(col_count);
                } else if declared_fields.len() != col_count {
                    return Err(reader.fail(format!(
                        "fields declares {} columns and the first record has {col_count}",
                        declared_fields.len()
                    )));
                }
                let mut opened = BbiWriterOptions {
                    fields: declared_fields.clone(),
                    ..clone_options(&options)
                };
                opened.kind = BbiKind::BigBed;
                writer = Some(
                    BbiWriter::create(&output.to_string_lossy(), opened)
                        .map_err(|e| reader.fail(e))?,
                );
                for name in declared_fields.keys().skip(3) {
                    values.insert(name.clone(), String::new());
                }
            } else if fields.len() != col_count {
                return Err(reader.fail(format!(
                    "bed record has {} columns and the first one had {col_count}; a bigbed \
                     stores one shape of record",
                    fields.len()
                )));
            }

            let start = parse_i64(fields[1], "chromStart").map_err(|e| reader.fail(e))?;
            let end = parse_i64(fields[2], "chromEnd").map_err(|e| reader.fail(e))?;
            // Overwritten in place rather than rebuilt: an `IndexMap` keeps an
            // entry where it was first inserted, so a record past the first
            // costs no allocation and no rehash.
            for (index, slot) in values.values_mut().enumerate() {
                slot.clear();
                slot.push_str(fields[index + 3]);
            }
            writer
                .as_mut()
                .expect("opened above")
                .write_entry(fields[0], start, end, &values)
                .map_err(|e| reader.guard(e))?;
        }
        Ok(())
    })();

    if let Err(error) = result {
        if let Some(writer) = &mut writer {
            writer.abandon();
        }
        return Err(error);
    }

    // A BED holding no record at all still produces a valid, empty bigBed, with
    // whatever columns were asked for or the standard bed3.
    let mut writer = match writer {
        Some(writer) => writer,
        None => {
            if declared_fields.is_empty() {
                declared_fields = default_bed_fields(3);
            }
            let mut opened = BbiWriterOptions {
                fields: declared_fields,
                ..clone_options(&options)
            };
            opened.kind = BbiKind::BigBed;
            BbiWriter::create(&output.to_string_lossy(), opened)?
        }
    };
    writer.close()?;
    tracker.done_report();

    Ok(ConvertResult {
        format: TextFormat::Bed,
        line_count,
        item_count: writer.entry_count(),
        skipped_count: writer.skipped_count(),
        clipped_count: writer.clipped_count(),
        chr_sizes: writer.chr_sizes(),
    })
}

/// `BbiWriterOptions` is not `Clone` — a `ChrMap` in it is, but the struct is
/// built once per writer — and the bigBed path opens its writer lazily, so it
/// needs the options twice.
fn clone_options(options: &BbiWriterOptions) -> BbiWriterOptions {
    BbiWriterOptions {
        kind: options.kind,
        chr_sizes: options.chr_sizes.clone(),
        fields: options.fields.clone(),
        items_per_slot: options.items_per_slot,
        block_size: options.block_size,
        compression_level: options.compression_level,
        parallel: options.parallel,
        section_policy: options.section_policy,
        cost_model: options.cost_model,
    }
}