R/ggplotScaleFreeFit.R

In twesleyb/Uezu2016: iBioID Proteomics from Uezu, Kanak, Bradshaw et al., 2016

#' ggplotScaleFreeFit
#'
#' evaluate the scale free fit of a graph - adapted from WGCNA function
#'
#' @param connectivity - degree centrality vector
#'
#' @return none
#'
#' @author Tyler W Bradshaw, \email{tyler.w.bradshaw@duke.edu}
#'
#' @references none
#'
#' @keywords none
#'
#' @examples
#' ggplotScaleFreePlot(connectivity)
ggplotScaleFreeFit <- function(connectivity, nBreaks = 10, truncated = FALSE,
                               removeFirst = FALSE, main = "", ...) {

  ## NOTE: adapted directly from WGCNA::scaleFreeFitIndex(k, nBreaks = 10, removeFirst = FALSE)
  # https://www.rdocumentation.org/packages/WGCNA/versions/1.69/topics/scaleFreeFitIndex

  suppressPackageStartupMessages({
    require(ggplot2)
  })
  k <- connectivity
  discretized.k <- cut(k, nBreaks)
  dk <- tapply(k, discretized.k, mean)
  p.dk <- as.vector(tapply(k, discretized.k, length) / length(k))
  breaks1 <- seq(from = min(k), to = max(k), length = nBreaks + 1)
  hist1 <- suppressWarnings(hist(k,
    breaks = breaks1, equidist = FALSE,
    plot = FALSE, right = TRUE
  )) # ...
  dk2 <- hist1$mids
  dk <- ifelse(is.na(dk), dk2, dk)
  dk <- ifelse(dk == 0, dk2, dk)
  p.dk <- ifelse(is.na(p.dk), 0, p.dk)
  log.dk <- as.vector(log10(dk))
  if (removeFirst) {
    p.dk <- p.dk[-1]
    log.dk <- log.dk[-1]
  }
  log.p.dk <- as.numeric(log10(p.dk + 1e-09))
  lm1 <- lm(log.p.dk ~ log.dk)
  pvalue <- normalp::lmp(lm1)
  title <- paste0(
    main, "R2 =", as.character(round(summary(lm1)$adj.r.squared, 2)),
    ", slope =", round(lm1$coefficients[[2]], 2)
  )
  OUTPUT <- data.frame(
    scaleFreeRsquared = round(summary(lm1)$adj.r.squared, 2),
    slope = round(lm1$coefficients[[2]], 2)
  )
  # generate ggplot
  df <- as.data.frame(cbind(log.dk, log.p.dk))
  plot <- ggplot(df, aes(x = log.dk, y = log.p.dk))
  plot <- plot + geom_point(size = 2)
  plot <- plot +  ggtitle(title)
  plot <- plot + geom_abline(intercept = coef(lm1)[1], 
			     slope = coef(lm1)[2], 
			     color = "black", 
			     linetype = "dashed")
  plot <- plot + labs(y = expression(Log[10](p(k))))
  plot <- plot + labs(x = expression(Log[10](k)))
  plot <- plot + 
	  theme(plot.title = element_text(hjust = 0.5, 
					  color = "black", 
					  size = 11, face = "bold"))
  plot <- plot + 
	  theme(axis.title.x = element_text(hjust = 0.5, 
					    color = "black", 
					    size = 11))
  plot <- plot + 
	  theme(axis.title.y = element_text(hjust = 0.5, 
					    color = "black", 
					    size = 11))
  plot <- plot + theme(panel.background = element_blank())
  plot <- plot + theme(axis.line.x=element_line())
  plot <- plot + theme(axis.line.y=element_line())
  return(plot)
}