Search
Question: Rsamtools fetch MD:Z field from bam files
0
gravatar for Julie Zhu
4 months ago by
Julie Zhu3.8k
United States
Julie Zhu3.8k 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 4 months ago • written 4 months ago by Julie Zhu3.8k
0
gravatar for Martin Morgan
4 months 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 4 months 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 4 months ago by Martin Morgan ♦♦ 20k • written 4 months ago by Julie Zhu3.8k

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

Martin,

Thanks so much!

Best regards,

Julie

ADD COMMENTlink written 4 months ago by Julie Zhu3.8k
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: 183 users visited in the last hour