muma (version 1.4)

volcano: Volcano plots

Description

This function is implemented in the unique function for univariate statistical analysis 'univariate'. This function generates volcano plots of each variable in the data set according to the class definition provided in the second column of the file.

Usage

volcano(file, plot.vol)

Arguments

file

a connection or a character string giving the name of the file containing the variables (matrix columns) to plot.

plot.vol

logical, default is 'TRUE'. Whether to visualize or not volcano plots generated.

Details

For details see ?univariate.

Examples

## The function is currently defined as
function (file, plot.vol) 
{
    pwdfile = paste(getwd(), "/Univariate/DataTable.csv", sep = "")
    file = pwdfile
    x <- read.csv(file, sep = ",", header = TRUE)
    x.x = x[, 3:ncol(x)]
    rownames(x.x) = x[, 2]
    k = matrix(x[, 1], ncol = 1)
    x.n = cbind(k, x.x)
    sorted = x.n[order(x.n[, 1]), ]
    sorted.x = as.matrix(sorted[, -1], ncol = ncol(sorted) - 
        1)
    g = c()
    for (i in 1:nrow(sorted)) {
        if (any(g == sorted[i, 1])) {
            g = g
        }
        else {
            g = matrix(c(g, sorted[i, 1]), ncol = 1)
        }
    }
    NoF = nrow(g)
    dirout.fc = paste(getwd(), "/Univariate/Fold_Changes/", sep = "")
    dir.create(dirout.fc)
    dirout.vol = paste(getwd(), "/Univariate/Volcano_Plots/", 
        sep = "")
    dir.create(dirout.vol)
    for (i in 1:NoF) {
        for (j in 1:NoF) {
            if (i < j) {
                ni = paste("r.", i, ".csv", sep = "")
                nj = paste("r.", j, ".csv", sep = "")
                pwdi = paste(getwd(), "/Univariate/Groups/", 
                  ni, sep = "")
                pwdj = paste(getwd(), "/Univariate/Groups/", 
                  nj, sep = "")
                pv = paste(getwd(), "/Univariate/Pvalues/Pvalues_", 
                  i, "vs", j, ".csv", sep = "")
                I = read.csv(pwdi, header = TRUE)
                J = read.csv(pwdj, header = TRUE)
                I = I[, -1]
                J = J[, -1]
                meanI = matrix(colMeans(I), ncol = ncol(I))
                meanJ = matrix(colMeans(J), ncol = ncol(J))
                MeanI = matrix(rep(NA, ncol(I)), nrow = 1)
                MeanJ = matrix(rep(NA, ncol(I)), nrow = 1)
                for (m in 1:ncol(I)) {
                  if (meanI[, m] < 0 | meanJ[, m] < 0) {
                    MeanI[, m] = 1
                    MeanJ[, m] = 1
                  }
                  else {
                    MeanI[, m] = meanI[, m]
                    MeanJ[, m] = meanJ[, m]
                  }
                }
                FC = matrix(MeanI/MeanJ, nrow = ncol(I))
                rownames(FC) = colnames(I)
                fc.csvfile = paste("Fold_Change_", i, "vs", j, 
                  ".csv", sep = "")
                write.csv(FC, paste(dirout.fc, fc.csvfile, sep = ""))
                PV = read.csv(pv, header = TRUE)
                PV = matrix(PV[, -1], ncol = 1)
                logfc = log2(FC)
                logpv = -log10(PV)
                colpv = matrix(rep(NA, nrow(PV)), ncol = 1)
                for (p in 1:nrow(PV)) {
                  if (logfc[p, ] < -0.3219281 | logfc[p, ] > 
                    0.2630344) {
                    if (logpv[p, ] > 1.30103) {
                      colpv[p, ] = "navy"
                    }
                    else {
                      colpv[p, ] = "dark grey"
                    }
                  }
                  else {
                    colpv[p, ] = "dark grey"
                  }
                }
                max.fc = 1.3 * (max(abs(logfc)))
                V = paste(dirout.vol, "VolcanoPlot_", i, "vs", 
                  j, ".pdf", sep = "")
                pospv = matrix(rep(NA, nrow(PV)), ncol = 1)
                for (p in 1:nrow(PV)) {
                  if (logfc[p, ] < 0) {
                    pospv[p, ] = 2
                  }
                  else {
                    pospv[p, ] = 4
                  }
                }
                pdf(V)
                plot(logfc, logpv, col = colpv, pch = 19, xlim = c(-max.fc, 
                  max.fc), xlab = "Log2 (Fold Change)", ylab = "Log10 (Pvalue)", 
                  main = paste("Volcano Plot ", i, " vs ", j, 
                    sep = ""), sub = "(Variables in Blue are significant (Pvalue<0.05) and showed Fold Changes >1.2 or <0.8)")
                text(logfc, logpv, labels = colnames(sorted.x), 
                  cex = 0.8, pos = pospv, col = colpv)
                axis(2, at = c(-1, 150), pos = c(-0.3219281, 
                  0), col = "blue", lwd = 0.3)
                axis(2, at = c(-1, 150), pos = c(0.2630344, 0), 
                  col = "blue", lwd = 0.3)
                axis(1, at = c(-150, 150), pos = c(1.30103, 0), 
                  col = "blue", lwd = 0.3)
                dev.off()
                if (plot.vol) {
                  dev.new()
                  plot(logfc, logpv, col = colpv, pch = 19, xlim = c(-max.fc, 
                    max.fc), xlab = "Log2 (Fold Change)", ylab = "Log10 (Pvalue)", 
                    main = paste("Volcano Plot ", i, " vs ", 
                      j, sep = ""), sub = "(Variables in Blue are significant (Pvalue<0.05) and showed Fold Changes >1.2 or <0.8)")
                  text(logfc, logpv, labels = colnames(sorted.x), 
                    cex = 0.8, pos = pospv, col = colpv)
                  axis(2, at = c(-1, 150), pos = c(-0.3219281, 
                    0), col = "blue", lwd = 0.3)
                  axis(2, at = c(-1, 150), pos = c(0.2630344, 
                    0), col = "blue", lwd = 0.3)
                  axis(1, at = c(-150, 150), pos = c(1.30103, 
                    0), col = "blue", lwd = 0.3)
                }
            }
        }
    }
  }