First off, I think you can simplify your function to:

ExtractGATCFragments <- function(chr) {
  gaps(matchPattern(chr, pattern='GATC'))
}

But to get to the result you might want, instead of the bsapply and
RangedData stuff, just do:

gr <- vmatchPDict("GATC", Hsapiens)

The 'gr' represents the ranges of GATC in the genome. It matches on both
strands, which is not useful for you since GATC palindromic, so it might be
slower than you want. This also assumes you're not interested in keeping
the genomic sequence around (it's big, and why?). But anyway, just do this
to get your result:

ans <- subset(gaps(gr), strand=="+") # may require updating your R/Bioc

But if you do want to use bsapply, you should at least coerce to GRanges,
not RangedData, which is not suitable for this type of data. Unfortunately,
GenomeData is a dinosaur and there is no direct coercion to GRanges, but
RangesList is a link to the past:

gr <- as(as(GATC.fragments, "RangesList"), "GRanges")

Really what we need for this use case is an lapply method for BSgenome,
then we would just coerce to List, which would result in a RangesList,
assuming the user function returned Ranges, and then coerce to GRanges.
And/or clean up vmatchPDict,BSgenome so that it has an option to match only
to the positive strand.

Michael


On Wed, May 7, 2014 at 8:05 AM, Ludo Pagie <ludo.pagie@gmail.com> wrote:

> Hi all,
>
> I'm getting an error when trying to coerce a GenomeData object to a
> RangedData object. Here's the code I have used together with some
> output (in particular the error msg when coercing):
>
>
> ######################################################################
> # import human genome sequence
> library(BSgenome.Hsapiens.UCSC.hg19)
>
> # Virtual digest of the entire genome in GATC fragments:
> # matchPattern to find GATC motifs, using bsapply
> ExtractGATCFragments <- function(chr) {
> # function for finding GATC sites in a chromosome and returning the
> fragments
>   # starting and ending with GATC
>   GATC.match <- matchPattern(chr, pattern='GATC')
>   start <- c(1, start(GATC.match)) # not sure what happens if
> chromosome starts with GATC
>   end <- c(end(GATC.match), length(chr))
>   Views(subject=unmasked(chr),start=start, end=end)
> }
> pm <- new('BSParams', X=Hsapiens, FUN=function(chr)
> ExtractGATCFragments(chr) )
> # create the GenomeData object:
> GATC.fragments <- bsapply(pm)
>
> # is it a genomeData??
> class(GATC.fragments)
> # [1] "GenomeData"
> # attr(,"package")
> # [1] "BSgenome"
>
> # coerce it to a RangedData:
> as(GATC.fragments, "RangedData")
> # Error in FUN(X[[1L]], ...) :
> #   no method or default for coercing “XStringViews” to “RangedData”
>
>
>
> > sessionInfo()
> R version 3.0.2 (2013-09-25)
> Platform: x86_64-unknown-linux-gnu (64-bit)
>
> locale:
>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
>  [3] LC_TIME=nl_NL.UTF-8        LC_COLLATE=en_US.UTF-8
>  [5] LC_MONETARY=nl_NL.UTF-8    LC_MESSAGES=en_US.UTF-8
>  [7] LC_PAPER=nl_NL.UTF-8       LC_NAME=C
>  [9] LC_ADDRESS=C               LC_TELEPHONE=C
> [11] LC_MEASUREMENT=nl_NL.UTF-8 LC_IDENTIFICATION=C
>
> attached base packages:
> [1] parallel  stats     graphics  grDevices utils     datasets  methods
> [8] base
>
> other attached packages:
> [1] BSgenome.Hsapiens.UCSC.hg19_1.3.19 BSgenome_1.30.0
> [3] Biostrings_2.30.1                  GenomicRanges_1.14.4
> [5] XVector_0.2.0                      IRanges_1.20.6
> [7] BiocGenerics_0.8.0
>
> loaded via a namespace (and not attached):
> [1] stats4_3.0.2 tools_3.0.2
> ######################################################################
>
> Can anybody point out to me why the object can not be coerced.
>
> Thanks, Ludo
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor@r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives:
> http://news.gmane.org/gmane.science.biology.informatics.conductor

	[[alternative HTML version deleted]]

