R/Roe.R

In zhupinglab/Aodav: Advanced omics data analysis and visualization

#' @title Roe
#' @description Calculate Roe score for each cell cluster, Ro/e, Chi-square test
#' @param data meta data of Seurat object, e.g. obj@meta.data
#' @param condition feature to facet cell source
#' @param cellType feature to annoate cell cluster
#' @param samples feature to group samples
#' @param ctrl the value of control group
#' @examples
#' \dontrun{
#' Roe(data, condition, cellType, samples)
#' }
#' @export
#' @import Seurat

Roe <- function(data, condition, cellType, samples, ctrl="Ctrl"){
  cellType.class <- unique(data[[cellType]][ data[[samples]] == ctrl])
  samples.class <- unique(data[[samples]])
  condition.class <- unique(data[[condition]])
  this.df <- c()
  # O/E
  for(cdt in condition.class){
    cdt.data <- data[ data[[condition]] == cdt & data[[cellType]] %in% cellType.class,]
    ctrl.total <- nrow(cdt.data[ cdt.data[[samples]] == ctrl,])
    for(ctyp in cellType.class){
      ctrl.ctyp <- nrow(cdt.data[ cdt.data[[samples]] == ctrl & cdt.data[[cellType]] == ctyp, ])
      ratio.expected <- ctrl.ctyp / ctrl.total
      for(ptt in samples.class){
        if(ptt != ctrl){
          ptt.total <- nrow(cdt.data[ cdt.data[[samples]] == ptt, ])
          ptt.ctyp <- nrow(cdt.data[ cdt.data[[samples]] == ptt & cdt.data[[cellType]] == ctyp, ])
          ptt.O2E <- ptt.ctyp / (ptt.total * ratio.expected)
          ppt.chisq.pvalue <- chisq.test(
            matrix(
              c(ctrl.ctyp,
                ctrl.total - ctrl.ctyp,
                ptt.ctyp,
                ptt.total - ptt.ctyp),
              nrow=2
            ))$p.value
          this.df <- rbind(this.df, c(cdt, ctyp, ptt, ptt.O2E, ppt.chisq.pvalue))
        }
      }
    }
  }
  this.df <- as.data.frame(this.df, stringsAsFactors=F)
  colnames(this.df) <- c("condition", "cellType", "samples", "O2E", "chisq.p")
  this.df$O2E <- as.numeric(this.df$O2E)
  this.df$chisq.p <- as.numeric(this.df$chisq.p)
  this.estimate <- c()
  for(cdt in condition.class){
    for(ctyp in cellType.class){
      this.O2E <- as.numeric(this.df$O2E[this.df$condition==cdt & this.df$cellType==ctyp])
      this.O2E <- this.O2E[!is.na(this.O2E)]
      this.mean <- mean(this.O2E)
      this.sem <- sd(this.O2E, na.rm=T)/sqrt(length(this.O2E))
      this.estimate <- rbind(this.estimate, c(cdt, ctyp, this.mean, this.sem))
    }
  }
  this.estimate <- as.data.frame(this.estimate, stringsAsFactors=F)
  colnames(this.estimate) <- c("condition", "cellType", "meanO2E", "sem")
  this.estimate$meanO2E <- as.numeric(this.estimate$meanO2E)
  this.estimate$sem <- as.numeric(this.estimate$sem)
  return(list("samples"=this.df, "group"=this.estimate))
}