Rsamtools fetch MD:Z field from bam files
2
0
Entering edit mode
Julie Zhu ★ 4.3k
@julie-zhu-3596
Last seen 5 months ago
United States

Hi, 

Is it possible to fetch the MD:Z field, e.g., 30G44 from bam files using Rsamtools or other Bioconductor packages? Thanks!

Best regards,

Julie

 

 

Rsamtools fetch MD:Z field from bam files • 1.1k views
ADD COMMENT
0
Entering edit mode
@martin-morgan-1513
Last seen 3 days ago
United States

Are ScanBamParam() arguments tag and tagFilter what you are looking for? See ?ScanBamParam.

ADD COMMENT
0
Entering edit mode

Martin,

Thanks for the quick response!

I tried to set the tag parameter, but NULL is returned. Here is the code and sessionInfo. 

pram1 <- ScanBamParam(
    tag=c( "MD", "XN", "AS"),
    what=c("rname", "strand", "pos", "qwidth",  "cigar", "mapq", "qual")
)
bam <- scanBam(bamPaths, index = bamIndex, pram = pram1)
str(bam[[1]][["tag"]])
NULL

I also tried with the example code in ScanBamParam with MD added to the tag, but for some reason, NULL is set in the MD field.

fl <- system.file("extdata", "ex1.bam", package="Rsamtools", mustWork=TRUE)
p4 <- ScanBamParam(tag=c("NM", "H1", "MD"), what="flag")
bam4 <- scanBam(fl, param=p4)     
str(bam4[[1]][["tag"]])
List of 3
 $ NM: int [1:3307] 0 0 0 5 0 1 0 1 0 1 ...
 $ H1: int [1:3307] 0 0 0 0 0 1 0 1 0 0 ...
 $ MD: NULL

Best regards,

Julie

sessionInfo()
R Under development (unstable) (2016-06-30 r70858)
Platform: x86_64-apple-darwin12.5.0 (64-bit)
Running under: OS X Mountain Lion 10.8.5
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 
[8] methods   base     

other attached packages:
[1] Rsamtools_1.27.11     Biostrings_2.43.2     XVector_0.15.0       
[4] GenomicRanges_1.27.18 GenomeInfoDb_1.11.6   IRanges_2.9.14       
[7] S4Vectors_0.13.5      BiocGenerics_0.21.2  

loaded via a namespace (and not attached):
[1] zlibbioc_1.21.0    tools_3.4.0        BiocParallel_1.9.4 bitops_1.0-6 
ADD REPLY
0
Entering edit mode

There is a typo in your call to ScanBamParam -- pram= instead of param=. For the system file, there is no record with the MD tag, so no way for R to know what type of vector to allocate for it, hence NULL.

ADD REPLY
0
Entering edit mode
Julie Zhu ★ 4.3k
@julie-zhu-3596
Last seen 5 months ago
United States

Martin,

Thanks so much!

Best regards,

Julie

ADD COMMENT

Login before adding your answer.

Traffic: 690 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