R/q.rank.R

In takakoizumi/QuEST: Constructing Quantitative Environment Sensor using Transcriptomes

##' Visualize R-squared value distribution in gene-environment interaction
##'
##' @title Visualize R-squared value distribution in gene-environment interaction
##' @param x A data matrix (row: samples, col: genes).
##' @param y A vector of an environment in which the samples were collected.
##' @param method A string to specify the method of regression for calculating R-squared values.
##' "linear" (default), "quadratic" or "cubic" regression model can be specified.
##' @param lower.thr The lower threshold of R-squared value to be included in QuEST model (default: 0).
##' @param n.gene The number of genes to be included in QuEST model (default: ncol(x)).
##' @param upper.xlim The upper limitation of x axis (i.e., the number of genes) in the resulted figure (default: ncol(x)).
##' @importFrom ggplot2 ggplot aes geom_line geom_area xlim ylim xlab ylab
##' geom_hline geom_text geom_vline theme element_rect element_blank
##' @importFrom stats lm
##' @return A rank order plot
##' @examples
##' data(Pinus)
##' train <- q.clean(Pinus$train)
##' target <- Pinus$target
##' train <- q.sort(train, target)
##' q.rank(train, target, lower.thr = 0.5)
##' @author Takahiko Koizumi
##' @export
q.rank <- function(x, y, method = "linear", lower.thr = 0, n.gene = ncol(x), upper.xlim = ncol(x)){
  R2 <- NULL
  degree <- switch(method,
                   "linear" = 1,
                   "quadratic" = 2,
                   "cubic" = 3,
                   stop("Select the <method> linear, quadratic, or cubic")
  )

  if(lower.thr > 0 & n.gene != ncol(x)){
    stop("Don't specify <lower.thr> and <n.gene> at a time")
  }
  if(lower.thr < 0 | lower.thr > 1){
    stop("<lower.thr> should be within the range of 0-1")
  }
  if(n.gene < 0){
    stop("<n.gene> should not be a negative value")
  }else if(n.gene > ncol(x)){
    stop(paste("<n.gene> must not exceed", ncol(x), sep = " "))
  }


  ## calculate R2 values
  result <- rep(NA, ncol(x))
  for(i in 1:ncol(x)){
    result[i] <- summary(lm(y ~ poly(x[, i], degree = degree, raw = TRUE)))$r.squared
  }

  h <- data.frame(
    rank = rep(NA, ncol(x)),
    R2 = result
  )

  ## sort genes in descending order of R2
  h <- h[order(h$R2, decreasing = TRUE), ]
  h$rank <- 1:nrow(h)

  ## set the number of genes to be used
  if(lower.thr != 0){
    n.gene <- length(result[result >= lower.thr])
  }

  if(lower.thr > 0){
    options(warn = -1)
    ggplot() +
      geom_line(data = h, aes(x = rank, y = R2), colour = "black") +
      geom_area(data = h, aes(x = rank, y = R2), fill = "gray") +
      geom_area(data = h[1:n.gene, ], aes(x = rank, y = R2), fill = "red") +
      xlim(0, upper.xlim) +
      ylim(0, 1) +
      xlab("Rank") +
      ylab("R-squared value") +
      geom_hline(aes(yintercept = lower.thr), colour = "black", linetype = "dotted") +
      geom_text(aes(0, lower.thr, label = format(round(lower.thr, 2), nsmall = 2), hjust = -0.2, vjust = 1.3)) +
      geom_vline(aes(xintercept = n.gene), colour = "black", linetype = "dotted") +
      geom_text(aes(n.gene, 0.03, label = n.gene, hjust = 1.3)) +
      theme(panel.background = element_rect(fill = "transparent", colour = "black"),
            panel.grid = element_blank(),
            strip.background = element_blank(),
            legend.position = "none")
  }else{
    R2val <- h[n.gene, "R2"]
    options(warn = -1)
    ggplot() +
      geom_line(data = h, aes(x = rank, y = R2), colour = "black") +
      geom_area(data = h, aes(x = rank, y = R2), fill = "gray") +
      geom_area(data = h[1:n.gene, ], aes(x = rank, y = R2), fill = "red") +
      xlim(0, upper.xlim) +
      ylim(0, 1) +
      xlab("Rank") +
      ylab("R-squared value") +
      geom_hline(aes(yintercept = R2val), colour = "black", linetype = "dotted") +
      geom_text(aes(0, R2val, label = format(round(R2val, 2), nsmall = 2), hjust = -0.2, vjust = 1.3)) +
      geom_vline(aes(xintercept = n.gene), colour = "black", linetype = "dotted") +
      geom_text(aes(n.gene, 0.03, label = n.gene, hjust = 1.3)) +
      theme(panel.background = element_rect(fill = "transparent", colour = "black"),
            panel.grid = element_blank(),
            strip.background = element_blank(),
            legend.position = "none")
  }
}