is the Rsamtools::pileup example broken?
1
2
Entering edit mode
@jeremy-leipzig-4924
Last seen 6.1 years ago
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   
pileup rsamtools • 1.2k views
ADD COMMENT
1
Entering edit mode

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

ADD REPLY
4
Entering edit mode
@martin-morgan-1513
Last seen 5 months ago
United States

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 COMMENT

Login before adding your answer.

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