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