I'm not really sure why Kasper rearranged
the IlluminaHumanMethylation450kmanifest package; it's more of a pain to
use now.  In the first release of minfi, there was a nice little example of
how to cast the source sequences, strand, and interrogated site from the
HM450k manifest into a GRanges object and go from there.  Perhaps another
project stole all of his time.  Anyways, the following should be useful for
people using either 'minfi' or 'methylumi'; I would suggest that, for
visualization and annotation purposes, a SummarizedExperiment is the
logical terminus for processed DNA methylation array data, as it segues
more naturally into BSseq, ChIPseq, and RNAseq integration.

In order to level the playing ground, I will construct a Ranges object
using the annotations at GEO which the TCGA ADF now references.

Since I'm refreshing both the 450k.db and probe packages (the latter has a
bug in the example script; don't use it), I will put something like the
following into the documentation.  The bottleneck for refreshing the .db
package is the slowness of build checks with db0-type chip-style schemas
and my slow learning of how to use a FeatureDb instead.  I strongly believe
that GRanges are the appropriate data structure for this type of annotation
data.  Personally I've been combining 450k results with RNAseq and ChIPseq
data to look at some (IMHO) rather interesting phenomena, and having the
flexibility of a GRanges-based object (a SummarizedExperiment, to be
specific) is enormously useful.  On the other hand, I haven't found that
segmentation-based approaches seem to be very useful for the 450k data, as
some fine-scale changes get lost when you compare to BS-seq data in similar
cell types.  In my experience, it makes more sense to define genomic
regions and sweep up changes in those, rather than segmenting DNA
methylation, but that's just one data point.  So, here goes.

require(GEOquery)
require(GenomicRanges)
require(GenomicFeatures)
require(BSgenome.Hsapiens.UCSC.hg19)

hm450k <- getGEO('GPL13534')

## Here's a gotcha that you might like to avoid:
## (one of the reasons I'm bothering with a 450k.db update is that I have
resolved it)
##
anyMissing(as.numeric(hm450k@dataTable@table$MAPINFO))
## TRUE
##
## exclude these (SNP probes)
noMap <- which(is.na(as.numeric(hm450k@dataTable@table$MAPINFO)))
hm450 <- with( hm450k@dataTable@table[ -noMap, ],
                GRanges( paste0('chr', CHR),
                         IRanges(as.numeric(MAPINFO), width=1),
                         probe=Name,
                         probeStrand=Strand,
                         sourceSeq=DNAStringSet(SourceSeq),
                        )
              )
names(hm450) <- values(hm450)$probe
values(hm450)$probeType = Rle(ifelse(substr(names(hm450),1,2)=='cg', 'CpG',
'CpH'))

Now we have something resembling annotation (but not really, for reasons
I'll demonstrate later).

## inspect it:
head(hm450, 1)
GRanges with 1 range and 4 elementMetadata cols:
             seqnames             ranges strand |       probe probeStrand
                                       sourceSeq probeType
                <Rle>          <IRanges>  <Rle> | <character> <character>
                                  <DNAStringSet>     <Rle>
  cg00035864     chrY [8553009, 8553009]      * |  cg00035864           F
AGACACTAGCAGTCTTGTCCACATAGACCCTTGAATTTATCTCAAATTCG       CpG

I propose that you always sanity check your work at each step.  Here's one
way I typically do it:

hm450.subset <- hm450[ sample(names(hm450), 10) ]
values(hm450.subset)$siteSeq <- getSeq(Hsapiens, resize(hm450.subset, 2))
show(hm450.subset)
## GRanges with 10 ranges and 5 elementMetadata cols:
##               seqnames                 ranges strand |        probe
probeStrand                                          ## sourceSeq probeType
       siteSeq
##                  <Rle>              <IRanges>  <Rle> |  <character>
<character>                                     ## <DNAStringSet>     <Rle>
<DNAStringSet>
##   ch.X.665616F     chrX [ 41078352,  41078352]      * | ch.X.665616F
      F CAGGCTTATGACAACTTAAGCTTGAGTGATCACTTACTAAGAGCAGTACT       CpH
      CA
##     cg19794706    chr14 [ 75643232,  75643232]      * |   cg19794706
      R GGAGATGGCAAGGACCAATCTGGGGCCGAGCAGGAACAAAAGCAGCAACG       CpG
      CG
##     cg15391651     chr1 [165132273, 165132273]      * |   cg15391651
      F CGGAGCCCTGAGTGTGCACAAAGCACCACTATGCCAGAGTGATGTTATCA       CpG
      CG


Suppose you wanted to mask probes whose interrogated site overlaps with
repeats in the genome.  This is easy:

library(parallel)
REPEATS <- do.call(c, mclapply(seqlevels(HM450.HG19), function(ch) {
  GRanges(ch, union(masks(Hsapiens[[ch]])$RM, masks(Hsapiens[[ch]])$TRF))
}))
rpts.450 <- subsetByOverlaps(hm450, REPEATS)
summary(rpts.450)
## Length   Class    Mode
##  78164 GRanges      S4

The example gets to the point where the Ranges/FeatureDb-based version is;
I'll post another in a while, using the distanceToNearest functionality to
see how far a given probe is from a TSS, either UCSC or FANTOM-annotated,
and assorted other things that use a TxDb object to compute quickly. I am
never going to put up a PDF with a "protocol" (the #$%@?!) that I'm not
1000% sure of ever again...

I strongly suggest that you take whatever container your data is in, and
turn it into a SummarizedExperiment with the hm450 object as the rowData.
 I believe you will thank me later on, to the point that I've written a
generic to do this for MethyLumiSet objects, and could do the same for
MethyLumiM and MethylSet objects if there is demand.  I would like to patch
Gviz to handle SummarizedExperiments soon, too.

One thing I have not done here is to mask SNPs, because without MAF
information (or, better, SNP information on a subject), I feel like it may
be swatting a fly with a Buick to mask all CpG/CpH overlaps with dbSNP
regardless of minor allele frequency.  You can use the Probe_SNPs
and Probe_SNPs_10 columns in the manifest if you want to mask in that
fashion against dbSNP build 131, but there are other things that I find
wanting in that approach (i.e. you should widen the interrogated sequences
appropriately).  Anyways...

Hope this helps.  Basically, if you want to look at the probes on the 450k
platform which fall into a given feature or class of features, you assemble
a GRanges object matching those data (perhaps ChIPseq peaks, perhaps UCSC
tracks, perhaps FeatureDb output) and use subsetByOverlaps with the
appropriate arguments to isolate the desired probes in those regions.  The
trick is finding feature data.

--t




On Fri, May 4, 2012 at 7:15 AM, Ed Siefker <ebs15242@gmail.com> wrote:

> I've run some arrays through minfi and have a table of differentially
> methylated positions.  I can get the probe sequence from
> IlluminaHumanMethylation450kmanifest, but it's not clear to me how to
> translate that into genes, promoters, chromosomes or anything
> biologically informative.
>
> What do people usually do at this point?
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor@r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives:
> http://news.gmane.org/gmane.science.biology.informatics.conductor
>



-- 
*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]]

