I am working on PAR-CLIP analysis using the wavClusteR package. When I try to load a bam file, I get the following error:

> filename<-("test.bam")
> Bam <- readSortedBam(filename = filename)
> Bam
> countTable <- getAllSub( Bam, minCov = 10 )
Loading required package: doMC
Considering substitutions, n = 1999, processing in 2 chunks
   chunk #: 1
   chunk #: 2
Error in { :
  task 1 failed - "solving row 247: range cannot be determined from the supplied arguments (too many NAs)"
In addition: Warning message:
executing %dopar% sequentially: no parallel backend registered

Any ideas on what could be causing this? I suspect this it has something to do with GRanges. Here is the portion of the file that is being noted in the error message:

> Bam[247]
GRanges object with 1 range and 2 metadata columns:
      seqnames               ranges strand |
         <Rle>            <IRanges>  <Rle> |
  [1]     chr1 [91853066, 91853115]      + |
                                                    qseq          MD
                                          <DNAStringSet> <character>
  [1] ACTGAGCTCGCCTTAGGACACCTGCGTTACCGTTTGACAGGTGTACCGCC        45A3
  -------
  seqinfo: 25 sequences from an unspecified genome; no seqlengths

This seems pretty normal, though. Suggestions?

Session info:

> sessionInfo()
R version 3.2.3 (2015-12-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS release 6.3 (Final)

locale:
[1] C

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

other attached packages:
[1] wavClusteR_2.4.1     Rsamtools_1.22.0     Biostrings_2.38.4
[4] XVector_0.10.0       GenomicRanges_1.22.4 GenomeInfoDb_1.6.3
[7] IRanges_2.4.8        S4Vectors_0.8.11     BiocGenerics_0.16.1

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.4                wmtsa_2.0-0
 [3] compiler_3.2.3             RColorBrewer_1.1-2
 [5] futile.logger_1.4.1        plyr_1.8.3
 [7] GenomicFeatures_1.22.13    bitops_1.0-6
 [9] futile.options_1.0.0       iterators_1.0.8
[11] tools_3.2.3                zlibbioc_1.16.0
[13] rpart_4.1-10               biomaRt_2.26.1
[15] mclust_5.2                 lattice_0.20-33
[17] RSQLite_1.0.0              gtable_0.2.0
[19] foreach_1.4.3              DBI_0.3.1
[21] gridExtra_2.2.1            cluster_2.0.3
[23] rtracklayer_1.30.4         stringr_1.0.0
[25] ifultools_2.0-1            ade4_1.7-4
[27] nnet_7.3-12                grid_3.2.3
[29] Biobase_2.30.0             AnnotationDbi_1.32.3
[31] survival_2.38-3            XML_3.98-1.4
[33] BiocParallel_1.4.3         foreign_0.8-66
[35] latticeExtra_0.6-28        Formula_1.2-1
[37] seqinr_3.1-3               ggplot2_2.1.0
[39] lambda.r_1.1.7             splus2R_1.2-0
[41] magrittr_1.5               splines_3.2.3
[43] Hmisc_3.17-3               MASS_7.3-45
[45] scales_0.4.0               codetools_0.2-14
[47] GenomicAlignments_1.6.3    SummarizedExperiment_1.0.2
[49] colorspace_1.2-6           stringi_1.0-1
[51] acepack_1.3-3.3            RCurl_1.95-4.8
[53] munsell_0.4.3

Hi Stephen,

to generate a substitution count table, wavClusteR splits a sorted Bam file by strand, flags all entries with an MD value other than perfect matches, and generates a GRangesList object where each slot contains the set of genomic regions that share the same MD. Elements of this list are then processed in chunks of size 1000 (using multiple cores if enabled). 

You sent me a chunk of the Bam file that was sufficient to generate the error you reported. This error is caused by the alignment:

$10^A38 
GRanges object with 1 range and 2 metadata columns:
      seqnames             ranges strand |                                               qseq          MD
         <Rle>          <IRanges>  <Rle> |                                     <DNAStringSet> <character>
  [1]    chr21 [9827478, 9827527]      + | CCGCAGGCGCGCAATTACCCACTCCCGACCCGGGGAGGTAGTGACGAAAA      10^A38
 

because the MD value 10^A38 (10 matches followed by deletion of one A and 38 matches) contains an indel. wavCluster does not support indels. How did you align you PAR-CLIP data? 

Let me know. 

Hope this helps,

Federico

Hi Stephen,

thank you for posting your issue on the support site. I realised the tag for the threads related to wavClusteR has been consistently misspelled to wavcluter. Is it possible to edit it?

Re. your issue, would it be possible for you to try out v2.5.2 in BioC devel and report back? As we discussed earlier today by email, the issue you reported has nothing to do with line 247 in your Bam file. Rather, it refers to the row 247 in the count table.

I will work on it and report back asap.

Let me know.

Federico

Hi Stephen,

did you have the chance to run your code with version 2.5.2 of the package?