Is there an efficient way to check an input BAM in R?

Question

I'm writing a function in R for an R package which takes as input a BAM.

my_func = function(my_bam){
    if (!file.exists(my_bam)){
        stop("Input a valid BAM file.")
    }

    ## do awesome stuff with BAM

}

Trying to write robust software, I would like to do more of a check whether this BAM is actually a BAM. Above, I'm just using file.exists().

I could also do something simple like "check whether the file has the extension *.bam" with

grepl('*.bam', my_bam)  ## TRUE for file.bam, FALSE for file.txt

However, is there a quick way to check that the input is actually a BAM? I would prefer to throw an error here, if the user inputs a *txt or *bedpe.

I don't think it's necessary to validate the BAM per se, whereby I believe most users rely upon Picard:

https://broadinstitute.github.io/picard/command-line-overview.html#ValidateSamFile

My first thought was Rsamtools for this operation:

https://bioconductor.org/packages/release/bioc/html/Rsamtools.html

However, it doesn't appear that there is a function here which checks the user inputs a BAM, e.g. BamFile() doesn't appear to throw an error

> library(Rsamtools)
> fake_bam = 'nonsense_string.bam'   ## a string I made up, not a real BAM
> BamFile(fake_bam)
class: BamFile 
path: nonsense_string.bam
index: NA
isOpen: FALSE 
yieldSize: NA 
obeyQname: FALSE 
asMates: FALSE 
qnamePrefixEnd: NA 
qnameSuffixStart: NA 

How could I "check" that the user is inputting a BAM file? Perhaps pipe to samtools (which doesn't feel very stable in terms of software)?

Answer 1

You won’t know whether it’s a valid BAM file until you attempt reading it in full; even Devon’s method (which I recommend) only does a superficial check — i.e. it checks whether (parts of) the BAM header can be read correctly. However, this doesn’t guarantee that the file isn’t truncated or contains invalid data.

This is of course a more general problem, but it’s perhaps a bit stranger for BAMs (and other biological formats), as it’s relatively frequent to end up with technically invalid files that were produced by existing software. The only recourse is actually trying to use the file, and to handle errors appropriately.

For a very quick check™ you could also pass the file through gunzip and test whether the “magic bytes” at the beginning of the file equal 'BAM\1':

bgzf_stream = gzfile(filename, 'r')
on.exit(close(bgzf_stream))
magic = readChar(bgzf_stream, 4)
stopifnot(identical(magic, 'BAM\1'))

Answer 2

Take your BamFile(fake_bam) and pass it through open(). That will return a value if fake_bam is a real BAM (or CRAM or maybe SAM) file and through an error otherwise. You can then wrap all of that in a trycatch():

trycatch(open(BamFile(fake_bam)), do something, do something else)

Alternatively, ditch the trycatch() and let open() return the error if an invalid file in entered.

Answer 3

Just expanding upon the answers by @DevonRyan and @KonradRudolph, here's one way to put this in an R package:

Either using open(BamFile() from Devon's answer, which needs Rsamtools:

my_func = function(my_bam){

    out = tryCatch(open(Rsamtools::BamFile(my_bam)), error = function(e) e)
    if (inherits(out, "error")){
        stop("Input a valid BAM file.")
    }

    ## do awesome stuff with BAM
}

Or from Konrad's answer

my_func = function(my_bam){

    magic = readChar(gzfile(filename, 'r'), 4)
    if (!identical(magic, 'BAM\1')){
        stop("Input a valid BAM file.")
    }

    ## do awesome stuff with BAM
}