I have been a regular user of the wateRmelon package since its inception, and have encountered new errors with wateRmelon_1.18.0. I was wondering if anyone has run into the same issue. Normalization using dasen() on the new Illumina EPIC chips (aka 850K) gives the following error:

dasen(melon.LumiSet)

Error in (function (od, vd)  :
  object and replacement value dimnames differ

This is a methyLumiSet object:

Object Information:
MethyLumiSet (storageMode: lockedEnvironment)
assayData: 866895 features, 1919 samples
  element names: betas, methylated, pvals, unmethylated
protocolData: none
phenoData
  sampleNames: 11 23 ... 1050 (1919 total)
  varLabels: PlateNumber Position ... barcodes (42 total)
  varMetadata: labelDescription
featureData
  featureNames: cg00000029 cg00000103 ... rs9839873 (866895 total)
  fvarLabels: Probe_ID DESIGN COLOR_CHANNEL
  fvarMetadata: labelDescription
experimentData: use 'experimentData(object)'
Annotation: IlluminaHumanMethylationEpic
Major Operation History:
            submitted            finished
1 2017-01-08 20:00:02 2017-01-08 21:14:54
2 2017-01-08 20:00:02 2017-01-08 21:14:54
3 2017-01-08 21:15:06 2017-01-08 21:22:30
                                                               command
1 NChannelSetToMethyLumiSet2(NChannelSet = dats, parallel = parallel,
2                                                    n = n, oob = oob)
3                                           Subset of 866895 features.

There have been a number of changes to key functions in methylation-related packages (e.g., ChAMP) to deal with the EPIC chip, so I thought that was the most likely contributing factor. So I went back to some 450K data from a few months ago that had been normalized using wateRmelon's dasen() without any previous issue, and I ran the same call on a different dataset set (same object name here though).

dasen(melon.LumiSet)

Error in rbind(deparse.level, ...) :
  numbers of columns of arguments do not match

Below is traceback() result on the previous object (the 850K data), after subsetting it a bit so what is returned is more readable (otherwise it just fills the output with the betas):

traceback()
11: stop("object and replacement value dimnames differ")
10: (function (od, vd)
    {
        if (is.null(vd))
            od <- seq_along(od)
        else if (!setequal(od, vd))
            stop("object and replacement value dimnames differ")
        od
    })(dots[[1L]][[2L]], dots[[2L]][[2L]])
9: mapply(FUN = f, ..., SIMPLIFY = FALSE)
8: Map(function(od, vd) {
       if (is.null(vd))
           od <- seq_along(od)
       else if (!setequal(od, vd))
           stop("object and replacement value dimnames differ")
       od
   }, dimnames(obj), dimnames(value))
7: .validate_assayDataElementReplace(obj, value)
6: assayDataElementReplace(object, "betas", value)
5: `betas<-`(`*tmp*`, value = c(0.765455893254262, 0.346322423811164,
   0.79320412999323, 0.84092723540146, 0.84896966594138, 0.834306200303351,
   0.901385362090579, 0.884037322605878, 0.60005993815556, 0.787555761543364,
   0.900060127125923, 0.1115827362559, 0.432910551304066, 0.806950808334054,
   0.603133083995342, 0.855765117635621, 0.864837304162418, 0.78520105886802,
   1.14472340361489, 0.438024138169323, 0.904903918466845, 0.856004612884625,
   0.872655767543147, 1.34053973279993, 0.294631710362047, 0.806198267564966,
   0.741435562805873, 0.24167156574564, 0.848020263964805, 0.596505864455023,
   0.838222781251156, 0.830319888734353, 0.884913868105098, 0.359941799245014,
   0.961472543633483, 0.770877200155365, 0.777847152847153, 0.350500377170557,
   0.842492462311558, 0.84092723540146, 0.831593420583592, 0.841105430183357,
   0.881447587354409, 0.884327469911387, 0.60005993815556, 0.751856705985146,
   0.874604430379747, 0.1115827362559, 0.421208057191832, 0.84480122324159,
   0.610154944639215, 0.889023552700683, 0.781185636249664, 0.837959643552187,
   1.14472340361489, 0.359941799245014, 0.905788423153693, 0.865380176353816,
   0.887914384320679, 1.26850888264193, 0.327938531960111, 0.788886796006781,
   0.762795732178754, 0.231658890102254, 0.88911126995473, 0.414400682843556,
   0.845163937483681, 0.822439379369493, 0.878624255514104, 0.393743740060147,
   0.961472543633483, 0.789270988132365, 0.743733794295592, 0.34234574114154,
   0.83549565436453, 0.84092723540146, 0.824683304972585, 0.810903884847547,
   0.892335766423358, 0.907196352979486, 0.60005993815556, 0.767401749414808,
   0.88911126995473, 0.1115827362559, 0.381207993361825, 0.806515825094553,
   0.904546732075583, 0.87263624109378, 0.788774573733275, 0.811194653299916,
   1.14472340361489, 0.421208057191832, 0.914182111200645, 0.893555153017447,
   0.9039458622315, 1.26850888264193, 0.284943181818182, 0.813448054848797,
   0.697977988226261, 0.24167156574564, 0.874163994502978, 0.531670084301663,
   0.841304430936269, 0.822439379369493, 0.886090748230536, 0.387614481944666,
   0.800086418050182, 0.770877200155365))
4: `betas<-`(`*tmp*`, value = c(0.765455893254262, 0.346322423811164,
   0.79320412999323, 0.84092723540146, 0.84896966594138, 0.834306200303351,
   0.901385362090579, 0.884037322605878, 0.60005993815556, 0.787555761543364,
   0.900060127125923, 0.1115827362559, 0.432910551304066, 0.806950808334054,
   0.603133083995342, 0.855765117635621, 0.864837304162418, 0.78520105886802,
   1.14472340361489, 0.438024138169323, 0.904903918466845, 0.856004612884625,
   0.872655767543147, 1.34053973279993, 0.294631710362047, 0.806198267564966,
   0.741435562805873, 0.24167156574564, 0.848020263964805, 0.596505864455023,
   0.838222781251156, 0.830319888734353, 0.884913868105098, 0.359941799245014,
   0.961472543633483, 0.770877200155365, 0.777847152847153, 0.350500377170557,
   0.842492462311558, 0.84092723540146, 0.831593420583592, 0.841105430183357,
   0.881447587354409, 0.884327469911387, 0.60005993815556, 0.751856705985146,
   0.874604430379747, 0.1115827362559, 0.421208057191832, 0.84480122324159,
   0.610154944639215, 0.889023552700683, 0.781185636249664, 0.837959643552187,
   1.14472340361489, 0.359941799245014, 0.905788423153693, 0.865380176353816,
   0.887914384320679, 1.26850888264193, 0.327938531960111, 0.788886796006781,
   0.762795732178754, 0.231658890102254, 0.88911126995473, 0.414400682843556,
   0.845163937483681, 0.822439379369493, 0.878624255514104, 0.393743740060147,
   0.961472543633483, 0.789270988132365, 0.743733794295592, 0.34234574114154,
   0.83549565436453, 0.84092723540146, 0.824683304972585, 0.810903884847547,
   0.892335766423358, 0.907196352979486, 0.60005993815556, 0.767401749414808,
   0.88911126995473, 0.1115827362559, 0.381207993361825, 0.806515825094553,
   0.904546732075583, 0.87263624109378, 0.788774573733275, 0.811194653299916,
   1.14472340361489, 0.421208057191832, 0.914182111200645, 0.893555153017447,
   0.9039458622315, 1.26850888264193, 0.284943181818182, 0.813448054848797,
   0.697977988226261, 0.24167156574564, 0.874163994502978, 0.531670084301663,
   0.841304430936269, 0.822439379369493, 0.886090748230536, 0.387614481944666,
   0.800086418050182, 0.770877200155365))
3: .local(mns, fudge, ...)
2: dasen(melon.LumiSet[35:70, 1:3])
1: dasen(melon.LumiSet[35:70, 1:3])

So one thing that stands out to me is that some betas (highlighted) have values outside [0,1]. Now the normalization quit due to error, so maybe these betas have not finished being processed.

Any thoughts?

 
BiocInstaller::biocValid()

* sessionInfo()

R version 3.3.2 (2016-10-31)
Platform: x86_64-redhat-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

locale:
 [1] LC_CTYPE=en_US.UTF-8          LC_NUMERIC=C
 [3] LC_TIME=en_US.UTF-8           LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=en_US.UTF-8       LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_US.UTF-8          LC_NAME=en_US.UTF-8
 [9] LC_ADDRESS=en_US.UTF-8        LC_TELEPHONE=en_US.UTF-8
[11] LC_MEASUREMENT=en_US.UTF-8    LC_IDENTIFICATION=en_US.UTF-8

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

other attached packages:
 [1] IlluminaHumanMethylationEPICanno.ilm10b2.hg19_0.6.0
 [2] stringr_1.1.0
 [3] xlsx_0.5.7
 [4] xlsxjars_0.6.1
 [5] rJava_0.9-8
 [6] wateRmelon_1.18.0
 [7] illuminaio_0.16.0
 [8] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.0
 [9] ROC_1.50.0
[10] lumi_2.26.3
[11] methylumi_2.20.0
[12] FDb.InfiniumMethylation.hg19_2.2.0
[13] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
[14] GenomicFeatures_1.26.2
[15] matrixStats_0.51.0
[16] ggplot2_2.2.1
[17] reshape2_1.4.2
[18] scales_0.4.1
[19] ChAMP_2.4.1
[20] IlluminaHumanMethylationEPICmanifest_0.3.0
[21] Illumina450ProbeVariants.db_1.10.0
[22] DMRcate_1.10.2
[23] DMRcatedata_1.10.1
[24] DSS_2.14.0
[25] bsseq_1.10.0
[26] FEM_3.2.0
[27] graph_1.52.0
[28] org.Hs.eg.db_3.4.0
[29] impute_1.48.0
[30] igraph_1.0.1
[31] corrplot_0.77
[32] marray_1.52.0
[33] limma_3.30.7
[34] Matrix_1.2-7.1
[35] AnnotationDbi_1.36.0
[36] ChAMPdata_2.2.0
[37] minfi_1.20.2
[38] bumphunter_1.14.0
[39] locfit_1.5-9.1
[40] iterators_1.0.8
[41] foreach_1.4.3
[42] Biostrings_2.42.1
[43] XVector_0.14.0
[44] SummarizedExperiment_1.4.0
[45] GenomicRanges_1.26.2
[46] GenomeInfoDb_1.10.2
[47] IRanges_2.8.1
[48] S4Vectors_0.12.1
[49] Biobase_2.34.0
[50] BiocGenerics_0.20.0
[51] BiocInstaller_1.24.0

loaded via a namespace (and not attached):
  [1] R.utils_2.5.0
  [2] RSQLite_1.1-1
  [3] htmlwidgets_0.8
  [4] grid_3.3.2
  [5] trimcluster_0.1-2
  [6] BiocParallel_1.8.1
  [7] munsell_0.4.3
  [8] codetools_0.2-15
  [9] preprocessCore_1.36.0
 [10] nleqslv_3.0.3
 [11] statmod_1.4.27
 [12] miniUI_0.1.1
 [13] colorspace_1.3-2
 [14] fastICA_1.2-0
 [15] knitr_1.15.1
 [16] robustbase_0.92-7
 [17] isva_1.8
 [18] biovizBase_1.22.0
 [19] diptest_0.75-7
 [20] R6_2.2.0
 [21] doParallel_1.0.10
 [22] flexmix_2.3-13
 [23] bitops_1.0-6
 [24] reshape_0.8.6
 [25] assertthat_0.1
 [26] nnet_7.3-12
 [27] gtable_0.2.0
 [28] affy_1.52.0
 [29] sva_3.22.0
 [30] ensembldb_1.6.2
 [31] genefilter_1.56.0
 [32] rtracklayer_1.34.1
 [33] lazyeval_0.2.0
 [34] acepack_1.4.1
 [35] GEOquery_2.40.0
 [36] dichromat_2.0-0
 [37] checkmate_1.8.2
 [38] yaml_2.1.14
 [39] backports_1.0.4
 [40] httpuv_1.3.3
 [41] qvalue_2.6.0
 [42] Hmisc_4.0-2
 [43] tools_3.3.2
 [44] nor1mix_1.2-2
 [45] affyio_1.44.0
 [46] RColorBrewer_1.1-2
 [47] DNAcopy_1.48.0
 [48] siggenes_1.48.0
 [49] Rcpp_0.12.8
 [50] plyr_1.8.4
 [51] base64enc_0.1-3
 [52] zlibbioc_1.20.0
 [53] purrr_0.2.2
 [54] RCurl_1.95-4.8
 [55] BiasedUrn_1.07
 [56] rpart_4.1-10
 [57] openssl_0.9.6
 [58] cluster_2.0.5
 [59] magrittr_1.5
 [60] data.table_1.10.0
 [61] colourpicker_0.3
 [62] mvtnorm_1.0-5
 [63] whisker_0.3-2
 [64] missMethyl_1.8.0
 [65] mime_0.5
 [66] xtable_1.8-2
 [67] RPMM_1.20
 [68] XML_3.98-1.5
 [69] mclust_5.2.1
 [70] gridExtra_2.2.1
 [71] biomaRt_2.30.0
 [72] tibble_1.2
 [73] KernSmooth_2.23-15
 [74] R.oo_1.21.0
 [75] htmltools_0.3.5
 [76] mgcv_1.8-15
 [77] Formula_1.2-1
 [78] tidyr_0.6.0
 [79] DBI_0.5-1
 [80] MASS_7.3-45
 [81] fpc_2.1-10
 [82] permute_0.9-4
 [83] quadprog_1.5-5
 [84] R.methodsS3_1.7.1
 [85] Gviz_1.18.1
 [86] RefFreeEWAS_2.0
 [87] GenomicAlignments_1.10.0
 [88] registry_0.3
 [89] IlluminaHumanMethylation450kmanifest_0.4.0
 [90] foreign_0.8-67
 [91] plotly_4.5.6.9000
 [92] annotate_1.52.1
 [93] rngtools_1.2.4
 [94] pkgmaker_0.22
 [95] multtest_2.30.0
 [96] beanplot_1.2
 [97] ruv_0.9.6
 [98] doRNG_1.6
 [99] VariantAnnotation_1.20.2
[100] digest_0.6.11
[101] base64_2.0
[102] htmlTable_1.8
[103] dendextend_1.3.0
[104] kernlab_0.9-25
[105] shiny_0.14.2
[106] Rsamtools_1.26.1
[107] gtools_3.5.0
[108] modeltools_0.2-21
[109] nlme_3.1-128
[110] jsonlite_1.2
[111] viridisLite_0.1.3
[112] BSgenome_1.42.0
[113] lattice_0.20-34
[114] httr_1.2.1
[115] DEoptimR_1.0-8
[116] survival_2.39-5
[117] GO.db_3.4.0
[118] interactiveDisplayBase_1.12.0
[119] shinythemes_1.1.1
[120] prabclus_2.2-6
[121] class_7.3-14
[122] stringi_1.1.2
[123] AnnotationHub_2.6.4
[124] latticeExtra_0.6-28
[125] memoise_1.0.0
[126] dplyr_0.5.0

* Out-of-date packages
        Package   LibPath                 Installed Built   ReposVer
bold    "bold"    "/home/share/R/library" "0.3.5"   "3.3.1" "0.4.0"
gdsfmt  "gdsfmt"  "/home/share/R/library" "1.10.0"  "3.3.1" "1.10.1"
rgl     "rgl"     "/home/share/R/library" "0.96.0"  "3.3.1" "0.97.0"
RSQLite "RSQLite" "/home/share/R/library" "1.1-1"   "3.3.1" "1.1-2"
tidyr   "tidyr"   "/home/share/R/library" "0.6.0"   "3.3.1" "0.6.1"
xml2    "xml2"    "/home/share/R/library" "1.0.0"   "3.3.1" "1.1.0"
        Repository
bold    "https://cran.rstudio.com/src/contrib"
gdsfmt  "https://bioconductor.org/packages/3.4/bioc/src/contrib"
rgl     "https://cran.rstudio.com/src/contrib"
RSQLite "https://cran.rstudio.com/src/contrib"
tidyr   "https://cran.rstudio.com/src/contrib"
xml2    "https://cran.rstudio.com/src/contrib"

update with biocLite()

* Packages too new for Bioconductor version '3.4'

       Version      LibPath
plotly "4.5.6.9000" "/home/share/R/library"
readxl "0.1.1.9000" "/home/share/R/library"

downgrade with biocLite(c("plotly", "readxl"))

Error: 6 package(s) out of date; 2 package(s) too new

Hello Wade,

This is a very interesting problem and I have not seen such an error before. I cannot seem to recreate this error using my current BioC install. And I do not think there is an issue with the code within wateRmelon as dasen has remained unchanged for a very long time (and has been working on EPIC arrays since migrating to them)

It appears to me that dasen is still correctly normalizing the data just at the end of the function when it appends the information to the resultant object it fails. My initial thought is that the error is originating from the Biobase package - which could mean that the 'betas<-' replacement methods packaged in methylumi for methylumiset objects are becoming redundant. But this doesn't really explain anything as I am also using Biobase_2.34.0. I will look into this further.

The betas > 1 is probably an artefact of using quantile normalization on a small subset of probes (such as the example in the trace-back) and the variation between the probes may be large enough to offset the normalized intensities, resulting in weird betas. This shouldn't happen when normalizing the full-set of probes.

As a solution you could try running dasen manually, and then doing the replacement manually:

# Generate sentrix positions (R0#C0#) from vector barcodes:
roco <- unlist(data.frame(strsplit(barcodes,  '_'),stringsAsFactors=FALSE)[2,])
norm <- dasen(mns = methylated(melon.LumiSet),
                    uns = unmethylated(melon.LumiSet),
                    onetwo = fData(melon.LumiSet)[,fot(melon.LumiSet)],
                    fudge = 100,
                    roco = roco, # Otherwise NULL
                    ret2 = TRUE # FALSE if you only want beta matrix
                    )
melon.LumiSet.norm <- melon.LumiSet
betas(melon.LumiSet.norm)        <- norm$beta
methylated(melon.LumiSet.norm)   <- norm$methylated
unmethylated(melon.LumiSet.norm) <- norm$unmethylated

Hope this is helpful.

Tyler

Hi Tyler,

Thanks for the quick reply and suggestion. I have tried that and the problem still exists when you attempt to replace any of the data in the slots using '<-'.

See below for a trivial case in the next to last line where I just take the data and add zero to it, and attempt to replace it back. The classes and dimnames all seem fine.

Wade

> # Generate sentrix positions (R0#C0#) from vector barcodes:
>
> roco <- unlist(data.frame(strsplit(pData(melon.LumiSet.pf)$barcodes,  '_'),stringsAsFactors=FALSE)[2,])
>
>
>
> # subsetting large object to reduce comptutational time just to see if it works
>
> norm <- dasen(mns = methylated(melon.LumiSet.pf[,1:100]),
+
+               uns = unmethylated(melon.LumiSet.pf[,1:100]),
+
+               onetwo = fData(melon.LumiSet.pf[,1:100])[,fot(melon.LumiSet.pf[,1:100])],
+
+               fudge = 100,
+
+               roco = roco[1:100], # Otherwise NULL
+
+               ret2 = TRUE # FALSE if you only want beta matrix
+
+ )
>
>
>
> melon.LumiSet.norm <- melon.LumiSet.pf[,1:100]
>
> betas(melon.LumiSet.norm)        <- norm$beta
Error in (function (od, vd)  :
  object and replacement value dimnames differ
>
>
>
> all(rownames(betas(melon.LumiSet.norm))==rownames(norm$beta))
[1] TRUE
>
> all(colnames(betas(melon.LumiSet.norm))==colnames(norm$beta))
[1] TRUE
>
>
>
> str(betas(melon.LumiSet.norm))
 num [1:863715, 1:100] 0.627 0.862 0.841 0.884 0.884 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:863715] "cg00000029" "cg00000103" "cg00000109" "cg00000155" ...
  ..$ : chr [1:100] "200223270003_R01C01" "200223270003_R02C01" "200223270003_R03C01" "200223270003_R04C01" ...
>
> str(norm$beta)
 num [1:863715, 1:100] 0.615 0.846 0.825 0.886 0.89 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:863715] "cg00000029" "cg00000103" "cg00000109" "cg00000155" ...
  ..$ : chr [1:100] "200223270003_R01C01" "200223270003_R02C01" "200223270003_R03C01" "200223270003_R04C01" ...
>
>
>
> betas(melon.LumiSet.norm)        <- betas(melon.LumiSet.norm)+0
Error in (function (od, vd)  :
  object and replacement value dimnames differ
>
>
>
> is(melon.LumiSet.norm)
[1] "MethyLumiSet"     "MethyLumi"        "methylData"       "eSet"
[5] "VersionedBiobase" "Versioned"
>

Interesting about the names producing that, but I don't think it is caused by NAs in my data:

anyis.na(colnames(betas(melon.LumiSet.norm))))
[1] FALSE
anyis.na(colnames(norm$beta)))
[1] FALSE

any(is.na(rownames(betas(melon.LumiSet.norm))))
[1] FALSE
any(is.na(rownames(norm$beta)))
[1] FALSE

I will start a fresh session and only load in the wateRmelon package, and try again.

Update: If I just call load(wateRmelon) I still get the error.

For now, I am going to just construct a new methyLumiSet object and proceed. Thanks for the time you spent looking into it.