Search
Question: Rsamtools fetch MD:Z field from bam files
0
gravatar for Julie Zhu
8 weeks ago by
Julie Zhu3.7k
United States
Julie Zhu3.7k wrote:

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

 

 

ADD COMMENTlink modified 8 weeks ago • written 8 weeks ago by Julie Zhu3.7k
0
gravatar for Martin Morgan
8 weeks ago by
Martin Morgan ♦♦ 20k
United States
Martin Morgan ♦♦ 20k wrote:

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

ADD COMMENTlink written 8 weeks ago by Martin Morgan ♦♦ 20k

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 REPLYlink modified 8 weeks ago by Martin Morgan ♦♦ 20k • written 8 weeks ago by Julie Zhu3.7k

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 REPLYlink written 8 weeks ago by Martin Morgan ♦♦ 20k
0
gravatar for Julie Zhu
8 weeks ago by
Julie Zhu3.7k
United States
Julie Zhu3.7k wrote:

Martin,

Thanks so much!

Best regards,

Julie

ADD COMMENTlink written 8 weeks ago by Julie Zhu3.7k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.2.0
Traffic: 239 users visited in the last hour