[Bioc-sig-seq] rtracklayer and coverage() from GenomeData

[Chris Seidel]

Hello,

I'm getting an error using export() to generate a wig file from a
coverage result.

In an earlier thread (Subject: Re: [Bioc-sig-seq] WIG from coverage, Sep
30, 2009 at 3:15 PM) there was an example of going from reads to
coverage to wig, and Michael Lawrence had suggested:

> The most convenient path would be to perform the analysis over all
> chromosomes and lanes at once. The extendReads() function should do
this for
> you, e.g:
>
> extendedReads <- extendReads(ypd, seqLen=150)
>
> The 'extendedReads' above is a GenomeDataList. You could then do something
> like:
>
> cov <- gdapply(extendedReads, coverage)
>
> Now 'cov' is another GenomeDataList. Just subset out one lane and export
> it:
>
> export(cov$wce, "wce.wig")

When I try this, I get an error:

# read in data and convert to genome data objects
aln <- readAligned(spath, pattern="s_1_export.txt", type="SolexaExport",
filter=filt)
ip <- as(aln,"GenomeData")
aln <- readAligned(spath, pattern="s_2_export.txt", type="SolexaExport",
filter=filt)
wce <- as(aln,"GenomeData")

# create GenomeDataList from ip and wce
cb <- GenomeDataList(list(ip=ip, wce=wce))

# extend reads
cb.ext <- extendReads(cb, seqLen=50)

# get coverage for both lanes across all chromosomes
cb.cov <- gdapply(cb.ext, coverage)

# export the ip lane to wig
export(cb.cov$ip, "cb_ip.wig")

Error in object[order(start(object)), ] :
    cannot combine 'i' row values from different spaces

However, it works correctly if I "export" one chromosome at a time:

export(cb.cov$ip$chr2L, "cb_ip.wig")

Is this a bug, or am I doing something wrong? Can export take a
subsetted GenomeDataList object? I even tried coercing the GenomeData
element to RangedData, but I got the same error while trying to export.

-Chris Seidel

> sessionInfo()
R version 2.10.0 Under development (unstable) (2009-09-03 r49555)
x86 <- 64-unknown-linux-gnu

other attached packages:
  [1] rtracklayer <- 1.5.17 RCurl <- 1.2-0        bitops <- 1.0-4.1    
chipseq <- 0.1.27
[5] ShortRead <- 1.3.36   lattice <- 0.17-25    BSgenome <- 1.13.14  
Biostrings <- 2.13.47
[9] IRanges <- 1.3.87

loaded via a namespace (and not attached):
  [1] Biobase <- 2.5.8 XML <- 2.6-0     grid <- 2.10.0   hwriter <- 1.1
  tools <- 2.10.0