R/DE.plot.R

In SoniaTC/NOISeq: Exploratory analysis and differential expression for RNA-seq data

DE.plot <- function (output, q = NULL, graphic = c("MD","expr","chrom","distr"),
                     pch = 20, cex = 0.5, col = 1, pch.sel = 1, cex.sel = 0.6, col.sel = 2,
                     log.scale = TRUE, chromosomes = NULL, join = FALSE,...) {
  
  mypar <- par(no.readonly = TRUE)
  
  if (class(output) != "Output") 
    stop("Error. Output argument must contain an object generated by noiseq or noiseqbio functions.\n")
  
  graphic <- match.arg(graphic)
  
  noiseqbio = "theta" %in% colnames(output@results[[1]])[1:4]
  
  
  ## MD plot
  if (graphic == "MD") {
    
    if (noiseqbio) {      
      
      M = output@results[[1]][,"log2FC"]
      D = abs(output@results[[1]][,1] - output@results[[1]][,2])+1
      names(M) = names(D) = rownames(output@results[[1]])
      
      plot(M, D, pch = pch, xlab = "M", ylab = "D", cex = cex, col = col, log = "y",...)
      
      if(!is.null(q)) {
        mySelection = rownames(degenes(output,q))
        points(M[mySelection], D[mySelection], col = col.sel, pch = pch.sel, cex = cex.sel)        
      }     
      
    } else {
      
      plot(output@results[[1]][,"M"], (1+output@results[[1]][,"D"]), pch = pch, 
           xlab = "M", ylab = "D", cex = cex, col = col, log = "y",...)
      
      if(!is.null(q)) {
        mySelection = rownames(degenes(output,q))
        points(output@results[[1]][mySelection,"M"], output@results[[1]][mySelection,"D"]+1, col = col.sel, pch = pch.sel, cex = cex.sel)
      }      
      
    }          
  } 
  
  
  ## Expression plot: Condition1 vs Condition2
  else if (graphic == "expr") { 
    
    data <- cbind(output@results[[1]][1],output@results[[1]][2])
    rownames(data) <- rownames(output@results[[1]])
    colnames(data) <- colnames(output@results[[1]][c(1,2)])
    
    if (log.scale) {
      
      k <- min(data, na.rm=TRUE)
      escala = logscaling(data, base = 2, k = k)
      
      plot(escala$data, pch = pch, cex = cex, col = col, yaxt = "n", xaxt = "n",...)
      axis(side = 1, at = escala$at, labels = escala$labels)
      axis(side = 2, at = escala$at, labels = escala$labels)
      
      if(!is.null(q)) {
        
        mySelection = rownames(degenes(output,q))  
        points(escala$data[mySelection,], col = col.sel, pch = pch.sel, cex = cex.sel)
      }
      
    } else {
      
      plot(data, pch = pch, cex = cex, col = col,...)
      
      if(!is.null(q)) {
        
        mySelection = rownames(degenes(output,q))  
        points(data[mySelection,], col = col.sel, pch = pch.sel, cex = cex.sel)
        
      }
    }
  }
  
  
  ## MANHATTAN PLOT
  else if (graphic == "chrom") {
    
    mydata <- data.frame(as.character(output@results[[1]][,"Chrom"]),
                         output@results[[1]][,c("GeneStart","GeneEnd")],output@results[[1]][,1:2])
    mydata <- na.omit(mydata)
    
    colnames(mydata) <- c("chr", "start", "end", colnames(mydata)[-c(1:3)])
    
    if (is.null(chromosomes)) {  # todos los cromosomas
      chromosomes <- unique(mydata$chr)
      chromosomes = sort(chromosomes)
    }
    
    print("REMEMBER. You are plotting these chromosomes and in this order:")
    print(chromosomes)
    
    # logarithmic scale
    if (log.scale) {
      if(min(mydata[,-c(1:3)], na.rm = TRUE) < 1) {
        kk <- -min(mydata[,-c(1:3)], na.rm = TRUE)+1
      } else { kk <- 0 }
      mydata[,-c(1:3)] <- log2(mydata[,-c(1:3)]+kk)
    }  
    
    # Selecting chromosomes and ordering positions
    ordenat = NULL
    for (cromo in chromosomes) {
      myselec = mydata[mydata[,"chr"] == cromo,]
      myselec = myselec[order(myselec[,"start"]),]
      ordenat = rbind(ordenat, myselec)
    }  
    
    sel.ord <- NULL
    
    if (!is.null(q)) {
      # up-regulated in first condition
      cond1 <- rownames(degenes(output,q,M="up"))
      # up-regulated in second condition
      cond2 <- rownames(degenes(output,q,M="down"))
      
      cond1 <- (rownames(ordenat) %in% cond1)*1
      cond2 <- (rownames(ordenat) %in% cond2)*1
      
      sel.ord  <- cbind(cond1, cond2)
      rownames(sel.ord) <- rownames(ordenat)
    }
    
    chr.long <- aggregate(ordenat[,"end"],
                          by = list(as.character(ordenat[,"chr"])), max)
    chr.long <- chr.long[match(chromosomes, chr.long[,1]),]
    
    
    if (join) { # si todos los cromosomas van en el mismo plot
      
      total.long <- sum(as.numeric(chr.long$x))
      
      chr.start <- cumsum(c(1,chr.long$x[-length(chr.long$x)]))
      names(chr.start) <- chr.long[,1]
      
      plot(c(1,total.long), c(-max(ordenat[,5]), max(ordenat[,4])),
           type = "n", xlab = "", ylab = "Expression data", xaxt = "n")
      
      axis(side = 1, at = chr.start, labels = chr.long[,1], font = 2)
      
      abline(h = 0, lty = 2, lwd = 0.5)    
      
      for (ch in chromosomes) {
        
        dat.chr <- ordenat[which(ordenat[,"chr"] == ch),]
        dat.chr[,c("start","end")] <- dat.chr[,c("start","end")] +
          chr.start[ch] - 1
        
        rect(xleft = dat.chr[,"start"], ybottom = 0,
             xright = dat.chr[,"end"], ytop = dat.chr[,4],
             col = "grey", border = NA)
        
        rect(xleft = dat.chr[,"start"], ybottom = -dat.chr[,5],
             xright = dat.chr[,"end"], ytop = 0,
             col = "grey", border = NA)
        
        if (!is.null(q)) {
          
          aux <- which(rownames(sel.ord) %in% rownames(dat.chr))          
          
          sel.chr1 <- dat.chr[,4]*sel.ord[aux,1]
          sel.chr2 <- -dat.chr[,5]*sel.ord[aux,2]
          
          rect(xleft = dat.chr[,"start"], ybottom = 0,
               xright = dat.chr[,"end"], ytop = sel.chr1,
               col = 2, border = NA)
          
          rect(xleft = dat.chr[,"start"], ybottom = sel.chr2,
               xright = dat.chr[,"end"], ytop = 0,
               col = 3, border = NA)        
        }
        
        segments(x0 = dat.chr[,"start"], y0 = 0, x1 = dat.chr[,"end"],
                 y1 = 0, col = 4, lwd = 0.5) # annotated genes
        
      }        
      
      text(c(1,1), 0.9*c(max(ordenat[,4]), -max(ordenat[,5])),
           colnames(mydata)[4:5], font = 3, adj = 0, col = 2:3)
      
      
      
      
      
      # a plot for each chromosome
    } else {    
      
      num.chr <- length(chromosomes)
      
      k <- 20
      
      long.prop <- round((chr.long$x / min(chr.long$x))*k, 0)
      
      while (max(long.prop)*num.chr > 500 | max(long.prop) > 50) {
        
        k <- k-1
        
        long.prop <- round((chr.long$x / min(chr.long$x))*k, 0)
        
      }
      
      forlayout <- matrix(0, num.chr, max(long.prop))
      
      #print(dim(forlayout))
      
      for (i in 1:num.chr) {
        
        forlayout[i, 1:long.prop[i]] <- i
        
      }
      
      layout(forlayout)
      
      if (num.chr > 1)	par(mar=c(2,4.1,0.1,0.1))
      
      miylim <- c(-max(ordenat[,5], na.rm=TRUE), max(ordenat[,4],na.rm=TRUE))
      
      for (i in 1:num.chr) {
        
        apintar <- ordenat[which(ordenat[,"chr"] ==  chromosomes[i]), 2:5]
        
        plot(c(1,chr.long[i,2]), miylim, type = "n",
             xlab = "", ylab = chromosomes[i], xaxt = "n",
             font.lab = 2, cex.lab = 1.3)
        
        rect(xleft = apintar[,"start"], ybottom = 0, xright = apintar[,"end"],
             ytop = apintar[,3], col = "grey", border = NA)
        
        rect(xleft = apintar[,"start"], ybottom = -apintar[,4],
             xright = apintar[,"end"], ytop = 0, col = "grey", border = NA)
        
        if (!is.null(q)) {
          
          aux <- which(rownames(sel.ord) %in% rownames(apintar))          
          asel <- sel.ord[aux,]*apintar[,3:4]
          
          rect(xleft = apintar[,"start"], ybottom = 0, xright = apintar[,"end"],
               ytop = asel[,1], col = 2, border = NA)
          
          rect(xleft = apintar[,"start"], ybottom = -asel[,2],
               xright = apintar[,"end"], ytop = 0, col = 3, border = NA)        
          
        }      
        
        abline(h = 0, lty = 2, lwd = 0.5)
        
        segments(x0 = apintar[,"start"], y0 = 0, x1 = apintar[,"end"],
                 y1 = 0, col = 4, lwd = 0.5) # annotated genes
        
        etiq <- mypretty(c(1,chr.long[i,2]), n = 10)
        
        axis(side = 1, at = etiq, labels = etiq)
        
        text(c(1,1), 0.9*miylim, colnames(mydata)[5:4],
             font = 3, adj = 0, col = 3:2)
      }
    }
  }
  
  
  
  ## DEG distribution across biotypes/chromosomes
  else if (graphic == "distr") {    
    
    if(!is.null(q)) {     # Computing DEG
      mySelection = rownames(degenes(output,q))
      detect = rownames(output@results[[1]]) %in% mySelection
    } else {
      stop("You must specify a valid value for q\n")
    }
    
    
    if ("Chrom" %in% colnames(output@results[[1]])) {
      
      numplot = 1
      
      infobio = output@results[[1]][,"Chrom"]
      
      genome <- 100*table(infobio)/sum(table(infobio))
      
      ordre <- order(genome, decreasing = TRUE)      
      
      perdet1 <- genome*table(infobio, detect)[names(genome),2] / table(infobio)[names(genome)]
      
      perdet2 <- 100*table(infobio, detect)[names(genome),2] / sum(table(infobio, detect)[,2])
      
      ceros <- rep(0, length(genome))
      
      biotable1 <- as.matrix(rbind(genome[ordre], perdet1[ordre], perdet2[ordre], ceros))
      rownames(biotable1) <- c("genome", "degVSgenome", "deg", "ceros")
      
      if (!is.null(chromosomes)) biotable1 = biotable1[,chromosomes]
      
      ymaxL1 <- ceiling(max(biotable1, na.rm = TRUE))      
      
    } else { numplot = 0 }
    
    
    
    if ("Biotype" %in% colnames(output@results[[1]])) {
      
      numplot = c(numplot, 1)
      
      infobio = output@results[[1]][,"Biotype"]
      
      genome <- 100*table(infobio)/sum(table(infobio))
      
      ordre <- order(genome, decreasing = TRUE)      
      
      perdet1 <- genome*table(infobio, detect)[names(genome),2] / table(infobio)[names(genome)]
      
      perdet2 <- 100*table(infobio, detect)[names(genome),2] / sum(table(infobio, detect)[,2])
      
      ceros <- rep(0, length(genome))
      
      biotable2 <- as.matrix(rbind(genome[ordre], perdet1[ordre], perdet2[ordre], ceros))
      rownames(biotable2) <- c("genome", "degVSgenome", "deg", "ceros")
      
      higher2 = which(biotable2[1,] > 2)
      lower2 = which(biotable2[1,] <= 2)
      
      if (length(higher2) > 0) {        
        ymaxL2 <- ceiling(max(biotable2[,higher2], na.rm = TRUE)) 
        if (length(lower2) > 0) {
          ymaxR2 <- ceiling(max(biotable2[,lower2], na.rm = TRUE))
          biotable2[,lower2] <- biotable2[,lower2]*ymaxL2/ymaxR2
        } else { ymaxR2 = ymaxL2 }
        
      } else { 
        ymaxR2 <- ceiling(max(biotable2[,lower2], na.rm = TRUE))
        ymaxL2 = ymaxR2               
      }     
      
      
      
    } else { numplot = c(numplot, 0) }
    
    
    
    
    # Plot
    if (sum(numplot) == 0) stop("Biotype or chromosome information is needed for this plot\n")
    
    if (sum(numplot) == 1) {  # 1 Plot 
      
      par(mar = c(10, 4, 2, 2))
      
      if (numplot[1] == 1) {
        
        barplot(biotable1[c(1,3),], main = "DEG distribution across chromosomes",
                xlab = NULL, ylab = "%features", axis.lty = 1, legend = FALSE,
                beside = TRUE, col = c("grey", 2), las = 2,
                ylim = c(0, ymaxL1), border = c("grey", 2))
        
        barplot(biotable1[c(2,4),], main = "DEG distribution across chromosomes",
                xlab = NULL, ylab = "%features", axis.lty = 1, legend = FALSE,
                beside = TRUE, col = c(2, 1), las = 2, density = 30,
                ylim = c(0, ymaxL1), border = 2, add = TRUE)
        
        legend(x = "topright", bty = "n", horiz = FALSE,
               fill = c("grey", 2, 2), density = c(NA,30,NA),
               border = c("grey", 2, 2),
               legend = c("%Chrom in genome", "%DEG in Chrom", "%Chrom in DEG"))
        
      }
      
      
      if (numplot[2] == 1) {
        
        barplot(biotable2[c(1,3),], main = "DEG distribution across biotypes",
                xlab = NULL, ylab = "%features", axis.lty = 1, legend = FALSE,
                beside = TRUE, col = c("grey", 4), las = 2,
                ylim = c(0, ymaxL2), border = c("grey", 4))
        
        barplot(biotable2[c(2,4),], main = "DEG distribution across biotypes",
                xlab = NULL, ylab = "%features", axis.lty = 1, legend = FALSE,
                beside = TRUE, col = c(4, 1), las = 2, density = 30,
                ylim = c(0, ymaxL2), border = 4, add = TRUE)
        
        axis(side=4, at = pretty(c(0,ymaxL2), n = 5), 
             labels = round(pretty(c(0,ymaxL2), n = 5)*ymaxR2/ymaxL2, 1))
        
        if (ymaxR2 != ymaxL2) {
          abline(v = 3*length(higher2) + 0.5, col = 3, lwd = 2, lty = 2)
        }    
        
        legend(x = "topright", bty = "n", horiz = FALSE, 
               fill = c("grey", 4, 4), density = c(NA,30,NA),
               border = c("grey", 4, 4),
               legend = c("%Biotype in genome", "%DEG in Biotype", "%Biotype in DEG"))        
      } 
    }
    
    
    if (sum(numplot) == 2) {  # 2 Plots 
      
      par(mar = c(10, 4, 2, 2), mfrow = c(1,2))
      
      # Chromosomes
      barplot(biotable1[c(1,3),], main = "DEG distribution across chromosomes",
              xlab = NULL, ylab = "%features", axis.lty = 1, legend = FALSE,
              beside = TRUE, col = c("grey", 2), las = 2,
              ylim = c(0, ymaxL1), border = c("grey", 2))
      
      barplot(biotable1[c(2,4),], main = "DEG distribution across chromosomes",
              xlab = NULL, ylab = "%features", axis.lty = 1, legend = FALSE,
              beside = TRUE, col = c(2, 1), las = 2, density = 30,
              ylim = c(0, ymaxL1), border = 2, add = TRUE)
      
      legend(x = "topright", bty = "n", horiz = FALSE,
             fill = c("grey", 2, 2), density = c(NA,30,NA), border = c("grey", 2, 2),
             legend = c("%Chrom in genome", "%DEG in Chrom", "%Chrom in DEG"))
      
      
      # Biotypes
      barplot(biotable2[c(1,3),], main = "DEG distribution across biotypes",
              xlab = NULL, ylab = "%features", axis.lty = 1, legend = FALSE,
              beside = TRUE, col = c("grey", 4), las = 2,
              ylim = c(0, ymaxL2), border = c("grey", 4))
      
      barplot(biotable2[c(2,4),], main = "DEG distribution across biotypes",
              xlab = NULL, ylab = "%features", axis.lty = 1, legend = FALSE,
              beside = TRUE, col = c(4, 1), las = 2, density = 30,
              ylim = c(0, ymaxL2), border = 4, add = TRUE)
      
      axis(side=4, at = pretty(c(0,ymaxL2), n = 5), 
           labels = round(pretty(c(0,ymaxL2), n = 5)*ymaxR2/ymaxL2, 1))
      
      if (ymaxR2 != ymaxL2) {
        abline(v = 3*length(higher2) + 0.5, col = 3, lwd = 2, lty = 2)
      }    
      
      legend(x = "topright", bty = "n", horiz = FALSE, 
             fill = c("grey", 4, 4), density = c(NA,30,NA),
             border = c("grey", 4, 4),
             legend = c("%Biotype in genome", "%DEG in Biotype", "%Biotype in DEG"))        
    }
  }

  par(mypar)

}




## Auxiliar function
mypretty <- function(x, n = 5) {
  
  mywidth <- diff(x)
  
  mybin <- ceiling(mywidth/(n-1))
  
  mylabels0 <- x[1] + mybin*(0:(n-1))
  
  ndig <- nchar(as.character(mylabels0))
  #print(ndig)
  
  mylabels <- mylabels0
  
  k <- 0
  
  for (i in 2:(n-1)) {
    
    mylabels[i] <- (round(mylabels0[i]/10^(ndig[i]-1),k))*10^(ndig[i]-1)
    
    while (mylabels[i] == mylabels[i-1]) {
      
      k <- k+1
      
      mylabels[i] <- (round(mylabels0[i]/10^(ndig[i]-1),k))*10^(ndig[i]-1)
      
    }
    
  }
  
  mylabels
}