Question: is the Rsamtools::pileup example broken?
2
gravatar for Jeremy Leipzig
22 months ago by
Jeremy Leipzig70 wrote:
BiocInstaller::biocLite(("RNAseqData.HNRNPC.bam.chr14"))
BiocInstaller::biocLite("Rsamtools")
library(Rsamtools)
library("RNAseqData.HNRNPC.bam.chr14")
fl <- RNAseqData.HNRNPC.bam.chr14_BAMFILES[1]
res <- pileup(fl)
head(res)
> res <- pileup(fl)
> 
  > head(res)
[1] seqnames   pos        strand     nucleotide count     
<0 rows> (or 0-length row.names)
> table(res$strand, res$nucleotide)

A C G T N = - +
  + 0 0 0 0 0 0 0 0
- 0 0 0 0 0 0 0 0
* 0 0 0 0 0 0 0 0

> sessionInfo()
R version 3.4.1 (2017-06-30)
Platform: x86_64-redhat-linux-gnu (64-bit)
Running under: Amazon Linux AMI 2017.09

Matrix products: default
BLAS/LAPACK: /efs/R/lib/libRblas.so

locale:
  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
[5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                 
[9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

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

other attached packages:
  [1] RNAseqData.HNRNPC.bam.chr14_0.12.0 Rsamtools_1.26.2                   Biostrings_2.42.1                 
[4] XVector_0.14.1                     GenomicRanges_1.26.4               GenomeInfoDb_1.10.3               
[7] IRanges_2.8.2                      S4Vectors_0.12.2                   BiocGenerics_0.20.0               

loaded via a namespace (and not attached):
  [1] pcaPP_1.9-72         Rcpp_0.12.11         DEoptimR_1.0-8       compiler_3.4.1       BiocInstaller_1.24.0 RColorBrewer_1.1-2  
[7] zlibbioc_1.20.0      bitops_1.0-6         tools_3.4.1          dotCall64_0.9-5      lattice_0.20-35      graph_1.52.0        
[13] mvtnorm_1.0-6        spam_2.1-2           hexbin_1.27.1        cluster_2.0.6        IDPmisc_1.1.17       fields_9.0          
[19] maps_3.2.0           grid_3.4.1           robustbase_0.92-8    Biobase_2.34.0       rrcov_1.4-3          R6_2.2.2            
[25] BiocParallel_1.8.2   latticeExtra_0.6-28  corpcor_1.6.9        matrixStats_0.52.2   flowFramePlus_0.1.0  MASS_7.3-48         
[31] assertthat_0.2.0     flowCore_1.40.6      KernSmooth_2.23-15   flowViz_1.38.0       RCurl_1.95-4.9   
rsamtools pileup • 388 views
ADD COMMENTlink modified 22 months ago by Martin Morgan ♦♦ 24k • written 22 months ago by Jeremy Leipzig70
1

same after update to Rsamtools_1.30.0 and RNAseqData.HNRNPC.bam.chr14_0.16.0 

ADD REPLYlink written 22 months ago by Jeremy Leipzig70
Answer: is the Rsamtools::pileup example broken?
4
gravatar for Martin Morgan
22 months ago by
Martin Morgan ♦♦ 24k
United States
Martin Morgan ♦♦ 24k wrote:

Nice question! the problem is that the base qualities of the reads in the example file are, for some reason, low quality, so they fail to pass the min_base_quality default of PileupParam(). So...

> fl <- RNAseqData.HNRNPC.bam.chr14_BAMFILES[1]
> param <- PileupParam(min_base_quality = 0L)
> res <- pileup(fl, pileupParam = param)
> head(res)
  seqnames      pos strand nucleotide count
1    chr14 19069583      +          T     1
2    chr14 19069584      +          G     1
3    chr14 19069585      +          A     1
4    chr14 19069586      +          G     1
5    chr14 19069587      +          A     1
6    chr14 19069588      +          A     1
> table(res$strand, res$nucleotide)

         A      C      G      T      N      =      -      +
  + 535270 569367 560276 539855      0      0    474      0
  - 537559 558364 571841 538627      0      0    465      0
  *      0      0      0      0      0      0      0      0

This commit introduced the change that caused the problem (changing the default min_base_quality to 13 from 10, to be consistent with samtools). I'll update the man page.

ADD COMMENTlink modified 22 months ago • written 22 months ago by Martin Morgan ♦♦ 24k
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 16.09
Traffic: 206 users visited in the last hour