RSamtools idxstatsBam not working with remote bam file
1
0
Entering edit mode
mja18 ▴ 40
@mja18-23865
Last seen 3.6 years ago
United States

idxstatsBam is failing when used on a remote (http) bam file. The corresponding .bai index file is definitely there on the server (from bash samtools idxstats works fine). The bam file and index also seem to be fine; if downloaded to my laptop, idxstatsBam works fine.

Other RSamtools functions, such as scanBam seem to find the remote just fine, as shown below.

Would really like to use idxstatsBam with http-hosted bam files for quick way to count all reads in a file. Any thoughts?


# does not find remote index
idxstatsBam("http://plantsmallrnagenes.science.psu.edu/ccm-b0.32/alignments/SRR3977128_ca.bam")
Error in idxstatsBam(bam, ...) : file.exists(index(file)) is not TRUE

# but scanBam can find and use the remote index
which <- IRangesList(Ccam0.32_scaffold111 = IRanges(801256, 801507))
what <- scanBamWhat()
param <- ScanBamParam(which=which, what=what)
scanBam("http://plantsmallrnagenes.science.psu.edu/ccm-b0.32/alignments/SRR3977128_ca.bam", 
        param = param)

sessionInfo()
R version 4.0.2 (2020-06-22)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.7

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] Rsamtools_2.4.0      Biostrings_2.56.0    XVector_0.28.0       GenomicRanges_1.40.0 GenomeInfoDb_1.24.2 
[6] IRanges_2.22.2       S4Vectors_0.26.1     BiocGenerics_0.34.0 

loaded via a namespace (and not attached):
 [1] crayon_1.4.1           digest_0.6.27          bitops_1.0-6           evaluate_0.14         
 [5] zlibbioc_1.34.0        rlang_0.4.10           rstudioapi_0.13        rmarkdown_2.7         
 [9] BiocParallel_1.22.0    tools_4.0.2            tinytex_0.30           RCurl_1.98-1.3        
[13] xfun_0.22              yaml_2.2.1             compiler_4.0.2         htmltools_0.5.1.1     
[17] knitr_1.31             GenomeInfoDbData_1.2.3
Rsamtools • 1.3k views
ADD COMMENT
1
Entering edit mode
@martin-morgan-1513
Last seen 4 months ago
United States

Can you tell me if the following works for you? In a new R session

library(Rsamtools)

setMethod(idxstatsBam, "BamFile",
          function(file, index=file, ...)
{
    if (!isOpen(file)) {
        open(file)
        on.exit(close(file))
    }
    result <- .Call(Rsamtools:::.idxstats_bamfile, Rsamtools:::.extptr(file))
    seqnames <- factor(result[[1]], levels=sortSeqlevels(unique(result[[1]])))
    o <- order(seqnames)
    data.frame(seqnames=seqnames[o], seqlength=result[[2]][o],
               mapped=result[[3]][o], unmapped=result[[4]][o])
})

and then running your own code

idxstatsBam("http://plantsmallrnagenes.science.psu.edu/ccm-b0.32/alignments/SRR3977128_ca.bam")

I will update Rsamtools along these lines...

ADD COMMENT
0
Entering edit mode

Thanks Martin. Yes, that works.

ADD REPLY

Login before adding your answer.

Traffic: 927 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6