R/2-3.RMT.R

In MetaNet: Network Analysis for Omics Data

# =========2.2RMT optimize=====

#' Get RMT threshold for a correlation matrix
#'
#' @param occor.r a corr object or a correlation matrix
#' @param min_threshold min_threshold
#' @param max_threshold max_threshold
#' @param step step
#' @param gif render a .gif file?
#' @param verbose verbose
#' @param out_dir output dir
#'
#' @return a r-threshold
#' @export
#' @references
#' J. Zhou, Y. Deng, FALSE. Luo, Z. He, Q. Tu, X. Zhi, (2010) Functional Molecular Ecological Networks, doi:10.1128/mBio.00169-10.
#' <https://matstat.org/content_en/RMT/RMThreshold_Intro.pdf>
#' @examples
#' \donttest{
#' data(otutab, package = "pcutils")
#' t(otutab) -> totu
#' c_net_calculate(totu) -> corr
#' rmt(corr)
#' # recommend: 0.69
#' c_net_build(corr, r_threshold = 0.69) -> co_net_rmt
#' }
RMT_threshold <- function(occor.r, out_dir, min_threshold = 0.5, max_threshold = 0.8,
                          step = 0.02, gif = FALSE, verbose = FALSE) {
  nwd <- getwd()
  on.exit(setwd(nwd))

  setwd(out_dir)
  if (inherits(occor.r, "corr")) occor.r <- occor.r$r
  if (!dir.exists("./RMT_temp")) dir.create("./RMT_temp")
  diag(occor.r) <- 0

  if (max_threshold >= max(abs(occor.r))) max_threshold <- (max(abs(occor.r)) - step)
  if (min_threshold >= max_threshold) min_threshold <- max_threshold - 10 * step

  thres_seq <- seq(min_threshold, max_threshold, step)

  res <- data.frame()
  for (i in seq_len(length(thres_seq))) {
    threshold <- thres_seq[i]
    if (!verbose) pcutils::dabiao(paste0("Calculating", i, ":  threshold =", signif(threshold, 3)), print = TRUE)
    corr_r1 <- occor.r
    corr_r1[abs(corr_r1) < threshold] <- 0
    # calculate eigenvalues
    rand.mat <- corr_r1
    eigenvalues <- eigen(rand.mat, only.values = TRUE)$values
    eigenvalues <- eigenvalues[order(eigenvalues)] / max(abs(eigenvalues))
    eigenvalues <- pcutils::remove.outliers(unique(eigenvalues))

    # get the NNDS
    { # uf <- rm.unfold.gauss(eigenvalues,pop.up = TRUE)
      dens <- density(eigenvalues, kernel = "gaussian")
      midpoints <- \(x)(x[-length(x)] + 0.5 * diff(x))
      scale.function <- approx(dens$x, dens$y, xout = midpoints(eigenvalues))
      ev.spacing <- diff(eigenvalues)
      ev.spacing <- ev.spacing * scale.function$y
      ev.spacing <- ev.spacing / mean(ev.spacing)
    }

    ev.spacing <- ev.spacing[ev.spacing <= 3]
    # test whether fit possion?
    p_ks_test <- ks.test(unique(ev.spacing), "pexp", 1)$p.value
    # get sse
    # sse = rm.sse(ev.spacing)
    sse <- get_sse(ev.spacing)
    log_sse <- log(sse)

    # maximum likelihood
    evs <- ev.spacing[ev.spacing != 0]
    N <- length(evs)
    log_LE <- -sum(evs) / N
    log_LW <- log(pi / 2) + sum(log(evs)) / N - 0.25 * pi * sum(evs^2) / N

    # save png
    {
      histo <- hist(ev.spacing, breaks = seq(min(ev.spacing), max(ev.spacing), len = 51), plot = FALSE)
      grDevices::png(paste0("RMT_temp/rmt_nnsd", i, ".png"), height = 600, width = 700, res = 130)
      nnsd_plot(
        histo = histo, title = "Eigenvalue spacing distribution (NNSD)", threshold = threshold,
        dis_GOE = log_LW, dis_possion = log_LE, p_ks_test = p_ks_test
      )
      grDevices::dev.off()
    }
    res <- rbind(res, data.frame(threshold, p_ks_test, log_sse, log_LW, log_LE))
  }
  message("The Intermediate files saved in ", out_dir, "/RMT_temp/.")
  # transfer to gif
  if (gif) {
    lib_ps("gifski", library = FALSE)
    gifski::gifski(paste0("RMT_temp/rmt_nnsd", seq_len(length(thres_seq)), ".png"),
      gif_file = "RMT_temp/rmt_nnsd.gif"
    )
  }
  r_threshold <- (res[which(res$log_LW == min(res$log_LW)), "threshold"] +
    res[which(res$log_LE == max(res$log_LE)), "threshold"]) / 2
  res <- list(res = res, r_threshold = r_threshold)
  class(res) <- c("rmt_res", class(res))
  return(res)
}

#' Plot a rmt_res
#'
#' @param x rmt_res
#' @param ... Additional arguments
#'
#' @return ggplot
#' @exportS3Method
#' @method plot rmt_res
plot.rmt_res <- function(x, ...) {
  threshold <- value <- variable <- xi <- y <- NULL
  res <- x$res
  linedf <- data.frame(
    variable = c("p_ks_test", "log_sse", "log_LW", "log_LE"),
    xi = c(
      res[which(res$p_ks_test == max(res$p_ks_test))[1], "threshold"],
      res[which(res$log_sse == min(res$log_sse))[1], "threshold"],
      res[which(res$log_LW == min(res$log_LW))[1], "threshold"],
      res[which(res$log_LE == max(res$log_LE))[1], "threshold"]
    ),
    x = max(res$threshold) - min(res$threshold),
    y = apply(res[, -1], 2, max)
  )

  reshape2::melt(res, "threshold") -> md

  # filter(threshold<0.77)%>%
  p <- ggplot(md, aes(threshold, value)) +
    geom_point(aes(col = variable)) +
    geom_line(aes(col = variable)) +
    scale_color_manual(values = get_cols(4, "col1")) +
    facet_wrap(. ~ variable, scales = "free_y") +
    theme_bw() +
    xlab(NULL) +
    geom_text(data = linedf, aes(x = xi - 0.02 * x, y = 0.5 * y, label = xi)) +
    geom_vline(data = linedf, aes(xintercept = xi), linetype = 2, col = "red") +
    theme(legend.position = "none")

  message(paste("recommend r_threshold: ", mean(linedf$xi)))
  return(p)
}

nnsd_plot <- \(histo = histo, title = title, threshold = threshold,
  dis_GOE = dis_GOE, dis_possion = dis_possion, p_ks_test = p_ks_test) {
  plot(histo, freq = FALSE, col = "#F4FCA1", main = title, font.main = 1, xlab = "eigenvalue spacing", ylab = "PDF of eigenvalue spacing")
  {
    actual.ymax <- par("yaxp")[2]
    x0 <- -log(actual.ymax * 0.98)
    possion_dis <- \(x)exp(-x)
    graphics::curve(possion_dis,
      from = max(x0, min(histo$breaks)),
      to = max(histo$breaks), n = 1001, add = TRUE, type = "l", lty = 1, col = "#EB34FF", lwd = 2
    )
  }
  {
    GOE <- function(x) pi / 2 * x * exp(-pi / 4 * x^2)
    graphics::curve(GOE,
      from = min(histo$breaks),
      to = max(histo$breaks), n = 1001, add = TRUE, type = "l",
      lty = 1, col = "blue", lwd = 2
    )
  }

  if ((!is.na(dis_GOE)) && (!is.na(dis_possion))) {
    graphics::mtext(side = 3, paste(
      "Distance to GOE =", signif(dis_GOE, 3),
      "\nDistance to Possion =", signif(dis_possion, 3), "; ks_test p.value for possion =", signif(p_ks_test, 3)
    ), col = "#878787", cex = 0.6)
  }

  if (!is.na(threshold)) graphics::mtext(side = 4, paste("threshold =", signif(threshold, 4)))

  graphics::legend("topright", inset = 0.05, c("Possion", "GOE"), col = c("#EB34FF", "blue"), lty = 1, lwd = 2, cex = 0.8)
}

trapez <- \(x, y){
  ind <- 2:length(x)
  as.double((x[ind] - x[ind - 1]) %*% (y[ind] + y[ind - 1])) / 2
}

get_sse <- \(ev.spacing){
  dens <- density(ev.spacing)
  N <- 20
  x <- seq(min(ev.spacing), max(ev.spacing), len = 1000)
  A <- exp(-min(ev.spacing)) - exp(-max(ev.spacing))
  xs <- numeric(N + 1)
  xs[1] <- min(ev.spacing)
  for (i in 1:N) xs[i + 1] <- -log(exp(-xs[i]) - A / N)
  area <- numeric(N)
  for (i in 1:N) {
    xsec <- x[(x > xs[i]) & (x < xs[i + 1])]
    xsec <- c(xs[i], xsec, xs[i + 1])
    ysec <- approx(dens$x, dens$y, xout = xsec)$y
    area[i] <- trapez(xsec, ysec)
  }
  sse <- sum((area[i] - A / N)^2)
  sse
}

#' Get RMT threshold for a correlation matrix roughly
#'
#' @export
#' @return recommend threshold
#' @rdname RMT_threshold
rmt <- function(occor.r, min_threshold = 0.5, max_threshold = 0.85, step = 0.01) {
  if (inherits(occor.r, "corr")) occor.r <- occor.r$r
  NNSD <- \(x)abs(diff(x))

  s <- seq(0, 3, 0.1)
  poisson_d <- exp(-s)
  nnsdpois <- density(NNSD(poisson_d))

  ps <- c()
  threshold <- c()

  for (i in seq(min_threshold, max_threshold, step)) {
    corr_r1 <- occor.r
    corr_r1[abs(corr_r1) < i] <- 0
    {
      eigen_res <- sort(eigen(corr_r1)$value)
      # spline to eigen_res
      check <- tryCatch(ssp <- smooth.spline(eigen_res, control.spar = list(low = 0, high = 3)),
        error = \(e) {
          TRUE
        }
      )
      if (rlang::is_true(check)) next
      nnsdw <- density(NNSD(ssp$y))
      chival <- sum((nnsdw$y - nnsdpois$y)^2 / nnsdpois$y / 1e3)
    }

    ps <- c(ps, chival)
    threshold <- c(threshold, i)
    if (((i * 100) %% 5 == 0)) {
      message(paste0("Calculating: ", i))
    }
  }

  res <- data.frame(threshold, ps)
  recommend_thres <- res[which.min(res[, 2]), 1]
  p <- ggplot(res, aes(threshold, ps)) +
    geom_point() +
    geom_vline(xintercept = recommend_thres, linetype = 2, col = "red") +
    geom_text(x = recommend_thres + 0.01, y = 0.5 * max(res$ps), label = recommend_thres) +
    theme_bw(base_size = 15)
  print(p)

  res1 <- res[(res$threshold < (recommend_thres + 0.05)) & (res$threshold > (recommend_thres - 0.05)), ]
  p <- ggplot(res1, aes(threshold, ps)) +
    geom_point() +
    geom_vline(xintercept = recommend_thres, linetype = 2, col = "red") +
    geom_text(x = recommend_thres + 0.01, y = 0.5 * max(res1$ps), label = recommend_thres) +
    theme_bw(base_size = 15)
  print(p)

  message("We recommend r-threshold: ", recommend_thres, ", you can calculate again in a smaller region")
  recommend_thres
}