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

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

ADD COMMENTlink written 11 months ago by Martin Morgan ♦♦ 21k

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

Martin,

Thanks so much!

Best regards,

Julie

ADD COMMENTlink written 11 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: 136 users visited in the last hour