R/runKappa.R

In xlucpu/MOVICS: Multi-Omics integration and VIsualization in Cancer Subtyping

#' @name runKappa
#' @title Run consistency evaluation using Kappa statistic
#' @description Calculate Kappa statistic to measure the consistency between two appraisements
#' @param subt1 A numeric vector to indicate the first appraisement. Order should be exactly the same with subt2 for each sample.
#' @param subt2 A numeric vector to indicate the second appraisement. Order should be exactly the same with subt1 for each sample.
#' @param subt1.lab A string value to indicate the label of the first subtype.
#' @param subt2.lab A string value to indicate the label of the second subtype.
#' @param fig.path A string value to indicate the output path for storing the consistency heatmap.
#' @param fig.name A string value to indicate the name of the consistency heatmap.
#' @param width A numeric value to indicate the width of output figure.
#' @param height A numeric value to indicate the height of output figure.
#' @details This function evaluates the consistency between two appraisements that targets to the same cohort.
#' For example, the NTP-predicted subtype amd PAM-predicted subtype of external cohort, or the current subtype and predicted subtype of discovery cohort.
#' Therefore, the arguments `subt1` and `subt2` can be the `clust` column of `clust.res` derived from `getMOIC()` with one specified algorithm
#' or `get\%algorithm_name\%` or `getConsensusMOIC()` with a list of multiple algorithms or `runNTP()` or `runPAM()`.
#' However, subtypes identified from different algorithm (i.e., `get\%algorithm_name1\%` and `get\%algorithm_name2\%`) can not be evaluated
#' because the subtype 1 identified from the first algorithm may not be the same subtype 1 from the second algorithm.
#' @return A figure of consistency heatmap (.pdf).
#' @export
#' @importFrom grDevices dev.copy2pdf
#' @examples # There is no example and please refer to vignette.
runKappa <- function(subt1     = NULL,
                     subt2     = NULL,
                     subt1.lab = NULL,
                     subt2.lab = NULL,
                     fig.path  = getwd(),
                     fig.name  = "constheatmap",
                     width     = 5,
                     height    = 5) {

  subt1 <- as.vector(as.character(subt1))
  subt2 <- as.vector(as.character(subt2))

  # check sample
  if(length(subt1) != length(subt2)) {
    stop("subtypes identified from different cohorts.")
  }

  # check subtype
  if(!identical(sort(unique(subt1)),sort(unique(subt2)))) {
    stop("subtypes fail to matched from two appraisements.")
  }

  # check argument
  if(is.null(subt1.lab) | is.null(subt2.lab)) {
    stop("label for subtype1 and subtype2 must be both indicated.")
  }

  # calculate consistency
  comb.subt <- data.frame(subt1 = paste0("CS", subt1),
                          subt2 = paste0("CS", subt2),
                          stringsAsFactors = F)

  tab_classify <- as.data.frame.array(table(comb.subt$subt1,comb.subt$subt2))

  # n.moic <- length(unique(moic.res$clust.res$clust))
  # if(n.moic < 5) {
  #   p.fisher <- fisher.test(tab_classify, workspace = 1e9)$p.value
  # } else {
  #   message("--using simulated p value in fisher's exact test.")
  #   p.fisher <- fisher.test(tab_classify, simulate.p.value = T)$p.value
  # }

  # calculate kappa statistic
  x <- table(comb.subt$subt1,comb.subt$subt2)
  nr <- nrow(x); nc <- ncol(x); N <- sum(x)
  Po <- sum(diag(x))/N; Pe <- sum(rowSums(x) * colSums(x)/N)/N
  kappa <- (Po - Pe)/(1 - Pe)
  seK0 <- sqrt(Pe/(N * (1 - Pe)))
  p.v <- 1 - pnorm(kappa/seK0)
  p.lab <- ifelse(p.v < 0.001, "P < 0.001", paste0("P = ", format(round(p.v,3), nsmall = 3)))

  # generate consistency table
  blue   <- "#204F8D"
  lblue  <- "#498EB9"
  dwhite <- "#B6D1E8"
  white  <- "#E6EAF7"

  par(bty="n", mgp = c(2,0.5,0), mar = c(4.1,4.1,4.1,2.1),tcl=-.25, font.main=3)
  par(xpd=NA)
  plot(c(0,ncol(tab_classify)),c(0,nrow(tab_classify)),
       col = "white",
       xlab = "",xaxt = "n",
       ylab = "",yaxt = "n")
  title(paste0("Consistency between ",subt1.lab," and ",subt2.lab,"\nKappa = ", format(round(kappa,3), nsmall = 3),
               "\n", p.lab),
        adj = 0, line = 0)

  # add y-axis
  axis(2, at = 0.5:(nrow(tab_classify)-0.5), labels = FALSE)
  text(y = 0.5:(nrow(tab_classify)-0.5),
       par("usr")[1],
       labels = rownames(tab_classify)[nrow(tab_classify):1],
       srt = 0, pos = 2, xpd = TRUE)
  mtext(paste0("Subtypes derived from ", subt1.lab), side = 2, line = 3)

  # add x-axis
  axis(1, at = 0.5:(ncol(tab_classify)-0.5), labels = FALSE)
  text(x = 0.5:(ncol(tab_classify)-0.5),
       par("usr")[1] - 0.2,
       labels = colnames(tab_classify),
       srt = 45, pos = 1, xpd = TRUE)
  mtext(paste0("Subtypes derived from ", subt2.lab), side = 1, line = 3)

  # generate colors
  input_matrix <- as.matrix(tab_classify)
  mat.max = max(input_matrix)
  unq.value <- unique(sort(as.vector(input_matrix)))
  rbPal <- colorRampPalette(c(white,dwhite,lblue,blue))
  col.vec <- rbPal(max(unq.value) + 1)
  col.mat <- matrix(NA,byrow = T,ncol = ncol(input_matrix),nrow = nrow(input_matrix))

  # fill matrix
  for (i in 1:nrow(col.mat)) {
    for (j in 1:ncol(col.mat)) {
      col.mat[i,j] <- col.vec[input_matrix[i,j] + 1]
    }
  }

  # generate heatmap
  x_size <- ncol(input_matrix)
  y_size <- nrow(input_matrix)

  my_xleft = rep(c(0:(x_size-1)),each = x_size)
  my_xright = my_xleft + 1
  my_ybottom = rep(c((y_size-1):0),y_size)
  my_ytop = my_ybottom + 1
  rect(xleft = my_xleft,
       ybottom = my_ybottom,
       xright = my_xright,
       ytop = my_ytop,
       col=col.mat,
       border = F)

  # fill count
  text(my_xleft + 0.5,my_ybottom + 0.5,input_matrix, cex = 1.3)

  # output to pdf
  outFig <- paste0(fig.name,".pdf")
  invisible(dev.copy2pdf(file = file.path(fig.path, outFig), width = width, height = height))
}