R/cellMap.R

In zhupinglab/Aodav: Advanced omics data analysis and visualization

#' @title cellMap
#' @description Map single cells to each cell clusters. This functions embeds 'scmap' tool and improves to enable multiple rounds of projection to increase the mapping accuracy of cells with similar experssion patterns.
#' @param ref Seurat object with defined clusters stored in ref@ident as a reference.
#' @param data Seurat object with cells to be mapped to the reference.
#' @param map list of clusters to map accordingly. Defualt is NULL to use all clusters in ref@ident.
#' @param pvalue p value to filter genes used for index construction.
#' @param adjP adjusted p value to filter genes used for index construction only if pvalue is NULL.
#' @param fc fold change of expression to filter genes used for index construction.
#' @param threshold similarity threshould for mapping cells to clusters.
#' @param method strategy for gene selection. Default is "pwDEG". In case few genes were selected, the "selectFeatures" function in scmap is also implemented to select 200 genes.
#' @examples
#' \dontrun{
#' round1 <- c("CD4-C1-CCR7", "CD4-C2-SELL", "CD4-C3-CD82", "CD4-C4-GZMH")
#' round2 <- c("CD4-C1-CCR7", "CD4-C2-SELL")
#' projection <- cellMap(ctrl.cd4T.obj, ctrl.cd4T.obj, map=list(round1, round2), threshold=0.3)
#' table(ctrl.cd4T.obj@ident)
#' table(projection$round1$scmap_cluster_labs)
#' table(projection$round2$scmap_cluster_labs)
#' }
#' @export
#' @import SingleCellExperiment
#' @import scmap

cellMap <- function(ref, data, map=NULL, pvalue=NULL, adjP=0.05, fc=0.25, threshold=0.7, method="pwDEG"){
  assignment <- list()
  if(is.null(map)){
    map <- list("roud1"=levels(ref@ident))
  }
  for(round in 1:length(map)){
    print(paste0("Projection round", round, " for cluster: ", paste0(map[[round]], collapse=","), " ..."))
    # subset data for each round
    this.ref <- SubsetData(ref, cells.use=row.names(ref@meta.data)[ ref@ident %in% map[[round]]])
    if( round == 1){
      this.data <- data
    } else {
      pre.cells <- assignment[[paste0("round", round - 1)]]$cells
      this.data <- SubsetData(data, cells.use=pre.cells[assignment[[paste0("round", round - 1)]]$scmap_cluster_labs %in% map[[round]]])
    }
    this.ref@meta.data$cellGroup <- this.ref@ident

    # single cell experiment object
    this.ref.sce <- SingleCellExperiment(
      assays = list("normcounts"=as.matrix(this.ref@raw.data[, row.names(this.ref@meta.data)])),
      colData = this.ref@meta.data
    )
    this.data.sce <- SingleCellExperiment(
      assays = list("normcounts"=as.matrix(this.data@raw.data[, row.names(this.data@meta.data)])),
      colData = this.data@meta.data
    )
    logcounts(this.ref.sce) <- log2(normcounts(this.ref.sce) + 1)
    logcounts(this.data.sce) <- log2(normcounts(this.data.sce) + 1)
    rowData(this.ref.sce)$feature_symbol <- rownames(this.ref.sce)
    rowData(this.data.sce)$feature_symbol <- rownames(this.data.sce)
    # find differentially expressed genes for index construction
    if(method == "pwDEG"){
      print(paste0("Find pwDEGs for clusters ..."))
      this.pwDEGs <- pwDEG(this.ref)
      if(is.null(pvalue)){
        print(paste0("Using adjusted p value <= ", adjP, " to filter genes"))
        this.pwDEGs <- unique(this.pwDEGs$gene[ this.pwDEGs$p_val_adj <= adjP & abs(this.pwDEGs$avg_logFC) >= fc])
      } else {
        print(paste0("Using p value <= ", pvalue, " to filter genes"))
        this.pwDEGs <- unique(this.pwDEGs$gene[ this.pwDEGs$p_val <= pvalue & abs(this.pwDEGs$avg_logFC) >= fc])
      }
      rowData(this.ref.sce)$scmap_features <- ifelse(row.names(this.ref.sce) %in% this.pwDEGs, TRUE, FALSE)
    } else {
      this.ref.sce <- selectFeatures(this.ref.sce, suppress_plot = TRUE, n_features = 200)
    }
    print(paste0(sum(rowData(this.ref.sce)$scmap_features == TRUE), " genes are selected to build index."))

    # construct index
    this.ref.sce <- indexCluster(this.ref.sce, cluster_col = "cellGroup")
    this.ref.idx <- metadata(this.ref.sce)$scmap_cluster_index
    this.ref.idx <- this.ref.idx[ rowSums(this.ref.idx) != 0,]
    print(paste0(nrow(this.ref.idx), " genes are retained in index."))

    # projection
    this.data.proj <- scmapCluster(
      projection = this.data.sce,
      threshold = threshold,
      index_list = list(ctrl = this.ref.idx)
    )
    this.data.proj$cells <- row.names(this.data@meta.data)
    assignment[[paste0("round", round)]] <- this.data.proj
  }

  assignment$merge$cells <- assignment$round1$cells
  assignment$merge$scmap_cluster_labs <- assignment$round1$scmap_cluster_labs
  for(round in 1:length(map)){
    if(round > 1){
      this.assign <- assignment[[paste0("round", round)]]
      cells.idx <- match(this.assign$cells, assignment$merge$cells)
      assignment$merge$scmap_cluster_labs[cells.idx] <- this.assign$scmap_cluster_labs
    }
  }
  print("Cell map done!")
  return(assignment)
}