R/functions_se.R

In szymczak-lab/harmonizeGeneExprData: Harmonization of gene expression data from GEO and ArrayExpress

#' Make SummarizedExperiment object
#'
#' Combine expression data, phenotype data, gene annotation and meta data into
#' SummarizedExperiment object including checks on samples and genes.
#'
#' @param pheno [data.frame] phenotype data
#' @param meta.data.l [list] meta data (as extracted by
#' \code{\link{extract_meta_data}})
#' @param assays [list] matrices of gene expression data, e.g. raw and voom
#' transformed counts for RNAseq (as prepared by
#' \code{\link{prepare_count_data}})) or raw and SCAN normalized expression
#' values for array data (as prepared by \code{\link{prepare_array_data}}))
#' @param anno [data.frame] gene annotation (e.g. gene identifier, chromosomal
#' position, gene names) with genes in rows and information in columns
#'
#' @return [SummarizedExperiment object] as generated by
#' \code{\link[SummarizedExperiment]{SummarizedExperiment}}
#' @export

make_se_object <- function(
  pheno,
  meta.data.l = NULL,
  assays,
  anno = NULL) {

  ## checks
  samples.a = colnames(assays[[1]])
  if (nrow(pheno) > length(samples.a)) {
    warning("some samples not in assays\n")
  }
  if (nrow(pheno) < length(samples.a)) {
    warning("some samples not in pheno\n")
  }
  samples.use = intersect(rownames(pheno),
                          samples.a)
  if (length(samples.use) == 0) {
    stop("no overlapping samples")
  }

  genes.a = rownames(assays[[1]])
  if (!is.null(anno)) {
    if (nrow(anno) < length(genes.a)) {
      warning("some genes not in anno")
    }
    genes.use = intersect(rownames(anno),
                          genes.a)
    if (length(genes.use) == 0) {
      stop("no overlapping genes")
    }

    anno = anno[genes.use, ]
  } else {
    genes.use = rownames(assays[[1]])
  }

  ## extract genes and samples
  assays = lapply(assays, function(x) {
    x[genes.use, samples.use]
  })

  se = SummarizedExperiment::SummarizedExperiment(
    assays = assays,
    rowData = anno,
    colData = pheno[samples.use, ],
    metadata = meta.data.l)

  return(se)
}