I'm sure this is simple but I haven't managed to do it and was wondering if someone could give me a simple version. I have a list of about 10,000 150bp genome co-ordinates in a bed file- which I can import into granges object if necessary. I would like to see which of a variety of tracks supplied in the UCSC browser overlap with my co-ordinates eg CpG islands, exons and other things. I can do this via the laborious download of individual tracks as a bed and then use bedtools closest but I was sure I could do this directly from the browser.

So I started to try to use rtracklayer package but it isn't as straightforward as I thought it might be.

So far I've got:

library(GenomicRanges)
library(rtracklayer)
library(Rsamtools)

BEDfile = "Mybedfile.bed"

#Load the annotation and reduce it
BED <- import.gff(BEDfile, format="BED", genome="hg19", asRangedData=F)

session <- browserSession("UCSC")

track(session, "BED") <- BED

......but what do I do from here to get info from the UCSC tracks that overlap my co-ordinates? Do I have to download the annotating track Im interested in and then do findOverlaps with my track?