R/genome-build.R

In misha: Toolkit for Analysis of Genomic Data

# Build a misha genome from a name.
#
# Public surface (exported):
#   gdb.build_genome()           - build from registry-resolved recipe
#   gdb.list_genomes()           - list resolvable genome names
#   gdb.genome_info()            - show resolved recipe without building
#   gdb.install_gff3_converter() - pre-install UCSC's gff3ToGenePred binary

# ---------------------------------------------------------------------------
# Internal constants
# ---------------------------------------------------------------------------

.MISHA_GENOME_SOURCES <- c("ucsc", "ucsc-hub", "ncbi", "s3", "manual", "local")

.UCSC_GOLDENPATH <- "https://hgdownload.soe.ucsc.edu/goldenPath"
.NCBI_DATASETS_API <- "https://api.ncbi.nlm.nih.gov/datasets/v2"

# NCBI backend chromosome naming options. See ?gdb.build_genome.
.NCBI_CHROM_NAMING_VALUES <- c("sequence_name", "ucsc", "accession")
.NCBI_DEFAULT_CHROM_NAMING <- "sequence_name"

# RefSeqLink columns we expose as gene annotations (used by UCSC backend).
# Schema: 19 columns as ncbiRefSeqLink.txt.gz ships them. The C++ importer
# requires the count to match exactly.
.UCSC_NCBI_REFSEQ_LINK_COLS <- c(
    "id", "status", "name", "product", "mrnaAcc", "protAcc", "locusLinkId",
    "omimId", "hgnc", "genbank", "pseudo", "gbkey", "source", "gene_biotype",
    "gene_synonym", "ncrna_class", "note", "description", "externalId"
)

# ---------------------------------------------------------------------------
# Registry parsing and resolution
# ---------------------------------------------------------------------------

# Normalize a raw registry entry to a recipe list with $source set.
# Accepts:
#   - character scalar  -> {source: local, path: <string>}  (legacy form)
#   - named list with $source
.normalize_recipe <- function(entry, name) {
    if (is.character(entry) && length(entry) == 1) {
        return(list(source = "local", path = entry))
    }
    if (!is.list(entry)) {
        stop(sprintf(
            "Genome '%s' has invalid registry entry: expected string or mapping, got %s",
            name, class(entry)[[1]]
        ), call. = FALSE)
    }
    if (is.null(entry$source) || !is.character(entry$source) || length(entry$source) != 1) {
        stop(sprintf("Genome '%s' is missing required 'source:' field", name), call. = FALSE)
    }
    if (!entry$source %in% .MISHA_GENOME_SOURCES) {
        stop(sprintf(
            "Genome '%s' has unknown source '%s'. Valid sources: %s",
            name, entry$source, paste(.MISHA_GENOME_SOURCES, collapse = ", ")
        ), call. = FALSE)
    }
    entry
}

.validate_recipe <- function(recipe, name) {
    src <- recipe$source
    miss <- function(field) {
        stop(sprintf("Genome '%s' (source: %s) is missing required field '%s'", name, src, field), call. = FALSE)
    }
    if (src == "ucsc") {
        if (is.null(recipe$assembly)) miss("assembly")
    } else if (src == "ucsc-hub") {
        if (is.null(recipe$accession)) miss("accession")
        if (!grepl("^GC[FA]_[0-9]+\\.[0-9]+$", recipe$accession)) {
            stop(sprintf(
                "Genome '%s': accession '%s' does not match GC[FA]_<digits>.<digits>",
                name, recipe$accession
            ), call. = FALSE)
        }
        if (!is.null(recipe$chrom_naming)) {
            if (!is.character(recipe$chrom_naming) ||
                length(recipe$chrom_naming) != 1L ||
                !nzchar(recipe$chrom_naming)) {
                stop(sprintf(
                    "Genome '%s': chrom_naming must be a non-empty single string.",
                    name
                ), call. = FALSE)
            }
        }
    } else if (src == "ncbi") {
        if (is.null(recipe$accession)) miss("accession")
        if (!grepl("^GC[FA]_[0-9]+\\.[0-9]+$", recipe$accession)) {
            stop(sprintf(
                "Genome '%s': accession '%s' does not match GC[FA]_<digits>.<digits>",
                name, recipe$accession
            ), call. = FALSE)
        }
        if (!is.null(recipe$chrom_naming) &&
            !recipe$chrom_naming %in% .NCBI_CHROM_NAMING_VALUES) {
            stop(sprintf(
                "Genome '%s': chrom_naming '%s' invalid. Valid values: %s",
                name, recipe$chrom_naming,
                paste(.NCBI_CHROM_NAMING_VALUES, collapse = ", ")
            ), call. = FALSE)
        }
    } else if (src == "s3") {
        if (is.null(recipe$assembly)) miss("assembly")
    } else if (src == "local") {
        if (is.null(recipe$path)) miss("path")
    } else if (src == "manual") {
        if (is.null(recipe$fasta)) miss("fasta")
    }
    invisible(recipe)
}

# Read a single registry YAML file. Returns a named list <name> -> recipe (raw,
# not yet validated).
.parse_genome_registry <- function(path) {
    if (!file.exists(path)) {
        stop(sprintf("Registry file does not exist: %s", path), call. = FALSE)
    }
    y <- tryCatch(
        yaml::read_yaml(path),
        error = function(e) stop(sprintf("Failed to parse YAML registry %s: %s", path, conditionMessage(e)), call. = FALSE)
    )
    if (is.null(y$genome)) {
        return(list())
    }
    if (!is.list(y$genome) || is.null(names(y$genome))) {
        stop(sprintf("Registry %s: 'genome:' must be a named mapping", path), call. = FALSE)
    }
    y$genome
}

# Walk up from getwd() to git root looking for a misha.yaml.
.find_project_misha_yaml <- function() {
    dir <- normalizePath(getwd(), mustWork = FALSE)
    while (TRUE) {
        candidate <- file.path(dir, "misha.yaml")
        if (file.exists(candidate)) {
            return(candidate)
        }
        if (file.exists(file.path(dir, ".git"))) {
            return(NULL)
        }
        parent <- dirname(dir)
        if (parent == dir) {
            return(NULL)
        }
        dir <- parent
    }
}

.builtin_registry_path <- function() {
    system.file("genomes.yaml", package = "misha")
}

# Resolve `name` through the chain. Returns a list:
#   list(recipe = <list>, resolved_from = <character>)
.resolve_genome <- function(name, registry = NULL) {
    if (!is.character(name) || length(name) != 1 || !nzchar(name)) {
        stop("name must be a non-empty string", call. = FALSE)
    }

    sources <- list()
    if (!is.null(registry)) {
        sources[[length(sources) + 1]] <- list(label = sprintf("registry arg (%s)", registry), path = registry)
    }
    opt <- getOption("misha.genome_registry")
    if (!is.null(opt)) {
        sources[[length(sources) + 1]] <- list(label = sprintf("getOption('misha.genome_registry') (%s)", opt), path = opt)
    }
    proj <- .find_project_misha_yaml()
    if (!is.null(proj)) {
        sources[[length(sources) + 1]] <- list(label = sprintf("project misha.yaml (%s)", proj), path = proj)
    }
    builtin <- .builtin_registry_path()
    if (nzchar(builtin)) {
        sources[[length(sources) + 1]] <- list(label = "built-in (inst/genomes.yaml)", path = builtin)
    }

    for (src in sources) {
        entries <- tryCatch(.parse_genome_registry(src$path), error = function(e) {
            stop(sprintf("Error reading %s: %s", src$label, conditionMessage(e)), call. = FALSE)
        })
        if (name %in% names(entries)) {
            recipe <- .normalize_recipe(entries[[name]], name)
            .validate_recipe(recipe, name)
            return(list(recipe = recipe, resolved_from = src$label))
        }
    }

    # Pattern fallback for UCSC mammal hub accessions.
    if (grepl("^GC[FA]_[0-9]+\\.[0-9]+$", name)) {
        recipe <- list(source = "ucsc-hub", accession = name)
        return(list(
            recipe = recipe,
            resolved_from = "pattern fallback (UCSC mammal hub accession)"
        ))
    }

    stop(sprintf(
        "Genome '%s' not found in any registry. Searched: %s. To define it, add an entry to a misha.yaml or use gdb.create() directly.",
        name, paste(vapply(sources, `[[`, character(1), "label"), collapse = "; ")
    ), call. = FALSE)
}

# ---------------------------------------------------------------------------
# Helpers for downloads + post-build annotation loading
# ---------------------------------------------------------------------------

# Download URL to a destination (binary mode). Returns dest path.
.download_to <- function(url, dest, verbose = TRUE) {
    if (verbose) message(sprintf("Downloading %s ...", url))
    utils::download.file(url, dest, mode = "wb", quiet = !verbose)
    dest
}

# Heal UCSC TSV escape artifacts in description/note fields:
#   - backslash + (CR? + LF)  -> space    (line continuation)
#   - backslash + tab          -> space    (intra-field tab escape)
#   - stray CR (carriage return) -> nothing  (Windows line endings R's
#     readLines() would otherwise treat as line terminators)
#
# Loads the whole file into memory; UCSC tables are small (a few MB) so
# this is fine in practice and lets us do the regex fixups on the full
# string without chunk-boundary worries.
.heal_ucsc_tsv_escapes <- function(in_path, out_path) {
    con_in <- if (grepl("\\.gz$", in_path)) gzfile(in_path, "rb") else file(in_path, "rb")
    on.exit(close(con_in), add = TRUE)
    chunks <- list()
    repeat {
        chunk <- readBin(con_in, raw(), n = 4L * 1024L * 1024L)
        if (!length(chunk)) break
        chunks[[length(chunks) + 1L]] <- chunk
    }
    s <- if (length(chunks)) rawToChar(unlist(chunks, use.names = FALSE)) else ""
    # Order matters: do CRLF before LF so we consume the whole sequence.
    s <- gsub("\\\\\r\n", " ", s, perl = TRUE)
    s <- gsub("\\\\\n", " ", s, perl = TRUE)
    s <- gsub("\\\\\t", " ", s, perl = TRUE)
    s <- gsub("\r", "", s, perl = TRUE, fixed = FALSE)
    con_out <- file(out_path, "wb")
    on.exit(close(con_out), add = TRUE)
    writeBin(charToRaw(s), con_out)
    out_path
}

# Normalize a UCSC TSV table to exactly N columns per row. UCSC sometimes
# ships rows with stray embedded tabs (e.g. in description fields) - these
# lines have NF != expected count and would crash the C++ importer's strict
# column-count check. Short rows are padded with empty strings; long rows
# have their trailing extras joined back into the last column with " ".
.normalize_ucsc_tsv <- function(in_path, out_path, n_cols) {
    con_in <- if (grepl("\\.gz$", in_path)) gzfile(in_path, "rt") else file(in_path, "rt")
    on.exit(close(con_in), add = TRUE)
    con_out <- file(out_path, "wt")
    on.exit(close(con_out), add = TRUE)
    repeat {
        lines <- readLines(con_in, n = 50000L, warn = FALSE)
        if (!length(lines)) break
        fields <- strsplit(lines, "\t", fixed = TRUE)
        out <- vapply(fields, function(x) {
            if (length(x) == n_cols) {
                paste(x, collapse = "\t")
            } else if (length(x) < n_cols) {
                paste(c(x, rep("", n_cols - length(x))), collapse = "\t")
            } else {
                # Join overflow into the last column so the row has exactly n_cols fields.
                paste(c(x[seq_len(n_cols - 1)], paste(x[n_cols:length(x)], collapse = " ")),
                    collapse = "\t"
                )
            }
        }, character(1))
        writeLines(out, con_out)
    }
    out_path
}

# Trim a genePred-format file to the classic 12-col layout the C++
# gintervals_import_genes expects. Two input shapes are handled:
#   - 16 cols: UCSC extended genePred with leading bin column (ncbiRefSeq,
#     refGene, knownGene from goldenPath/database/). Take cols 2-13.
#   - 15 cols: extended genePred without bin (gff3ToGenePred output, NCBI).
#     Take cols 1-12.
# In both cases the resulting 12 cols are:
#   [name, chrom, strand, txStart, txEnd, cdsStart, cdsEnd, exonCount,
#    exonStarts, exonEnds, score, name2]
# Cols 11-12 occupy the C++ importer's PROTEINID/ALIGNID slots - read but
# unused, so populating with score/name2 is harmless.
#
# Accepts .txt or .txt.gz, writes a .txt with 12-col content. Streams the
# input so memory stays bounded.
.normalize_ucsc_genepred <- function(in_path, out_path) {
    con_in <- if (grepl("\\.gz$", in_path)) gzfile(in_path, "rt") else file(in_path, "rt")
    on.exit(close(con_in), add = TRUE)
    con_out <- file(out_path, "wt")
    on.exit(close(con_out), add = TRUE)

    chunk_size <- 50000L
    take_range <- NULL
    repeat {
        lines <- readLines(con_in, n = chunk_size, warn = FALSE)
        if (!length(lines)) break
        fields <- strsplit(lines, "\t", fixed = TRUE)
        nf <- vapply(fields, length, integer(1))
        if (is.null(take_range)) {
            shape <- nf[[1]]
            take_range <- if (shape == 16L) {
                2:13
            } else if (shape == 15L) {
                1:12
            } else if (shape == 12L) {
                1:12
            } else {
                stop(sprintf(
                    "genePred file %s: unsupported column count %d (expected 12, 15, or 16)",
                    in_path, shape
                ), call. = FALSE)
            }
        }
        bad <- nf < max(take_range)
        if (any(bad)) {
            stop(sprintf(
                "genePred file %s: row %d has only %d fields, need at least %d",
                in_path, which(bad)[1], nf[which(bad)[1]], max(take_range)
            ), call. = FALSE)
        }
        trimmed <- vapply(
            fields, function(x) paste(x[take_range], collapse = "\t"),
            character(1)
        )
        writeLines(trimmed, con_out)
    }
    out_path
}

# Decompress a .gz file (file -> file with .gz removed). Returns the new path.
.gunzip_to_file <- function(gz_path, out_path = sub("\\.gz$", "", gz_path)) {
    con_in <- gzfile(gz_path, "rb")
    on.exit(close(con_in), add = TRUE)
    con_out <- file(out_path, "wb")
    on.exit(close(con_out), add = TRUE)
    repeat {
        chunk <- readBin(con_in, raw(), n = 1024 * 1024)
        if (length(chunk) == 0) break
        writeBin(chunk, con_out)
    }
    out_path
}

# Parse UCSC rmsk.txt(.gz) - 17 columns. Returns a data.frame with intervals
# columns plus name/class/family. Strand is normalized to numeric (1/-1/0).
.parse_ucsc_rmsk <- function(file) {
    cols <- c(
        "bin", "swScore", "milliDiv", "milliDel", "milliIns",
        "genoName", "genoStart", "genoEnd", "genoLeft", "strand",
        "repName", "repClass", "repFamily", "repStart", "repEnd",
        "repLeft", "id"
    )
    df <- utils::read.table(file,
        sep = "\t", header = FALSE, col.names = cols,
        quote = "", comment.char = "", stringsAsFactors = FALSE,
        na.strings = character(0)
    )
    data.frame(
        chrom = df$genoName,
        start = df$genoStart,
        end = df$genoEnd,
        strand = ifelse(df$strand == "+", 1L, ifelse(df$strand == "-", -1L, 0L)),
        name = df$repName,
        class = df$repClass,
        family = df$repFamily,
        stringsAsFactors = FALSE
    )
}

.parse_ucsc_cpg_island <- function(file) {
    cols <- c(
        "bin", "chrom", "chromStart", "chromEnd", "name",
        "length", "cpgNum", "gcNum", "perCpg", "perGc", "obsExp"
    )
    df <- utils::read.table(file,
        sep = "\t", header = FALSE, col.names = cols,
        quote = "", comment.char = "", stringsAsFactors = FALSE,
        na.strings = character(0)
    )
    data.frame(
        chrom = df$chrom,
        start = df$chromStart,
        end = df$chromEnd,
        name = df$name,
        length = df$length,
        cpgNum = df$cpgNum,
        perCpg = df$perCpg,
        perGc = df$perGc,
        obsExp = df$obsExp,
        stringsAsFactors = FALSE
    )
}

.parse_ucsc_cytoband <- function(file) {
    cols <- c("chrom", "chromStart", "chromEnd", "name", "gieStain")
    df <- utils::read.table(file,
        sep = "\t", header = FALSE, col.names = cols,
        quote = "", comment.char = "", stringsAsFactors = FALSE,
        na.strings = character(0)
    )
    data.frame(
        chrom = df$chrom,
        start = df$chromStart,
        end = df$chromEnd,
        name = df$name,
        stain = df$gieStain,
        stringsAsFactors = FALSE
    )
}

# Write genome_info.yaml - a record of where this groot came from.
.write_genome_info <- function(groot, name, recipe, sets, files = list()) {
    info <- list(
        name = name,
        source = recipe$source,
        downloaded_at = format(Sys.time(), "%Y-%m-%dT%H:%M:%SZ", tz = "UTC"),
        misha_version = as.character(utils::packageVersion("misha")),
        sets = as.list(sets),
        recipe = recipe,
        files = files
    )
    .gwith_umask(yaml::write_yaml(info, file.path(groot, "genome_info.yaml")))
    invisible(NULL)
}


.ncbi_datasets_zip_url <- function(accession, include) {
    stopifnot(is.character(include), length(include) >= 1L)
    sprintf(
        "%s/genome/accession/%s/download?include_annotation_type=%s",
        .NCBI_DATASETS_API, accession, paste(include, collapse = ",")
    )
}

# Fetch NCBI Datasets /dataset_report for an accession. Returns the parsed
# top-level list (yaml::yaml.load handles JSON, which is a YAML subset, and
# is already in Imports). Throws on network/parse failure; callers wrap in
# tryCatch.
.ncbi_dataset_report <- function(accession, timeout = 30) {
    url <- sprintf(
        "%s/genome/accession/%s/dataset_report",
        .NCBI_DATASETS_API, accession
    )
    h <- curl::new_handle()
    curl::handle_setopt(h, timeout = timeout, useragent = "misha")
    resp <- curl::curl_fetch_memory(url, handle = h)
    if (resp$status_code >= 400) {
        stop(sprintf("dataset_report HTTP %d for %s", resp$status_code, accession),
            call. = FALSE
        )
    }
    yaml::yaml.load(rawToChar(resp$content))
}

# Pure: extract annotation provenance + organism info from a parsed
# dataset_report list. has_annotation is TRUE iff the assembly has any
# annotation provider on file (NCBI RefSeq or community submitter).
.ncbi_parse_annotation_info <- function(report) {
    reports <- report$reports %||% list()
    rep <- if (length(reports)) reports[[1L]] else list()
    ai <- rep$annotation_info %||% list()
    org <- rep$organism %||% list()
    asm <- rep$assembly_info %||% list()
    provider <- as.character(ai$provider %||% "")
    list(
        has_annotation  = nzchar(provider),
        provider        = provider,
        annotation_name = as.character(ai$name %||% ""),
        organism_name   = as.character(org$organism_name %||% ""),
        organism_tax_id = if (is.null(org$tax_id)) NA_integer_ else as.integer(org$tax_id),
        assembly_name   = as.character(asm$assembly_name %||% "")
    )
}

# Pure: decide what to do with the requested `sets` given the parsed
# annotation info. Trims 'genes' if the assembly has no annotation, and
# emits an actionable warning naming the accession (and any RefSeq-companion
# hint the caller resolved).
.ncbi_resolve_sets_with_preflight <- function(sets, info, accession, hint = "") {
    out <- list(sets = sets, warnings = character(0))
    if ("genes" %in% sets && !info$has_annotation) {
        msg <- sprintf(
            "NCBI accession %s has no annotation (annotation_info empty). 'genes' will be skipped.",
            accession
        )
        if (nzchar(hint)) {
            msg <- paste0(msg, "\n  Annotated alternative: ", hint)
        }
        out$warnings <- c(out$warnings, msg)
        out$sets <- setdiff(sets, "genes")
    }
    out
}

# Look up a RefSeq-annotated assembly for the same taxon and return a short
# hint string, e.g. "GCF_000003025.6 (NCBI Sus scrofa Annotation Release 106)".
# Returns "" on any failure or if the suggestion would be the same accession.
# Best-effort: callers must tolerate "".
.ncbi_suggest_annotated_alternative <- function(tax_id, current_accession,
                                                timeout = 15) {
    if (is.na(tax_id)) {
        return("")
    }
    url <- sprintf(
        "%s/genome/taxon/%d/dataset_report?filters.has_annotation=true&page_size=1",
        .NCBI_DATASETS_API, as.integer(tax_id)
    )
    h <- curl::new_handle()
    curl::handle_setopt(h, timeout = timeout, useragent = "misha")
    resp <- tryCatch(curl::curl_fetch_memory(url, handle = h), error = function(e) NULL)
    if (is.null(resp) || resp$status_code >= 400) {
        return("")
    }
    j <- tryCatch(yaml::yaml.load(rawToChar(resp$content)),
        error = function(e) NULL
    )
    reports <- j$reports %||% list()
    if (!length(reports)) {
        return("")
    }
    rep <- reports[[1L]]
    suggested <- as.character(rep$accession %||% "")
    if (!nzchar(suggested) || identical(suggested, current_accession)) {
        return("")
    }
    nm <- (rep$annotation_info %||% list())$name %||%
        (rep$annotation_info %||% list())$provider %||% ""
    if (nzchar(nm)) {
        sprintf("%s (%s)", suggested, nm)
    } else {
        suggested
    }
}

# Parse NCBI Datasets sequence_report.jsonl into a data.frame.
# Returns: data.frame(refseqAccession, genbankAccession, chrName, sequenceName,
# role, length).
.parse_ncbi_sequence_report <- function(path) {
    if (!file.exists(path)) {
        stop(sprintf("Sequence report not found: %s", path), call. = FALSE)
    }
    lines <- readLines(path, warn = FALSE)
    lines <- lines[nzchar(lines)]
    parsed <- lapply(lines, function(l) {
        # Tiny JSON parser via yaml (yaml is a JSON superset; already a dep).
        yaml::yaml.load(l)
    })
    pull <- function(field, default = NA_character_) {
        vapply(parsed, function(r) {
            v <- r[[field]]
            if (is.null(v) || !length(v)) default else as.character(v)[[1]]
        }, character(1))
    }
    pull_int <- function(field) {
        vapply(parsed, function(r) {
            v <- r[[field]]
            if (is.null(v) || !length(v)) NA_real_ else as.numeric(v)[[1]]
        }, numeric(1))
    }
    data.frame(
        refseqAccession  = pull("refseqAccession"),
        genbankAccession = pull("genbankAccession"),
        chrName          = pull("chrName"),
        sequenceName     = pull("sequenceName"),
        role             = pull("role"),
        length           = pull_int("length"),
        stringsAsFactors = FALSE
    )
}

# UCSC-style canonical name from a sequence_report row.
# Assembled molecules: chr1, chrX, chrM (UCSC uses chrM not chrMT).
# Unplaced: chrUn_<accession_underscore_dot_to_v>, e.g. NW_026256937.1 -> chrUn_NW026256937v1.
# Unlocalized: chr<chrom>_<acc>_random.
.ncbi_to_ucsc_name <- function(refseq_acc, chr_name, role) {
    # Build the UCSC-friendly suffix from an accession: drop underscore + dot,
    # encode version as v<N>: "NW_026256937.1" -> "NW026256937v1".
    encode_acc <- function(acc) {
        s <- gsub("_", "", acc, fixed = TRUE)
        s <- sub("\\.([0-9]+)$", "v\\1", s)
        s
    }
    if (role == "assembled-molecule") {
        if (chr_name %in% c("MT", "M")) {
            return("chrM")
        }
        return(paste0("chr", chr_name))
    }
    if (role == "unlocalized-scaffold") {
        return(paste0("chr", chr_name, "_", encode_acc(refseq_acc), "_random"))
    }
    paste0("chrUn_", encode_acc(refseq_acc))
}

# Synthesize a UCSC-shaped chromAlias data.frame from a parsed NCBI
# sequence_report. Five columns chosen so the existing column-detect +
# match_by_length passes in gdb.install_intervals can resolve any common
# groot naming convention (chr1.., NC_*, GenBank, bare chr name, HAL hybrid).
# Empty cells (e.g. assemblies without GenBank twin) stay as "" -- alias
# detection is bp-weighted and tolerates partial columns.
.ncbi_seqrep_to_alias_df <- function(seqrep) {
    sequence_name <- ifelse(seqrep$role == "assembled-molecule",
        seqrep$chrName,
        seqrep$sequenceName
    )
    ucsc <- mapply(.ncbi_to_ucsc_name,
        seqrep$refseqAccession, seqrep$chrName, seqrep$role,
        USE.NAMES = FALSE
    )
    data.frame(
        accession = seqrep$refseqAccession,
        genbank = seqrep$genbankAccession,
        sequence_name = sequence_name,
        chr_name = seqrep$chrName,
        ucsc = ucsc,
        stringsAsFactors = FALSE
    )
}

# Extract the <assembly_name> suffix for `accession` from an NCBI FTP parent-
# directory listing. NCBI Datasets /dataset_report returns {} for some older
# / suppressed accessions (e.g. GCF_000001635.26 GRCm38.p6) so the
# assembly_name field is missing; the FTP listing is the canonical fallback.
.ncbi_ftp_assembly_name_from_dir <- function(accession, listing) {
    pat <- sprintf("%s_([^/\"<> ]+)/", accession)
    m <- regmatches(listing, regexpr(pat, listing, perl = TRUE))
    if (!length(m)) {
        return("")
    }
    sub(pat, "\\1", m[[1L]])
}

# NCBI FTP genomes-all directory for an assembly. Accession's nine-digit
# numeric portion is split into three triplets:
#   GCF_000001405.40 -> GCF/000/001/405/GCF_000001405.40_GRCh38.p14
.ncbi_ftp_assembly_dir <- function(accession, assembly_name) {
    prefix <- substr(accession, 1, 3)
    digits <- sub("^GC[AF]_", "", accession)
    digits <- sub("\\..*$", "", digits)
    parts <- substring(digits, c(1, 4, 7), c(3, 6, 9))
    sprintf(
        "https://ftp.ncbi.nlm.nih.gov/genomes/all/%s/%s/%s/%s/%s_%s",
        prefix, parts[1], parts[2], parts[3], accession, assembly_name
    )
}

# Build a named character vector: original FASTA/GFF id -> target canonical
# chrom name, given the desired chrom_naming and a parsed sequence report.
# Always indexed by the refseqAccession used in the FASTA/GFF (the field NCBI
# uses as the seqid).
.build_ncbi_rename_map <- function(seqrep, chrom_naming) {
    chrom_naming <- match.arg(chrom_naming, .NCBI_CHROM_NAMING_VALUES)
    targets <- if (chrom_naming == "accession") {
        seqrep$refseqAccession
    } else if (chrom_naming == "sequence_name") {
        # Use NCBI chrName for assembled molecules; refseq accession for
        # everything else (chrName="Un" would collide).
        ifelse(seqrep$role == "assembled-molecule",
            seqrep$chrName,
            seqrep$refseqAccession
        )
    } else { # ucsc
        mapply(.ncbi_to_ucsc_name,
            seqrep$refseqAccession, seqrep$chrName, seqrep$role,
            USE.NAMES = FALSE
        )
    }
    if (anyDuplicated(targets)) {
        dups <- unique(targets[duplicated(targets)])
        stop(sprintf(
            "chrom_naming='%s' produced duplicate names: %s",
            chrom_naming, paste(utils::head(dups, 5), collapse = ", ")
        ), call. = FALSE)
    }
    setNames(targets, seqrep$refseqAccession)
}

# Stream-rewrite FASTA: replace headers ">acc ..." with ">new_name". Works
# on plain or .gz input; output is plain.
.rename_fasta_headers <- function(in_path, out_path, rename_map, verbose = TRUE) {
    con_in <- if (grepl("\\.gz$", in_path)) gzfile(in_path, "rt") else file(in_path, "rt")
    on.exit(close(con_in), add = TRUE)
    con_out <- file(out_path, "wt")
    on.exit(close(con_out), add = TRUE)
    n_renamed <- 0L
    n_unmapped <- 0L
    repeat {
        lines <- readLines(con_in, n = 100000L, warn = FALSE)
        if (!length(lines)) break
        is_header <- startsWith(lines, ">")
        if (any(is_header)) {
            headers <- lines[is_header]
            ids <- sub("^>([^[:space:]]+).*$", "\\1", headers, perl = TRUE)
            new <- rename_map[ids]
            # Treat both NA (key absent) and empty (key present but no target)
            # as "unmapped, keep original". UCSC chromAlias has gaps -- e.g.
            # Bos mutus's MT row has refseq=NC_006380.3 but genbank="". A
            # naive rename would rewrite the header to ">", which misha's
            # FASTA loader then surfaces as the default chrom name "contig".
            unmapped <- is.na(new) | !nzchar(new)
            new[unmapped] <- ids[unmapped]
            n_renamed <- n_renamed + sum(!unmapped)
            n_unmapped <- n_unmapped + sum(unmapped)
            lines[is_header] <- paste0(">", new)
        }
        writeLines(lines, con_out)
    }
    if (verbose) {
        message(sprintf(
            "  Renamed %d FASTA contigs (%d unmapped, kept original id).",
            n_renamed, n_unmapped
        ))
    }
    invisible(out_path)
}

# Stream-rewrite GFF3: replace seqid (col 1) using rename_map. Comment lines
# preserved unchanged. Lines whose seqid isn't in the map are kept as-is.
.rename_gff3_seqids <- function(in_path, out_path, rename_map, verbose = TRUE) {
    con_in <- if (grepl("\\.gz$", in_path)) gzfile(in_path, "rt") else file(in_path, "rt")
    on.exit(close(con_in), add = TRUE)
    con_out <- file(out_path, "wt")
    on.exit(close(con_out), add = TRUE)
    n_renamed <- 0L
    n_unmapped <- 0L
    repeat {
        lines <- readLines(con_in, n = 100000L, warn = FALSE)
        if (!length(lines)) break
        is_data <- nzchar(lines) & !startsWith(lines, "#")
        if (any(is_data)) {
            data_lines <- lines[is_data]
            tab_pos <- regexpr("\t", data_lines, fixed = TRUE)
            seqids <- ifelse(tab_pos > 0, substring(data_lines, 1, tab_pos - 1), data_lines)
            rest <- ifelse(tab_pos > 0, substring(data_lines, tab_pos), "")
            new <- rename_map[seqids]
            unmapped <- is.na(new)
            new[unmapped] <- seqids[unmapped]
            n_renamed <- n_renamed + sum(!unmapped)
            n_unmapped <- n_unmapped + sum(unmapped)
            lines[is_data] <- paste0(new, rest)
        }
        writeLines(lines, con_out)
    }
    if (verbose) {
        message(sprintf(
            "  Rewrote %d GFF3 records' seqids (%d unmapped passed through).",
            n_renamed, n_unmapped
        ))
    }
    invisible(out_path)
}

# ---------------------------------------------------------------------------
# Public API
# ---------------------------------------------------------------------------

#' Build a misha genome database from a name
#'
#' Builds a misha genomic database for a named assembly. Resolves the name
#' through the registry chain (or pattern-fallback for \code{GC[FA]_*}
#' accessions), downloads the FASTA, calls \code{\link{gdb.create}} to build
#' the seq-only groot, then dispatches to \code{\link{gdb.install_intervals}}
#' for the requested sets.
#'
#' For details on resolution, sources, sets, and chromosome-alias handling,
#' see \code{\link{gdb.install_intervals}}.
#'
#' @param name Genome name (registry key, alias, or \code{GC[FA]_*} accession).
#' @param path Output directory; must not exist.
#' @param registry Optional path to an explicit registry YAML.
#' @param sets Subset of \code{c("genes", "rmsk", "cgi", "cytoband")}.
#'   Empty vector \code{character(0)} = sequence-only build.
#' @param prefix Character scalar prepended to set names (see
#'   \code{\link{gdb.install_intervals}}).
#' @param gene_sets Named character vector mapping the four
#'   \code{gintervals.import_genes()} roles to on-disk set names; \code{NA} skips
#'   a role.
#' @param gtf_priority Character vector ordering GTF source preference.
#' @param min_coverage Minimum fraction of groot chroms that must appear in a
#'   chromAlias column for that column to be picked as canonical (forwarded to
#'   \code{\link{gdb.install_intervals}}). Default \code{1.0} (strict). Lower
#'   to e.g. \code{0.99} when a column has small gaps -- typical when a target
#'   column doesn't span every contig (e.g. UCSC's \code{genbank} column has
#'   no value for the mitochondrion in many hubs, leaving 1 stray chrom).
#'   Honored only by the \code{ucsc-hub} backend; supplying a non-default value
#'   for any other source is an error (raised before any download).
#' @param match_by_length Forwarded to \code{\link{gdb.install_intervals}}.
#'   When \code{TRUE} (default), complements column-based canonical detection
#'   with a per-row length match for alias rows the chosen column couldn't
#'   cover, and switches asset translation to a cross-column per-row lookup
#'   so GFFs in any naming scheme import cleanly. Set \code{FALSE} for the
#'   stricter single-column-only behavior.
#' @param target_chroms Optional character vector of chrom names the resulting
#'   groot should align to (typically the output of \code{halStats
#'   --sequenceStats}, the chrom names in a HAL file you intend to liftover
#'   against). When supplied, misha auto-picks the chromAlias column whose
#'   values cover \code{target_chroms} best and uses that column as the
#'   canonical naming, instead of \code{chrom_naming}. Honored only by the
#'   \code{ucsc-hub} backend; supplying it with any other source is an error
#'   (raised before any download).
#' @param target_lengths Optional numeric vector aligned with
#'   \code{target_chroms} (typically the second field of
#'   \code{halStats --sequenceStats}). When supplied alongside
#'   \code{target_chroms} and with \code{match_by_length = TRUE}, this is
#'   the strong-guarantee path: misha force-aligns the hub FASTA to
#'   \code{target_chroms}, placing every target on its chromAlias row by
#'   name match across columns or unique-on-both-sides length pairing. If
#'   any target can't be placed, the build errors (in the pre-flight,
#'   before the multi-GB FASTA download). On success the resulting groot's
#'   chrom names are exactly \code{target_chroms} (alias rows not in
#'   \code{target_chroms} keep their original FASTA-header accession).
#'   Honored only by the \code{ucsc-hub} backend.
#' @param chrom_naming Optional override for the recipe's \code{chrom_naming}.
#'   Selects which name space the canonical chrom names should come from. For
#'   \code{ucsc-hub}: any chromAlias column (\code{"ucsc"}, \code{"genbank"},
#'   \code{"refseq"}, \code{"ncbi"}), plus the friendly aliases
#'   \code{"sequence_name"} (= \code{"assembly"}) and \code{"accession"} (keep
#'   the FASTA's source column). For \code{ncbi}: \code{"sequence_name"} (default),
#'   \code{"ucsc"}, or \code{"accession"}. \code{NULL} (default) keeps whatever
#'   the recipe specifies.
#' @param format \code{"indexed"} or \code{"per-chromosome"}; \code{NULL} =>
#'   \code{getOption("gmulticontig.indexed_format", TRUE)}.
#' @param verbose If \code{TRUE}, prints progress.
#' @return None (invisible \code{NULL}). The installed gene-derived sets
#'   (\code{tss}, \code{exons}, \code{utr3}, \code{utr5}) carry a \code{name}
#'   column (transcript/RNA accession) and a \code{geneName} column (gene symbol
#'   from the source annotation; blank when the source has no symbol).
#'
#' @seealso \code{\link{gdb.install_intervals}}, \code{\link{gdb.create}},
#'   \code{\link{gdb.list_genomes}}, \code{\link{gdb.genome_info}}.
#'
#' @examples
#' \dontrun{
#' gdb.build_genome("hg38", path = "~/genomes/hg38")
#' gdb.build_genome("GCA_004023825.1",
#'     path   = "~/genomes/arctic_fox",
#'     prefix = "intervs.global."
#' )
#' # Match HAL/Cactus canonical names (GenBank accessions like JH880237.1):
#' gdb.build_genome("GCF_000298355.1",
#'     path         = "~/genomes/Bos_mutus",
#'     chrom_naming = "genbank",
#'     prefix       = "intervs.global."
#' )
#' }
#'
#' @export
gdb.build_genome <- function(name,
                             path = name,
                             registry = NULL,
                             sets = c("genes", "rmsk", "cgi", "cytoband"),
                             prefix = "",
                             gene_sets = c(
                                 tss = "tss", exons = "exons",
                                 utr3 = "utr3", utr5 = "utr5"
                             ),
                             gtf_priority = c(
                                 "ncbiRefSeq", "bestRefSeq",
                                 "ensGene", "augustus", "xenoRefGene"
                             ),
                             chrom_naming = NULL,
                             target_chroms = NULL,
                             target_lengths = NULL,
                             min_coverage = 1.0,
                             match_by_length = TRUE,
                             format = NULL,
                             verbose = TRUE) {
    if (!is.numeric(min_coverage) || length(min_coverage) != 1L ||
        min_coverage <= 0 || min_coverage > 1) {
        stop("`min_coverage` must be a single number in (0, 1].", call. = FALSE)
    }
    if (!is.logical(match_by_length) || length(match_by_length) != 1L) {
        stop("`match_by_length` must be a single logical.", call. = FALSE)
    }
    if (!is.null(target_chroms)) {
        if (!is.character(target_chroms) || !length(target_chroms) ||
            anyNA(target_chroms)) {
            stop("`target_chroms` must be a non-empty character vector with no NAs.",
                call. = FALSE
            )
        }
    }
    if (!is.null(target_lengths)) {
        if (is.null(target_chroms)) {
            stop("`target_lengths` requires `target_chroms`.", call. = FALSE)
        }
        if (!is.numeric(target_lengths) || anyNA(target_lengths) ||
            any(target_lengths <= 0)) {
            stop("`target_lengths` must be a positive numeric vector with no NAs.",
                call. = FALSE
            )
        }
        if (length(target_lengths) != length(target_chroms)) {
            stop(sprintf(
                "`target_lengths` (%d) must align with `target_chroms` (%d).",
                length(target_lengths), length(target_chroms)
            ), call. = FALSE)
        }
    }
    if (file.exists(path)) {
        stop(sprintf(
            "Output path '%s' already exists; refusing to overwrite. Choose a fresh path.",
            path
        ), call. = FALSE)
    }
    if (length(sets)) {
        sets <- match.arg(sets,
            choices = c("genes", "rmsk", "cgi", "cytoband"),
            several.ok = TRUE
        )
    }

    res <- .resolve_genome(name, registry = registry)
    recipe <- res$recipe
    if (!is.null(chrom_naming)) {
        recipe$chrom_naming <- chrom_naming
        .validate_recipe(recipe, name)
    }
    if (verbose) {
        message(sprintf(
            "Resolved '%s' from %s -> source=%s",
            name, res$resolved_from, recipe$source
        ))
    }

    if (recipe$source != "ucsc-hub") {
        if (!is.null(target_chroms)) {
            stop(sprintf(
                "`target_chroms` is honored only for ucsc-hub sources; got source='%s'.",
                recipe$source
            ), call. = FALSE)
        }
        if (!is.null(target_lengths)) {
            stop(sprintf(
                "`target_lengths` is honored only for ucsc-hub sources; got source='%s'.",
                recipe$source
            ), call. = FALSE)
        }
        if (!isTRUE(all.equal(min_coverage, 1.0))) {
            stop(sprintf(
                "`min_coverage` is honored only for ucsc-hub sources; got source='%s'.",
                recipe$source
            ), call. = FALSE)
        }
    }

    # From here on, any failure must clean up `path` if we created it. The
    # gate is placed AFTER the file.exists(path) guard above, so a stale
    # pre-existing directory is never unlinked.
    success <- FALSE
    on.exit(
        if (!success && dir.exists(path)) unlink(path, recursive = TRUE),
        add = TRUE
    )

    prefetched_alias <- NULL
    if (recipe$source == "ucsc-hub" && length(sets)) {
        pf_workdir <- tempfile("misha_hub_preflight_")
        dir.create(pf_workdir, recursive = TRUE)
        on.exit(unlink(pf_workdir, recursive = TRUE), add = TRUE)
        prefetched_alias <- .hub_preflight_coverage(
            accession       = recipe$accession,
            target_chroms   = target_chroms,
            target_lengths  = target_lengths,
            chrom_naming    = recipe$chrom_naming,
            min_coverage    = min_coverage,
            match_by_length = match_by_length,
            workdir         = pf_workdir,
            verbose         = verbose
        )
    }

    seq_info <- .build_seq(recipe, path,
        target_chroms    = target_chroms,
        target_lengths   = target_lengths,
        format           = format,
        prefetched_alias = prefetched_alias,
        match_by_length  = match_by_length,
        verbose          = verbose
    )
    gdb.init(path, rescan = TRUE)

    if (length(sets)) {
        gdb.install_intervals(
            groot              = path,
            source             = recipe,
            sets               = sets,
            prefix             = prefix,
            gene_sets          = gene_sets,
            gtf_priority       = gtf_priority,
            overwrite          = FALSE,
            registry           = NULL,
            target_chroms      = target_chroms,
            target_lengths     = target_lengths,
            min_coverage       = min_coverage,
            match_by_length    = match_by_length,
            prefetched_alias   = prefetched_alias,
            verbose            = verbose,
            .from_build_genome = TRUE
        )
    }

    .write_genome_info(path, name, recipe, sets, files = seq_info$files_record)
    success <- TRUE
    invisible(NULL)
}


#' List resolvable genome names
#'
#' Returns a data frame describing every genome resolvable from the active
#' registry chain (see \code{\link{gdb.build_genome}} for the chain order).
#'
#' @param registry Optional path to an explicit registry YAML, overriding the
#'   resolution chain.
#' @return A data frame with columns:
#'   \itemize{
#'     \item \code{name} - registry key.
#'     \item \code{source} - backend (\code{ucsc}, \code{ncbi}, \code{s3},
#'       \code{local}, \code{manual}).
#'     \item \code{detail} - assembly / accession / path.
#'     \item \code{resolved_from} - which registry the entry came from.
#'   }
#'
#' @seealso \code{\link{gdb.build_genome}}, \code{\link{gdb.genome_info}}.
#'
#' @examples
#' gdb.list_genomes()
#'
#' @export
gdb.list_genomes <- function(registry = NULL) {
    sources <- list()
    if (!is.null(registry)) {
        sources[[length(sources) + 1]] <- list(label = sprintf("registry arg (%s)", registry), path = registry)
    }
    opt <- getOption("misha.genome_registry")
    if (!is.null(opt)) {
        sources[[length(sources) + 1]] <- list(label = sprintf("getOption (%s)", opt), path = opt)
    }
    proj <- .find_project_misha_yaml()
    if (!is.null(proj)) {
        sources[[length(sources) + 1]] <- list(label = sprintf("project (%s)", proj), path = proj)
    }
    builtin <- .builtin_registry_path()
    if (nzchar(builtin)) {
        sources[[length(sources) + 1]] <- list(label = "built-in", path = builtin)
    }

    rows <- list()
    seen <- character(0)
    for (src in sources) {
        entries <- .parse_genome_registry(src$path)
        for (nm in names(entries)) {
            if (nm %in% seen) next
            seen <- c(seen, nm)
            recipe <- tryCatch(.normalize_recipe(entries[[nm]], nm), error = function(e) NULL)
            if (is.null(recipe)) next
            detail <- recipe$assembly %||% recipe$accession %||% recipe$path %||%
                (if (!is.null(recipe$fasta)) paste(head(recipe$fasta, 1), collapse = ",") else NA_character_)
            rows[[length(rows) + 1]] <- data.frame(
                name = nm,
                source = recipe$source,
                detail = detail,
                resolved_from = src$label,
                stringsAsFactors = FALSE
            )
        }
    }
    if (!length(rows)) {
        return(data.frame(name = character(), source = character(), detail = character(), resolved_from = character(), stringsAsFactors = FALSE))
    }
    do.call(rbind, rows)
}


#' Inspect a resolved genome recipe without building
#'
#' Resolves \code{name} through the registry chain and returns the recipe (a
#' list) along with the source it was resolved from. Useful for previewing
#' what \code{\link{gdb.build_genome}} would do.
#'
#' @param name Genome name.
#' @param registry Optional path to an explicit registry YAML.
#' @return A list with components \code{recipe} (the resolved recipe) and
#'   \code{resolved_from} (the registry source).
#'
#' @seealso \code{\link{gdb.build_genome}}, \code{\link{gdb.list_genomes}}.
#'
#' @examples
#' gdb.genome_info("hg38")
#' gdb.genome_info("GCF_009806435.1")
#'
#' @export
gdb.genome_info <- function(name, registry = NULL) {
    .resolve_genome(name, registry = registry)
}


#' Pre-install UCSC's gff3ToGenePred binary
#'
#' Downloads UCSC's \code{gff3ToGenePred} static binary (~25 MB) into
#' \code{tools::R_user_dir("misha", "cache")/bin/}, verifies its SHA256, and
#' makes it executable. Used by \code{\link{gdb.build_genome}} when the
#' \code{ncbi} backend (or the \code{manual} backend with
#' \code{genes_format: gff3}) is invoked. Calling it directly is useful in CI
#' or in non-interactive scripts where the consent prompt would otherwise
#' fail.
#'
#' Override the binary location by setting environment variable
#' \code{MISHA_GFF3_TO_GENEPRED} to a binary you provide (for example, one
#' installed via \code{conda install -c bioconda ucsc-gff3togenepred}). This
#' is the recommended workaround on systems whose glibc is older than the one
#' UCSC's prebuilt binary requires.
#'
#' @param force If \code{TRUE}, skip the consent prompt and re-download even
#'   if the binary is already cached.
#' @return The cache path of the installed binary (invisibly).
#'
#' @examples
#' \dontrun{
#' gdb.install_gff3_converter()
#' Sys.setenv(MISHA_GFF3_TO_GENEPRED = "/path/to/your/gff3ToGenePred")
#' }
#'
#' @export
gdb.install_gff3_converter <- function(force = FALSE) {
    invisible(.install_gff3_converter(force = force))
}

#' Pre-install UCSC's gtfToGenePred binary
#'
#' Mirrors \code{\link{gdb.install_gff3_converter}}. Required for the
#' \code{ucsc-hub} backend's \code{genes} set (UCSC mammal hubs ship GTFs).
#'
#' Override the binary location by setting environment variable
#' \code{MISHA_GTF_TO_GENEPRED}.
#'
#' @param force If \code{TRUE}, skip consent prompt and re-download even if cached.
#' @return The cache path (invisibly).
#' @examples
#' \dontrun{
#' gdb.install_gtf_converter()
#' Sys.setenv(MISHA_GTF_TO_GENEPRED = "/path/to/your/gtfToGenePred")
#' }
#' @export
gdb.install_gtf_converter <- function(force = FALSE) {
    invisible(.install_gtf_converter(force = force))
}

#' Install interval sets onto an existing groot
#'
#' Given an existing groot and a source recipe (or registry name, or accession),
#' fetches the relevant annotation files and installs interval sets - one or
#' more of \code{genes / rmsk / cgi / cytoband}.
#'
#' Decoupled from \code{\link{gdb.build_genome}} so that:
#' \itemize{
#'   \item users with a private FASTA build can layer canonical annotations onto it;
#'   \item failed installs can be resumed without re-fetching the FASTA;
#'   \item the same groot can host annotations from multiple sources under
#'         different prefixes (e.g. \code{intervs.global.}, \code{intervs.repeats.}).
#' }
#'
#' @param groot Path to a misha groot. \code{NULL} uses the active groot.
#' @param source Either a registry name, a recipe \code{list}, or a bare
#'   \code{GC[FA]_<digits>.<digits>} accession.
#' @param sets Subset of \code{c("genes", "rmsk", "cgi", "cytoband")}.
#' @param prefix Character scalar prepended verbatim to each set name. Include
#'   the trailing dot if you want one (e.g. \code{"intervs.global."}).
#' @param gene_sets Named character vector mapping
#'   \code{c("tss", "exons", "utr3", "utr5")} to the on-disk set name. \code{NA}
#'   value skips that role.
#' @param gtf_priority Character vector ordering GTF source preference for
#'   sources that ship multiple GTFs (currently \code{ucsc-hub}). First found wins.
#' @param overwrite If \code{FALSE} (default), error on existing target sets.
#'   If \code{TRUE}, remove existing sets before saving.
#' @param registry Optional path to a registry YAML; overrides the resolution chain.
#' @param target_chroms Optional character vector to pin the canonical column
#'   to (e.g. chrom names from \code{halStats --sequenceStats} for a HAL you
#'   intend to liftover against). When \code{NULL} (default) misha uses the
#'   groot's own chrom names (i.e. picks the alias column matching whatever
#'   is currently in the database). When supplied, misha picks the alias
#'   column matching \code{target_chroms} instead, and switches detection to
#'   count-weighted coverage (bp weighting requires lengths, which target
#'   chrom lists typically don't carry).
#' @param target_lengths Optional numeric vector aligned with
#'   \code{target_chroms}. Only honored when this call originates from
#'   \code{\link{gdb.build_genome}} (the groot was just force-aligned to
#'   \code{target_chroms}); standalone calls ignore it and use the strict
#'   column gate. When honored, canonical is set to a synthetic
#'   \code{".target_chroms"} column populated by name match + unique-length
#'   pairing, so \code{chrom_aliases.tsv} writes \code{target_chroms} as
#'   canonical with all other chromAlias columns as aliases.
#' @param min_coverage Minimum fraction that must be covered by a chromAlias
#'   column for it to be picked as the canonical mapping. Default \code{1.0}
#'   (strict). On the groot side this is bp-weighted (fraction of genome
#'   basepairs covered) - a long-tail of small unmapped contigs (e.g. a 16 kb
#'   mitochondrion missing from UCSC's \code{genbank} column out of a 3 Gb
#'   genome) costs ~0.0005% rather than ~1/N. On the source-file side
#'   (asset chroms read from a GTF/GFF) the metric is the count-weighted
#'   fraction of distinct names. Unmapped contigs receive no annotations.
#' @param match_by_length If \code{TRUE} (default), complement the
#'   column-based canonical detection with a per-row length-based fill:
#'   alias rows whose chosen column is empty are paired with a groot chrom
#'   of the same length, but only when the length is unique on both sides
#'   (ambiguous lengths are skipped, never guessed). Asset translation also
#'   switches to a per-row cross-column lookup, so a GFF in any naming
#'   scheme (or mixed schemes) imports cleanly. Currently honored only by
#'   the \code{ucsc-hub} backend (which ships per-contig lengths in
#'   \code{<acc>.chrom.sizes.txt}); other backends are unaffected. Set
#'   \code{FALSE} for the stricter single-column-only behavior.
#' @param force If \code{FALSE} (default), any requested set that the source
#'   doesn't provide raises an error and aborts the call before touching the
#'   groot. If \code{TRUE}, missing sets are demoted to a single summary
#'   warning and the available sets are installed.
#' @param verbose If \code{TRUE}, prints progress.
#' @param prefetched_alias Optional pre-fetched chromAlias bundle (the
#'   return value of \code{.hub_preflight_coverage}). When supplied the
#'   ucsc-hub fetcher reuses it instead of re-downloading. Internal; users
#'   never set this directly.
#' @param .from_build_genome Internal flag; when \code{TRUE},
#'   \code{gdb.build_genome} signals that it has already rescanned the groot
#'   and we can skip the entry rescan. Users never set this directly.
#' @return Invisible \code{NULL}. Side effects: writes \code{.interv} files under
#'   \code{<groot>/tracks/}, extends \code{<groot>/chrom_aliases.tsv}, appends to
#'   \code{<groot>/genome_info.yaml}, and re-initializes the active groot.
#'
#'   The gene-derived sets (\code{tss}, \code{exons}, \code{utr3}, \code{utr5})
#'   carry a \code{name} column (transcript/RNA accession) and a \code{geneName}
#'   column (gene symbol from the source annotation; blank when the source has
#'   no symbol).
#'
#' @seealso \code{\link{gdb.build_genome}}, \code{\link{gdb.install_gtf_converter}}.
#'
#' @examples
#' \dontrun{
#' # Standalone install on an existing groot.
#' gdb.install_intervals(
#'     groot  = "/genomes/arctic_fox",
#'     source = "GCA_004023825.1",
#'     prefix = "intervs.global."
#' )
#'
#' # Layered: private FASTA groot + intervals from a UCSC hub assembly.
#' gdb.install_intervals(
#'     groot  = "/genomes/my_private",
#'     source = list(source = "ucsc-hub", accession = "GCF_009806435.1"),
#'     sets   = c("genes", "rmsk")
#' )
#' }
#' @export
gdb.install_intervals <- function(groot,
                                  source,
                                  sets = c("genes", "rmsk", "cgi", "cytoband"),
                                  prefix = "",
                                  gene_sets = c(
                                      tss = "tss", exons = "exons",
                                      utr3 = "utr3", utr5 = "utr5"
                                  ),
                                  gtf_priority = c(
                                      "ncbiRefSeq", "bestRefSeq",
                                      "ensGene", "augustus", "xenoRefGene"
                                  ),
                                  overwrite = FALSE,
                                  registry = NULL,
                                  target_chroms = NULL,
                                  target_lengths = NULL,
                                  min_coverage = 1.0,
                                  match_by_length = TRUE,
                                  force = FALSE,
                                  verbose = TRUE,
                                  prefetched_alias = NULL,
                                  .from_build_genome = FALSE) {
    sets <- match.arg(sets,
        choices = c("genes", "rmsk", "cgi", "cytoband"),
        several.ok = TRUE
    )
    if (!is.null(target_chroms) &&
        (!is.character(target_chroms) || !length(target_chroms) || anyNA(target_chroms))) {
        stop("`target_chroms` must be a non-empty character vector with no NAs.",
            call. = FALSE
        )
    }
    if (!is.null(target_lengths)) {
        if (is.null(target_chroms)) {
            stop("`target_lengths` requires `target_chroms`.", call. = FALSE)
        }
        if (!is.numeric(target_lengths) || anyNA(target_lengths) ||
            any(target_lengths <= 0)) {
            stop("`target_lengths` must be a positive numeric vector with no NAs.",
                call. = FALSE
            )
        }
        if (length(target_lengths) != length(target_chroms)) {
            stop(sprintf(
                "`target_lengths` (%d) must align with `target_chroms` (%d).",
                length(target_lengths), length(target_chroms)
            ), call. = FALSE)
        }
    }
    if (!is.numeric(min_coverage) || min_coverage <= 0 || min_coverage > 1) {
        stop("`min_coverage` must be in (0, 1].", call. = FALSE)
    }
    if (!is.logical(match_by_length) || length(match_by_length) != 1L) {
        stop("`match_by_length` must be a single logical.", call. = FALSE)
    }
    if (!is.null(groot)) {
        if (!dir.exists(groot) ||
            !file.exists(file.path(groot, "chrom_sizes.txt"))) {
            stop(sprintf(
                "'%s' is not a misha groot (no chrom_sizes.txt). ",
                groot
            ), call. = FALSE)
        }
        # gdb.build_genome already rescanned this groot at line ~882; skip the
        # redundant rescan to avoid an extra chrom_sizes/ALLGENOME read (slow
        # on NFS for many-contig assemblies).
        if (!.from_build_genome) {
            gdb.init(groot, rescan = TRUE)
        }
    } else if (!exists("GROOT", envir = .misha)) {
        stop("No active groot and no `groot` argument supplied.", call. = FALSE)
    }
    groot <- get("GROOT", envir = .misha)

    # Resolve source: list -> use as recipe; string -> registry/pattern.
    recipe <- if (is.list(source)) {
        .normalize_recipe(source, "<arg>")
        .validate_recipe(source, "<arg>")
        source
    } else if (is.character(source) && length(source) == 1L) {
        res <- .resolve_genome(source, registry = registry)
        if (verbose) {
            message(sprintf(
                "Resolved source '%s' from %s -> source=%s",
                source, res$resolved_from, res$recipe$source
            ))
        }
        res$recipe
    } else {
        stop("`source` must be a length-1 character or a recipe list.", call. = FALSE)
    }

    if (recipe$source %in% c("local", "s3")) {
        stop(sprintf("source '%s' has no fetchable intervals.", recipe$source),
            call. = FALSE
        )
    }

    workdir <- tempfile("misha_install_intervals_")
    dir.create(workdir, recursive = TRUE)
    on.exit(unlink(workdir, recursive = TRUE), add = TRUE)

    fetcher <- switch(recipe$source,
        ucsc = .ucsc_fetch_assets,
        `ucsc-hub` = function(r, s, w, v) {
            .hub_fetch_assets(r, s, w, gtf_priority,
                prefetched_alias = prefetched_alias, verbose = v
            )
        },
        ncbi = .ncbi_fetch_assets,
        manual = .manual_fetch_assets,
        stop(sprintf("No fetcher for source '%s'", recipe$source), call. = FALSE)
    )
    assets <- fetcher(recipe, sets, workdir, verbose)

    # Detect requested sets the source couldn't provide. Fetchers warn-and-skip
    # individually, which historically meant gdb.install_intervals returned
    # cleanly even when some requested sets weren't installed. Surface that as
    # an error (or a single summary warning under force=TRUE).
    missing_sets <- sets[vapply(sets, function(s) is.null(assets[[s]]), logical(1))]
    if (length(missing_sets) > 0L) {
        msg <- sprintf(
            "Requested set(s) not available from source '%s': %s",
            recipe$source, paste(missing_sets, collapse = ", ")
        )
        if (isTRUE(force)) {
            warning(msg, call. = FALSE)
        } else {
            stop(sprintf(
                "%s\nPass force = TRUE to demote this to a warning and install the available sets.",
                msg
            ), call. = FALSE)
        }
    }
    installed_sets <- setdiff(sets, missing_sets)

    # chromAlias: detect groot column and source columns; build translator closure.
    # Coverage is bp-weighted on the groot side -- a long-tail of small unmapped
    # contigs (e.g. a 16 kb mitochondrion missing from UCSC's genbank column)
    # then costs ~0.0005% instead of ~1/N, which is the meaningful figure for
    # whole-genome work.
    allg <- get("ALLGENOME", envir = .misha)[[1]]
    groot_chroms <- as.character(allg$chrom)
    groot_lengths <- as.numeric(allg$end - allg$start)
    alias_df <- assets$chrom_alias$df
    # `target_chroms` lets the caller pin the canonical column to whatever
    # external naming they're aligning to (e.g. HAL/halStats output) without
    # having to hand-pick a chromAlias column. When NULL (default) we fall
    # back to "what's currently in the groot", which is the typical case.
    detect_chroms <- target_chroms %||% groot_chroms
    detect_lengths <- if (is.null(target_chroms)) groot_lengths else NULL
    # Force-align canonical to target_chroms when gdb.build_genome ran the
    # same force-align in .build_seq_ucsc_hub (so the FASTA's headers are now
    # target_chroms). Inject a synthetic ".target_chroms" column carrying the
    # per-row assignment; downstream this becomes the canonical column,
    # which is later filled for non-target rows by the standard
    # match_by_length groot-side length-fill. Gated on .from_build_genome so
    # standalone gdb.install_intervals (where the groot wasn't built with
    # this branch) keeps strict gate semantics.
    force_align <- .from_build_genome && !is.null(alias_df) &&
        !is.null(target_chroms) && !is.null(target_lengths) &&
        isTRUE(match_by_length) && !is.null(assets$chrom_alias$row_lengths)
    if (force_align) {
        alias_df[[".target_chroms"]] <- .assign_target_chroms_per_row(
            alias_df, target_chroms, target_lengths,
            assets$chrom_alias$row_lengths
        )
        groot_col <- structure(".target_chroms",
            overlap = 1.0, bp_weighted = FALSE
        )
        if (verbose) {
            message("  Force-align canonical = target_chroms (synthetic .target_chroms column).")
        }
    } else {
        # When match_by_length is on, the pre-rescue column pick is just a
        # starting point for the length-fill / length-override / name-override
        # passes; `min_coverage` should apply to the post-rescue canonical
        # below, not to the single-column score. Hybrid HAL conventions
        # (e.g. Phylo447: assembled = UCSC names, unplaced = bare GenBank
        # accessions) routinely have no single column above 99% but the
        # rescues take final coverage to 100%. Picking at `min_coverage = 0`
        # always returns the best-scoring column.
        select_min <- if (isTRUE(match_by_length)) 0 else min_coverage
        groot_col <- if (!is.null(alias_df)) {
            .detect_alias_column(alias_df, detect_chroms,
                min_coverage = select_min, chrom_lengths = detect_lengths
            )
        } else {
            NA_character_
        }
        if (!is.null(alias_df) && is.na(groot_col)) {
            label <- if (is.null(target_chroms)) "groot" else "target_chroms"
            .coverage_gate(alias_df, detect_chroms, detect_lengths,
                min_coverage = min_coverage, label = label
            )
        }
    }
    if (!is.null(alias_df)) {
        groot_col_chr <- as.character(groot_col)
        # When match_by_length, complement the chosen column with a per-row
        # length-match fill. The chosen column may have empty cells (e.g. MT
        # row's genbank cell is blank in many UCSC hubs); length matching
        # provides a canonical for those rows by pairing unique-on-both-sides
        # lengths. The combined per-row canonical lives in a virtual
        # ".canonical" column that downstream code (chrom_aliases.tsv writer,
        # asset translator) treats as the groot column.
        canonical_col <- groot_col_chr
        if (match_by_length) {
            base_canonical <- alias_df[[groot_col_chr]]
            # Treat NA as empty for fill purposes.
            base_canonical[is.na(base_canonical)] <- ""
            row_lengths <- assets$chrom_alias$row_lengths
            n_filled <- 0L
            n_overridden <- 0L
            if (!is.null(row_lengths)) {
                filled <- .length_match_fill(
                    base_canonical, row_lengths,
                    groot_chroms, groot_lengths
                )
                # Account empty/NA after fill as still-empty.
                filled[is.na(filled)] <- ""
                n_filled <- sum(nzchar(filled) & !nzchar(base_canonical))
                # Second pass: rows where canonical is non-empty but doesn't
                # appear in the groot (different naming convention for the
                # same physical contig -- e.g. canonical="chrM" but the groot
                # uses the GenBank accession "AY172581.1") get overridden by
                # unique-on-both-sides length pairing. Rat MT is the
                # canonical example: chromAlias and report agree on "chrM"
                # but the groot follows a HAL that used GenBank accessions.
                pre_override <- filled
                filled <- .length_match_override(
                    filled, row_lengths, groot_chroms, groot_lengths
                )
                n_overridden <- sum(filled != pre_override)
                base_canonical <- filled
            }
            # Third pass: name-based fallback. For any row whose canonical
            # still doesn't appear in the groot, scan the row's other
            # columns for a value that does -- and use it. Catches rows
            # where the chosen canonical column was empty/different but
            # some other naming-scheme column (e.g. GenBank-Accn) carries
            # the exact groot name (typical: bare-accession contigs in
            # HAL-built grooots).
            pre_name <- base_canonical
            base_canonical <- .name_match_override(
                base_canonical, alias_df, groot_col_chr, groot_chroms
            )
            n_name_overridden <- sum(base_canonical != pre_name)
            alias_df[[".canonical"]] <- base_canonical
            canonical_col <- ".canonical"
            if (verbose && n_filled > 0L) {
                message(sprintf(
                    "  Length-matched %d alias row(s) the '%s' column couldn't cover.",
                    n_filled, groot_col_chr
                ))
            }
            if (verbose && n_overridden > 0L) {
                message(sprintf(
                    "  Length-overrode %d alias row(s) where canonical didn't match groot but length paired uniquely.",
                    n_overridden
                ))
            }
            if (verbose && n_name_overridden > 0L) {
                message(sprintf(
                    "  Name-overrode %d alias row(s) by picking a non-canonical column whose value matched the groot.",
                    n_name_overridden
                ))
            }
            if (verbose && match_by_length && is.null(row_lengths)) {
                message("  match_by_length=TRUE but no per-row lengths fetched; falling back to column-only.")
            }
        }
        canonical_vals <- alias_df[[canonical_col]]
        unmapped <- groot_chroms[!groot_chroms %in% canonical_vals]
        # Post-rescue gate: when match_by_length=TRUE, min_coverage applies
        # to the final canonical column (after all rescue passes), not to
        # the single-column score above. The pre-rescue gate is skipped
        # precisely so the rescues get a chance; this is where we re-check.
        if (isTRUE(match_by_length) && length(unmapped)) {
            unmapped_bp <- sum(groot_lengths[!groot_chroms %in% canonical_vals])
            final_cov <- 1 - unmapped_bp / sum(groot_lengths)
            if (final_cov < min_coverage) {
                stop(sprintf(
                    "After alias resolution (column-detect + length-fill + length-override + name-override), %d groot contigs (%s bp, %.4f%% of genome) remain unmapped; final coverage %.4f%% < min_coverage %.4f%%.\n  %s\nHint: lower `min_coverage` (e.g. min_coverage = %.2f) if you accept the unmapped fraction.",
                    length(unmapped),
                    format(unmapped_bp, big.mark = ","),
                    100 * unmapped_bp / sum(groot_lengths),
                    100 * final_cov,
                    100 * min_coverage,
                    .diagnose_unmapped_chroms(
                        unmapped, alias_df, canonical_col,
                        groot_chroms, groot_lengths,
                        assets$chrom_alias$row_lengths, groot
                    ),
                    floor(final_cov * 100) / 100
                ), call. = FALSE)
            }
        }
        if (verbose && length(unmapped)) {
            unmapped_bp <- sum(groot_lengths[!groot_chroms %in% canonical_vals])
            message(sprintf(
                "  %d groot contigs (%s bp, %.4f%% of genome) remain unmapped and will receive no annotations.\n  %s",
                length(unmapped),
                format(unmapped_bp, big.mark = ","),
                100 * unmapped_bp / sum(groot_lengths),
                .diagnose_unmapped_chroms(
                    unmapped,
                    alias_df,
                    canonical_col,
                    groot_chroms,
                    groot_lengths,
                    assets$chrom_alias$row_lengths,
                    groot
                )
            ))
        }
        .merge_chrom_aliases_tsv(groot, alias_df, canonical_col)
        gdb.init(groot, rescan = TRUE) # reload CHROM_ALIAS
    }

    # When match_by_length is on, build the cross-column reverse index once
    # so we don't rebuild it per asset.
    rev_idx <- if (!is.null(alias_df) && match_by_length) {
        .build_alias_rev_index(alias_df, canonical_col)
    } else {
        NULL
    }

    # Helper: per-asset translator.
    # When match_by_length is on, translate per-row across all alias columns
    # so a GFF in any naming scheme (or mixed schemes) imports cleanly. Else
    # fall back to single-column detection at `min_coverage`.
    make_translator <- function(asset_chroms, asset_label) {
        if (is.null(alias_df)) {
            return(NULL)
        }
        if (match_by_length) {
            return(function(rows, chrom_col) {
                .translate_chroms_per_row(rows, chrom_col, rev_idx)
            })
        }
        src_col <- .detect_alias_column(alias_df, unique(asset_chroms),
            min_coverage = min_coverage
        )
        if (is.na(src_col)) {
            scores <- attr(src_col, "scores")
            stop(sprintf(
                "chromAlias has no column with %.0f%% coverage of distinct chroms in %s.\nPer-column overlap counts: %s\nFirst 5 unmapped chroms: %s",
                100 * min_coverage,
                asset_label,
                paste(sprintf("%s=%d/%d", names(scores), scores, length(unique(asset_chroms))),
                    collapse = ", "
                ),
                paste(
                    utils::head(setdiff(
                        unique(asset_chroms),
                        unlist(alias_df, use.names = FALSE)
                    ), 5L),
                    collapse = ", "
                )
            ), call. = FALSE)
        }
        # Strip attributes for stable indexing into alias_df.
        src_col_chr <- as.character(src_col)
        function(rows, chrom_col) {
            .translate_chroms(rows, chrom_col, alias_df, src_col_chr, canonical_col)
        }
    }

    # Genes.
    if ("genes" %in% sets && !is.null(assets$genes)) {
        # Sample first ~100 chroms from genePred/GTF/GFF for translator detection.
        chroms <- .sample_chroms_from_file(assets$genes$file, assets$genes$format)
        translator <- make_translator(chroms, "genes file")
        asset <- c(assets$genes, list(translate = translator))
        .install_genes_set(asset,
            prefix = prefix, gene_sets = gene_sets,
            overwrite = overwrite, verbose = verbose
        )
    }
    # rmsk / cgi / cytoband: identical pipeline (parse -> translate -> install).
    # Genes is structurally different (streamed, translator passed as asset).
    simple_set_specs <- list(
        rmsk = list(
            parser = function(a) {
                if (a$format == "rmsk-out") .parse_rm_out(a$file, verbose) else .parse_ucsc_rmsk(a$file)
            },
            installer = .install_rmsk_set
        ),
        cgi = list(
            parser = function(a) .parse_ucsc_cpg_island(a$file),
            installer = .install_cgi_set
        ),
        cytoband = list(
            parser = function(a) .parse_ucsc_cytoband(a$file),
            installer = .install_cytoband_set
        )
    )
    for (key in names(simple_set_specs)) {
        if (!(key %in% sets) || is.null(assets[[key]])) next
        spec <- simple_set_specs[[key]]
        df <- spec$parser(assets[[key]])
        if (!is.null(alias_df)) {
            # match_by_length translates per-row via rev_idx and ignores
            # asset_chroms; skip the unique() (rmsk df can have ~7M rows).
            asset_chroms <- if (match_by_length) NULL else unique(df$chrom)
            translator <- make_translator(asset_chroms, key)
            df <- translator(df, "chrom")
            # Drop rows whose chrom didn't translate to a groot name (NA
            # from rev_idx misses, "" from rows whose canonical was unset
            # by the 3-pass resolution). Mirror the genes path: the C++
            # importers (and downstream gintervals.save) treat empty
            # chrom as a hard error.
            keep <- !is.na(df$chrom) & nzchar(df$chrom)
            if (any(!keep)) {
                if (verbose) {
                    message(sprintf(
                        "  %s: %d rows dropped (chrom didn't translate to a groot contig).",
                        key, sum(!keep)
                    ))
                }
                df <- df[keep, , drop = FALSE]
            }
        }
        spec$installer(df, prefix = prefix, overwrite = overwrite, verbose = verbose)
    }

    # Provenance: append to genome_info.yaml. Record only the sets that
    # actually got installed (under force=TRUE some requested sets may have
    # been missing and skipped).
    .append_tracks_to_genome_info(groot, recipe, installed_sets, prefix)

    # Final reload + summary.
    gdb.init(groot, rescan = TRUE)
    if (verbose) .install_intervals_summary(groot, recipe, installed_sets, prefix)
    invisible(NULL)
}

# Sample distinct chroms from a genePred/GTF/GFF for translator detection.
# Reads up to 50,000 lines.
.sample_chroms_from_file <- function(file, format) {
    con <- if (grepl("\\.gz$", file)) gzfile(file, "rt") else file(file, "rt")
    on.exit(close(con), add = TRUE)
    chroms <- character(0)
    chunk_size <- 50000L
    repeat {
        lines <- readLines(con, n = chunk_size, warn = FALSE)
        if (!length(lines)) break
        # genePred: chrom in column 2; GTF/GFF: column 1.
        col_idx <- if (format == "genepred") 2L else 1L
        # Skip comment lines.
        lines <- lines[!startsWith(lines, "#") & nzchar(lines)]
        if (!length(lines)) next
        f <- strsplit(lines, "\t", fixed = TRUE)
        chroms <- unique(c(
            chroms,
            vapply(
                f, function(x) if (length(x) >= col_idx) x[[col_idx]] else NA_character_,
                character(1)
            )
        ))
        if (length(chroms) > 5000L) break # plenty for detection
    }
    chroms[!is.na(chroms) & nzchar(chroms)]
}

.install_intervals_summary <- function(groot, recipe, sets, prefix) {
    message("\ngdb.install_intervals: completed")
    message(sprintf("  groot:   %s", groot))
    message(sprintf(
        "  source:  %s%s", recipe$source,
        if (!is.null(recipe$accession)) {
            sprintf(" (%s)", recipe$accession)
        } else if (!is.null(recipe$assembly)) {
            sprintf(" (%s)", recipe$assembly)
        } else {
            ""
        }
    ))
    message(sprintf("  prefix:  %s", if (nzchar(prefix)) sprintf("\"%s\"", prefix) else "(none)"))
    message(sprintf("  sets:    %s", paste(sets, collapse = ", ")))
}

.append_tracks_to_genome_info <- function(groot, recipe, sets, prefix) {
    info_path <- file.path(groot, "genome_info.yaml")
    info <- if (file.exists(info_path)) yaml::read_yaml(info_path) else list()
    if (is.null(info$tracks)) info$tracks <- list()
    ts <- format(Sys.time(), "%Y-%m-%dT%H:%M:%SZ", tz = "UTC")
    for (s in sets) {
        info$tracks[[length(info$tracks) + 1L]] <- list(
            set          = paste0(prefix, s),
            source       = recipe$source,
            installed_at = ts
        )
    }
    .gwith_umask(yaml::write_yaml(info, info_path))
    invisible(NULL)
}