Plot QQ-plot using ggplot2

ggQQplot.R

ggQQplot <- function(...,
                     colVec = c("black", "orange", "red", "purple"),
                     sizeVec = rep(0.3, length(colVec)),
                     legendLabel = NULL,
                     title = "",
                     openRange = F,
                     sampling = 1,
                     signList = NULL) {
  if (!is.list(..1)) {
    pList <- list(...)
  } else {
    pList <- ..1
  }
  
  if (is.null(legendLabel)) {
    legendLabel = colVec
  }
  
  if (length(pList) > length(colVec)) {
    stop("pList needs to be shorter than colVec!!!")
  }
  
  # Check signList if provided
  if (!is.null(signList)) {
    if (length(signList) != length(pList)) {
      stop("signList must have the same length as pList!!!")
    }
    # Validate that each element in signList has same length as corresponding pList element
    for (i in 1:length(pList)) {
      if (length(signList[[i]]) != length(pList[[i]])) {
        stop(paste0(
          "signList[[",
          i,
          "]] must have the same length as pList[[",
          i,
          "]]!!!"
        ))
      }
    }
  }
  
  if (openRange) {
    pList <- lapply(pList, function(x)
      x[x > 0 | x < 1])
    # Also filter signList if provided
    if (!is.null(signList)) {
      signList <- lapply(1:length(pList), function(i) {
        keep <- pList[[i]] > 0 | pList[[i]] < 1
        signList[[i]][keep]
      })
    }
  }
  
  if (sampling < 1) {
    sampleIndices <- lapply(pList, function(x)
      sample(1:length(x), length(x) * sampling))
    pList <- lapply(1:length(pList), function(i)
      pList[[i]][sampleIndices[[i]]])
    # Also sample signList if provided
    if (!is.null(signList)) {
      signList <- lapply(1:length(signList), function(i)
        signList[[i]][sampleIndices[[i]]])
    }
  }
  
  # Determine axis labels based on whether signs are provided
  xLabel <- ifelse(
    !is.null(signList),
    "Empirical signed -log<sub>10</sub> (<i>P</i>)",
    "Empirical -log<sub>10</sub> (<i>P</i>)"
  )
  yLabel <- ifelse(
    !is.null(signList),
    "Observed signed -log<sub>10</sub> (<i>P</i>)",
    "Observed -log<sub>10</sub> (<i>P</i>)"
  )
  
  g <- ggplot() +
    ggtheme_tp(base_size = 10) +
    ggplot2::labs(title = title, x = xLabel, y = yLabel) +
    ggplot2::geom_abline(
      slope = 1,
      intercept = 0,
      col = "red",
      lty = 2
    ) +
    ggplot2::theme(
      axis.title.x = ggtext::element_markdown(size = 10),
      axis.title.y = ggtext::element_markdown(size = 10)
    )
  
  # Add horizontal line at 0 for signed plots
  if (!is.null(signList)) {
    g <- g + ggplot2::geom_hline(yintercept = 0,
                                 col = "gray50",
                                 lty = 2)
  }
  
  for (i in 1:length(pList)) {
    n <- length(pList[[i]])
    
    # Calculate empirical distribution
    xValue <- -log10(ppoints(n))
    
    # Calculate observed values
    yValue <- pList[[i]]
    yValue[yValue == 0] <- min(yValue[yValue > 0])
    yValue <- -log10(yValue)
    
    # Apply signs if provided
    if (!is.null(signList)) {
      sortOrder <- order(abs(yValue), decreasing = TRUE)
      yValue <- yValue[sortOrder]
      signs <- signList[[i]][sortOrder]
      yValue <- yValue * signs
      xValue <- xValue * signs
    } else {
      yValue <- sort(yValue, decreasing = TRUE)
    }
    
    # Create a proper data frame
    plot_data <- data.frame(x = xValue, y = yValue, group = legendLabel[i])
    
    g <- g + ggrastr::geom_point_rast(
      data = plot_data,
      mapping = aes(x = x, y = y, color = group),
      size = sizeVec[i]
    )
  }
  
  g <- g +
    ggplot2::scale_color_manual(limits = legendLabel,
                                values = colVec,
                                name = NULL)
  
  return(g)
}