

There are a couple of how too's about this but briefly it

Say you go a bed file of reads  of "whatever you want to get the
coverage of" ... the bed file is a dataframe with start, end ,and
chromosome info:
Otherwise you can use ShortRead to get the RangeData object of reads...
whatever you a a rangeList or rangeData object of what you are
counting..


Turn the bed file, (as a dataframe) into a Ranges Data object: You can
use a RangeList if you want ...


rl<-RangedData(IRanges(start=bed[,"start"],end=bed[,"end"]),space=bed[,"chr"])
________________________________________

So this has chromosomes in the space section..

You the want coverage on a set of "other regions".... say the exons,
chromosomes, Cpg islands ... whatever:

Eg for exons:
library(GenomicFeatures.Mmusculus.UCSC.mm9)
library(org.Mm.eg.db)

genes<-geneMouse()
the.exons<-exons_deprecated(genes)

"the.exons" is not a ranged data object: (a RangeList can also be used)
of UCSC exons or name you own..

the.exons<- RangesList(
chr10 =
IRanges(start=starts[starts[,"chrom"]=="chr10","txStart"],end=starts[starts[,"chrom"]=="chr10","txStart"]),
chr11  =
IRanges(start=starts[starts[,"chrom"]=="chr11","txStart"],end=starts[starts[,"chrom"]=="chr11","txStart"]),
chr12  =
IRanges(start=starts[starts[,"chrom"]=="chr12","txStart"],end=starts[starts[,"chrom"]=="chr12","txStart"]),
....)

so in summary the.exons is a rangeData or Rangelist objuct of "the
regions you want the coverage on" (al las  start and end of 
lapply(ranges(RangedData),coverage, start = 1, end = 1000) 


________________________________________-

## do coverage calculation

library(BSgenome.Mmusculus.UCSC.mm9)
mouse.chromlens<-seqlengths(Mmusculus)

Get a rleList objust of the coverage over the genome:

cov<-coverage(rl,width=mouse.chromlens) ## list names are the chromosome
names/spcae names in rl...you don't have to use "width"
##

cov<-cov[names(the.exons)] # make sure list in the same order can this not checked in the  line below please check the chromosome names are consistentin rl and the.exons

you now have the coverge on all the chromosomes.... but for other specified regions use:


regionViews <- RleViewsList(rleList = cov, rangesList = ranges(the.exons))  ###


## your done!
____________________________________________

##dumP the reults as a dataframe::
the.counts<-{}
order.chromos<-names(the.exons)
for (i in 1:length(order.chromos)){
  chromo<-order.chromos[i]
  the.counts<-rbind(the.counts,cbind(chromo,start(the.exons[chromo]),end(the.exons[chromo]), viewSums(regionViews[[chromo]]) ))
}
colnames(the.counts) <- c("chr","start","end","count")

Cheers
Paul





-----Original Message-----
From: Jane Landolin <landolin@gmail.com>
To: bioc-sig-sequencing@r-project.org
Subject: [Bioc-sig-seq] Coverage for RangeList?
Date: Mon, 19 Apr 2010 16:40:54 -0700


Hi all,

I just joined the list, and was hoping to get some advice about the
best way to work with the various classes that deal with genomic range
info. I am a postdoc in the Celniker Lab at Lawrence Berkeley National
Lab, working on the modENCODE project.

So my question has to do with using various functions that operate on
the IRanges level, while keeping my data in RangedData or RangeList
objects.  I should be able to use rdapply or lapply, but it is not
very clear to me when I should use one over another.

If I want to calculate the chromosome-by-chromosome coverage for a
RangeData object, I can do something like:

lapply(ranges(RangedData),coverage, start = 1, end = 1000)

but how can I do it for varying starts and ends?

Alternatively, is it better to write a function and call it using
rdapply?  What have you guys been doing for this scenario?


Thanks in advance,
Jane Landolin

_______________________________________________
Bioc-sig-sequencing mailing list
Bioc-sig-sequencing@r-project.org
https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing

	[[alternative HTML version deleted]]

