R/GlowwormSubpopulation.R

In Hannahglover/Glowworm: Mapping genes to single cell RNAseq

#' @title Generate ranked Glowworm scores for  specific subpopulation/s
#' @param GlowwormScaleOutput Output from GlowwormScale
#' @param TargetPopulations A vector of populations for ranking
#' @return This function takes the output from GlowwormGlobal and generates a rank score specific for the population/s of interest \cr \cr Exports a list containing:
#' \itemize{
#'     \item MetaData - A copy of the metadata
#'     \item GlowwormScaleOutput - A scaled data frame of n genes by n populations in metadata
#'     \item SumStats - A data frame containing Total enrichment score (Abundance), SD enrichment score (Specificity) and
#'     the population specific rank score.}
#'     \cr \cr Output list can be used to generate a stacked bar plot via GlowwormOutputStacked.
#' @export

GlowwormSubpopulation = function(GlowwormScaleOutput, TargetPopulations){
  OutputList = list()
  OutputList[["MetaData"]] = GlowwormScaleOutput[["MetaData"]]
  GeneData = GlowwormScaleOutput[["GlowwormScaleOutput"]]
  OutputList[["GlowwormScaleOutput"]] = GeneData
  GeneData$Pop = NULL
  GeneData$Broad = NULL
  GeneData_t = as.data.frame(t(GeneData))
  colnames(GeneData_t) = gsub("\\^.*", "", colnames(GeneData_t))

  Population = GeneData_t %>% dplyr::select(all_of(TargetPopulations))
  Population$PopulationSum = rowSums(Population)

  SumStats = as.data.frame(matrix(ncol = 1, nrow = dim(GeneData_t)[1]))
  row.names(SumStats) = row.names(GeneData_t)
  SumStats$Total = sqrt(rowSums(GeneData_t))
  SumStats$SD = apply(GeneData_t[,-ncol(GeneData_t)], 1, sd)
  SumStats$PopulationSum = Population$PopulationSum
  SumStats$RankScore = (SumStats$Total * SumStats$SD) * (SumStats$PopulationSum/SumStats$Total)
  SumStats = data.frame(SumStats[order(-SumStats$RankScore),])
  SumStats = subset(SumStats, SumStats$RankScore > 0)
  SumStats$V1 = NULL
  OutputList[["SumStats"]] = SumStats
  if(dim(SumStats)[1] == 0){stop("No genes enriched in specified subpopulation")}
  return(OutputList)
}