R/object_access.R

In NCBI-Hackathons/RClusterCT: Classifier for Single-cell RNA-seq Using Cell Clusters

#' An example Seurat object
#'
#' @return a Seurat object populated with data
#' from the [pbmc_matrix_small] scRNA-seq dataset, additionally
#' annotated with cluster assignments.
#'
#' @importFrom SeuratObject CreateSeuratObject CreateDimReducObject VariableFeatures
#' @export
so_pbmc <- function() {
    x <- pbmc_example_data()
    so <- SeuratObject::CreateSeuratObject(x$mat,
        meta.data = x$metadata
    )
    umap_dr <- SeuratObject::CreateDimReducObject(
        embeddings = x$umap,
        key = "umap_",
        assay = "RNA"
    )
    if (is_seurat_v5()) {
        so <- SeuratObject::SetAssayData(
            so,
            "data",
            SeuratObject::LayerData(so, layer = "counts")
        )
    } else {
        so <- SeuratObject::SetAssayData(
            so,
            "data",
            SeuratObject::GetAssayData(so, slot = "counts")
        )
    }
    so[["umap"]] <- umap_dr
    SeuratObject::VariableFeatures(so) <- x$vargenes
    so
}

#' An example SingleCellExperiment object
#'
#' @return a SingleCellExperiment object populated with data
#' from the [pbmc_matrix_small] scRNA-seq dataset, additionally
#' annotated with cluster assignments.
#'
#' @export
sce_pbmc <- function() {
    x <- pbmc_example_data()

    cols_to_keep <- c(
        "orig.ident",
        "nCount_RNA",
        "nFeature_RNA",
        "percent.mt",
        "RNA_snn_res.0.5",
        "classified"
    )

    md <- x$metadata[, cols_to_keep]
    # rename to more sce-like names
    colnames(md) <- c(
        "cell_source",
        "sum",
        "detected",
        "subsets_Mito_percent",
        "clusters",
        "cell_type"
    )
    SingleCellExperiment::SingleCellExperiment(
        list(
            counts = x$mat,
            logcounts = x$mat
        ),
        colData = md,
        reducedDims = list(
            UMAP = x$umap
        )
    )
}

pbmc_example_data <- function() {
    mat <- clustifyr::pbmc_matrix_small
    md <- clustifyr::pbmc_meta
    umap_cols <- c("UMAP_1", "UMAP_2")
    umap <- as.matrix(md[, umap_cols])
    md <- md[, setdiff(colnames(md), umap_cols)]
    vargenes <- clustifyr::pbmc_vargenes

    list(
        mat = mat,
        metadata = md,
        umap = umap,
        vargenes = vargenes
    )
}

#' Function to access object data
#' @return expression matrix, with genes as row names,
#' and cell types as column names
#' @export
object_data <- function(object, ...) {
    UseMethod("object_data", object)
}

#' @rdname object_data
#' @param object SingleCellExperiment or Seurat object
#' @param slot data to access
#' @param n_genes number of genes limit for Seurat variable genes, by default 1000,
#'   set to 0 to use all variable genes (generally not recommended)
#' @param ... additional arguments
#' @examples
#' so <- so_pbmc()
#' mat <- object_data(
#'     object = so,
#'     slot = "data"
#' )
#' mat[1:3, 1:3]
#' @export
object_data.Seurat <- function(
    object,
    slot,
    n_genes = 1000,
    ...) {
    if (slot == "data") {
        temp <- get_seurat_matrix(object, ...)
        return(temp)
    } else if (slot == "meta.data") {
        return(object@meta.data)
    } else if (slot == "var.genes") {
        vars <- SeuratObject::VariableFeatures(object)

        if (is.null(vars) || length(vars) <= 1) {
            message("variable genes not found, please manually specify with query_genes argument")
        }
        if ((length(vars) > n_genes) & (n_genes > 0)) {
            vars <- vars[seq_len(n_genes)]
        }

        return(vars)
    } else {
        stop(slot, " access method not implemented")
    }
}

#' @importFrom utils packageVersion
is_seurat_v5 <- function() {
    utils::packageVersion("SeuratObject") >= "5.0.0"
}

extract_v5_matrix <- function(x, ...) {
    ob_layers <- SeuratObject::Layers(x)
    if ("data" %in% ob_layers) {
        res <- SeuratObject::LayerData(x, layer = "data", ...)
    } else if ("counts" %in% ob_layers) {
        message("Unable to find 'data' layer, using 'count' layer instead")
        res <- SeuratObject::LayerData(x, layer = "counts", ...)
    } else {
        da <- DefaultAssay(x)
        stop(
            "\nUnable to find data or count layer in ", da, " Assay of SeuratObject\n",
            "Extracting data from V5 objects with multiple count\n",
            "or data layers is not supported"
        )
    }
    res
}

extract_v4_matrix <- function(x) {
    res <- SeuratObject::GetAssayData(x, layer = "data")

    if (length(res) == 0) {
        message("Unable to find 'data' slot, using 'count' slot instead")
        res <- SeuratObject::GetAssayData(x, layer = "count")
    }

    res
}

get_seurat_matrix <- function(x, warn = TRUE) {
    ob_assay <- SeuratObject::DefaultAssay(x)
    if (warn && ob_assay != "RNA") {
        warning(
            "Default assay of input Seurat object is ", ob_assay, "\n",
            "Data will be used from this assay rather than RNA"
        )
    }

    if (is_seurat_v5()) {
        res <- extract_v5_matrix(x)
    } else {
        res <- extract_v4_matrix(x)
    }
    res
}

#' @rdname object_data
#' @param object object after tsne or umap projections
#'  and clustering
#' @param slot data to access
#' @param ... additional arguments
#' @importFrom SingleCellExperiment logcounts colData
#' @examples
#' sce <- sce_pbmc()
#' mat <- object_data(
#'     object = sce,
#'     slot = "data"
#' )
#' mat[1:3, 1:3]
#' @export
object_data.SingleCellExperiment <- function(
    object,
    slot,
    ...) {
    if (slot == "data") {
        return(SingleCellExperiment::logcounts(object))
    } else if (slot == "meta.data") {
        return(as.data.frame(SingleCellExperiment::colData(object)))
    } else {
        stop(slot, " access method not implemented")
    }
}

#' Function to write metadata to object
#' @return object with newly inserted metadata columns
#' @export
write_meta <- function(object, ...) {
    UseMethod("write_meta", object)
}

#' @rdname write_meta
#' @param object object after tsne or umap projections
#'  and clustering
#' @param meta new metadata dataframe
#' @param ... additional arguments
#' @examples
#' so <- so_pbmc()
#' obj <- write_meta(
#'     object = so,
#'     meta = seurat_meta(so)
#' )
#' @export
write_meta.Seurat <- function(
    object,
    meta,
    ...) {
    object_new <- object
    object_new@meta.data <- meta
    object_new
}

#' @rdname write_meta
#' @param object object after tsne or umap projections
#'  and clustering
#' @param meta new metadata dataframe
#' @param ... additional arguments
#' @importFrom SingleCellExperiment colData
#' @importFrom S4Vectors DataFrame
#' @importFrom SummarizedExperiment colData<-
#' @examples
#' sce <- sce_pbmc()
#' obj <- write_meta(
#'     object = sce,
#'     meta = object_data(sce, "meta.data")
#' )
#' @export
write_meta.SingleCellExperiment <- function(
    object,
    meta,
    ...) {
    colData(object) <- S4Vectors::DataFrame(meta)
    object
}

#' Function to convert labelled seurat object to avg expression matrix
#' @return reference expression matrix, with genes as row names,
#' and cell types as column names
#' @examples
#' so <- so_pbmc()
#' ref <- seurat_ref(so, cluster_col = "seurat_clusters")
#' @export
seurat_ref <- function(seurat_object, ...) {
    UseMethod("seurat_ref", seurat_object)
}

#' @rdname seurat_ref
#' @param seurat_object seurat_object after tsne or umap projections
#'  and clustering
#' @param cluster_col column name where classified cluster names
#'  are stored in  seurat meta data, cannot be "rn"
#' @param var_genes_only whether to keep only var_genes in the final
#' matrix output, could also look up genes used for PCA
#' @param assay_name any additional assay data, such as ADT, to include.
#'  If more than 1, pass a vector of names
#' @param method whether to take mean (default) or median
#' @param subclusterpower whether to get multiple averages per
#'  original cluster
#' @param if_log input data is natural log,
#' averaging will be done on unlogged data
#' @param ... additional arguments
#' @export
seurat_ref.Seurat <- function(
    seurat_object,
    cluster_col = "classified",
    var_genes_only = FALSE,
    assay_name = NULL,
    method = "mean",
    subclusterpower = 0,
    if_log = TRUE,
    ...) {
    if (is(seurat_object, "Seurat")) {
        temp_mat <- object_data(seurat_object, "data")

        if (is.logical(var_genes_only) && var_genes_only) {
            temp_mat <- temp_mat[object_data(seurat_object, "var.genes"), ]
        } else if (var_genes_only == "PCA") {
            temp_mat <-
                temp_mat[rownames(object_data(seurat_object, "pca")), ]
        }

        if (!is.null(assay_name)) {
            og_assay <- SeuratObject::DefaultAssay(seurat_object)
            assay_name <- setdiff(assay_name, og_assay)
            temp_mat <- temp_mat[0, ]
            for (element in assay_name) {
                SeuratObject::DefaultAssay(seurat_object) <- element
                temp_mat2 <- object_data(seurat_object, "data", warn = FALSE)
                temp_mat <- rbind(temp_mat, as.matrix(temp_mat2))
            }
            SeuratObject::DefaultAssay(seurat_object) <- og_assay
        }
    } else {
        stop("Input is not a compatible Seurat object")
    }

    temp_res <- average_clusters(
        temp_mat,
        object_data(seurat_object, "meta.data"),
        cluster_col = cluster_col,
        method = method,
        subclusterpower = subclusterpower,
        if_log = if_log
    )

    temp_res
}

#' Function to convert labelled seurat object to fully prepared metadata
#' @return dataframe of metadata, including dimension reduction plotting info
#' @examples
#' so <- so_pbmc()
#' m <- seurat_meta(so)
#' @export
seurat_meta <- function(seurat_object, ...) {
    UseMethod("seurat_meta", seurat_object)
}

#' @rdname seurat_meta
#' @param seurat_object seurat_object after tsne or
#'  umap projections and clustering
#' @param dr dimension reduction method
#' @param ... additional arguments
#' @export
seurat_meta.Seurat <- function(
    seurat_object,
    dr = "umap",
    ...) {
    dr2 <- dr

    mdata <- object_data(seurat_object, "meta.data")
    temp_col_id <- get_unique_column(mdata, "rn")

    temp_dr <-
        tryCatch(
            as.data.frame(seurat_object@reductions[[dr2]]@cell.embeddings),
            error = function(e) {
                message("cannot find dr info")
                return(NA)
            }
        )
    if (!is.data.frame(temp_dr)) {
        return(mdata)
    } else {
        temp_dr <- tibble::rownames_to_column(temp_dr, temp_col_id)
        temp_meta <- tibble::rownames_to_column(mdata, temp_col_id)
        temp <- dplyr::left_join(temp_meta, temp_dr, by = temp_col_id)
        if (tibble::has_rownames(temp)) {
            temp <- tibble::remove_rownames(temp)
        }
        return(tibble::column_to_rownames(temp, temp_col_id))
    }
}

#' Function to convert labelled object to avg expression matrix
#' @return reference expression matrix, with genes as row names,
#'  and cell types as column names
#' @export
object_ref <- function(input, ...) {
    UseMethod("object_ref", input)
}

#' @rdname object_ref
#' @param input object after tsne or umap projections and clustering
#' @param cluster_col column name where classified cluster names
#' are stored in seurat meta data, cannot be "rn"
#' @param var_genes_only whether to keep only var.genes in the
#' final matrix output, could also look up genes used for PCA
#' @param assay_name any additional assay data, such as ADT, to include.
#'  If more than 1, pass a vector of names
#' @param method whether to take mean (default) or median
#' @param lookuptable if not supplied, will look
#' in built-in table for object parsing
#' @param if_log input data is natural log,
#' averaging will be done on unlogged data
#' @param ... additional arguments
#' @examples
#' so <- so_pbmc()
#' object_ref(
#'     so,
#'     cluster_col = "seurat_clusters"
#' )
#' @export
object_ref.default <- function(
    input,
    cluster_col = NULL,
    var_genes_only = FALSE,
    assay_name = NULL,
    method = "mean",
    lookuptable = NULL,
    if_log = TRUE,
    ...) {
    if (!is(input, "seurat")) {
        input_original <- input
        temp <- parse_loc_object(
            input,
            type = class(input),
            expr_loc = NULL,
            meta_loc = NULL,
            var_loc = NULL,
            cluster_col = cluster_col,
            lookuptable = lookuptable
        )
        if (!(is.null(temp[["expr"]]))) {
            message("recognized object type - ", class(input))
        }
        input <- temp[["expr"]]
        metadata <- temp[["meta"]]
        query_genes <- temp[["var"]]
        if (is.null(cluster_col)) {
            cluster_col <- temp[["col"]]
        }
    }

    temp_mat <- input
    if (is.logical(var_genes_only) && var_genes_only) {
        temp_mat <- temp_mat[query_genes, ]
    }

    temp_res <- average_clusters(
        temp_mat,
        metadata,
        cluster_col = cluster_col,
        method = method,
        if_log = if_log
    )

    temp_res
}

#' @rdname object_ref
#' @export
object_ref.Seurat <- function(
    input,
    cluster_col = NULL,
    var_genes_only = FALSE,
    assay_name = NULL,
    method = "mean",
    lookuptable = NULL,
    if_log = TRUE,
    ...) {
    temp_mat <- object_data(input, "data")
    metadata <- object_data(input, "meta.data")
    query_genes <- object_data(input, "var.genes")
    if (is.null(cluster_col)) {
        message("please indicate metadata column containing cell identities")
    }

    if (is.logical(var_genes_only) && var_genes_only) {
        temp_mat <- temp_mat[query_genes, ]
    }

    temp_res <- average_clusters(
        temp_mat,
        metadata,
        cluster_col = cluster_col,
        method = method,
        if_log = if_log
    )

    temp_res
}

#' @rdname object_ref
#' @export
object_ref.SingleCellExperiment <- function(
    input,
    cluster_col = NULL,
    var_genes_only = FALSE,
    assay_name = NULL,
    method = "mean",
    lookuptable = NULL,
    if_log = TRUE,
    ...) {
    temp_mat <- object_data(input, "data")
    metadata <- object_data(input, "meta.data")
    if (is.null(cluster_col)) {
        message("please indicate metadata column containing cell identities")
    }

    temp_res <- average_clusters(
        temp_mat,
        metadata,
        cluster_col = cluster_col,
        method = method,
        if_log = if_log
    )

    temp_res
}