R/scIntegrate.R

In ncrna/SCKIT: Yet Another Single-Cell Analysis Toolkit (Yeskit)

# Set a default value if the left value is null @param lhs Default value to use
# @param rhs The value to use if lhs is null @author Hadley Wickham @references
# https://adv-r.hadley.nz/functions.html#missing-arguments
#' @export

`%||%` <- function(lhs, rhs) {
    if (!is.null(lhs)) {
        lhs
    } else {
        rhs
    }
}

#' Integrate seurat objects into one
#' @title scIntegrate
#' @param object.list A list of seurat objects
#' @param object.names An array of seurat object names
#' @param nVariable Number of features to select as top variable features
#' @param nPC Number of PCs to use
#' @param res Resolution parameter to set. Default res=0.5
#' @param batch.rm Remove batch effect with 'seurat' or 'harmony'. 
#'   Default batch.rm='harmony'
#' @return Seurat object.

#' @author rstatistics
#' @export

scIntegrate <- function(object.list = NULL, object.names = NULL, 
                        nVariable = 2000, nPC = 30, res = 0.5, 
                        batch.rm = "harmony") {
  # Merge each seurat object
  meta.list <- list()
  meta.table <- object.list[[1]]@meta.data
  for (i in 2:length(object.list)) {
    meta.table <- dplyr::bind_rows(meta.table, object.list[[i]]@meta.data)
  }
  meta.table[is.na(meta.table)] <- 0

  object <- merge(object.list[[1]], y = object.list[2:length(object.list)], 
    add.cell.ids = as.character(object.names), project = "seurat")
  rownames(meta.table) <- rownames(object@meta.data)
  object@meta.data <- meta.table

  object <- Seurat::NormalizeData(object, verbose = TRUE)
  object <- Seurat::FindVariableFeatures(object, selection.method = "vst", 
    nfeatures = nVariable)
  object <- Seurat::ScaleData(object, vars.to.regress = c("nCount_RNA", 
    "percent.mito"), verbose = TRUE)
  object <- Seurat::RunPCA(object, pc.genes = Seurat::VariableFeatures(object),
    npcs = nPC, verbose = TRUE)

  object <- harmony::RunHarmony(object = object, group.by.vars = "sample", 
    assay.use = Seurat::DefaultAssay(object), plot_convergence = TRUE)

  if (batch.rm == "harmony") {
    object <- Seurat::FindNeighbors(object, reduction = "harmony", 
      dims = seq_len(nPC))
    object <- Seurat::FindClusters(object, resolution = res)
    object <- Seurat::RunUMAP(object, reduction = "harmony", 
      dims = seq_len(nPC))
    object <- Seurat::RunTSNE(object, reduction = "harmony", 
      dims = seq_len(nPC), do.fast = TRUE)
  } else if (batch.rm == "seurat") {
    object <- Seurat::FindNeighbors(object, reduction = "pca", 
      dims = seq_len(nPC))
    object <- Seurat::FindClusters(object, resolution = res)
    object <- Seurat::RunUMAP(object, reduction = "pca", dims = seq_len(nPC))
    object <- Seurat::RunTSNE(object, reduction = "pca", dims = seq_len(nPC))
  } else {
    stop("batch.rm must be 'harmony' or 'seurat'!")
  }

  cols <- NA
  if (length(levels(Seurat::Idents(object))) <= 36) {
    cols = c("#1660A7", "#FF6A00", "#219418", "#CD0C18", "#814BB2", "#794339",
            "#DC59B6", "#CC79A7", "#FF0000", "#11B3C6", "#AFB400", "#00FFFF", 
            "#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", 
            "#D55E00", "#D55E00", "#CC79A7", "#00AFBB", "#E69F00", "#009E73", 
            "#56B4E9", "#F0E442", "#0072B2", "#D55E00", "#CC79A7", "#4477AA", 
            "#EE6677", "#228833", "#CCBB44", "#66CCEE", "#AA3377", "#BBBBBB")
    object@misc$cols = cols
  }
  # Add reduction coordinates into Meta.data slot
  reductions <- intersect(c("pca", "tsne", "umap"), names(object))
  for (reduction in reductions) {
    meta_ids <- gsub("coord", toupper(reduction), c("coord_1", "coord_2"))
    coord <- Seurat::Embeddings(object = object, reduction = reduction)[, c(1, 2)]
    object <- Seurat::AddMetaData(object = object, metadata = coord, 
      col.name = meta_ids)
  }
  return(object)
}