fix: transparently handle regular-gzip (bcftools -O z) VCFs (#32)

DESCRIPTION

@@ -3,7 +3,7 @@ Type: Package
 Title: Annotation of Genetic Variants
 Description: Annotate variants, compute amino acid coding changes,
         predict coding outcomes.
-Version: 1.59.0
+Version: 1.59.2
 Authors@R: c(
         person("Valerie", "Oberchain", role="aut"),
         person("Martin", "Morgan", role="aut"),

R/methods-readVcf.R

@@ -55,6 +55,42 @@ setMethod(readVcf, c(file="character", param="missing"),
              row.names=row.names, ...)
 })
 
+.is_bgzf <- function(path)
+{
+    ## BGZF blocks carry a specific extra-field signature in the gzip header:
+    ##   byte  3 (FLG)  has the FEXTRA bit (0x04) set
+    ##   bytes 13-14 (first extra sub-field ID) equal 'B','C' (0x42, 0x43)
+    ## Regular gzip files are missing that sub-field even when FLG has FEXTRA.
+    bytes <- tryCatch(readBin(path, "raw", n = 18L), error = function(e) raw())
+    if (length(bytes) < 18L) return(FALSE)
+    ## magic: 0x1f 0x8b
+    if (!identical(bytes[1:2], as.raw(c(0x1f, 0x8b)))) return(FALSE)
+    ## FLG byte (4th) must have FEXTRA (0x04)
+    if (bitwAnd(as.integer(bytes[4L]), 0x04L) == 0L) return(FALSE)
+    ## SI1='B'(0x42), SI2='C'(0x43) at bytes 13-14 (1-based)
+    identical(bytes[13:14], as.raw(c(0x42, 0x43)))
+}
+
+.ensure_bgzf <- function(path)
+{
+    ## If the file is regular gzip (not BGZF), re-compress to a BGZF tempfile.
+    ## Returns the (possibly new) path; the caller is responsible for cleanup
+    ## via .bgzf_tempfiles if needed.
+    if (!grepl("\\.gz$", path, ignore.case = TRUE)) return(path)
+    if (.is_bgzf(path)) return(path)
+    message("Note: '", basename(path), "' appears to be regular gzip ",
+            "(not BGZF). Re-compressing to BGZF for htslib compatibility.")
+    tmp_vcf <- tempfile(fileext = ".vcf")
+    tmp_bgz <- paste0(tmp_vcf, ".bgz")
+    con_in  <- gzcon(file(path, "rb"))
+    txt <- readLines(con_in)
+    close(con_in)
+    writeLines(txt, tmp_vcf)
+    bgzip(tmp_vcf, tmp_bgz, overwrite = TRUE)
+    unlink(tmp_vcf)
+    tmp_bgz
+}
+
 .checkFile <- function(x)
 {
     if (1L != length(x)) 
@@ -63,6 +99,10 @@ setMethod(readVcf, c(file="character", param="missing"),
     if (grepl("\\.tbi$", x))
         return(TabixFile(sub("\\.tbi", "", x)))
 
+    ## Transparently handle regular-gzip files (e.g. from bcftools -O z)
+    ## by re-compressing to BGZF, which htslib/scanBcfHeader requires.
+    x <- .ensure_bgzf(x)
+
     ## Attempt to create TabixFile
     tryCatch(x <- TabixFile(x), error=function(e) return(x))
 

R/methods-scanVcfHeader.R

@@ -1,15 +1,155 @@
+### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+### Quote-aware VCF header line parser
+###
+### scanBcfHeader() (htslib) splits structured header fields on '=' without
+### respecting double-quoted values, so a Description or CommandLineOptions
+### containing '=' gets silently truncated.  We re-parse the raw header text
+### here and patch each DataFrame back to the correct values.
+###
+
+## Parse a single structured header body  "<ID=foo,Description="a=b",..."
+## into a named character vector, correctly handling "=" inside double quotes.
+.parseVcfHeaderBody <- function(body) {
+    ## Strip surrounding < >
+    body <- sub("^<", "", sub(">$", "", body))
+
+    ## Walk character by character, splitting on ',' outside quotes,
+    ## then splitting each token on the FIRST '=' outside quotes.
+    chars  <- strsplit(body, "")[[1]]
+    n      <- length(chars)
+    in_q   <- FALSE
+    tokens <- character(0)
+    start  <- 1L
+
+    for (i in seq_len(n)) {
+        ch <- chars[i]
+        if (ch == '"') {
+            in_q <- !in_q
+        } else if (ch == ',' && !in_q) {
+            tokens <- c(tokens, paste(chars[start:(i - 1L)], collapse = ""))
+            start  <- i + 1L
+        }
+    }
+    tokens <- c(tokens, paste(chars[start:n], collapse = ""))
+
+    ## For each token split on the first '=' that is NOT inside quotes
+    result <- character(0)
+    for (tok in tokens) {
+        tc   <- strsplit(tok, "")[[1]]
+        tn   <- length(tc)
+        tq   <- FALSE
+        split_at <- NA_integer_
+        for (j in seq_len(tn)) {
+            if (tc[j] == '"') {
+                tq <- !tq
+            } else if (tc[j] == '=' && !tq) {
+                split_at <- j
+                break
+            }
+        }
+        if (!is.na(split_at)) {
+            key <- paste(tc[seq_len(split_at - 1L)], collapse = "")
+            val <- paste(tc[(split_at + 1L):tn],     collapse = "")
+            ## strip surrounding quotes from value
+            val <- sub('^"(.*)"$', "\\1", val)
+            result[key] <- val
+        }
+        ## tokens without '=' (e.g. bare flags) are skipped — consistent with
+        ## what scanBcfHeader produces
+    }
+    result
+}
+
+## Parse all ##TYPE=<...> meta-lines from raw header text.
+## Returns a named list of data.frames (one per TYPE), each row being one
+## record, columns being the union of keys seen.  Row names are the ID field.
+.parseRawVcfHeader <- function(raw_lines) {
+    ## only structured lines:  ##KEY=<...>
+    struct <- grep("^##[A-Za-z0-9_]+=<", raw_lines, value = TRUE)
+    if (!length(struct))
+        return(list())
+
+    ## extract TYPE and BODY
+    type  <- sub("^##([A-Za-z0-9_]+)=<.*$", "\\1", struct)
+    body  <- sub("^##[A-Za-z0-9_]+=(<.*>)$", "\\1", struct)
+
+    ## parse each line into a named character vector
+    parsed <- lapply(body, .parseVcfHeaderBody)
+
+    ## group by TYPE
+    types_unique <- unique(type)
+    result <- vector("list", length(types_unique))
+    names(result) <- types_unique
+
+    for (tp in types_unique) {
+        rows <- parsed[type == tp]
+        ## union of all keys
+        all_keys <- unique(unlist(lapply(rows, names)))
+        mat <- matrix(NA_character_, nrow = length(rows), ncol = length(all_keys),
+                      dimnames = list(NULL, all_keys))
+        for (i in seq_along(rows)) {
+            k <- names(rows[[i]])
+            mat[i, k] <- rows[[i]][k]
+        }
+        df <- as.data.frame(mat, stringsAsFactors = FALSE)
+        rownames(df) <- if ("ID" %in% colnames(df)) df$ID else seq_len(nrow(df))
+        result[[tp]] <- df
+    }
+    result
+}
+
+## Patch the DataFrames produced by scanBcfHeader with correctly-parsed values
+## from the raw header lines.
+.patchVcfHeader <- function(hdr_list, raw_lines) {
+    parsed <- .parseRawVcfHeader(raw_lines)
+    if (!length(parsed))
+        return(hdr_list)
+
+    for (tp in names(parsed)) {
+        ref_df  <- parsed[[tp]]           # correctly-parsed data.frame
+        curr_df <- hdr_list[[tp]]         # possibly-truncated DataFrame (or NULL)
+
+        if (is.null(curr_df) || !is(curr_df, "DataFrame"))
+            next
+
+        ## For each column that exists in both, overwrite with re-parsed values
+        ## matched by row name (ID).
+        common_ids <- intersect(rownames(curr_df), rownames(ref_df))
+        if (!length(common_ids))
+            next
+
+        for (col in intersect(colnames(curr_df), colnames(ref_df))) {
+            curr_df[common_ids, col] <- ref_df[common_ids, col]
+        }
+        hdr_list[[tp]] <- curr_df
+    }
+    hdr_list
+}
+
+### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+### scanVcfHeader methods
+###
+
 setMethod(scanVcfHeader, "missing",
     function(file, ...)
 {
     VCFHeader()
 })
 
 setMethod(scanVcfHeader, "character",
-    function(file, ...) 
+    function(file, ...)
 {
     if (length(file)) {
-        hdr <- scanBcfHeader(file[[1]], ...)[[1]]
-        VCFHeader(hdr$Reference, hdr$Sample, hdr$Header)
+        f1       <- .ensure_bgzf(file[[1]])
+        hdr      <- scanBcfHeader(f1, ...)[[1]]
+        ## Read raw header lines; use the BGZF path for gzip, plain path otherwise
+        raw_lines <- if (grepl("\\.gz$", f1, ignore.case = TRUE))
+                         readLines(gzcon(file(f1, "rb")))
+                     else
+                         readLines(file[[1]])
+        raw_lines <- raw_lines[startsWith(raw_lines, "##")]
+        patched  <- .patchVcfHeader(hdr$Header, raw_lines)
+        VCFHeader(hdr$Reference, hdr$Sample, patched)
     } else {
         VCFHeader()
     }
@@ -18,5 +158,14 @@ setMethod(scanVcfHeader, "character",
 setMethod(scanVcfHeader, "TabixFile",
     function(file, ...)
 {
-    scanVcfHeader(path(file), ...)
+    if (isOpen(file)) {
+        ## already open: fall back to path-based method which reads raw lines
+        scanVcfHeader(path(file), ...)
+    } else {
+        hdr       <- scanBcfHeader(path(file), ...)[[1]]
+        raw_lines <- headerTabix(file)$header
+        raw_lines <- raw_lines[startsWith(raw_lines, "##")]
+        patched   <- .patchVcfHeader(hdr$Header, raw_lines)
+        VCFHeader(hdr$Reference, hdr$Sample, patched)
+    }
 })

inst/unitTests/test_readVcf-methods.R

@@ -258,3 +258,33 @@ test_buffer_realloc <- function()
     target <- ".GGGGGGGGG"
     checkIdentical(target, alt(vcf)[[1]])
 }
+
+test_regular_gzip_readVcf <- function()
+{
+    ## Issue #32: regular gzip (bcftools -O z) should be transparently
+    ## re-compressed to BGZF so scanBcfHeader does not fail.
+    fl_bgzf <- system.file("extdata", "chr7-sub.vcf.gz",
+                           package = "VariantAnnotation")
+
+    ## Build a regular gzip copy in a tempfile
+    tmp_vcf <- tempfile(fileext = ".vcf")
+    tmp_gz  <- paste0(tmp_vcf, ".gz")
+    con_in  <- gzcon(file(fl_bgzf, "rb"))
+    txt     <- readLines(con_in)
+    close(con_in)
+    writeLines(txt, tmp_vcf)
+    con_out <- gzfile(tmp_gz, "wb")
+    writeLines(txt, con_out)
+    close(con_out)
+    on.exit({ unlink(tmp_vcf); unlink(tmp_gz) }, add = TRUE)
+
+    ## .is_bgzf should return FALSE for the regular gzip file
+    checkTrue(!VariantAnnotation:::.is_bgzf(tmp_gz))
+    ## .is_bgzf should return TRUE for the original BGZF file
+    checkTrue(VariantAnnotation:::.is_bgzf(fl_bgzf))
+
+    ## readVcf should succeed (transparently re-compresses to BGZF)
+    suppressMessages(vcf <- readVcf(tmp_gz, genome = "hg19"))
+    checkTrue(is(vcf, "CollapsedVCF"))
+    checkTrue(nrow(vcf) > 0L)
+}

man/PolyPhenDb-class.Rd

@@ -22,32 +22,27 @@
 \section{Methods}{
   In the code below, \code{x} is a \code{PolyPhenDb} object.
   \describe{
-    \item{}{
-      \code{metadata(x)}:
+    \item{\code{metadata(x)}}{
       Returns \code{x}'s metadata in a data frame.
     }
-    \item{}{
-      \code{columns(x)}:
+    \item{\code{columns(x)}}{
       Returns the names of the \code{columns} that can be used to subset the
       data columns. For column descriptions see \code{?PolyPhenDbColumns}.
     }
-    \item{}{
-      \code{keys(x)}:
+    \item{\code{keys(x)}}{
       Returns the names of the \code{keys} that can be used to subset the
       data rows. The \code{keys} values are the rsid's.
     }
-    \item{}{
-      \code{select(x, keys = NULL, columns = NULL, ...)}:
-      Returns a subset of data defined by the character vectors \code{keys} 
+    \item{\code{select(x, keys = NULL, columns = NULL, ...)}}{
+      Returns a subset of data defined by the character vectors \code{keys}
       and \code{columns}. If no \code{keys} are supplied, all rows are
       returned. If no \code{columns} are supplied, all columns
       are returned. See \code{?PolyPhenDbColumns} for column descriptions.
-    } 
-    \item{}{
-      \code{duplicateRSID(x)}:
-      Returns a named list of duplicate rsid groups. The names are the 
-      \code{keys}, the list elements are the rsid's that have been 
-      reported as having identical chromosome position and alleles and 
+    }
+    \item{\code{duplicateRSID(x)}}{
+      Returns a named list of duplicate rsid groups. The names are the
+      \code{keys}, the list elements are the rsid's that have been
+      reported as having identical chromosome position and alleles and
       therefore translating into the same amino acid residue substitution.
     }
   }

man/SIFTDb-class.Rd

@@ -21,23 +21,19 @@
 \section{Methods}{
   In the code below, \code{x} is a \code{SIFTDb} object.
   \describe{
-    \item{}{
-      \code{metadata(x)}:
+    \item{\code{metadata(x)}}{
       Returns \code{x}'s metadata in a data frame.
     }
-    \item{}{
-      \code{columns(x)}:
+    \item{\code{columns(x)}}{
       Returns the names of the \code{columns} that can be used to subset the
       data columns.
     }
-    \item{}{
-      \code{keys(x)}:
+    \item{\code{keys(x)}}{
       Returns the names of the \code{keys} that can be used to subset the
       data rows. The \code{keys} values are the rsid's.
     }
-    \item{}{
-      \code{select(x, keys = NULL, columns = NULL, ...)}:
-      Returns a subset of data defined by the character vectors \code{keys} 
+    \item{\code{select(x, keys = NULL, columns = NULL, ...)}}{
+      Returns a subset of data defined by the character vectors \code{keys}
       and \code{columns}. If no \code{keys} are supplied, all rows are
       returned. If no \code{columns} are supplied, all columns
       are returned. For column descriptions see \code{?SIFTDbColumns}.