diff --git a/R/CNV.R b/R/CNV.R index 4a9bacb..413acb9 100644 --- a/R/CNV.R +++ b/R/CNV.R @@ -61,7 +61,7 @@ getGLADProfiles<-function(rnb.set,refbased=TRUE){ annot <- annotation(rnb.set) target <- rnb.set@target - if(target %in% c("probes27","probes450","probesEPIC", "probesEPICv2")){ + if(target %in% c("probes27","probes450","probesEPIC", "probesEPICv2", "probesMSA")){ target <- "probes450" } cnv.reference.data<-rnb.get.cnv.annotations(target) diff --git a/R/RnBSet-class.R b/R/RnBSet-class.R index 3eca4a1..73a3927 100644 --- a/R/RnBSet-class.R +++ b/R/RnBSet-class.R @@ -1177,6 +1177,8 @@ setMethod("mergeSamples", signature(object = "RnBSet"), platform<-"EPIC" }else if (object@target=="probesEPICv2"){ platform<-"EPICv2" + }else if (object@target=="probesMSA"){ + platform<-"MSA" }else if (object@target=="probes450"){ platform<-"450k" }else if(object@target=="probes27"){ diff --git a/R/RnBeadRawSet-class.R b/R/RnBeadRawSet-class.R index 459304c..102a0d0 100644 --- a/R/RnBeadRawSet-class.R +++ b/R/RnBeadRawSet-class.R @@ -173,7 +173,8 @@ RnBeadRawSet<-function( beta.offset=100, summarize.bead.counts=TRUE, summarize.regions=TRUE, - region.types = rnb.region.types.for.analysis(ifelse(platform=="MMBC", "mm10", ifelse(target=="probesEPICv2", "hg38", rnb.getOption("assembly")))), ## TODO: Implement assembly selection + region.types = rnb.region.types.for.analysis(ifelse(platform=="MMBC", "mm10", + ifelse(target=="probesEPICv2"|| target=="probesMSA", "hg38", rnb.getOption("assembly")))), ## TODO: Implement assembly selection useff=rnb.getOption("disk.dump.big.matrices"), ffcleanup=FALSE){ @@ -252,6 +253,9 @@ RnBeadRawSet<-function( }else if (platform =="EPICv2") { target <- "probesEPICv2" assembly <- "hg38" + }else if (platform =="MSA") { + target <- "probesMSA" + assembly <- "hg38" }else if (platform =="450k") { target <- "probes450" assembly <- ifelse(genome.assembly == "hg19", "hg19", "hg38") @@ -422,6 +426,8 @@ setAs("MethyLumiSet", "RnBeadRawSet", platform="EPIC" }else if(annotation(from)=="IlluminaMethylationEPICv2"){ platform="EPICv2" + }else if(annotation(from)=="IlluminaMethylationScreeningArray"){ + platform="MSA" }else if(annotation(from)=="IlluminaHumanMethylation450k"){ platform="450k" }else if(annotation(from)=="IlluminaHumanMethylation27k"){ @@ -507,6 +513,9 @@ setAs("RnBeadRawSet","MethyLumiSet", }else if(from@target == "probesEPICv2"){ probeIDs<-rnb.get.annotation("controlsEPICv2", assembly = "hg38")[,"Target"] probeIDs<-paste(probeIDs, unlist(sapply(table(probeIDs)[unique(probeIDs)], seq, from=1 )), sep=".") + }else if(from@target == "probesMSA"){ + probeIDs<-rnb.get.annotation("controlsMSA", assembly = "hg38")[,"Target"] + probeIDs<-paste(probeIDs, unlist(sapply(table(probeIDs)[unique(probeIDs)], seq, from=1 )), sep=".") }else if(from@target == "probes450"){ probeIDs<-rnb.get.annotation("controls450")[,"Target"] probeIDs<-paste(probeIDs, unlist(sapply(table(probeIDs)[unique(probeIDs)], seq, from=1 )), sep=".") @@ -521,6 +530,8 @@ setAs("RnBeadRawSet","MethyLumiSet", annotation(mset@QC) <- "IlluminaMethylationEPIC" }else if(from@target == "probesEPICv2"){ annotation(mset@QC) <- "IlluminaMethylationEPICv2" + }else if(from@target == "probesMSA"){ + annotation(mset@QC) <- "IlluminaMethylationScreeningArray" }else if(from@target == "probes450"){ annotation(mset@QC) <- "IlluminaHumanMethylation450k" }else if(from@target == "probes27"){ @@ -532,6 +543,8 @@ setAs("RnBeadRawSet","MethyLumiSet", annotation(mset) <- "IlluminaMethylationEPIC" }else if(from@target == "probesEPICv2"){ annotation(mset) <- "IlluminaMethylationEPICv2" + }else if(from@target == "probesMSA"){ + annotation(mset) <- "IlluminaMethylationScreeningArray" }else if(from@target == "probes450"){ annotation(mset) <- "IlluminaHumanMethylation450k" }else if(from@target == "probes27"){ @@ -558,7 +571,10 @@ setAs("RGChannelSet", "RnBeadRawSet", function(from, to) { } else if (assay.name == "IlluminaHumanMethylationEPICv2") { assay.name <- "probesEPICv2" platform.name <- "EPICv2" - } else if (assay.name == "IlluminaHumanMethylation450k") { + } else if (assay.name == "IlluminaHumanMethylationScreeningArray") { + assay.name <- "probesMSA" + platform.name <- "MSA" + } else if (assay.name == "IlluminaHumanMethylation450k") { assay.name <- "probes450" platform.name <- "450k" } else if (assay.name == "IlluminaHumanMethylation27k") { @@ -569,7 +585,8 @@ setAs("RGChannelSet", "RnBeadRawSet", function(from, to) { } ## Use RnBeads' mapping from probe IDs to addresses - probes.all <- rnb.get.annotation(assay.name, ifelse(assay.name=="probesEPICv2", "hg38", rnb.getOption("assembly"))) ## TODO: Impove genome build selection + is.EPICv2.or.MSA <- ifelse(assay.name == "probesEPICv2" || assay.name == "probesMSA", TRUE, FALSE) + probes.all <- rnb.get.annotation(assay.name, ifelse(is.EPICv2.or.MSA, "hg38", rnb.getOption("assembly"))) ## TODO: Impove genome build selection probes.all <- lapply(probes.all, function(x) { result <- as.data.frame(mcols(x)[, c("Design", "Color", "AddressA", "AddressB")]) rownames(result) <- names(x) @@ -577,7 +594,7 @@ setAs("RGChannelSet", "RnBeadRawSet", function(from, to) { } ) probes.all <- do.call(rbind, unname(probes.all)) - controls.all <- rnb.get.annotation(sub("^probes", "controls", assay.name), ifelse(assay.name=="probesEPICv2", "hg38", rnb.getOption("assembly"))) ## TODO: Impove genome build selection + controls.all <- rnb.get.annotation(sub("^probes", "controls", assay.name), ifelse(is.EPICv2.or.MSA, "hg38", rnb.getOption("assembly"))) ## TODO: Impove genome build selection controls.all <- controls.all[, "ID"] ## Extract data on signals @@ -674,7 +691,8 @@ setAs("RnBeadRawSet", "RGChannelSet", function(from, to){ # probes.all <- do.call(rbind, unname(probes.all)) probes.all <- annotation(from)[,c("Design", "Color", "AddressA", "AddressB")] genome.assembly<-rnb.getOption("assembly") - controls.all <- rnb.get.annotation(sub("^probes", "controls", assay.name), ifelse(assay.name=="probesEPICv2", "hg38", genome.assembly)) + is.EPICv2.or.MSA<-ifelse(assay.name == "probesEPICv2" || assay.name == "probesMSA", TRUE, FALSE) + controls.all <- rnb.get.annotation(sub("^probes", "controls", assay.name), ifelse(is.EPICv2.or.MSA, "hg38", genome.assembly)) controls.all <- controls.all[, "ID"] # Obtain methylated and unmethylated intensities @@ -747,6 +765,8 @@ setAs("RnBeadRawSet", "RGChannelSet", function(from, to){ anno <- "IlluminaHumanMethylationEPIC" }else if(assay.name %in% "probesEPICv2"){ anno <- "IlluminaHumanMethylationEPICv2" + }else if(assay.name %in% "probesMSA"){ + anno <- "IlluminaHumanMethylationScreeningArray" }else if(assay.name %in% "probes450"){ anno <- "IlluminaHumanMethylation450k" }else if(assay.name %in% "probes27"){ @@ -1305,6 +1325,8 @@ intensities.by.color<-function(raw.set, if (add.controls) { if (raw.set@target == "probesEPICv2") { ncd <- rnb.get.annotation("controlsEPICv2", assembly = "hg38") + } else if (raw.set@target == "probesMSA") { + ncd <- rnb.get.annotation("controlsMSA", assembly = "hg38") } else { ncd <- rnb.get.annotation(ifelse(raw.set@target == "probesEPIC", "controlsEPIC", "controls450")) } diff --git a/R/RnBeadSet-class.R b/R/RnBeadSet-class.R index 4709208..4b8bf8f 100644 --- a/R/RnBeadSet-class.R +++ b/R/RnBeadSet-class.R @@ -136,6 +136,7 @@ setMethod("initialize", "RnBeadSet", .Object@qc<-qc genome.assembly<-rnb.getOption("assembly") + is.EPICv2.or.MSA<-ifelse(target == "probesEPICv2" || target == "probesMSA", TRUE, FALSE) callNextMethod(.Object, pheno=pheno, @@ -143,7 +144,7 @@ setMethod("initialize", "RnBeadSet", meth.sites=meth.sites, covg.sites=covg.sites, status=status, - assembly=ifelse(target=="probesMMBC", "mm10", ifelse(target=="probesEPICv2", "hg38", genome.assembly)), + assembly=ifelse(target=="probesMMBC", "mm10", ifelse(is.EPICv2.or.MSA, "hg38", genome.assembly)), target=target ) @@ -246,6 +247,9 @@ RnBeadSet<-function( } else if (platform == "EPICv2") { target <- "probesEPICv2" assembly <- "hg38" + } else if (platform == "MSA") { + target <- "probesMSA" + assembly <- "hg38" } else if (platform == "450k") { target <- "probes450" assembly <- ifelse(genome.assembly == "hg19", "hg19", "hg38") diff --git a/R/annotations.R b/R/annotations.R index a5a5b30..14b3227 100644 --- a/R/annotations.R +++ b/R/annotations.R @@ -98,7 +98,24 @@ EPIC.CONTROL.TARGETS <- c( "staining" = "STAINING", "target removal" = "TARGET REMOVAL") -EPICv2.CONTROL.TARGETS <- c( ## TODO: Not validated for EPIC v2 (Most likely it's compatible) +EPICv2.CONTROL.TARGETS <- c( + "bisulfite conversion I" = "BISULFITE CONVERSION I", + "bisulfite conversion II" = "BISULFITE CONVERSION II", + "extension" = "EXTENSION", + "hybridization" = "HYBRIDIZATION", + "negative control" = "NEGATIVE", + "non-polymorphic" = "NON-POLYMORPHIC", + "norm A" = "NORM_A", + "norm C" = "NORM_C", + "norm G" = "NORM_G", + "norm T" = "NORM_T", + "restoration" = "RESTORATION", + "specificity I" = "SPECIFICITY I", + "specificity II" = "SPECIFICITY II", + "staining" = "STAINING", + "target removal" = "TARGET REMOVAL") + +MSA.CONTROL.TARGETS <- c( ## TODO: Not validated for MSA (Most likely it's compatible) "bisulfite conversion I" = "BISULFITE CONVERSION I", "bisulfite conversion II" = "BISULFITE CONVERSION II", "extension" = "EXTENSION", @@ -766,7 +783,7 @@ rnb.get.assemblies <- function() { #' rnb.get.annotation("promoters") #' } rnb.get.annotation <- function(type = "CpG", assembly = "hg19") { - if (type == "probesEPICv2" || type == "controlsEPICv2") { + if (type == "probesEPICv2" || type == "controlsEPICv2" || type == "probesMSA" || type == "controlsMSA") { assembly = "hg38" } @@ -1364,7 +1381,7 @@ rnb.get.chromosomes <- function(assembly = "hg19") { #' Extracts all control probe types in the HumanMethylation450 assay. #' #' @param target A singleton of type \code{character}, specifying the microarray platform. -#' \code{"probesEPICv2"},\code{"probesEPIC"},\code{"probes450"} and \code{"probes27"} correspond to MethylationEPICv2, +#' \code{"probesEPICv2"},\code{"probesMSA"},\code{"probesEPIC"},\code{"probes450"} and \code{"probes27"} correspond to MethylationEPICv2, #' Methylation EPIC, HumanMethylation450, and HumanMethylation27 microarrays respectively. #' #' @return \code{character} vector of control targets. @@ -1380,6 +1397,8 @@ rnb.infinium.control.targets <- function(target="probes450") { return(EPIC.CONTROL.TARGETS) }else if(target=="probesEPICv2"){ return(EPICv2.CONTROL.TARGETS) + }else if(target=="probesMSA"){ + return(MSA.CONTROL.TARGETS) }else if(target=="probes450"){ return(HM450.CONTROL.TARGETS) }else if(target=="probes27"){ diff --git a/R/batch.R b/R/batch.R index e4c1f46..43b9213 100644 --- a/R/batch.R +++ b/R/batch.R @@ -707,7 +707,6 @@ rnb.section.dreduction.internal <- function(report, pcoordinates, sample.phenoty } else { report <- rnb.add.section(report, "Low-dimensional Representation", stext) } - ## Mention that only some sites and/or regions are selected if (!is.null(attr(pcoordinates, "selected"))) { selected <- attr(pcoordinates, "selected") diff --git a/R/batch.quality.R b/R/batch.quality.R index d0a8235..1362eee 100644 --- a/R/batch.quality.R +++ b/R/batch.quality.R @@ -59,6 +59,11 @@ rnb.execute.batch.qc <- function(rnb.set, pcoordinates, permutations = NULL) { id.col<-"ID" type.col<-"Target" ctrls<-"controlsEPICv2" + } else if(rnb.set@target=="probesMSA"){ + CONTROL.TYPES <- rnb.infinium.control.targets("probesMSA")[c(1:4, 6, 11:14)] + id.col<-"ID" + type.col<-"Target" + ctrls<-"controlsMSA" } else if(rnb.set@target=="probes450"){ CONTROL.TYPES <- rnb.infinium.control.targets("probes450")[c(1:4, 6, 11:14)] id.col<-"ID" diff --git a/R/combining.R b/R/combining.R index 299a574..3c3dc2c 100644 --- a/R/combining.R +++ b/R/combining.R @@ -122,6 +122,7 @@ rnb.combine.pheno <- function(dataset1, dataset2) { #' @author Yassen Assenov #' @export rnb.combine.arrays <- function(dataset1, dataset2, type="common") { + ## TODO: Update for MSA if (!inherits(dataset1, "RnBeadSet")) { stop("Invalid value for dataset1") } @@ -132,7 +133,7 @@ rnb.combine.arrays <- function(dataset1, dataset2, type="common") { stop("Incompatible assemblies") } i <- c(dataset1@target, dataset2@target) - common.platform <- c("probesEPICv2" = "EPICv2", "probesEPIC" = "EPIC", "probes450" = "450k", "probes27" = "27k") ## TODO: Validate for EPICv2 (on EPICv1 ?) + common.platform <- c("probesEPICv2" = "EPICv2", "probesEPIC" = "EPIC", "probes450" = "450k", "probes27" = "27k") if (!(i[1] %in% names(common.platform))) { stop("Unsupported platform for dataset1") } diff --git a/R/controlPlots.R b/R/controlPlots.R index d6eaeee..e2578fe 100644 --- a/R/controlPlots.R +++ b/R/controlPlots.R @@ -61,6 +61,8 @@ rnb.plot.control.boxplot <- function( meta <- rnb.get.annotation("controlsEPIC", assembly = "hg19") }else if(rnb.set@target=="probesEPICv2"){ meta <- rnb.get.annotation("controlsEPICv2", assembly = "hg38") + }else if(rnb.set@target=="probesMSA"){ + meta <- rnb.get.annotation("controlsMSA", assembly = "hg38") }else if(rnb.set@target=="probes450"){ meta <- rnb.get.annotation("controls450", assembly = "hg19") }else if(rnb.set@target=="probes27"){ @@ -73,6 +75,8 @@ rnb.plot.control.boxplot <- function( types<-rnb.infinium.control.targets(rnb.set@target)[c(14,4,3,15,1:2,12:13,6,11)] }else if(rnb.set@target=="probesEPICv2"){ types<-rnb.infinium.control.targets(rnb.set@target)[c(14,4,3,15,1:2,12:13,6,11)] + }else if(rnb.set@target=="probesMSA"){ + types<-rnb.infinium.control.targets(rnb.set@target)[c(14,4,3,15,1:2,12:13,6,11)] }else if(rnb.set@target=="probes450"){ types<-rnb.infinium.control.targets(rnb.set@target)[c(13,4,3,14,1:2,11:12,6)] }else if(rnb.set@target=="probes27"){ @@ -85,12 +89,12 @@ rnb.plot.control.boxplot <- function( warning("Unoptimized probe type, plotting performance may be decreased") } - if(rnb.set@target=="probes450" || rnb.set@target=="probesEPIC" || rnb.set@target=="probesMMBC" || rnb.set@target=="probesEPICv2"){ + if(rnb.set@target=="probes450" || rnb.set@target=="probesEPIC" || rnb.set@target=="probesMMBC" || rnb.set@target=="probesEPICv2" || rnb.set@target=="probesMSA"){ rownames(meta)<-meta[["ID"]] ### TODO: Remove the following passage ### for testing purposes only! - if(rnb.set@target=="probesEPIC" || rnb.set@target=="probesEPICv2"){ - meta<-rnb.update.controlsEPIC.enrich(meta) ## TODO: Not validated for EPICv2 + if(rnb.set@target=="probesEPIC" || rnb.set@target=="probesEPICv2" || rnb.set@target=="probesMSA"){ + meta<-rnb.update.controlsEPIC.enrich(meta) ## TODO: Not validated for EPICv2 & MSA } meta <- meta[type == meta[["Target"]], ] ids<-as.character(meta[["ID"]]) @@ -125,7 +129,7 @@ rnb.plot.control.boxplot <- function( scales<-lapply(qc(rnb.set), get.unified.scale) - if(rnb.set@target=="probes450" || rnb.set@target=="probesEPIC" || rnb.set@target=="probesMMBC" || rnb.set@target=="probesEPICv2"){ + if(rnb.set@target=="probes450" || rnb.set@target=="probesEPIC" || rnb.set@target=="probesMMBC" || rnb.set@target=="probesEPICv2" || rnb.set@target=="probesMSA"){ ## Shorten the words describing probe's expected intensity INTENSITIES <- c("Background" = "Bgnd", "High" = "High", "Low" = "Low", "Medium" = "Med") levels(meta[, "Expected Intensity"]) <- INTENSITIES[levels(meta[, "Expected Intensity"])] @@ -220,14 +224,14 @@ rnb.plot.negative.boxplot<- function( ...) { - if(rnb.set@target=="probesEPIC" || rnb.set@target=="probesEPICv2" || rnb.set@target=="probesMMBC" ){ + if(rnb.set@target=="probesEPIC" || rnb.set@target=="probesEPICv2" || rnb.set@target=="probesMMBC" || rnb.set@target=="probesMSA"){ genome.assembly<-rnb.getOption("assembly") - meta <- rnb.get.annotation(gsub("probes", "controls", rnb.set@target), assembly=c("probesEPIC"=genome.assembly, "probesMMBC"="mm10", "probesEPICv2"="hg38")[rnb.set@target]) + meta <- rnb.get.annotation(gsub("probes", "controls", rnb.set@target), assembly=c("probesEPIC"=genome.assembly, "probesMMBC"="mm10", "probesEPICv2"="hg38", "probesMSA"="hg38")[rnb.set@target]) ## Extract intensities of the control probes ### TODO: Remove the following passage ### for testing purposes only! - if(rnb.set@target=="probesEPIC" || rnb.set@target=="probesEPICv2"){ - meta<-rnb.update.controlsEPIC.enrich(meta) ## TODO: Not validated for EPICv2 (duplicate) + if(rnb.set@target=="probesEPIC" || rnb.set@target=="probesEPICv2" || rnb.set@target=="probesMSA"){ + meta<-rnb.update.controlsEPIC.enrich(meta) ## TODO: Not validated for EPICv2 & MSA (duplicate) } meta <- meta["NEGATIVE" == meta[["Target"]], ] ids<-as.character(meta[["ID"]]) @@ -367,9 +371,9 @@ rnb.plot.control.barplot<-function( { - if(rnb.set@target=="probesEPIC" || rnb.set@target=="probesEPICv2" || rnb.set@target=="probesMMBC"){ + if(rnb.set@target=="probesEPIC" || rnb.set@target=="probesEPICv2" || rnb.set@target=="probesMMBC" || rnb.set@target=="probesMSA"){ genome.assembly<-rnb.getOption("assembly") ## TODO: Improve genome build selection - control.meta.data <- rnb.get.annotation(gsub("probes", "controls", rnb.set@target), assembly=c("probesEPIC"=genome.assembly, "probesEPICv2"="hg38", "probesMMBC"="mm10")[rnb.set@target]) + control.meta.data <- rnb.get.annotation(gsub("probes", "controls", rnb.set@target), assembly=c("probesEPIC"=genome.assembly, "probesEPICv2"="hg38", "probesMSA"="hg38", "probesMMBC"="mm10")[rnb.set@target]) ### TODO: Remove the following passage ### for testing purposes only! #if(rnb.set@target=="probesEPIC"){ @@ -438,7 +442,7 @@ rnb.plot.control.barplot<-function( ## get meta information - if(rnb.set@target=="probes450" || rnb.set@target=="probesEPIC" || rnb.set@target=="probesMMBC" || rnb.set@target=="probesEPICv2"){ + if(rnb.set@target=="probes450" || rnb.set@target=="probesEPIC" || rnb.set@target=="probesMMBC" || rnb.set@target=="probesEPICv2" || rnb.set@target=="probesMSA"){ meta<-subset(control.meta.data, ID==id) }else if(rnb.set@target=="probes27"){ meta<-subset(control.meta.data, Address==id) @@ -465,7 +469,7 @@ rnb.plot.control.barplot<-function( ### plot green channel - if(rnb.set@target=="probes450" || rnb.set@target=="probesEPIC" || rnb.set@target=="probesMMBC" || rnb.set@target=="probesEPICv2"){ + if(rnb.set@target=="probes450" || rnb.set@target=="probesEPIC" || rnb.set@target=="probesMMBC" || rnb.set@target=="probesEPICv2" || rnb.set@target=="probesMSA"){ main_txt_grn<-paste(probe, meta[,"Description"], "green channel", if(meta[, "Evaluate Green"]=="+") meta[, "Expected Intensity"] else "Background", sep=": ") main_txt_red<-paste(probe, meta[,"Description"], "red channel",if(meta[, "Evaluate Red"]=="+") meta[, "Expected Intensity"] else "Background", sep=": ") }else{ @@ -531,7 +535,7 @@ rnb.plot.control.barplot<-function( #' @author Yassen Assenov #' @noRd rnb.get.snp.matrix <- function(dataset, threshold.nas = 1) { - if (dataset@target %in% c("probes450", "probesEPIC", "probesMMBC", "probesEPICv2")) { + if (dataset@target %in% c("probes450", "probesEPIC", "probesMMBC", "probesEPICv2", "probesMSA")) { result <- meth(dataset, row.names=TRUE) result <- result[grep("^rs", rownames(result)), , drop = FALSE] } else if (dataset@target == "probes27") { diff --git a/R/dataImport.R b/R/dataImport.R index a2117ef..db40165 100644 --- a/R/dataImport.R +++ b/R/dataImport.R @@ -24,6 +24,10 @@ INTENSITY.SUMMARIZATION.INFO<-list( "typeIred"=list(Design="I", Color="Red", Msource="Red", Usource="Red", Maddress="AddressB", Uaddress="AddressA"), "typeIgrn"=list(Design="I", Color="Grn", Msource="Grn", Usource="Grn", Maddress="AddressB", Uaddress="AddressA"), "typeII"=list(Design="II", Color="Both", Msource="Grn", Usource="Red", Maddress="AddressA", Uaddress="AddressA")), + probesMSA=list( + "typeIred"=list(Design="I", Color="Red", Msource="Red", Usource="Red", Maddress="AddressB", Uaddress="AddressA"), + "typeIgrn"=list(Design="I", Color="Grn", Msource="Grn", Usource="Grn", Maddress="AddressB", Uaddress="AddressA"), + "typeII"=list(Design="II", Color="Both", Msource="Grn", Usource="Red", Maddress="AddressA", Uaddress="AddressA")), probes450=list( "typeIred"=list(Design="I", Color="Red", Msource="Red", Usource="Red", Maddress="AddressB", Uaddress="AddressA"), "typeIgrn"=list(Design="I", Color="Grn", Msource="Grn", Usource="Grn", Maddress="AddressB", Uaddress="AddressA"), @@ -283,7 +287,7 @@ read.data.dir<-function(dir, ifelse(nrow(beta.table)<30000L,"27k","450k")) }else{ platform.dict <- list("probesEPIC" = "EPIC", "probesEPICv2" = "EPICv2", - "probes27" = "27k", "probes450" = "450k") + "probesMSA" = "MSA", "probes27" = "27k", "probes450" = "450k") platform <- platform.dict$platform rm(platform.dict) } @@ -696,10 +700,11 @@ read.idat.files <- function(base.dir, "probes450"="HumanMethylation450", "probesEPIC"="MethylationEPIC", "probesEPICv2"="MethylationEPICv2", + "probesMSA"="MethylationScreeningArray", "probesMMBC"="MouseMethylationBeadChip") rnb.info(paste("Detected platform:", txt[platform])) - if (platform == "probesEPICv2" && genome.assembly != "hg38") { - rnb.info(paste0("MethylationEPICv2 is not supported for this session's genome assembly: ", genome.assembly, ". Changing genome assembly to: hg38")) + if ( (platform == "probesEPICv2" | platform == "probesMSA") && genome.assembly != "hg38") { + rnb.info(paste0(txt[platform], " is not supported for this session's genome assembly: ", genome.assembly, ". Changing genome assembly to: hg38")) rnb.options(assembly = "hg38") ## EPICv2 is only annotated in the RnBeads.hg38 package genome.assembly<-rnb.getOption("assembly") } @@ -712,6 +717,7 @@ read.idat.files <- function(base.dir, "probes450"=ifelse(genome.assembly == "hg19", "hg19", "hg38"), "probesEPIC"=ifelse(genome.assembly == "hg19", "hg19", "hg38"), "probesEPICv2"="hg38", + "probesMSA"="hg38", "probesMMBC"="mm10") annot_gr<-rnb.get.annotation(platform, genome[platform]) @@ -730,6 +736,11 @@ read.idat.files <- function(base.dir, ctrls.target.col<-"Target" neg.ctrl.indexes<-which(annot.ctrls[["Target"]]=="NEGATIVE") }else if(platform=="probesEPICv2"){ + id.col<-"ID" + ctrls.address.col<-"ID" + ctrls.target.col<-"Target" + neg.ctrl.indexes<-which(annot.ctrls[["Target"]]=="NEGATIVE") + }else if(platform=="probesMSA"){ id.col<-"ID" ctrls.address.col<-"ID" ctrls.target.col<-"Target" @@ -900,6 +911,9 @@ read.idat.files <- function(base.dir, }else if(platform %in% "probesEPICv2"){ rnb.platform<-"EPICv2" assembly<-"hg38" + }else if(platform %in% "probesMSA"){ + rnb.platform<-"MSA" + assembly<-"hg38" }else{ rnb.platform<-"MMBC" assembly<-"mm10" @@ -909,11 +923,10 @@ read.idat.files <- function(base.dir, sample.sheet<-data.frame(barcodes=barcode) } - # saveRDS(M, "/Users/baris.kalem/Code/RnBeads_Project/Kaur_CompareEPICv1_EPICv2/duplicated_probes_bug/M_good.RDS") ### solve the problem of duplicated probes ### in each pair select those that have a lower detection p-value - if(rnb.platform=="MMBC" || rnb.platform=="EPICv2"){ + if(rnb.platform=="MMBC" || rnb.platform=="EPICv2" || rnb.platform=="MSA"){ probe_names<-annot[["Name"]] dup_probe_names<-unique(probe_names[duplicated(probe_names)]) @@ -1798,8 +1811,11 @@ rnb.detect.infinium.platform <- function(idat.fnames){ if (all(file.sizes>10000000)) { return("probesEPIC") } + if (all(file.sizes>5000000)) { + return("probesMSA") ## idat file size is ~ 5MB + } if (all(file.sizes<5000000)) { - return("probesMMBC") + return("probesMMBC") ## idat file size is ~ 4.7MB } if (all(file.sizes<10000000)) { return("probes450") diff --git a/R/differentialMethylation.R b/R/differentialMethylation.R index ab3ae65..29abb1f 100644 --- a/R/differentialMethylation.R +++ b/R/differentialMethylation.R @@ -180,7 +180,7 @@ rowOneSampleTP <- function(X,mu=0,alternative="two.sided"){ #' limmaP #' -#' applies hierarchical modeling anlalogous to differential expression employed in the \code{limma} package and returns +#' applies hierarchical modeling analogous to differential expression employed in the \code{limma} package and returns #' p-values for differential methylation #' @author Fabian Mueller #' @param X Matrix on which the test is performed for every row @@ -1582,7 +1582,7 @@ rnb.section.diffMeth.introduction <- function(diffmeth,report){ #include information on the p-value method site.test.method <- get.site.test.method(diffmeth) txt <- c( - "In the following anlyses, p-values on the site level were computed using the ",site.test.method," method. " + "In the following analyses, p-values on the site level were computed using the ",site.test.method," method. " ) if (site.test.method == "limma"){ txt <- c(txt, @@ -2059,9 +2059,6 @@ rnb.section.diffMeth.site <- function(rnbSet,diffmeth,report,gzTable=FALSE){ sectionText <- paste(sectionText,"",sep="") rnb.add.paragraph(report, sectionText) - # nv_diff <- dmt ## New nv dmt df - # saveRDS(nv_diff, "/Users/baris.kalem/Code/RnBeads_Project/nv_probes_EPICv2/debug/nv_diff.RDS") - logger.info("Enriching nv-probes comparison table with HGNC symbols. see: GeneSymbol column in the finished table.") rnb.require("biomaRt") mart <- useMart("ENSEMBL_MART_ENSEMBL", dataset = "hsapiens_gene_ensembl",host="https://feb2023.archive.ensembl.org") diff --git a/R/exportGEO.R b/R/exportGEO.R index b9a689e..3d6820d 100644 --- a/R/exportGEO.R +++ b/R/exportGEO.R @@ -45,15 +45,15 @@ prepareGEOSampleInfoTemplate<-function(rnb.set, pheno.cols=1:ncol(pheno(rnb.set) }else{ sample.descr<-NULL } - - platform <- c("HumanMethylation27", "HumanMethylation450", "HumanMethylationEPIC", "HumanMethylationEPICv2")[match(rnb.set@target, c('probes27','probes450','probesEPIC', 'probesEPICv2'))] + platform <- c("HumanMethylation27", "HumanMethylation450", "HumanMethylationEPIC", "HumanMethylationEPICv2", "HumanMethylationScreeningArray")[match(rnb.set@target, c('probes27','probes450','probesEPIC', 'probesEPICv2', 'probesMSA'))] chipVersion<-"unknown" organism <- c("Homo sapiens", "Mus muscullus", "Rattus norvegicus")[which(!is.na(pmatch( c('hg','mm','rn'), assembly(rnb.set))))] platformID<-c( "HumanMethylation27"="GPL8490", "HumanMethylation450"="GPL13534", "HumanMethylationEPIC"="GPL23976", ## Alternative: GPL21145 - "HumanMethylationEPICv2"="GPL33022")[platform] + "HumanMethylationEPICv2"="GPL33022", + "HumanMethylationScreeningArray"="GPL34394")[platform] templateTitle <- c( "Sample_title", diff --git a/R/loading.R b/R/loading.R index 3366af9..4008c4f 100644 --- a/R/loading.R +++ b/R/loading.R @@ -423,6 +423,7 @@ rnb.section.import<-function(report, object, data.source, data.type=rnb.getOptio "probes450"="HumanMethylation450", "probesEPIC"="MethylationEPIC", "probesEPICv2"="MethylationEPICv2", + "probesMSA"="MethylationScreeningArray", "probesMMBC"="MouseMethylationBeadChip") txt[platform] }, @@ -584,7 +585,7 @@ rnb.step.import <- function(data.source, data.type = rnb.getOption("import.defau logger.status(c("Loaded data from", d.source)) ## Perform sex prediction - if (rnb.getOption("import.sex.prediction") && inherits(object, c("RnBeadRawSet","RnBiseqSet")) && object@target %in% c("probes450","probesEPIC","probesEPICv2","CpG")) { + if (rnb.getOption("import.sex.prediction") && inherits(object, c("RnBeadRawSet","RnBiseqSet")) && object@target %in% c("probes450","probesEPIC","probesEPICv2","probesMSA","CpG")) { object <- rnb.execute.sex.prediction(object) if (is.null(object@inferred.covariates$sex)) { object@inferred.covariates$sex <- FALSE diff --git a/R/main.R b/R/main.R index 44cfb3a..8d8345e 100644 --- a/R/main.R +++ b/R/main.R @@ -1287,7 +1287,7 @@ rnb.run.inference <- function(rnb.set, dir.reports, } if (inherits(rnb.set,"RnBSet") && rnb.getOption("inference.age.prediction") - && !(rnb.set@target == "probesEPIC" | rnb.set@target == "probesEPICv2")){ ## There is no age prediction for EPICv1 & EPICv2 + && !(rnb.set@target == "probesEPIC" | rnb.set@target == "probesEPICv2" | rnb.set@target == "probesMSA")){ ## There is no age prediction for EPICv1, EPICv2 & MSA ph <- pheno(rnb.set) ages <- ph$predicted_ages if(is.null(ages)){ diff --git a/R/normalization.R b/R/normalization.R index 9a52e27..920867c 100644 --- a/R/normalization.R +++ b/R/normalization.R @@ -157,6 +157,8 @@ rnb.execute.normalization<-function( disable.method(FALSE, 'not supported for Infinium 27k') } else if (bgcorr.method == "methylumi.noob" && object@target == "probesEPICv2") { disable.method(FALSE, 'methylumi.noob is not supported for MethylationEPICv2') + } else if (bgcorr.method == "methylumi.noob" && object@target == "probesMSA") { + disable.method(FALSE, 'methylumi.noob is not supported for MethylationScreeningArray') } else if (bgcorr.method == "methylumi.noob" && object@target == "probesEPIC") { disable.method(FALSE, 'methylumi.noob is not supported for MethylationEPIC') } else if (bgcorr.method == "methylumi.noob" && object@target == "probesMMBC") { @@ -228,6 +230,9 @@ rnb.execute.normalization<-function( if (inherits(object, "RnBeadSet") && object@target == "probesEPICv2" && !(method %in% accepted)) { disable.method(TRUE, 'not supported for HumanMethylationEPICv2') } + if (inherits(object, "RnBeadSet") && object@target == "probesMSA" && !(method %in% accepted)) { + disable.method(TRUE, 'not supported for HumanMethylationScreeningArray') + } accepted<-setdiff(.rnb.options[["accepted"]][["normalization.method"]], c("illumina", "minfi.funnorm", "swan", "wm.tost","wm.fuks", "wm.swan")) if (inherits(object, "RnBeadSet") && object@target == "probesMMBC" && !(method %in% accepted)) { disable.method(TRUE, 'not supported for Mouse Methylation Bead Chip') diff --git a/R/normalizationMethods.R b/R/normalizationMethods.R index bbf5d11..4bf70db 100644 --- a/R/normalizationMethods.R +++ b/R/normalizationMethods.R @@ -512,6 +512,12 @@ get.platform.tokens<-function(platform){ dict$bg_token<-"NEGATIVE" dict$id_col<-"ID" dict$trg_col<-"Target" + } else if(platform=="probesMSA"){ + dict$norm_token["Cy3"] <- 'NORM_(C|G)' + dict$norm_token["Cy5"] <- 'NORM_(A|T)' + dict$bg_token<-"NEGATIVE" + dict$id_col<-"ID" + dict$trg_col<-"Target" } else if(platform=="probesMMBC"){ dict$norm_token["Cy3"] <- 'NORM_(C|G)' dict$norm_token["Cy5"] <- 'NORM_(A|T)' diff --git a/R/options.R b/R/options.R index 44c4337..86a2062 100644 --- a/R/options.R +++ b/R/options.R @@ -539,7 +539,7 @@ rnb.is.option <- function(txt) { #' Skip the check of the loaded RnBSet object after loading. Helps with keeping the memory profile down} #' \item{\bold{\code{import.idat.platform}}\code{ = NULL}}{ #' Character specifying the Infinium platform that is uses. Has to be one of \code{'probes27'}, -#' \code{'probes450'}, \code{'probesEPIC'} or \code{'probesEPICv2'}. If \code{'auto'}, the platform is automatically detected +#' \code{'probes450'}, \code{'probesEPIC'}, \code{'probesEPICv2'} or \code{'probesMSA'}. If \code{'auto'}, the platform is automatically detected #' from the IDAT file names.} #' \item{\bold{\code{import.sex.prediction}}\code{ = TRUE}}{ #' Flag indicating if sex prediction is to be performed. Sex prediction is supported for Infinium 450k, EPIC, EPICv2 diff --git a/R/qualityControl.R b/R/qualityControl.R index c13b9a5..e7b7e51 100644 --- a/R/qualityControl.R +++ b/R/qualityControl.R @@ -369,6 +369,8 @@ add.qc.boxplots<-function(report, object){ ctypes<-rnb.infinium.control.targets(object@target)[c(14,4,3,15,1:2,12:13,6,11)] }else if(object@target=="probesEPICv2"){ ctypes<-rnb.infinium.control.targets(object@target)[c(14,4,3,15,1:2,12:13,6,11)] + }else if(object@target=="probesMSA"){ + ctypes<-rnb.infinium.control.targets(object@target)[c(14,4,3,15,1:2,12:13,6,11)] }else if(object@target=="probes450"){ ctypes<-rnb.infinium.control.targets(object@target)[c(13,4,3,14,1:2,11:12,6)] }else if(object@target=="probes27"){ @@ -399,6 +401,9 @@ add.qc.barplots<-function(report, object, sample.batch.size=50){ }else if(object@target=="probesEPICv2"){ cmd <- rnb.get.annotation("controlsEPICv2", assembly = "hg38") ctypes<-unique(cmd$Target)[unique(cmd$Target) %in% rnb.infinium.control.targets("probesEPICv2")[c(14,4,3,15,1:2,12:13,6,11)]] + }else if(object@target=="probesMSA"){ + cmd <- rnb.get.annotation("controlsMSA", assembly = "hg38") + ctypes<-unique(cmd$Target)[unique(cmd$Target) %in% rnb.infinium.control.targets("probesMSA")[c(14,4,3,15,1:2,12:13,6,11)]] ## TODO: Check control targets of the MSA platform }else if(object@target=="probes450"){ cmd <- rnb.get.annotation("controls450", assembly = rnb.getOption("assembly")) ctypes<-unique(cmd$Target)[unique(cmd$Target) %in% rnb.infinium.control.targets("probes450")[c(13,4,14,3,1:2,11:12,6)]] @@ -425,7 +430,8 @@ add.qc.barplots<-function(report, object, sample.batch.size=50){ cplots<-lapply(ctypes, function(type){ - if(object@target=="probes450" || object@target=="probesEPIC" || object@target=="probesMMBC" || object@target=="probesEPICv2"){ + if(object@target=="probes450" || object@target=="probesEPIC" || object@target == "probesMMBC" || + object@target == "probesEPICv2" || object@target == "probesMSA"){ cmdt <- cmd[cmd[["Target"]] == type, ] pn<-paste(type, 1:(dim(cmdt)[1]), sep=".") }else if(object@target=="probes27"){ @@ -439,7 +445,8 @@ add.qc.barplots<-function(report, object, sample.batch.size=50){ report=report, writeToFile=TRUE, numeric.names=TRUE, width=8, height=6, low.png=100, high.png=300, verbose=TRUE, name.prefix=portions[portion.id]) - if(object@target=="probes450" || object@target=="probesEPIC" || object@target=="probesMMBC" || object@target=="probesEPICv2"){ + if(object@target=="probes450" || object@target=="probesEPIC" || object@target == "probesMMBC" || + object@target == "probesEPICv2" || object@target == "probesMSA"){ names(plots)<-paste(type, 1:(dim(cmdt)[1])) }else if(object@target=="probes27"){ names(plots)<-as.character(cmdt$Name) @@ -461,7 +468,8 @@ add.qc.barplots<-function(report, object, sample.batch.size=50){ names(sn[[1]])<-portions - if(object@target=="probes450" || object@target=="probesEPIC" || object@target == "probesMMBC" || object@target == "probesEPICv2"){ + if(object@target=="probes450" || object@target=="probesEPIC" || object@target == "probesMMBC" || + object@target == "probesEPICv2" || object@target == "probesMSA"){ names(sn[[2]])<-1:length(plot.names) }else if(object@target=="probes27"){ names(sn[[2]])<-match(plot.names,cmd$Name[cmd$Type %in% rnb.infinium.control.targets("probes27")[c(10,3,2,11,1,9,6)]]) diff --git a/R/readGEO.R b/R/readGEO.R index f931a1c..da9981b 100644 --- a/R/readGEO.R +++ b/R/readGEO.R @@ -12,6 +12,7 @@ GEO.PLATFORMS <- c( "GPL8490" = "probes27", "GPL13534" = "probes450", "GPL16304" = "probes450", + "GPL34394" = "probesMSA", "GPL21145" = "probesEPIC", "GPL33022" = "probesEPICv2") ## https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL33022 @@ -326,7 +327,6 @@ rnb.geo.parse.series.matrix <- function(fname, verbose) { #' @author Yassen Assenov, modified by Baris Kalem #' @export rnb.read.geo <- function(accession = NULL, verbose = logger.isinitialized(), destdir = tempdir()) { - ## TODO: Make it compatible with EPIC v2 if (verbose) { rnb.logger.start("Loading GEO Data Series") } @@ -363,6 +363,6 @@ rnb.read.geo <- function(accession = NULL, verbose = logger.isinitialized(), des ## Parse the series matrix file result <- rnb.geo.parse.series.matrix(fname, verbose) - x <- c("probes27" = "27k", "probes450" = "450k", "probesEPIC" = "EPIC", "probesEPICv2" = "EPICv2")[result[[3]]] + x <- c("probes27" = "27k", "probes450" = "450k", "probesEPIC" = "EPIC", "probesEPICv2" = "EPICv2", "probesMSA" = "MSA")[result[[3]]] RnBeadSet(pheno = result[[1]], betas = result[[2]], platform = x) } diff --git a/R/sesame.R b/R/sesame.R index 92e8375..24c19d8 100644 --- a/R/sesame.R +++ b/R/sesame.R @@ -50,6 +50,8 @@ rnb.execute.pOOBAH <- function(raw.set, anno.table = NULL, pval.thresh = 0.05, v platform = "EPIC" }else if(raw.set@target == "probesEPICv2"){ platform = "EPICv2" + }else if(raw.set@target == "probesMSA"){ + platform = "MSA" }else{ stop("Invalid value for platform") } diff --git a/R/sex.R b/R/sex.R index fc1a280..0e8c229 100644 --- a/R/sex.R +++ b/R/sex.R @@ -57,6 +57,8 @@ rnb.get.XY.shifts <- function(rnb.set, signal.type = "raw") { probes.bad <- lapply(probes.bad, function(x) { which((mcols(x)[, "SNPs 3 Alternative"] != 0)) }) } else if(target == 'probesEPICv2'){ probes.bad <- lapply(probes.bad, function(x) { which((mcols(x)[, "SNPs 3 Alternative"] != 0)) }) + } else if(target == 'probesMSA'){ + probes.bad <- lapply(probes.bad, function(x) { which((mcols(x)[, "SNPs 3 Alternative"] != 0)) }) } else if (target == 'probes450'){ probes.bad <- lapply(probes.bad, function(x) { which((mcols(x)[, "SNPs 3"] != 0)) }) } @@ -237,7 +239,7 @@ rnb.execute.sex.prediction <- function(rnb.set) { stop("invalid value for rnb.set") } if (inherits(rnb.set, "RnBeadRawSet")) { - if (rnb.set@target != "probes450" && rnb.set@target != "probesEPIC" && rnb.set@target != "probesEPICv2") { + if (rnb.set@target != "probes450" && rnb.set@target != "probesEPIC" && rnb.set@target != "probesEPICv2" && rnb.set@target != "probesMSA") { stop("unsupported platform") } shifts <- rnb.get.XY.shifts(rnb.set) diff --git a/inst/extdata/options.txt b/inst/extdata/options.txt index 489909f..5ad9a25 100644 --- a/inst/extdata/options.txt +++ b/inst/extdata/options.txt @@ -28,7 +28,7 @@ import.bed.test logical no no TRUE import.bed.test.only logical no no FALSE import.skip.object.check logical no no FALSE import.idat.chunk.size integer no yes (1, NULL -import.idat.platform character no no (probes27,probes450,probesEPIC,probesEPICv2,probesMMBC,auto) auto +import.idat.platform character no no (probes27,probes450,probesEPIC,probesEPICv2,probesMSA,probesMMBC,auto) auto import.sex.prediction logical no no TRUE import.dpval.method character no no (controls,pOOBAH) controls qc logical no no TRUE