Add CRAM support

DESCRIPTION

@@ -10,14 +10,15 @@ 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(
     person("Martin", "Morgan", role = "aut"),
     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",

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,15 +344,15 @@ 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
 
     o bamPath and bamIndicies applied to BamViews returns named
       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.
 

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"))

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"))
 

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,21 +83,35 @@ 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 <-
     function(con, ...)
 {
     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")
 })

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

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()