On Fri, Oct 9, 2009 at 4:55 PM, Dykema, Karl <[email protected]> wrote:
> Hello All,
>
> I am attempting to calculate the mean number of reads for every exon that
> was represented on our custom Nimblegen capture chip. The capture regions
> don't line up exactly with the exon coordinates so I identified all exons
> that have an overlap, about 9500. I have previously done some crude sampling
> of the pileup output and my results were a lot better than what I seem to be
> calculating now. Currently I am running coverage() on my alignedRead object
> and then looping through the capture regions something like this:
>
> curCov <- coverage(lane1)
> for(q in c(1:22,"X","Y")){
> capturedExonsOnChrom <-
> capturedExons[which(capturedExons$chrom==q),]
> for(i in 1:length(capturedExonsOnChrom)) {
> exonCoordinates <-
> capturedExonsOnChrom[i]$exonStart:capturedExonsOnChrom[i]$exonEnd
> capturedExonMeanCoverage[] <-
> mean(as.vector(curCov[[q]])[exonCoordinates])
> }
> }
>
> First off, I am assuming that there is a more efficient method to do this.
> Is that the case? Also, I have also been using GenomeGraphs to create plots
> of the coverage and those plots also seem to show plenty of reads where my
> calculations above indicate none. I am hoping there is an error in my
> current method because the previous coverage numbers were quite encouraging.
> Thanks in advance for your help.
>
>
Hi, Karl.
I put together a little vignette that includes, among other things, a little
vignette that includes a short example of calculation of given a capture
method (Agilent, but the idea is the same). The reads are read in using the
ShortRead package and the capture regions are in the form of a simple bed
file. The vignette is here:
http://watson.nci.nih.gov/~sdavis/sequencing_example.pdf
Feel free to take a look at it for a slightly more efficient way of doing
things. I'm not guaranteeing total correctness in the code and not all
parts will apply to nimblegen, I suppose, but it might get you started.
Sean
[[alternative HTML version deleted]]
_______________________________________________
Bioc-sig-sequencing mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing
