R/makeSummarizedExperimentFromExpressionSet.R

In Bioconductor/SummarizedExperiment: A container (S4 class) for matrix-like assays

##
## makeSummarizedExperimentFromExpressionSet

## coercion

.from_rowRanges_to_FeatureData <- function(from)
{
    if (is(from, "GRanges")) {
        fd <- .from_GRanges_to_FeatureData(from)
    } else if (is(from, "GRangesList")) {
        fd <- .from_GRangesList_to_FeatureData(from)
    } else {
        stop("class ", sQuote(class(from)),
             " is not a supported type for rowRanges coercion")
    }
    featureNames(fd) <- names(from)
    fd
}

.from_GRanges_to_FeatureData <- function(from)
{
    data <- as.data.frame(from)

    ## the first mcols are automatically included in the data.frame from
    ## as.data.frame, the secondary mcols holds the metadata for the first
    ## metadata columns.
    metaData <- mcols(mcols(from, use.names=FALSE), use.names=FALSE)
    if (is.null(metaData)) {
        metaData <- as.data.frame(matrix(ncol=0, nrow=NCOL(data)))
    } else {
        metaData <- as.data.frame(metaData)
    }
    AnnotatedDataFrame(data, metaData)
}
.from_GRangesList_to_FeatureData <- function(from)
{
    data <- as.data.frame(mcols(from, use.names=FALSE))

    ## the first mcols are automatically included in the data.frame from
    ## as.data.frame, the secondary mcols holds the metadata for the first
    ## metadata columns.
    metaData <- mcols(mcols(from, use.names=FALSE), use.names=FALSE)
    if (is.null(metaData)) {
        metaData <- as.data.frame(matrix(ncol=0, nrow=NCOL(data)))
    } else {
        metaData <- as.data.frame(metaData)
    }
    AnnotatedDataFrame(data, metaData)
}

.from_AnnotatedDataFrame_to_DataFrame <- function(from)
{
    df <- DataFrame(pData(from), row.names=rownames(from))
    mcols(df) <- DataFrame(varMetadata(from))
    df
}

.from_DataFrame_to_AnnotatedDataFrame <- function(df)
{
    data <- as(df, "data.frame")
    metaData <- mcols(df, use.names=FALSE)
    if (is.null(metaData)) {
        metaData <- as.data.frame(matrix(ncol=0, nrow=NCOL(data)))
    } else {
        metaData <- as(metaData, "data.frame")
    }
    AnnotatedDataFrame(data, metaData)
}

## If the ExpressionSet has featureData with range information make
## GRanges out of that, otherwise make an empty GRangesList with names
## from the featureNames
naiveRangeMapper <- function(from)
{
    nms <- featureNames(from)
    res <- tryCatch({
        makeGRangesFromDataFrame(pData(featureData(from)),
                                 keep.extra.columns = TRUE)
    }, error = function(e) {
        res <- relist(GRanges(), vector("list", length=length(nms)))
        mcols(res) <- .from_AnnotatedDataFrame_to_DataFrame(featureData(from))
        res
    })
    names(res) <- nms
    res
}

# Simple ProbeId to Range mapper
# Probes with multiple ranges are dropped
# The sign of the chromosome location is assumed to contain the strand
# information
probeRangeMapper <- function(from)
{
    annotation <- annotation(from)
    if (identical(annotation, character(0))) {
        return(naiveRangeMapper(from))
    }
    if (requireNamespace("annotate", quietly = TRUE)) {
        annotationPackage <- annotate::annPkgName(annotation)
        test <- require(annotationPackage, character.only = TRUE,
                        quietly = TRUE)
        if (test) {
            db <- get(annotationPackage, envir = asNamespace(annotationPackage))
            pid <- featureNames(from)
            locs <- AnnotationDbi::select(
                db, pid, columns = c("CHR", "CHRLOC", "CHRLOCEND"))
            locs <- na.omit(locs)
            dups <- duplicated(locs$PROBEID)
            if (any(dups)) {
                locs <- locs[!dups, , drop = FALSE]
            }
            strand <- ifelse(locs$CHRLOC > 0, "+", "-")
            res <- GRanges(seqnames = locs$CHR,
                           ranges = IRanges(abs(locs$CHRLOC),
                                            abs(locs$CHRLOCEND)),
                           strand = strand)
            names(res) <- locs$PROBEID

            if (NROW(res) < length(pid)) {
                warning(length(pid) - NROW(res),
                        " probes could not be mapped.", call. = FALSE)
            }
            res
        } else {
            stop("Failed to load ", sQuote(annotationPackage), " package",
                 call. = FALSE)
        }
    } else {
        stop("Failed to load annotate package", call. = FALSE)
    }
}

# Simple ProbeId to Gene mapper
# Is there a way to get the txDb given the annotation package?
geneRangeMapper <- function(txDbPackage, key = "ENTREZID")
{
    function(from) {
        annotation <- annotation(from)
        if (identical(annotation, character(0))) {
            return(naiveRangeMapper(from))
        }
        if (requireNamespace("annotate", quietly = TRUE)) {
            annotationPackage <- annotate::annPkgName(annotation)
            test <- require(annotationPackage, character.only = TRUE,
                            quietly = TRUE)
            if (test) {
                db <- get(annotationPackage,
                          envir = asNamespace(annotationPackage))
                pid <- featureNames(from)
                probeIdToGeneId <-
                    AnnotationDbi::mapIds(db, pid, key, "PROBEID")
                geneIdToProbeId <-
                    setNames(names(probeIdToGeneId), probeIdToGeneId)

                if (requireNamespace(txDbPackage, quietly = TRUE)) {
                    txDb <- get(txDbPackage, envir = asNamespace(txDbPackage))
                    genes <- GenomicFeatures::genes(txDb)
                    probesWithAMatch <-
                        probeIdToGeneId[probeIdToGeneId %in% names(genes)]
                    res <- genes[probesWithAMatch]
                    names(res) <- geneIdToProbeId[names(res)]
                    if (NROW(res) < length(pid)) {
                        warning(length(pid) - NROW(res),
                                " probes could not be mapped.", call. = FALSE)
                    }
                    res
                } else {
                    stop("Failed to load ", sQuote(txDbPackage), " package",
                         call. = FALSE)
                }

            } else {
                stop("Failed to load ", sQuote(annotationPackage), " package",
                     call. = FALSE)
            }
        } else {
            stop("Failed to load annotate package", call. = FALSE)
        }
    }
}

makeSummarizedExperimentFromExpressionSet <-
    function(from, mapFun = naiveRangeMapper, ...)
{
    mapFun <- match.fun(mapFun)
    rowRanges <- mapFun(from, ...)
    matches <- match(names(rowRanges),
                     featureNames(from),
                     nomatch = 0)
    from <- from[matches, drop = FALSE]
    assays <- as.list(assayData(from))
    colData <- .from_AnnotatedDataFrame_to_DataFrame(phenoData(from))
    metadata <- SimpleList(
        experimentData = experimentData(from),
        annotation = annotation(from),
        protocolData = protocolData(from)
    )

    SummarizedExperiment(
        assays = assays,
        rowRanges = rowRanges,
        colData = colData,
        metadata = metadata
    )
}

setAs("ExpressionSet", "RangedSummarizedExperiment", function(from)
{
    makeSummarizedExperimentFromExpressionSet(from)
})

setAs("ExpressionSet", "SummarizedExperiment", function(from)
{
    as(makeSummarizedExperimentFromExpressionSet(from), "SummarizedExperiment")
})

setAs("RangedSummarizedExperiment", "ExpressionSet",
      function(from)
{
    assayData <- list2env(as.list(assays(from)))

    numAssays <- length(assayData)

    if (numAssays == 0) {
        assayData$exprs <- new("matrix")
    } else if (!"exprs" %in% ls(assayData)) {
        ## if there isn't an exprs assay we need to pick one as exprs,
        ## so rename the first element exprs and issue a warning.
        exprs <- ls(assayData)[[1]]
        warning("No assay named ", sQuote("exprs"), " found, renaming ",
                exprs, " to ", sQuote("exprs"), ".")
        assayData[["exprs"]] <- assayData[[exprs]]
        rm(list=exprs, envir=assayData)
    }
    lockEnvironment(assayData, bindings = TRUE)

    featureData <- .from_rowRanges_to_FeatureData(rowRanges(from))
    phenoData <- .from_DataFrame_to_AnnotatedDataFrame(colData(from))

    metadata <- metadata(from)

    experimentData <- if (!is.null(metadata$experimentData)) {
        metadata$experimentData
    } else {
        MIAME()
    }

    annotation <- if (!is.null(metadata$annotation)) {
        metadata$annotation
    } else {
        character()
    }

    protocolData <- if (!is.null(metadata$protocolData)) {
        metadata$protocolData
    } else {
        annotatedDataFrameFrom(assayData, byrow=FALSE)
    }

    ExpressionSet(assayData,
                  phenoData = phenoData,
                  featureData = featureData,
                  experimentData = experimentData,
                  annotation = annotation,
                  protocolData = protocolData
                  )
})

setAs("SummarizedExperiment", "ExpressionSet", function(from)
    as(as(from, "RangedSummarizedExperiment"), "ExpressionSet")
)