Hello,
Triggered by the paper by Andrew Jaffe et al. on bump hunting to identify DMRs, I tried to analyze one of our Illumina 450K datasets using the bump hunting methodology implemented in the package charm (dmrFinder).
I found the lines of codes posted by Tim Triche on the BioC mailing list an excellent way to get started (http://comments.gmane.org/gmane.science.biology.informatics.conductor/44365). However, I now do have some practical questions/problems and would appreciate some feedback/pointers.

- If SWAN-normalized data is used as input for dmrFinder, is within and between sample normalization still required? I assume it is not, but if still so, what would be the recommended method to get started for these 450K arrays (loess resp. quantile or sqn)?
- Using default settings, dmrFinder didn't return any significant dmr's in my dataset, except when the minimum number of probes in a region was set at 1 (but this is equal to analyzing the CpGs one-by-one and not bump hunting). I was just wondering; did other people using dmrFinder on 450K arrays observe this as well? Or is this specific for my dataset because the effect size in my study is too small?
- What exactly is the value of the cutoff used to identify a dmr reflecting? The help pages state that it is a t-statistic cutoff used to identify probes as being in a DMR, so am I correct it is equal to 1-p.value (since the default cutoff is 0.995)?
- Finally, I am not able to produce some nice plots but don't understand what is going wrong. Anyone else an idea?

Thanks,
Guido


> sample.data
class: GenomicMethylSet
dim: 485512 20
exptData(0):
assays(2): Meth Unmeth
rownames(485512): cg13869341 cg14008030 ... cg08265308 cg14273923
rowData metadata column names(0):
colnames(20): 7766130001_R01C01 7766130001_R02C01 ... 7766130037_R04C02 7766130037_R05C02
colData names(13): Sample_ID Sample_Plate ... Basename filenames
Annotation
  array: IlluminaHumanMethylation450k
  annotation: ilmn.v1.2
Preprocessing
  Method: SWAN (based on a MethylSet preprocesses as 'Raw (no normalization or bg correction)'
  minfi version: 1.5.2
  Manifest version: 0.4.0
>


# When minimum number of probes per drm set at 2.
> results <- dmrFinder(eset=NULL, groups=groups, p=(assays(sample.data, withDimnames=F)$Meth)/(assays(sample.data, withDimnames=F)$Meth+assays(sample.data, withDimnames=F)$Unmeth),chr=as.character(seqnames(sample.data)),
pos=start(sample.data), pns=rownames(sample.data), withinSampleNorm="none", betweenSampleNorm="none", removeIf=expression(nprobes<2), package='IlluminaHumanMethylation')
Computing group medians and SDs for 2 groups:

1

2

Done.

Smoothing:

====================================================================================================================================================================
Done.
Finding DMRs for each pairwise comparison.

old-young
.....................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................0 DMR candidates found using cutoff=0.995.

Done
> head(results$tabs[[1]])
[1] chr        start      end        p1         p2         regionName indexStart indexEnd   nprobes    area       ttarea     diff       maxdiff
<0 rows> (or 0-length row.names)
>
> cpg.cur = read.delim("http://rafalab.jhsph.edu/CGI/model-based-cpg-islands-hg19.txt",as.is=TRUE)
>
> dmrPlot(dmr=results, which.table=1, which.plot=1:5, legend.size=1,
+ all.lines=FALSE, all.points=FALSE, colors.l=c("blue","red"),
+ colors.p=c("blue","black"), outpath=".", cpg.islands=cpg.cur, Genome=Hsapiens)
Smoothing:
====================================================================================================================================================================
Done.
Plotting regions in table 1
Plotting region 1
Error in plot.window(...) : need finite 'xlim' values
In addition: Warning messages:
1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf
3: In min(pos[Index]) : no non-missing arguments to min; returning Inf
4: In max(pos[Index]) : no non-missing arguments to max; returning -Inf



### when minimum number of probes set at 1
> results <- dmrFinder(eset=NULL, groups=groups, p=(assays(sample.data, withDimnames=F)$Meth)/(assays(sample.data, withDimnames=F)$Meth+assays(sample.data, withDimnames=F)$Unmeth),chr=as.character(seqnames(sample.data)), pos=start(sample.data), pns=rownames(sample.data), withinSampleNorm="none", betweenSampleNorm="none", removeIf=expression(nprobes<1), package='IlluminaHumanMethylation') ##
Computing group medians and SDs for 2 groups:

1

2

Done.

Smoothing:

====================================================================================================================================================================
Done.
Finding DMRs for each pairwise comparison.

old-young
.....................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................7482 DMR candidates found using cutoff=0.995.

Done
>
> head(results$tabs[[1]])
         chr     start       end        p1        p2 regionName indexStart indexEnd nprobes      area   ttarea       diff    maxdiff
159352 chr14  95046686  95046686 0.5991461 0.9343884 cg08210706     147284   147284       1 0.3352423 78.42597 -0.3352423 -0.3352423
349254 chr12 113516445 113516445 0.9081579 0.5608301 cg19353052     117535   117535       1 0.3473278 75.86928  0.3473278  0.3473278
183581 chr12   4699232   4699232 0.4798595 0.8823826 cg09581911     101344   101344       1 0.4025231 65.85157 -0.4025231 -0.4025231
459811  chr1  22425642  22425642 0.4948227 0.9094913 cg26348696      10893    10893       1 0.4146686 64.13913 -0.4146686 -0.4146686
385470  chr7   1097952   1097952 0.5278402 0.9391450 cg21655171     414247   414247       1 0.4113047 61.48937 -0.4113047 -0.4113047
217517 chr19  18273012  18273012 0.9377911 0.5412174 cg11594160     233921   233921       1 0.3965737 57.28228  0.3965737  0.3965737
>
> dmrPlot(dmr=results, which.table=1, which.plot=1:5, legend.size=1,
+ all.lines=FALSE, all.points=FALSE, colors.l=c("blue","red"),
+ colors.p=c("blue","black"), outpath=".", cpg.islands=cpg.cur, Genome=Hsapiens)
Smoothing:
====================================================================================================================================================================
Done.
Plotting regions in table 1
Plotting region 1
Error in matplot(pos[Index], M[Index, showpts], col = tcols, cex = 0.6,  :
  'x' and 'y' must have same number of rows
>
>
> sessionInfo()
R Under development (unstable) (2012-11-21 r61136)
Platform: x86_64-w64-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252    LC_MONETARY=English_United States.1252 LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] IlluminaHumanMethylation450kannotation.ilmn.v1.2_0.1.3 IlluminaHumanMethylation450kmanifest_0.4.0
 [3] minfi_1.5.2                                            Biostrings_2.27.7
 [5] GenomicRanges_1.11.15                                  IRanges_1.17.15
 [7] reshape_0.8.4                                          plyr_1.7.1
 [9] lattice_0.20-10                                        charm_2.5.6
[11] genefilter_1.41.0                                      RColorBrewer_1.0-5
[13] fields_6.7                                             spam_0.29-2
[15] SQN_1.0.5                                              nor1mix_1.1-3
[17] mclust_4.0                                             Biobase_2.19.0
[19] BiocGenerics_0.5.1

loaded via a namespace (and not attached):
[1] affxparser_1.31.1     affyio_1.27.1         annotate_1.37.2       AnnotationDbi_1.21.7  beanplot_1.1          BiocInstaller_1.9.4   bit_1.1-9
 [8] BSgenome_1.27.1       codetools_0.2-8       corpcor_1.6.4         DBI_0.2-5             ff_2.2-10             foreach_1.4.0         grid_2.16.0
[15] gtools_2.7.0          illuminaio_0.1.3      iterators_1.0.6       limma_3.15.4          MASS_7.3-22           matrixStats_0.6.2     multtest_2.15.0
[22] oligo_1.23.0          oligoClasses_1.21.5   preprocessCore_1.21.0 R.methodsS3_1.4.2     RSQLite_0.11.2        siggenes_1.33.0       splines_2.16.0
[29] stats4_2.16.0         survival_2.36-14      sva_3.5.0             tools_2.16.0          XML_3.95-0.1          xtable_1.7-0          zlibbioc_1.5.0
>


	[[alternative HTML version deleted]]

