First off, my apologies about the BSgenome requirement.  I realized
(...slowly...) that rtracklayer requires the appropriate build of BSgenome
to correctly assign the chromosome lengths when requested.  It's either
that or assume the user has connectivity to UCSC, or store them statically
(which seems like a bad idea, but perhaps it isn't).  I attempted to upload
an updated version of the package but due to missed communications with
Marc, it never made it onto the BioC repository.  Mea culpa.  That said...

The following is for a bunch of hits from an experiment, but you could just
swap in something like

if(!require(BSgenome.Hsapiens.UCSC.hg19))
biocLite('BSgenome.Hsapiens.UCSC.hg19')
if(!require(FDb.InfiniumMethylation.hg19))
biocLite('FDb.InfiniumMethylation.hg19')

library(FDb.InfiniumMethylation.hg19)
hm450 <- get450k()
hm450.CpGs <- split(hm450, hm450$probeType)$cg

strand(hm450.CpGs) <- '*' ## always worth double checking



and use hm450.CpGs in the following code where I am using varByAgeProbes...


library(Homo.sapiens)
txsByGene <- transcriptsBy(TxDb.Hsapiens.UCSC.hg19.knownGene, 'gene')
names(txsByGene) <- mget(names(txsByGene),
                               org.Hs.egSYMBOL, ifnotfound=NA)
genicHits <- subsetByOverlaps(varByAgeProbes, txsByGene)

genicOLs <- findOverlaps(varByAgeProbes, txsByGene)
genesHit <- sort(table(names(txsByGene)[subjectHits(genicOLs)]))

promotersByGene <- flank(txsByGene, 1500) ## adjust the number to suit your
definitions

promoterHits <- subsetByOverlaps(varByAgeProbes, promotersByGene)

promoterOLs <- findOverlaps(varByAgeProbes, promotersByGene)
promotersHit <-
sort(table(names(promotersByGene)[subjectHits(promoterOLs)]))

promoterGenicHits <- intersect( promoterHits, genicHits )  ## probes that
hit both regions

promotersAndGenes <- reduce(unlist(c(reduce(promotersByGene),
reduce(txsByGene))))

intergenicHits <- varByAgeProbes[ which(!varByAgeProbes
%over% promotersAndGenes) ]

## trust but verify:
length(intersect(intergenicHits, genicHits))
## [1] 0
length(intersect(intergenicHits, promoterHits))
## [1] 0

So what, right?  I mean, you could have figured this out from the
manifest...maybe (i.e., assuming you didn't use the GEO manifest with its
Excel-gene names, and assuming there are no recently updated or rare
transcripts of interest in your data).  But suppose you have a bunch of
ChIP-seq or DNAse-seq results, or footprints of transcription factor
binding sites, and you want to see which of those are hit.  Or H3K27me3
peak calls, or CTCF peaks.  Let's use the latter.  (There's a very specific
reason that I am looking at CTCF locations in this particular analysis)

require(rtracklayer)
CTCF.peaks <- import("CTCF_peaks.bed", genome='hg19')

subsetByOverlaps( intergenicHits, CTCF.peaks )

You *could* do the above with just the manifest, but it would be a colossal
pain in the... posterior.  Now consider if you were looking at
DNAse-inaccessible PU.1 sites, or HiC interacting domains, or the like.

Anyways, I'll patch up the issues with FDb.InfiniumMethylation.hg18/19, try
to make it a bit faster, and re-upload the packages over the weekend (or
before, but I worry that can't be done).  I hope the above clarifies my
intent with the package.  An alternative is to use minfi's mapToGenome()
function, FWIW.





On Tue, Jul 2, 2013 at 2:14 AM, Simone <enomis.bioc@gmail.com> wrote:

> Dear Tim,
>
> > You will need to use toggleProbes() to get all probe locations (this is
> an
> > undesirable side effect of the db0 package infrastructure, which is
> geared
> > towards expression arrays where probes that map to multiple accessions
> are a
> > "bad thing", and that is the reason I deprecated the 27k.db and 450k.db
> > packages).
>
> I understand. However, also with toggleProbes() I have got a problem.
> I tried for example:
>
> > x <- toggleProbes(IlluminaHumanMethylation450kPATH2PROBE, "all")
>
> Error in (function (classes, fdef, mtable)  :
>   unable to find an inherited method for function ‘toggleProbes’ for
> signature ‘"AnnDbBimap"’
>
> But for example the following works
>
> > x <- toggleProbes(IlluminaHumanMethylation450kENSEMBL2PROBE, "all")
>
> However, I would need the location information, so the column
> "transcript_location" (e.g. 5'UTR, TSS200, Body, etc.) and it seems
> that I cannot get it like this?
>
> > The FDb.InfiniumMethylation packages have the locations of the targeted
> CpGs
> > (and CpHs, and SNPs) for the corresponding build (just for the record,
> the
> > genome build targeted by 450k.db is hg19/GRch37), but they are returned
> > using get450k() or get27k() as a GRanges object.
>
> Currently, this doesn't work for me neither (see my previous post).
> Probably I'm doing something wrong ...
>
> > Essentially, the mistake I made in compiling the 27k.db and 450k.db
> packages
> > was to use an expression-centric (or gene- and transcript-centric)
> > infrastructure to represent probes targeting individual loci.  It's more
> > useful to think of the 450k array, in particular, as a SNP array, with
> > functions for a given target locus inferred after discovery.  It's also
> more
> > useful to think of the chip as a SNP array when considering
> differentially
> > methylated regions, which could be loosely regarded as the methylation
> > equivalent of copy number variations.  Last but not least, you can in
> fact
> > use the array to estimate copy number along the genome, as for
> determining
> > sex from chrX copy number.
>
> Thank you for your hints and explanations!
> But at the moment I am still stucked at a much simpler step, I would
> just need an easy way to distinguish probes located in promoters from
> those located in gene bodies (and to discard those which are annotated
> to both).
>
> Best wishes,
> Simone
>



-- 
*A model is a lie that helps you see the truth.*
*
*
Howard Skipper<http://cancerres.aacrjournals.org/content/31/9/1173.full.pdf>

	[[alternative HTML version deleted]]

