R/rsbplot.R

In cmcouto-silva/LD: Linkage Disequilibrium Analysis

rsbplot <- function (rsb.data, file.name, file.type = 'png', plot.pval = TRUE, ylim.scan = 2, pch = 16, cex = 0.3, 
                     cex.lab = 1.25, main = NA, cex.main = 1.5, cex.axis = 1) {
  
  substrRight <- function(x, n) substr(x, nchar(x)-n+1, nchar(x))
  
  if (!(grepl("Rsb", colnames(rsb.data)[3]))) {
    stop("Unrecognized column name for Rsb: the rsb.data is not a matrix generated by the ies2rsb() function...")
  }
  
  if (colnames(rsb.data)[4] == "-log10(p-value) [bilateral]") {
    yleg = expression("-" * log[10] * "[1" ~ "-" ~ "2" ~ "|" ~ 
                        Phi[scriptstyle(italic(Rsb))] ~ "-" ~ 0.5 * "|]")
  } else if (colnames(rsb.data)[4] == "-log10(p-value) [unilateral]") {
    yleg = expression("-" * log[10] * "[1" ~ "-" ~ Phi[scriptstyle(italic(Rsb))] * "|]")
  } else {
    warning("Unrecognized column name for p-values: plot.pval has been turned off")
    plot.pval = FALSE
  }
  
  lst_chrm <- unique(rsb.data$CHR)
  nbr_chrm <- length(lst_chrm)
  cum <- rep(0, (nbr_chrm + 1))
  
  if (nbr_chrm > 1) {
    for (i in 1:nbr_chrm) {
      cum[i + 1] <- tail(rsb.data$POSITION[rsb.data$CHR == lst_chrm[i]], n = 1)
    }
  }
  
  cum <- cumsum(cum)
  pos <- rep(0, length(rsb.data$CHR))
  pos_labels <- rep(0, nbr_chrm)
  
  for (i in 1:nbr_chrm) {
    pos[rsb.data$CHR == lst_chrm[i]] <- rsb.data$POSITION[rsb.data$CHR == lst_chrm[i]] + cum[which(lst_chrm == lst_chrm[i])]
    pos_labels[i] <- (cum[i] + cum[i + 1])/2
  }
  
  col_chr <- rsb.data$CHR
  if (nbr_chrm > 1) {
    
    # Plot
    filename <- ifelse(test = substrRight(file.name, 4L) == paste0(".", file.type), 
                       yes = paste0(file.name, ".rsbplot"), no = paste0(file.name, paste0(".score.", file.type)))
    
    Cairo::Cairo(file = filename, type = file.type, width = 1200, height = 600, units = "px", pointsize = 14, dpi = 72)
    
    par(mar = c(5, 5, 4, 2) + 0.1)
    plot(pos, rsb.data[, 3], pch = pch, cex = cex, las = 1, col = col_chr, xaxt = "n", xlab = "Chromosome", 
         ylab = expression(italic(Rsb)), cex.lab = cex.lab, main = main, cex.main = cex.main, cex.axis = cex.axis)
    axis(1, at = pos_labels, labels = lst_chrm, las = 1, cex.lab = cex.lab, cex.axis = 0.5)
    abline(h = c(-ylim.scan, ylim.scan), lty = 2)
    dev.off()
    
    if (plot.pval) {
      
      # Plot
      filename <- ifelse(test = substrRight(file.name, 4L) == paste0(".", file.type), 
                         yes = paste0(file.name, ".rsbplot"), no = paste0(file.name, paste0(".log.", file.type)))
      
      Cairo::Cairo(file = filename, type = file.type, width = 1200, height = 600, units = "px", pointsize = 14, dpi = 72)
      
      par(mar = c(5, 5, 4, 2) + 0.1)
      plot(pos, rsb.data[, 4], pch = pch, cex = cex, las = 1, col = col_chr, xaxt = "n", xlab = "Chromosome", 
           ylab = expression(-log[10] ~ "(" * italic(p) * "-value)"), cex.lab = cex.lab, main = main, 
           cex.main = cex.main, cex.axis = cex.axis)
      axis(1, at = pos_labels, labels = lst_chrm, las = 1, cex.lab = cex.lab, cex.axis = 0.5)
      abline(h = c(-ylim.scan, ylim.scan), lty = 2)
      dev.off()
      
    }
    
  } else {
    
    if (max(pos) < 1000) {
      scale <- 1
      unit = "(bp)"
    }
    else if (max(pos) < 1e+06) {
      scale <- 1000
      unit = "(kb)"
    }
    else if (max(pos) < 1e+09) {
      scale <- 1e+06
      unit = "(Mb)"
    }
    else {
      scale <- 1e+09
      unit = "(Gb)"
    }
    
    # Plot
    filename <- ifelse(test = substrRight(file.name, 4L) == paste0(".", file.type), 
                       yes = paste0(file.name, ".rsbplot"), no = paste0(file.name, paste0(".score.", file.type)))
    
    Cairo::Cairo(file = filename, type = file.type, width = 1200, height = 600, units = "px", pointsize = 14, dpi = 72)
    
    par(mar = c(5, 5, 4, 2) + 0.1)
    plot(pos/scale, rsb.data[, 3], pch = pch, col = col_chr, 
         xlab = paste("Position", unit), ylab = expression(italic(Rsb)), 
         cex.lab = cex.lab, main = main, cex.main = cex.main, cex.axis = cex.axis)
    axis(1, at = pos_labels, labels = lst_chrm, las = 1, cex.lab = cex.lab, cex.axis = 0.5)
    abline(h = c(-ylim.scan, ylim.scan), lty = 2)
    dev.off()
    
    if (plot.pval) {
      
      # Plot
      filename <- ifelse(test = substrRight(file.name, 4L) == paste0(".", file.type), 
                         yes = paste0(file.name, ".rsbplot"), no = paste0(file.name, paste0(".log.", file.type)))
      
      Cairo::Cairo(file = filename, type = file.type, width = 1200, height = 600, units = "px", pointsize = 14, dpi = 72)
      
      par(mar = c(5, 5, 4, 2) + 0.1)
      plot(pos/scale, rsb.data[, 4], pch = pch, col = col_chr, 
           xlab = paste("Position", unit), ylab = expression(-log[10] ~ "(" * italic(p) * "-value)"), 
           cex.lab = cex.lab, main = main, cex.main = cex.main, cex.axis = cex.axis)
      axis(1, at = pos_labels, labels = lst_chrm, las = 1, cex.lab = cex.lab, cex.axis = 0.5)
      abline(h = c(-ylim.scan, ylim.scan), lty = 2)
      dev.off()
      
    }
  }
  
}