Hi, I am using minfi to analyze an Illumina 450K dataset, but due to an experimental design that comprises more than 2 grps I am a lost on how to properly go for 'bumphunting'... In addition, I do have 3 queries regarding minfi's workflow. See below. Any suggestions would be appreciated! Thanks, Guido Design of my study: parallel design, before + after measurement for each subject (baseline vs treatment, = paired), young and aged subjects. We included blood cells from 10 subjects in our study. Subjects were divided into 2 grps, "young" (5x) and "old" (5x). At baseline (before treatment) samles were taken (10x), and after 48hr of treatment (10x). In total thus 20 samples were analysed on the 450k array. Each treatment sample can be compared with it's respective baseline (paired analysis). Of interest are: 1) difference between young vs old at baseline (unpaired analysis, 10 vs 10), 2) effect of treatment in young (paired analysis, 5 vs 5), 3) effect of treatment in old (paired analysis, 5 vs 5), 4) difference in response due to age, thus CpGs that respond differently in 3 vs 2 ["delta of the delta's"]. To identify the individual regulated CpGs for each comparison, I used limma (on the M-values) using the multilevel model (thus taking into account the paired structure using the duplicateCorrelation() function). However, considering the design of my experiment, how to optimally analyse for enriched regions using the bumphunter() function? Since I have multiple groups + paired structure, the design matrix doesn't contain a column (coefficient) that corresponds to each of the contrasts of interest (1-4). This is also returned as warning by the bumphunter() function, and referred is to the vignette [which one I don't know, but I guess it's from the library bumphunter]? ** Thus how to setup the design matrix that all 4 contrasts are analysed with bumphunter? Code: > norm.swan MethylSet (storageMode: lockedEnvironment) assayData: 485512 features, 20 samples element names: Meth, Unmeth phenoData sampleNames: 7766130001_R01C01 7766130001_R02C01 ... 7766130037_R05C02 (20 total) varLabels: Sample_ID Sample_Plate ... filenames (15 total) varMetadata: labelDescription Annotation array: IlluminaHumanMethylation450k annotation: ilmn12.hg19 Preprocessing Method: SWAN (based on a MethylSet preprocesses as 'Raw (no normalization or bg correction)' minfi version: 1.10.1 Manifest version: 0.4.0 > > # convert MethylSet (norm.swan) into GenomicRatioSet through RatioSet > norm.swan2 <- ratioConvert(norm.swan, what = "both", keepCN = TRUE) > norm.swan2 <- mapToGenome(norm.swan2) Loading required package: IlluminaHumanMethylation450kanno.ilmn12.hg19 > > # experimental factors > treatment <- pd$TrtmntGrps > subject <- as.factor(pd$SubjectID) > > design <- model.matrix(~0+ treatment + subject) > > design treatmentoldctrl treatmentoldWY treatmentyoungctrl treatmentyoungWY subject31 subject34 subject35 1 1 0 0 0 0 0 0 2 0 0 0 1 0 0 0 3 0 1 0 0 0 0 0 4 1 0 0 0 0 0 0 5 1 0 0 0 0 0 0 6 0 0 0 1 0 0 1 7 0 1 0 0 0 0 0 8 0 0 1 0 0 0 0 9 0 0 0 1 0 0 0 10 1 0 0 0 0 0 0 11 0 1 0 0 0 0 0 12 1 0 0 0 0 0 0 13 0 0 1 0 0 0 0 14 0 0 0 1 0 1 0 15 0 0 1 0 0 1 0 16 0 1 0 0 0 0 0 17 0 1 0 0 0 0 0 18 0 0 1 0 0 0 1 19 0 0 1 0 1 0 0 20 0 0 0 1 1 0 0 subject43 subject52 subject62 subject64 subject65 subject66 1 0 1 0 0 0 0 2 1 0 0 0 0 0 3 <<snip>> > #sample (nonsense) bumphunter run to check everything worked > bumps <- bumphunter(norm.swan2, design = design, coef=1, B = 1, type = "Beta", cutoff = 0.10) [bumphunterEngine] Using a single core (backend: doSEQ, version: 1.4.2). [bumphunterEngine] Computing coefficients. [bumphunterEngine] Performing 1 permutations. [bumphunterEngine] Computing marginal permutation p-values. [bumphunterEngine] cutoff: 0.1 [bumphunterEngine] Finding regions. [bumphunterEngine] Found 205534 bumps. [bumphunterEngine] Computing regions for each permutation. Loading required package: rngtools Loading required package: pkgmaker Loading required package: registry Attaching package: 'pkgmaker' The following object is masked from 'package:IRanges': new2 [bumphunterEngine] Estimating p-values and FWER. Warning message: In bumphunterEngine(getMethSignal(object, type), design = design, : The use of the permutation test (B>0), is not recommended with multiple covariates, (ncol(design)>2). See vignette for more information. > I also have three small queries on the workflow: ** I normalized my data using SWAN. To use these data as input for bumphunter(), I had to do a 2-step conversion [ratioConvert() followed by mapToGenome()]. Is this the way to do this, or is there another (direct) way? > norm.swan2 <- ratioConvert(norm.swan, what = "both", keepCN = TRUE) > norm.swan2 <- mapToGenome(norm.swan2) ** To confirm, in the bumphunter() function, when 'type=beta' a 'cutoff=0.10' means that the absolute difference between the beta values of the 2 grps should at least be 0.10? Thus at least a difference of 10% in methylation status of one of the CpGs in the region (bump)? ** What is the current way of annotating the 450k data? Making use of "FDb.InfiniumMethylation.hg19" from Tim, "IlluminaHumanMethylation450kanno.ilmn12.hg19" from Kasper, or an allternative method? > sessionInfo() R Under development (unstable) (2013-11-19 r64265) 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 [4] LC_NUMERIC=C LC_TIME=English_United States.1252 attached base packages: [1] parallel stats graphics grDevices utils datasets methods base other attached packages: [1] doRNG_1.6 rngtools_1.2.4 [3] pkgmaker_0.22 registry_0.2 [5] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.2.1 IlluminaHumanMethylation450kmanifest_0.4.0 [7] minfi_1.10.1 bumphunter_1.4.1 [9] locfit_1.5-9.1 iterators_1.0.7 [11] foreach_1.4.2 Biostrings_2.32.0 [13] XVector_0.4.0 GenomicRanges_1.16.3 [15] GenomeInfoDb_1.0.2 IRanges_1.22.6 [17] lattice_0.20-29 Biobase_2.24.0 [19] BiocGenerics_0.10.0 loaded via a namespace (and not attached): [1] annotate_1.42.0 AnnotationDbi_1.26.0 base64_1.1 beanplot_1.1 codetools_0.2-8 compiler_3.1.0 [7] DBI_0.2-7 digest_0.6.4 genefilter_1.46.1 grid_3.1.0 illuminaio_0.6.0 limma_3.20.2 [13] MASS_7.3-33 matrixStats_0.8.14 mclust_4.3 multtest_2.20.0 nlme_3.1-117 nor1mix_1.1-4 [19] plyr_1.8.1 preprocessCore_1.26.1 R.methodsS3_1.6.1 RColorBrewer_1.0-5 Rcpp_0.11.1 reshape_0.8.5 [25] RSQLite_0.11.4 siggenes_1.38.0 splines_3.1.0 stats4_3.1.0 stringr_0.6.2 survival_2.37-7 [31] tools_3.1.0 XML_3.98-1.1 xtable_1.7-3 zlibbioc_1.10.0 > --------------------------------------------------------- Guido Hooiveld, PhD Nutrition, Metabolism & Genomics Group Division of Human Nutrition Wageningen University Biotechnion, Bomenweg 2 NL-6703 HD Wageningen the Netherlands tel: (+)31 317 485788 fax: (+)31 317 483342 email: guido.hooiveld@wur.nl internet: http://nutrigene.4t.com http://scholar.google.com/citations?user=qFHaMnoAAAAJ http://www.researcherid.com/rid/F-4912-2010 [[alternative HTML version deleted]]