On Tue, Jun 30, 2009 at 8:32 AM, Simon Anders <anders@ebi.ac.uk> wrote:

> Hi Michael et al.
>
> Assume I have a list of AlignedRead objects, with data from several Solexa
> lanes. I would like to get a coverage vector over all the lanes.
>
> As 'coverage' takes only one AlignedRead object, I have two possibilities:
>
> (a) Concatenate the AlignedRead objects to a single big one. As far as I
> can see, 'rbind2' is not implemented for AlignedRead, and 'append' seemed
> very slow to me. It is probably not posible to do this without making a copy
> of all the data in memory.
>

I think this is what combineLanes() in the chipseq package does, and it
seems to work fine for us. Note that it takes a GenomeDataList, rather than
a list of AlignedReads, but ShortRead makes it easy to get a GenomeDataList,
if I remember right.

Of course, optimization is welcome.


> (b) Calculate the coverage for each AlignedRead object separately and add
> up the GenomeData objects. The `+` operator is not supported for GenomeData
> objects but it is for Rle objects. So I need to loop through the
> chromosomes.
>

> I've now written this short function for the purpose:
>
> sumUpCoverage <- function( lanes, seqLens )
> {
>   res <- NULL
>   for( i in 1:length(lanes) ) {
>      cvg <- coverage( lanes, width=seqLens )
>      if( is.null( res ) )
>         res <- cvg
>      else {
>         stopifnot( all( names(res) == names(cvg) ) )
>         for( seq in names(res) )
>            res[[seq]] <- res[[seq]] + cvg[[seq]]
>      }
>   }
>   res
> }
>
> This does not seem very elegant (and it takes 9 minutes, which, however,
> might be ok for 29 mio reads). Any idea how to do it better? (The use of
> 'for' instead of 'sapply' is on purpose: I hope to save memory that way.)
>
> And would it make sense to overload the `+` operator for GenomeData
> objects? If so, could I suggest adding this to BSgenome?
>
> Cheers
>  Simon
>
>
>

	[[alternative HTML version deleted]]

