R/methods-ExpressionSetIllumina.R

In markdunning/beadarray: Quality assessment and low-level analysis for Illumina BeadArray data

setMethod("initialize", "ExpressionSetIllumina",
          function(.Object,
                   assayData = assayDataNew(exprs=exprs,se.exprs=se.exprs, nObservations=nObservations, Detection=Detection, storage.mode="list"),

                   phenoData = new("AnnotatedDataFrame"),
                   exprs=new("matrix"),
                   se.exprs=new("matrix"),
                   nObservations=new("matrix"),
                   Detection=new("matrix"),
                   annotation = character(),
                   featureData = new("AnnotatedDataFrame"),
                   experimentData = new("MIAME")
  )
 {
            .Object<-callNextMethod(.Object,
                           assayData = assayData,
                           phenoData = phenoData,
                           experimentData = experimentData,
                           annotation = annotation,
                           featureData = featureData
			   )
            
            .Object
          })


setMethod("[", "ExpressionSetIllumina", function(x, i, j, ..., drop = FALSE) {
  if (missing(drop))
    drop <- FALSE
  
  if (!missing(j)) {
    phenoData(x) <- phenoData(x)[j,, ..., drop = drop]
    protocolData(x) <- protocolData(x)[j,, ..., drop = drop]
  }
  
  if (!missing(i))
    featureData(x) <- featureData(x)[i,,..., drop=drop]
  ## assayData; implemented here to avoid function call
  orig <- assayData(x)
  ###I took this code from the eSet methods in Biobase to allow for empty se.exprs, nObservations, Detection
  storage.mode <- Biobase:::assayDataStorageMode(orig)
  
  assayData(x) <-
    switch(storage.mode,
           environment =,
           lockedEnvironment = {
             aData <- new.env(parent=emptyenv())
             if (missing(i))                     # j must be present
               for(nm in ls(orig)) {
                 if(nrow(orig[[nm]])>0)  aData[[nm]] <- orig[[nm]][, j, ..., drop = drop]
                 else aData[[nm]] <- orig[[nm]]
               }
             
             else {                              # j may or may not be present
               if (missing(j))
                 for(nm in ls(orig)){
                   if(nrow(orig[[nm]])>0) aData[[nm]] <- orig[[nm]][i,, ..., drop = drop]
                   else aData[[nm]] <- orig[[nm]]
                 }
               else
                 for(nm in ls(orig)){ 
                   if(nrow(orig[[nm]])>0) aData[[nm]] <-orig[[nm]][i, j, ..., drop = drop]
                   else aData[[nm]] <- orig[[nm]]
                 }
             }
             if ("lockedEnvironment" == storage.mode) assayDataEnvLock(aData)
             aData
           },
           
           list = {
             if (missing(i))                     # j must be present
               lapply(orig, function(obj) {
                 if(nrow(obj)>0) obj[, j, ..., drop = drop]
                 else obj
               })
             
             else {                              # j may or may not be present
               if (missing(j))
                 lapply(orig, function(obj){ 
                   if(nrow(obj)>0) obj[i,, ..., drop = drop]
                   else obj
                 })
               
               
               else
                 lapply(orig, function(obj){
                   if(nrow(obj)>0) obj[i, j, ..., drop = drop]
                   else obj
                 }) 
             }
           }
    )
  
  
  
  
  
  if(!is.null(x@QC) && !missing(j)) x@QC<-x@QC[j,]
  
  if(!is.null(x@channelData) && missing(j)) x@channelData<-x@channelData
  
  
  x
  
})



setAs("ExpressionSet","ExpressionSetIllumina",
      function(from)
      {
        
        to <- new("ExpressionSetIllumina")
        to@assayData = assayDataNew(exprs=exprs(from),storage.mode="list")
        to@channelData <- list()
        to@channelData[[1]] <- rep("G", length(sampleNames(from)))
        
        phenoData(to) <- phenoData(from)  
        featureData(to) <- featureData(from)
        newanno <- switch(annotation(from), 
                                 GPL6947="Humanv3", 
                                 GPL10558="Humanv4", 
                                 GPL6887="Mousev2", 
                                 GPL6102="Humanv2")
        
          
        if(!is.null(newanno)) {
          annotation(to) <- newanno
          to <- addFeatureData(to)
        }
        to
      })

## we turn off the validty check here compared to the version in Biobase aimed
## at the standard ExpressionSet class.
## This may be a bad idea, but it solves the problem reported here:
## https://support.bioconductor.org/p/90046
setReplaceMethod("exprs", signature(object="ExpressionSetIllumina",value="matrix"),
                 function(object, value) assayDataElementReplace(object, "exprs", value, validate = FALSE))

setAs("ExpressionSetIllumina", "GRanges",
      function(from)
      {
        
        
        annoName <- annotation(from)
        
        annoLoaded <- require(paste("illumina", annoName, ".db",sep=""), character.only=TRUE)
        
        if(annoLoaded){
          
          
          mapEnv <-  as.name(paste("illumina", annoName, "GENOMICLOCATION",sep=""))
          fn <- featureNames(from)
          fn <- fn[which(fn %in% mappedkeys(eval(mapEnv)))]
          
          locs <- mget(fn,eval(mapEnv),ifnotfound=NA)
          
          locs <- lapply(locs, function(x) gsub(" ", ",", x,fixed=T))
          
          asLocMatrix <- function(str){
            x<- do.call("rbind",sapply(strsplit(as.character(str), ",",fixed=T)[[1]], function(x) as.vector(strsplit(x, ":",fixed=T))))
          }
          
          locMat <- lapply(locs, asLocMatrix)
          
          rn <- rep(names(locs), unlist(lapply(locMat, nrow)))
          
          locMat <- do.call("rbind", locMat)
          rng <- GRanges(locMat[,1], IRanges(as.numeric(locMat[,2]), as.numeric(locMat[,3]),names=rn),strand=locMat[,4])
          #mcols(rng) <- df[match(names(rng), rownames(df)),]
          
          mcols(rng) <- data.frame(fData(from)[rn,], exprs(from)[rn,])
          
          sort(rng)
          
        }
        
      }
      
)


setValidity("ExpressionSetIllumina", function(object) {
  assayDataValidMembers(assayData(object), c("exprs", "se.exprs", "nObservations"))
})

setMethod("dim", "ExpressionSetIllumina", function(x) {

	nFeatures = nrow(fData(x))
	nSamps = length(sampleNames(x))
	nChannels = length(channelNames(x))

    c("Features"=nFeatures, "Samples"=nSamps, "Channels"=nChannels)
 } )



setMethod("exprs", signature(object="ExpressionSetIllumina"), function(object) assayDataElement(object, "exprs"))

#setGeneric("exprs<-", function(object, value) standardGeneric("exprs<-"))

setReplaceMethod("exprs", signature(object="ExpressionSetIllumina",value="matrix"), function(object, value){
	assayDataElementReplace(object, "exprs", value)
})

setMethod("se.exprs", signature(object="ExpressionSetIllumina"), function(object) assayDataElement(object, "se.exprs"))

#setGeneric("se.exprs<-", function(object, value) standardGeneric("se.exprs<-"))

setReplaceMethod("se.exprs", signature(object="ExpressionSetIllumina",value="matrix"), function(object, value){
	assayDataElementReplace(object, "se.exprs", value)
})

setGeneric("nObservations", function(object) standardGeneric("nObservations"))

setMethod("nObservations", signature(object="ExpressionSetIllumina"), function(object) assayDataElement(object, "nObservations"))

setGeneric("nObservations<-", function(object, value) standardGeneric("nObservations<-"))

setReplaceMethod("nObservations", signature(object="ExpressionSetIllumina",value="matrix"), function(object, value){
	assayDataElementReplace(object, "nObservations", value)
})

setGeneric("Detection", function(object) standardGeneric("Detection"))

setMethod("Detection", signature(object="ExpressionSetIllumina"), function(object) assayDataElement(object, "Detection"))

setGeneric("Detection<-", function(object, value) standardGeneric("Detection<-"))

setReplaceMethod("Detection", signature(object="ExpressionSetIllumina",value="matrix"), function(object, value){
	assayDataElementReplace(object, "Detection", value)
})

setMethod("show", signature(object="ExpressionSetIllumina"), function(object) {
  callNextMethod(object)
  
  cat("QC Information\n")
  cat(" Available Slots:  ")
  cat(names(object@QC))
  nms=selectSome(colnames(object@QC@data))
  cat("\n  QC Items:", paste(nms, collapse=", "))
  nms=selectSome(sampleNames(object@QC))
  cat("\n  sampleNames:", paste(nms, collapse=", "))
  cat("\n")
})


setGeneric("qcData", function(object) standardGeneric("qcData"))

setMethod("qcData", signature(object="ExpressionSetIllumina"), function(object) object@QC@data)

setGeneric("SampleGroup", function(object) standardGeneric("SampleGroup"))

setMethod("SampleGroup", signature(object = "ExpressionSetIllumina"), function(object) object@SampleGroup)
setGeneric("SampleGroup<-", function(object, value) standardGeneric("SampleGroup<-"))

setReplaceMethod("SampleGroup",
                 signature=signature(
                   object="ExpressionSetIllumina",
                   value="character"),
                 function(object, value) {
                   object@SampleGroup <- value
                   object
                 })

#setGeneric("exprs<-", function(object, value) standardGeneric("exprs<-"))

#setReplaceMethod("exprs", "ExpressionSetIllumina", function(object, value){
#	assayDataElementReplace(object, "exprs", value)
#})

#setReplaceMethod("exprs", c("ExpressionSetIllumina", "matrix"), function(object, value) {
#  assayDataElementReplace(object, "exprs", value)
#})

#setReplaceMethod("se.exprs", c("ExpressionSetIllumina", "matrix"), function(object, value) {
#  assayDataElementReplace(object, "se.exprs", value)
#})

.mergeAssayData<-function(x, y, newdimnames) {
  # this is derived from assayData combine method
  # differences:
  # - allows different number of reporters/features
  # - will merge data from identical column names into 1 column ie rbind()) 
  # - only works on 2-dimensional assayData elements



  combineElement <- function(x, y) {
    outarr<-array(NA,dim=c(length(newdimnames[[1]]),length(newdimnames[[2]])),newdimnames)
    mode(outarr)<-mode(x)
    outarr[rownames(y),colnames(y)]<-y
    outarr[rownames(x),colnames(x)]<-x
    outarr
  }


  storage.mode <- storageMode(x)
  nmfunc <- assayDataElementNames

  if (storageMode(y) != storage.mode)
    stop(paste("assayData must have same storage, but are ",
               storage.mode, ", ", storageMode(y), sep=""))
  if (length(nmfunc(x)) != length(nmfunc(y)))
    stop("assayData have different numbers of elements:\n\t",
         paste(nmfunc(x), collapse=" "), "\n\t",
         paste(nmfunc(y), collapse=" "))
  if (!all(nmfunc(x) == nmfunc(y)))
    stop(paste("assayData have different element names:",
               paste(nmfunc(x), collapse=" "),
               paste(nmfunc(y), collapse=" "), sep="\n\t"))
               
  for (nm in nmfunc(x)) {
    x<-assayDataElementReplace(x,nm,combineElement(assayDataElement(x,nm),assayDataElement(y,nm)))
  }
  x
}

.mergePhenodata<-function(x , y, samples) {
  variables<-union(colnames(pData(x)),colnames(pData(y)))
  outarr<-array(data=NA,dim=c(length(samples),length(variables)),dimnames=list(samples,variables))
  outarr[sampleNames(y),colnames(pData(y))]<-as.matrix(pData(y))
  outarr[sampleNames(x),colnames(pData(x))]<-as.matrix(pData(x))
  pd<-data.frame(outarr)
  vardescs<-union(colnames(varMetadata(x)),colnames(varMetadata(y)))
  outarr<-array(data=NA,dim=c(length(variables),length(vardescs)),dimnames=list(variables,vardescs))
  outarr[colnames(pData(y)),colnames(varMetadata(y))]<-as.matrix(varMetadata(y))
  outarr[colnames(pData(x)),colnames(varMetadata(x))]<-as.matrix(varMetadata(x))
  vd<-data.frame(outarr)
  new("AnnotatedDataFrame", data=pd, varMetadata=vd)
}


#setMethod("combine", signature(x="ExpressionSetIllumina",y="ExpressionSetIllumina"), function(x, y, ...) {
setMethod("combine", signature(x="ExpressionSetIllumina",y="ExpressionSetIllumina"), function(x, y) {


  if (class(x) != class(y))
    stop(paste("objects must be the same class, but are ",
               class(x), ", ", class(y), sep=""))
  newdimnames<-list(union(featureNames(x),featureNames(y)),union(colnames(exprs(x)),colnames(exprs(y))))
  x <- .mergeAssayData(x, y, newdimnames)
  # a bit of a hack to only keep the union, and discard double entries

  newsamplenames = union(sampleNames(x), sampleNames(y))
 	
  phenoData(x) <- .mergePhenodata(x, y, newsamplenames)

  experimentData(x) <- combine(experimentData(x),experimentData(y))

  protocolData(x) <- combine(protocolData(x), protocolData(y))

    
  ## annotation -- constant
  if (any(annotation(x) != annotation(y))) {
    warning("objects have different annotations: ",
         annotation(x), ", ", annotation(y))
    annotation(x)<-unique(c(annotation(x),annotation(y)))
  }


  ##Preserve the channel names of the resulting object

  x@channelData[[1]] = c(x@channelData[[1]],y@channelData[[1]])

  x@QC = combine(x@QC,y@QC)
  
  x

})