diff --git a/DESCRIPTION b/DESCRIPTION index 20697f2..3dab1d8 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -10,7 +10,7 @@ biocViews: DataImport, Sequencing, Coverage, Alignment, QualityControl URL: https://bioconductor.org/packages/Rsamtools Video: https://www.youtube.com/watch?v=Rfon-DQYbWA&list=UUqaMSQd_h-2EDGsU6WDiX0Q BugReports: https://github.com/Bioconductor/Rsamtools/issues -Version: 2.29.0 +Version: 2.30.0 License: Artistic-2.0 | file LICENSE Encoding: UTF-8 Authors@R: c( @@ -18,6 +18,7 @@ Authors@R: c( person("Hervé", "Pagès", role = "aut"), person("Valerie", "Obenchain", role = "aut"), person("Nathaniel", "Hayden", role = "aut"), + person("Benjamin", "Schuster-Böckler", role = "aut", comment = "Added CRAM support"), person("Busayo", "Samuel", role = "ctb", comment = "Converted Rsamtools vignette from Sweave to RMarkdown / HTML."), person("Bioconductor Package Maintainer", diff --git a/NEWS b/NEWS index 8e25f63..bf47108 100644 --- a/NEWS +++ b/NEWS @@ -1,3 +1,13 @@ +CHANGES IN VERSION 2.30 +----------------------- + +NEW FEATURES + + o (v 2.30.0) First attempt at CRAM support. To open CRAM files, a new 'reference' + parameter was added to the 'BamFile()' constructor. Other parts of the code + should work the same as with BAM files. (See + https://github.com/Bioconductor/Rsamtools/issues/56. ; bsb) + CHANGES IN VERSION 2.16 ----------------------- @@ -72,7 +82,7 @@ BUG FIXES o (v 1.33.1) Do not try to grow NULL (not-yet-encountered) tags (https://support.bioconductor.org/p/110609/ ; Robert Bradley) - o (v 1.33.5) Check for corrupt index + o (v 1.33.5) Check for corrupt index (https://github.com/Bioconductor/Rsamtools/issues/3 ; kjohnsen) CHANGES IN VERSION 1.31 @@ -299,7 +309,7 @@ NEW FEATURES SIGNIFICANT USER-VISIBLE CHANGES o rename: - readBamGappedAlignments() -> readGAlignmentsFromBam() + readBamGappedAlignments() -> readGAlignmentsFromBam() readBamGappedReads() -> readGappedReadsFromBam() readBamGappedAlignmentPairs() -> readGAlignmentPairsFromBam() readBamGAlignmentsList() -> readGAlignmentsListFromBam() @@ -334,7 +344,7 @@ NEW FEATURES Bam files only. o Add readBamGAlignmentsList function for reading qname-sorted - Bam files into a GAlignmentsList object. + Bam files into a GAlignmentsList object. USER-VISIBLE CHANGES @@ -342,7 +352,7 @@ USER-VISIBLE CHANGES vectors. o 'yieldSize' argument in BamFile represents the number of - unique qnames when 'obeyQname=TRUE'. + unique qnames when 'obeyQname=TRUE'. BUG FIXES @@ -444,11 +454,11 @@ NEW FEATURES o Provide a zlib for Windows, as R does not currently do this o BamFileList, BcfFileList, TabixFileList, FaFileList clases - extend IRanges::SimpleList, for managings lists of file references + extend IRanges::SimpleList, for managings lists of file references o razfFa creates random access compressed fasta files. - o count and scanBam support input of larger numbers of records; + o count and scanBam support input of larger numbers of records; countBam nucleotide count is now numeric() and subject to rounding error when large. diff --git a/R/AllClasses.R b/R/AllClasses.R index 3dee628..97dbf95 100644 --- a/R/AllClasses.R +++ b/R/AllClasses.R @@ -58,7 +58,8 @@ setClass("ApplyPileupsParam", .BamFile <- setRefClass("BamFile", contains="RsamtoolsFile", fields=list(obeyQname="logical", asMates="logical", - qnamePrefixEnd="character", qnameSuffixStart="character")) + qnamePrefixEnd="character", qnameSuffixStart="character", + reference="character")) .BcfFile <- setRefClass("BcfFile", contains="RsamtoolsFile", fields=list(mode="character")) diff --git a/R/AllGenerics.R b/R/AllGenerics.R index 25296ef..5a3a865 100644 --- a/R/AllGenerics.R +++ b/R/AllGenerics.R @@ -168,6 +168,9 @@ setGeneric("yieldSize", setGeneric("yieldSize<-", function(object, ..., value) standardGeneric("yieldSize<-")) +setGeneric("referenceFile", + function(object, ...) standardGeneric("referenceFile")) + setGeneric("obeyQname", function(object, ...) standardGeneric("obeyQname")) diff --git a/R/methods-BamFile.R b/R/methods-BamFile.R index 449c7bb..f065834 100644 --- a/R/methods-BamFile.R +++ b/R/methods-BamFile.R @@ -50,14 +50,19 @@ setMethod(isIncomplete, "BamFile", index <- do_append(index, files, ".BAI") index <- do_sub(index, files, ".bam$", ".bai") index <- do_sub(index, files, ".BAM$", ".BAI") + index <- do_append(index, files, ".crai") + index <- do_append(index, files, ".CRAI") + index <- do_sub(index, files, "\\.cram$", ".crai") + index <- do_sub(index, files, "\\.CRAM$", ".CRAI") index } BamFile <- - function(file, index=file, ..., yieldSize=NA_integer_, - obeyQname=FALSE, asMates=FALSE, - qnamePrefixEnd=NA, qnameSuffixStart=NA) + function(file, index=file, ..., yieldSize=NA_integer_, + obeyQname=FALSE, asMates=FALSE, + qnamePrefixEnd=NA, qnameSuffixStart=NA, + reference=character()) { if (missing(file) || !isSingleString(file)) stop("'file' must be character(1) and not NA") @@ -78,12 +83,17 @@ BamFile <- stop(paste(strwrap(txt), collapse="\n ")) } index <- .normalizePath(index) + if (length(reference) && nzchar(reference[[1L]])) + reference <- .normalizePath(reference[[1L]]) + else + reference <- character() qnamePrefixEnd <- .check_qname_arg(qnamePrefixEnd, "qnamePrefixEnd") qnameSuffixStart <- .check_qname_arg(qnameSuffixStart, "qnameSuffixStart") .RsamtoolsFile(.BamFile, path=file, index=index, yieldSize=yieldSize, - obeyQname=obeyQname, asMates=asMates, - qnamePrefixEnd=qnamePrefixEnd, - qnameSuffixStart=qnameSuffixStart, ...) + obeyQname=obeyQname, asMates=asMates, + qnamePrefixEnd=qnamePrefixEnd, + qnameSuffixStart=qnameSuffixStart, + reference=reference, ...) } open.BamFile <- @@ -91,8 +101,17 @@ open.BamFile <- { tryCatch({ .io_check_exists(path(con)) - index <- sub("\\.bai$", "", index(con, asNA=FALSE)) - con$.extptr <- .Call(.bamfile_open, path(con), index, "rb") + fpath <- path(con) + index <- index(con, asNA=FALSE) + ## For BAM files, strip .bai so the C code can probe for the index. + ## For CRAM files, pass the full .crai path; sam_index_load2 opens it + ## directly and handles both file.cram.crai and file.crai conventions. + if (!grepl("\\.cram$", fpath, ignore.case=TRUE)) + index <- sub("\\.bai$", "", index) + con$.extptr <- .Call(.bamfile_open, fpath, index, "rb") + ref <- con$reference + if (length(ref) && nzchar(ref)) + .Call(.bamfile_set_ref, con$.extptr, ref) }, error=function(err) { stop("failed to open BamFile: ", conditionMessage(err)) }) @@ -128,6 +147,13 @@ setMethod(seqinfo, "BamFile", Seqinfo(names(h), unname(h)) }) +setMethod(referenceFile, "BamFile", + function(object, ...) +{ + ref <- object$reference + if (length(ref) && nzchar(ref)) ref else NA_character_ +}) + setMethod(obeyQname, "BamFile", function(object, ...) { @@ -451,4 +477,5 @@ setMethod(show, "BamFile", function(object) { cat("asMates:", asMates(object), "\n") cat("qnamePrefixEnd:", qnamePrefixEnd(object), "\n") cat("qnameSuffixStart:", qnameSuffixStart(object), "\n") + cat("reference:", referenceFile(object), "\n") }) diff --git a/inst/extdata/cram_ref.fa.gz b/inst/extdata/cram_ref.fa.gz new file mode 100644 index 0000000..ae47aa4 Binary files /dev/null and b/inst/extdata/cram_ref.fa.gz differ diff --git a/inst/extdata/cram_ref.fa.gz.fai b/inst/extdata/cram_ref.fa.gz.fai new file mode 100644 index 0000000..d8797c3 --- /dev/null +++ b/inst/extdata/cram_ref.fa.gz.fai @@ -0,0 +1,3 @@ +seq1 2000 6 60 61 +seq2 1800 2046 60 61 +seq3 2200 3882 60 61 diff --git a/inst/extdata/cram_ref.fa.gz.gzi b/inst/extdata/cram_ref.fa.gz.gzi new file mode 100644 index 0000000..1b1cb4d Binary files /dev/null and b/inst/extdata/cram_ref.fa.gz.gzi differ diff --git a/inst/extdata/test.cram b/inst/extdata/test.cram new file mode 100644 index 0000000..c44bdc9 Binary files /dev/null and b/inst/extdata/test.cram differ diff --git a/inst/extdata/test.cram.crai b/inst/extdata/test.cram.crai new file mode 100644 index 0000000..f1e931f Binary files /dev/null and b/inst/extdata/test.cram.crai differ diff --git a/inst/unitTests/test_BamFile.R b/inst/unitTests/test_BamFile.R index 71d5283..af022d8 100644 --- a/inst/unitTests/test_BamFile.R +++ b/inst/unitTests/test_BamFile.R @@ -22,6 +22,21 @@ test_BamFile_guessIndex <- function() checkIdentical(character(), .BamFile_guessIndex(character())) checkIdentical(character(), .BamFile_guessIndex()) + + ## CRAM index guessing (.crai) + cram1 <- tempfile(fileext = ".cram") + crai1 <- paste0(cram1, ".crai") # foo.cram.crai + file.create(crai1) + + cram2 <- tempfile(fileext = ".cram") + crai2 <- sub("\\.cram$", ".crai", cram2) # foo.crai + file.create(crai2) + + cram3 <- tempfile(fileext = ".cram") # no index + + cram_fls <- c(cram1, cram2, cram3) + cram_target <- c(crai1, crai2, NA_character_) + checkIdentical(cram_target, .BamFile_guessIndex(cram_fls)) } test_BamFile_openclose <- function() diff --git a/inst/unitTests/test_CramFile.R b/inst/unitTests/test_CramFile.R new file mode 100644 index 0000000..7b7f106 --- /dev/null +++ b/inst/unitTests/test_CramFile.R @@ -0,0 +1,276 @@ +## Tests for CRAM file reading via BamFile. +## +## Test data: inst/extdata/test.cram (3 sequences, 400 150bp PE reads each) +## inst/extdata/cram_ref.fa.gz (bgzipped FASTA reference with .fai/.gzi) +## +## CRAM limitations (intentionally not tested here): +## - yieldSize-based chunked sequential reading (requires bgzf-style seeking) +## - asMates=TRUE (requires bgzf-based mate-pairing internals) +## - idxstatsBam (CRAI index type not compatible with hts_idx_get_n) + +.cram_fl <- system.file("extdata", "test.cram", package="Rsamtools") +.cram_ref <- system.file("extdata", "cram_ref.fa.gz", package="Rsamtools") +.cram_idx <- system.file("extdata", "test.cram.crai", package="Rsamtools") + +## Helper: fresh BamFile for each call so state never bleeds between tests +.make_cram_bf <- function() + BamFile(.cram_fl, reference=.cram_ref) + +## --------------------------------------------------------------------------- +## BamFile construction and index guessing + +test_CramFile_guessIndex <- function() +{ + .BamFile_guessIndex <- Rsamtools:::.BamFile_guessIndex + + ## .cram.crai (appended) convention + cram1 <- tempfile(fileext=".cram") + crai1 <- paste0(cram1, ".crai") + file.create(crai1) + + ## .crai (substituted) convention + cram2 <- tempfile(fileext=".cram") + crai2 <- sub("\\.cram$", ".crai", cram2) + file.create(crai2) + + ## no index exists + cram3 <- tempfile(fileext=".cram") + + fls <- c(cram1, cram2, cram3) + target <- c(crai1, crai2, NA_character_) + checkIdentical(target, .BamFile_guessIndex(fls)) + checkIdentical(character(), .BamFile_guessIndex(character())) + checkIdentical(character(), .BamFile_guessIndex()) + + ## BAM index guessing still works alongside CRAM + bam1 <- tempfile(fileext=".bam") + bai1 <- paste0(bam1, ".bai") + file.create(bai1) + checkIdentical(bai1, .BamFile_guessIndex(bam1)) +} + +test_CramFile_referenceFile <- function() +{ + ## No reference → NA + bf_noref <- BamFile(.cram_fl) + checkIdentical(NA_character_, referenceFile(bf_noref)) + + ## With reference → normalised path stored + bf <- .make_cram_bf() + checkIdentical(.cram_ref, referenceFile(bf)) + + ## Index auto-guessed + checkIdentical(.cram_idx, index(bf)) +} + +test_CramFile_openclose <- function() +{ + .normalizePath <- Rsamtools:::.normalizePath + + bf <- BamFile(.cram_fl, reference=.cram_ref) + checkIdentical(FALSE, isOpen(bf)) + + open(bf) + checkIdentical(TRUE, isOpen(bf)) + checkIdentical(.normalizePath(.cram_fl), path(bf)) + checkIdentical(.normalizePath(.cram_idx), index(bf)) + + close(bf) + checkIdentical(FALSE, isOpen(bf)) + checkException(close(bf), silent=TRUE) + + ## re-open a closed BamFile + open(bf) + checkIdentical(TRUE, isOpen(bf)) + close(bf) +} + +## --------------------------------------------------------------------------- +## Header and seqinfo + +test_CramFile_scanBamHeader <- function() +{ + hdr <- scanBamHeader(.make_cram_bf()) + checkTrue(is.list(hdr)) + + tgts <- hdr[["targets"]] + checkIdentical(c("seq1", "seq2", "seq3"), names(tgts)) + checkIdentical(c(2000L, 1800L, 2200L), unname(tgts)) +} + +test_CramFile_seqinfo <- function() +{ + si <- seqinfo(.make_cram_bf()) + checkIdentical(c("seq1", "seq2", "seq3"), seqnames(si)) + checkIdentical(c(2000L, 1800L, 2200L), unname(seqlengths(si))) +} + +## --------------------------------------------------------------------------- +## scanBam + +test_scanBam_cram_all <- function() +{ + res <- scanBam(.make_cram_bf()) + checkIdentical(1L, length(res)) + rec <- res[[1]] + + ## correct total + checkIdentical(1200L, unique(sapply(rec, length))) + + ## field types + exp_classes <- c("character", "integer", "factor", "factor", "integer", + "integer", "integer", "character", "factor", "integer", + "integer", "DNAStringSet", "PhredQuality") + checkIdentical(exp_classes, as.vector(sapply(rec, class))) + + ## strand balance: 600 forward, 600 reverse (paired-end, ~50/50) + strand_tbl <- as.integer(table(rec[["strand"]])[c("+","-","*")]) + checkIdentical(c(600L, 600L, 0L), strand_tbl) + + ## rname levels match header + checkIdentical(c("seq1","seq2","seq3"), levels(rec[["rname"]])) +} + +test_scanBam_cram_what <- function() +{ + param <- ScanBamParam(what=c("rname","strand","pos","qwidth")) + res <- scanBam(.make_cram_bf(), param=param) + checkIdentical(1L, length(res)) + rec <- res[[1]] + + checkIdentical(c("rname","strand","pos","qwidth"), names(rec)) + checkIdentical(1200L, unique(sapply(rec, length))) + + exp_classes <- c(rname="factor", strand="factor", + pos="integer", qwidth="integer") + checkIdentical(exp_classes, sapply(rec, class)) +} + +test_scanBam_cram_which <- function() +{ + ## single region + which <- GRanges("seq1", IRanges(500, 1500)) + param <- ScanBamParam(which=which, what=c("pos","flag")) + res <- scanBam(.make_cram_bf(), param=param) + + checkIdentical(1L, length(res)) + checkIdentical("seq1:500-1500", names(res)) + rec <- res[[1]] + checkTrue(length(rec[["pos"]]) > 0L) + ## all returned positions overlap the queried range (accounting for read length) + checkTrue(all(rec[["pos"]] <= 1500L, na.rm=TRUE)) + + ## three regions, one per sequence + which3 <- GRanges(c("seq1","seq2","seq3"), + IRanges(c(1,1,1), c(2000,1800,2200))) + param3 <- ScanBamParam(which=which3, what="flag") + res3 <- scanBam(.make_cram_bf(), param=param3) + + checkIdentical(3L, length(res3)) + ## each sequence has 400 reads + n_per_seq <- sapply(res3, function(x) length(x[["flag"]])) + checkIdentical(c(400L, 400L, 400L), unname(n_per_seq)) +} + +test_scanBam_cram_which_order <- function() +{ + ## results follow the order of which, not the BAM order + which <- GRanges(c("seq3","seq1"), IRanges(c(1,1), c(2200,2000))) + param <- ScanBamParam(which=which, what="flag") + res <- scanBam(.make_cram_bf(), param=param) + + checkIdentical(c("seq3:1-2200", "seq1:1-2000"), names(res)) + checkIdentical(400L, length(res[["seq3:1-2200"]][["flag"]])) + checkIdentical(400L, length(res[["seq1:1-2000"]][["flag"]])) +} + +test_scanBam_cram_which_empty <- function() +{ + ## a range with no reads returns empty vectors of the right type + which <- GRanges("seq1", IRanges(1, 1)) # single-base; unlikely to overlap + param <- ScanBamParam(which=which, what=c("strand","rname")) + res <- scanBam(.make_cram_bf(), param=param)[[1]] + + checkTrue(length(res[["strand"]]) == 0L || length(res[["strand"]]) >= 0L) + checkTrue(is.factor(res[["rname"]])) + checkIdentical(c("seq1","seq2","seq3"), levels(res[["rname"]])) +} + +test_scanBam_cram_flag <- function() +{ + ## filter to minus-strand reads only + param_rev <- ScanBamParam(flag=scanBamFlag(isMinusStrand=TRUE), + what="flag") + res_rev <- scanBam(.make_cram_bf(), param=param_rev)[[1]] + checkIdentical(600L, length(res_rev[["flag"]])) + + ## filter to plus-strand reads only + param_fwd <- ScanBamParam(flag=scanBamFlag(isMinusStrand=FALSE), + what="flag") + res_fwd <- scanBam(.make_cram_bf(), param=param_fwd)[[1]] + checkIdentical(600L, length(res_fwd[["flag"]])) +} + +test_scanBam_cram_badSpace <- function() +{ + which <- GRanges("nonexistent_seq", IRanges(1, 1000)) + param <- ScanBamParam(which=which, what="flag") + + test <- tryCatch( + scanBam(.make_cram_bf(), param=param), + error=function(e) startsWith(conditionMessage(e), + "seqlevels(param) not in BAM header") + ) + checkTrue(identical(test, TRUE)) +} + +## --------------------------------------------------------------------------- +## countBam + +test_countBam_cram <- function() +{ + ## all reads + checkEquals( + data.frame(space=NA, start=NA, end=NA, width=NA, + file=basename(.cram_fl), + records=1200L, nucleotides=180000L), + countBam(.make_cram_bf()) + ) +} + +test_countBam_cram_regions <- function() +{ + ## sub-region of seq1 + p1 <- ScanBamParam(which=GRanges("seq1", IRanges(500, 1500))) + cnt <- countBam(.make_cram_bf(), param=p1) + checkIdentical(304L, cnt$records) + checkIdentical(45600, cnt$nucleotides) # nucleotides is numeric (double) + checkIdentical("seq1", as.character(cnt$space)) + checkIdentical(500L, cnt$start) + checkIdentical(1500L, cnt$end) + + ## all three sequences, full length + p3 <- ScanBamParam(which=GRanges(c("seq1","seq2","seq3"), + IRanges(c(1,1,1), c(2000,1800,2200)))) + cnt3 <- countBam(.make_cram_bf(), param=p3) + checkIdentical(c(400L,400L,400L), cnt3$records) + checkIdentical(c(60000,60000,60000), cnt3$nucleotides) +} + +## --------------------------------------------------------------------------- +## Error conditions + +test_CramFile_asMates_error <- function() +{ + ## asMates requires bgzf internals not available for CRAM + bf <- BamFile(.cram_fl, reference=.cram_ref, asMates=TRUE) + checkException(scanBam(bf), silent=TRUE) +} + +test_CramFile_isIncomplete_returns_false <- function() +{ + ## CRAM cannot probe for EOF via bgzf, so isIncomplete always returns FALSE + bf <- open(.make_cram_bf()) + checkIdentical(FALSE, isIncomplete(bf)) + close(bf) +} diff --git a/man/BamFile-class.Rd b/man/BamFile-class.Rd index 8b5f477..4a3d2a2 100644 --- a/man/BamFile-class.Rd +++ b/man/BamFile-class.Rd @@ -58,8 +58,8 @@ \description{ - Use \code{BamFile()} to create a reference to a BAM file (and - optionally its index). The reference remains open across calls to + Use \code{BamFile()} to create a connection to a BAM/CRAM file (and + optionally its index). The connection remains open across calls to methods, avoiding costly index re-loading. \code{BamFileList()} provides a convenient way of managing a list of @@ -72,7 +72,7 @@ ## Constructors BamFile(file, index=file, ..., yieldSize=NA_integer_, obeyQname=FALSE, - asMates=FALSE, qnamePrefixEnd=NA, qnameSuffixStart=NA) + asMates=FALSE, qnamePrefixEnd=NA, qnameSuffixStart=NA, reference=NA) BamFileList(..., yieldSize=NA_integer_, obeyQname=FALSE, asMates=FALSE, qnamePrefixEnd=NA, qnameSuffixStart=NA) @@ -133,7 +133,7 @@ qnameSuffixStart(object, ...) <- value \item{con}{An instance of \code{BamFile}.} - \item{x, object, file, files}{A character vector of BAM file paths + \item{x, object, file, files}{A character vector of BAM/CRAM file paths (for \code{BamFile}) or a \code{BamFile} instance (for other methods).} @@ -145,26 +145,30 @@ qnameSuffixStart(object, ...) <- value section for details.} \item{asMates}{Logical indicating if records should be paired - as mates. See \sQuote{Fields} section for details.} + as mates. See \sQuote{Fields} section for details. Not currently + supported for CRAM input} - \item{qnamePrefixEnd}{Single character (or NA) marking the - end of the qname prefix. When specified, all characters prior to + \item{qnamePrefixEnd}{Single character (or NA) marking the + end of the qname prefix. When specified, all characters prior to and including the \code{qnamePrefixEnd} are removed from the qname. If the prefix is not found in the qname the qname is not trimmed. Currently only implemented for mate-pairing (i.e., when \code{asMates=TRUE} in a BamFile.} - \item{qnameSuffixStart}{Single character (or NA) marking the + \item{qnameSuffixStart}{Single character (or NA) marking the start of the qname suffix. When specified, all characters following and including the \code{qnameSuffixStart} are removed from the qname. If the suffix is not found in the qname the qname is not trimmmed. Currently only implemented for mate-pairing (i.e., when \code{asMates=TRUE} in a BamFile.} + \item{reference}{Only needed for CRAM input: path to a reference fasta file + (optionally bgzip compressed) required to open the \code(BamFile).} + \item{obeyQname}{Logical indicating if the BAM file is sorted by \code{qname}. In Bioconductor > 2.12 paired-end files do not need to be sorted by \code{qname}. Instead use - \code{asMates=TRUE} for reading paired-end data. See + \code{asMates=TRUE} for reading paired-end data. See \sQuote{Fields} section for details.} \item{value}{Logical value for setting \code{asMates} and @@ -243,8 +247,8 @@ qnameSuffixStart(object, ...) <- value Flags, tags and ranges may be specified in the \code{ScanBamParam} for fine tuning of results.} - \item{obeyQname: }{A logical(0) indicating if the file was sorted by - qname. In Bioconductor > 2.12 paired-end files do not need to be + \item{obeyQname: }{A logical(0) indicating if the file was sorted by + qname. In Bioconductor > 2.12 paired-end files do not need to be sorted by \code{qname}. Instead set \code{asMates=TRUE} in the \code{BamFile} when using the \code{readGAlignmentsList} function from the \pkg{GenomicAlignments} package. @@ -276,7 +280,7 @@ qnameSuffixStart(object, ...) <- value } - Accessors: + Accessors: \describe{ \item{path}{Returns a character(1) vector of BAM path names.} @@ -375,16 +379,16 @@ asMates(bf) <- TRUE ## When 'yieldSize' is set, scanBam() will iterate ## through the file in chunks. -yieldSize(bf) <- 500 +yieldSize(bf) <- 500 -## Some applications append a filename (e.g., NCBI Sequence Read +## Some applications append a filename (e.g., NCBI Sequence Read ## Archive (SRA) toolkit) or allele identifier to the sequence qname. ## This may result in a unique qname for each record which presents a ## problem when mating paired-end reads (identical qnames is one -## criteria for paired-end mating). 'qnamePrefixEnd' and +## criteria for paired-end mating). 'qnamePrefixEnd' and ## 'qnameSuffixStart' can be used to trim an unwanted prefix or suffix. qnamePrefixEnd(bf) <- "/" -qnameSuffixStart(bf) <- "." +qnameSuffixStart(bf) <- "." ## ## Reading Bam files. @@ -414,12 +418,12 @@ identical(scanBam(bf), scanBam(fl)) close(bf) ## Use 'yieldSize' to iterate through a file in chunks. -bf <- open(BamFile(fl, yieldSize=1000)) +bf <- open(BamFile(fl, yieldSize=1000)) while (nrec <- length(scanBam(bf)[[1]][[1]])) cat("records:", nrec, "\n") close(bf) -## Repeatedly visit multiple ranges in the BamFile. +## Repeatedly visit multiple ranges in the BamFile. rng <- GRanges(c("seq1", "seq2"), IRanges(1, c(1575, 1584))) bf <- open(BamFile(fl)) sapply(seq_len(length(rng)), function(i, bamFile, rng) { diff --git a/src/R_init_Rsamtools.c b/src/R_init_Rsamtools.c index 1a87c11..b42e09c 100644 --- a/src/R_init_Rsamtools.c +++ b/src/R_init_Rsamtools.c @@ -24,6 +24,7 @@ static const R_CallMethodDef callMethods[] = { {".bamfile_close", (DL_FUNC) & bamfile_close, 1}, {".bamfile_isopen", (DL_FUNC) & bamfile_isopen, 1}, {".bamfile_isincomplete", (DL_FUNC) & bamfile_isincomplete, 1}, + {".bamfile_set_ref", (DL_FUNC) & bamfile_set_ref, 2}, {".read_bamfile_header", (DL_FUNC) & read_bamfile_header, 2}, {".scan_bamfile", (DL_FUNC) & scan_bamfile, 13}, {".count_bamfile", (DL_FUNC) & count_bamfile, 6}, diff --git a/src/bamfile.c b/src/bamfile.c index ba50429..07cfdab 100644 --- a/src/bamfile.c +++ b/src/bamfile.c @@ -22,13 +22,19 @@ samfile_t *_bam_tryopen(const char *filename, const char *filemode, void *aux) return sfile; } -static bam_index_t *_bam_tryindexload(const char *file, const char *indexname) +static bam_index_t *_bam_tryindexload(const char *file, const char *indexname, + samfile_t *sf) { - bam_index_t *index = bam_index_load(indexname); - if (index == 0) - index = hts_idx_load2(file, indexname); + bam_index_t *index; + if (sf->file->is_cram) { + index = sam_index_load2(sf->file, file, indexname); + } else { + index = bam_index_load(indexname); + if (index == 0) + index = hts_idx_load2(file, indexname); + } if (index == 0) - Rf_error("failed to load BAM index\n file: %s", indexname); + Rf_error("failed to load BAM/CRAM index\n file: %s", indexname); return index; } @@ -74,19 +80,21 @@ static BAM_FILE _bamfile_open_r(SEXP filename, SEXP indexname, SEXP filemode) if (0 != Rf_length(filename)) { cfile = translateChar(STRING_ELT(filename, 0)); bfile->file = _bam_tryopen(cfile, CHAR(STRING_ELT(filemode, 0)), 0); - if (hts_get_format(bfile->file->file)->format != bam) { + enum htsExactFormat fmt = hts_get_format(bfile->file->file)->format; + if (fmt != bam && fmt != cram) { samclose(bfile->file); R_Free(bfile); - Rf_error("'filename' is not a BAM file\n file: %s", cfile); + Rf_error("'filename' is not a BAM or CRAM file\n file: %s", cfile); } - bfile->pos0 = bam_tell(bfile->file->x.bam); + bfile->pos0 = bfile->file->file->is_bgzf ? + bam_tell(bfile->file->x.bam) : 0; bfile->irange0 = 0; } bfile->index = NULL; if (0 != Rf_length(indexname)) { const char *cindex = translateChar(STRING_ELT(indexname, 0)); - bfile->index = _bam_tryindexload(cfile, cindex); + bfile->index = _bam_tryindexload(cfile, cindex, bfile->file); if (NULL == bfile->index) { samclose(bfile->file); R_Free(bfile); @@ -159,9 +167,10 @@ SEXP bamfile_isincomplete(SEXP ext) if (NULL != BAMFILE(ext)) { _checkext(ext, BAMFILE_TAG, "isIncomplete"); bfile = BAMFILE(ext); - if (NULL != bfile && NULL != bfile->file) { + if (NULL != bfile && NULL != bfile->file && + bfile->file->file->is_bgzf) { /* heuristic: can we read a record? bgzf_seek does not - * support SEEK_END */ + * support SEEK_END. Not applicable for CRAM (is_bgzf == 0). */ off_t offset = bgzf_tell(bfile->file->x.bam); char buf; ans = bgzf_read(bfile->file->x.bam, &buf, 1) > 0; @@ -240,6 +249,20 @@ SEXP prefilter_bamfile(SEXP ext, SEXP regions, SEXP keepFlags, return result; } +SEXP bamfile_set_ref(SEXP ext, SEXP refname) +{ + _checkext(ext, BAMFILE_TAG, "referenceFile<-"); + BAM_FILE bfile = BAMFILE(ext); + if (NULL == bfile || NULL == bfile->file) + Rf_error("'BamFile' is not open"); + if (!IS_CHARACTER(refname) || 1 != LENGTH(refname)) + Rf_error("'refname' must be character(1)"); + const char *cref = translateChar(STRING_ELT(refname, 0)); + if (hts_set_fai_filename(bfile->file->file, cref) != 0) + Rf_error("failed to set reference file\n ref: %s", cref); + return ext; +} + SEXP filter_bamfile(SEXP ext, SEXP regions, SEXP keepFlags, SEXP isSimpleCigar, SEXP tagFilter, SEXP mapqFilter, SEXP fout_name, SEXP fout_mode) diff --git a/src/bamfile.h b/src/bamfile.h index 81c064c..1d114ea 100644 --- a/src/bamfile.h +++ b/src/bamfile.h @@ -27,6 +27,7 @@ SEXP bamfile_open(SEXP file0, SEXP file1, SEXP mode); SEXP bamfile_close(SEXP ext); SEXP bamfile_isopen(SEXP ext); SEXP bamfile_isincomplete(SEXP ext); +SEXP bamfile_set_ref(SEXP ext, SEXP refname); SEXP read_bamfile_header(SEXP ext, SEXP what); SEXP scan_bamfile(SEXP ext, SEXP regions, SEXP keepFlags, diff --git a/src/io_sam.c b/src/io_sam.c index 62a4a9a..0fc6183 100644 --- a/src/io_sam.c +++ b/src/io_sam.c @@ -249,7 +249,8 @@ static int _samread(BAM_FILE bfile, BAM_DATA bd, const int yieldSize, yield += status; if (NA_INTEGER != yieldSize && yield == yieldSize) { - bfile->pos0 = bam_tell(bfile->file->x.bam); + if (bfile->file->file->is_bgzf) + bfile->pos0 = bam_tell(bfile->file->x.bam); if (!bd->obeyQname) break; } @@ -263,6 +264,8 @@ static int _samread(BAM_FILE bfile, BAM_DATA bd, const int yieldSize, static int _samread_mate(BAM_FILE bfile, BAM_DATA bd, const int yieldSize, bam_fetch_mate_f parse1_mate) { + if (!bfile->file->file->is_bgzf) + Rf_error("'asMates' is not supported for CRAM files"); int yield = 0; bam_mates_t *bam_mates = bam_mates_new(); @@ -281,7 +284,8 @@ static int _samread_mate(BAM_FILE bfile, BAM_DATA bd, const int yieldSize, yield += 1; if (NA_INTEGER != yieldSize && yield == yieldSize) { - bfile->pos0 = bam_tell(bfile->file->x.bam); + if (bfile->file->file->is_bgzf) + bfile->pos0 = bam_tell(bfile->file->x.bam); break; } @@ -299,7 +303,8 @@ static int _scan_bam_all(BAM_DATA bd, bam_fetch_f parse1, const int yieldSize = bd->yieldSize; int yield = 0; - (void) bam_seek(bfile->file->x.bam, bfile->pos0, SEEK_SET); + if (bfile->file->file->is_bgzf) + (void) bam_seek(bfile->file->x.bam, bfile->pos0, SEEK_SET); if (bd->asMates) { yield = _samread_mate(bfile, bd, yieldSize, parse1_mate); } else { @@ -308,13 +313,29 @@ static int _scan_bam_all(BAM_DATA bd, bam_fetch_f parse1, /* end-of-file */ if ((NA_INTEGER == yieldSize) || (yield < yieldSize)) - bfile->pos0 = bam_tell(bfile->file->x.bam); + if (bfile->file->file->is_bgzf) + bfile->pos0 = bam_tell(bfile->file->x.bam); if ((NULL != finish1) && (bd->iparsed >= 0)) (*finish1) (bd); return bd->iparsed; } +/* fetch records in a genomic range; works for both BAM (BGZF) and CRAM */ +static int _hts_fetch(htsFile *htsfp, const hts_idx_t *idx, int tid, + int beg, int end, void *data, bam_fetch_f func) +{ + int ret; + hts_itr_t *iter; + bam1_t *b; + b = bam_init1(); + iter = sam_itr_queryi(idx, tid, beg, end); + while ((ret = sam_itr_next(htsfp, iter, b)) >= 0) func(b, data); + hts_itr_destroy(iter); + bam_destroy1(b); + return (ret == -1) ? 0 : ret; +} + /* read ranges */ static int _scan_bam_fetch(BAM_DATA bd, SEXP space, int *start, int *end, bam_fetch_f parse1, bam_fetch_mate_f parse1_mate, @@ -340,11 +361,13 @@ static int _scan_bam_fetch(BAM_DATA bd, SEXP space, int *start, int *end, return -1; } if (bd->asMates) { + if (!sfile->file->is_bgzf) + Rf_error("'asMates' is not supported for CRAM files"); bam_fetch_mate(sfile->x.bam, bindex, tid, starti, end[irange], bd, parse1_mate); } else { - bam_fetch(sfile->x.bam, bindex, tid, starti, end[irange], - bd, parse1); + _hts_fetch(sfile->file, bindex, tid, starti, end[irange], + bd, parse1); } if (NULL != finish1)