R/population.pairing_CellComm.R

In edroaldo/fusca: Framework for Unified Single-Cell Analysis

#' Compute cell-cell interaction network.
#'
#' Analyze the mean expression data generated by the computeValue method of
#' CellRouter to determine co-expression patterns of ligands and receptors.
#'
#' @param mean.expr list; the statistics of the expressed genes as calculated by
#' the computeValue method of CellRouter.
#' @param ligands vector; list of genes annotated as ligands.
#' @param receptors vector; list of genes annotated as receptors.
#' @param threshold numeric; threshold to select the genes in the population
#' according to their mean expression.
#' @param pairs the pairs;
#'
#' @return data frame; interactions with celltype1, celltype2, pair, ligand,
#' and receptor.
#'
#' @export
population.pairing <- function(mean.expr, ligands, receptors, threshold, pairs){
  interactions <- list()
  for(ct1 in names(mean.expr)){
    for(ct2 in names(mean.expr)){
      ct1.expressed <- mean.expr[[ct1]][c(ligands, receptors),]
      ct1.expressed <- ct1.expressed[which(ct1.expressed$p > threshold),]

      ct2.expressed <- mean.expr[[ct2]][c(ligands, receptors),]
      ct2.expressed <- ct2.expressed[which(ct2.expressed$p > threshold),]

      combs <- expand.grid(rownames(ct1.expressed), rownames(ct2.expressed))
      combs <- paste(combs$Var1, combs$Var2, sep='_')
      combs <- intersect(combs, pairs$Pair.Name)
      if(length(combs) > 0){
        ligs <- sapply(strsplit(combs, split='_', fixed=TRUE), function(x) (x[1]))
        recs <- sapply(strsplit(combs, split='_', fixed=TRUE), function(x) (x[2]))
        df <- data.frame(celltype1=rep(ct1, times=length(combs)),
                         celltype2=rep(ct2, times=length(combs)),
                         pair=combs, ligand=ligs, receptor=recs,
                         stringsAsFactors = FALSE)
        interactions[[paste(ct1,ct2,sep='_')]] <- df
      }

    }
  }
  interactions <- dplyr::bind_rows(interactions)#rbind.fill(interactions)

  return(interactions)
}