On Tue, Oct 26, 2010 at 5:02 PM, Sheena Scroggins <
sheena.scroggins@gmail.com> wrote:

> Hello,
>
> I'm trying to find a way to use rtracklayer to give me an output of the
> conservation scores for a genomic region. For example: I want to be able to
> define a range of data like the RangedData does:
>
> >  strand <- c("+", "+", "+")
> >  start <- c(17563082, 21920042, 61191807)
> >  end <- c(17605866, 21956374, 61290544)
> > rdTable <- RangedData(ranges = IRanges(start = start, end = end), strand
> =
> strand, space = chr)
> > rdTable
> RangedData with 3 rows and 1 value column across 2 spaces
>        space               ranges |      strand
>  <character>            <IRanges> | <character>
> 1        chr1 [17563082, 17605866] |           +
> 2        chr1 [61191807, 61290544] |           +
> 3        chr2 [21920042, 21956374] |           +
>
> But then combine it with rtracklayer so that I can find the conservation
> scores of the data ranges:
> >  s2 = browserSession("UCSC")
> > ct = track(s2, "cons46way")
> > ct
> UCSC track 'Primate Cons'
> UCSCData with 9974 rows and 1 value column across 1 space
>           space               ranges   |     score
>     <character>            <IRanges>   | <numeric>
> 1          chr21 [33031597, 33031597]   |  0.560087
> 2          chr21 [33031598, 33031598]   |  0.560087
> 3          chr21 [33031599, 33031599]   |  0.435717
> 4          chr21 [33031600, 33031600]   |  0.435717
> 5          chr21 [33031601, 33031601]   |  0.560087
> 6          chr21 [33031602, 33031602]   |  0.560087
> 7          chr21 [33031603, 33031603]   |  0.435717
> 8          chr21 [33031604, 33031604]   |  0.435717
> 9          chr21 [33031605, 33031605]   |  0.560087
> ...          ...                  ... ...       ...
> 9966       chr21 [33041562, 33041562]   | -0.266457
> 9967       chr21 [33041563, 33041563]   |  0.433850
> 9968       chr21 [33041564, 33041564]   |  0.507567
> 9969       chr21 [33041565, 33041565]   |  0.655000
> 9970       chr21 [33041566, 33041566]   | -0.155882
> 9971       chr21 [33041567, 33041567]   | -0.340173
> 9972       chr21 [33041568, 33041568]   |  0.507567
> 9973       chr21 [33041569, 33041569]   | -1.740790
> 9974       chr21 [33041570, 33041570]   | -1.925080
> >
>
> So that the output will be the different regions in my RangedData with
> their
> respective conservation scores.


Currently, rtracklayer uses the UCSC Table Browser via the web, which is
generally too limited and inefficient for genome-wide queries. It is
feasible to obtain the conservation scores for a small number of queries:

query <- UCSCTableQuery(s2, "cons46way")
range(query) <- ranges(rdTable[1,])
cons <- track(query)

And repeat the last two lines for each range.

In the long run, the best method is to download the conservation scores
track from their FTP site.   You could just mget the WIG file for each
chromosome:
http://hgdownload.cse.ucsc.edu/goldenPath/hg19/phastCons46way/vertebrate/.

Then use the UCSC bigwig utilities to convert the merged data into a bigwig
file. This is essentially an indexed database of the conservation scores,
and you can query it with rtracklayer using code like:

cons <- import("cons.bw", selection = ranges(rdTable))

That will give you a RangedData with a column containing the conservation
scores for those regions in rdTable.



> It would be REALLY great if you could help
> me also find a way to put a threshold on the conservation scores, or a way
> to filter them afterwards.
>
>
Filtering afterwards is straight-forward, e.g.:

filteredCons <- subset(cons, score > 2)

Michael

Thanks!
> Sheena
>
>        [[alternative HTML version deleted]]
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor@stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives:
> http://news.gmane.org/gmane.science.biology.informatics.conductor
>

	[[alternative HTML version deleted]]

