R/Plot_plsda_stat.R

In 13479776/statTarget: Statistical Analysis of Molecular Profiles

Plot_plsda_stat <- function(pcx, pcy, scaling, Labels) {
    pwd.score = paste(getwd(), "/PLS_DA_", scaling, "/PLSDA_Scores_", scaling, ".csv", sep = "")
    score = read.csv(pwd.score, header = TRUE)
    score.x = score[, -1]
    rownames(score.x) = score[, 1]
    pwd.load = paste(getwd(), "/PLS_DA_", scaling, "/PLSDA_Loadings_", scaling, ".csv", sep = "")
    loading = read.csv(pwd.load, header = TRUE)
    loading.x = loading[, -1]
    rownames(loading.x) = loading[, 1]
    pwd.p = paste(getwd(), "/PLS_DA_", scaling, "/PLSDA_P_", scaling, ".csv", sep = "")
    p = read.csv(pwd.p, header = TRUE)
    p.x = matrix(p[, -1], ncol = 1)
    pwd.ptot = paste(getwd(), "/PLS_DA_", scaling, "/PLSDA_Ptot_", scaling, ".csv", sep = "")
    p = read.csv(pwd.ptot, header = TRUE)
    slink = paste(getwd(), "/DataPretreatment", "/slink.csv", sep = "")
    slink = read.csv(slink, header = TRUE)
    pvar.a = p.x[pcx, ]/p
    pvar.b = p.x[pcy, ]/p
    pvar.ai = round(pvar.a * 100, 1)
    pvar.bi = round(pvar.b * 100, 1)
    cum = pvar.ai + pvar.bi
    xlab = paste("Component", pcx, "(", pvar.ai, "%)", sep = "")
    ylab = paste("Component", pcy, "(", pvar.bi, "%)", sep = "")
    
    max.pc1 = 1.3 * (max(abs(score.x[, pcx])))
    max.pc2 = 1.3 * (max(abs(score.x[, pcy])))
    lim = c()
    if (max.pc1 > max.pc2) {
        lim = c(-max.pc1, max.pc1)
    } else {
        lim = c(-max.pc2, max.pc2)
    }
    
    
    Epaxis <- dataEllipse_sT(score.x[, pcx], score.x[, pcy], levels = c(0.95), add = FALSE, draw = FALSE, 
        col = "black", lwd = 0.4, plot.points = FALSE, center.cex = 0.2)
    
    if (1.1 * max(Epaxis[, "x"]) < max(lim) & 1.1 * max(Epaxis[, "y"]) < max(lim)) {
        lim = lim
    } else {
        lim = c(1.2 * min(as.numeric(Epaxis)), 1.2 * max(as.numeric(Epaxis)))
    }
    
    
    pwdK = paste(getwd(), "/scaleData_", scaling, "/class.csv", sep = "")
    k = read.csv(pwdK, header = TRUE)
    k.s = k[order(k[, 2]), ]
    tutticolors = matrix(c(1, 2, 3, 4, 5, 6, 7, 8, "rosybrown4", "green4", "navy", "purple2", "orange", 
        "pink", "chocolate2", "coral3", "khaki3", "thistle", "turquoise3", "palegreen1", "moccasin", 
        "olivedrab3", "azure4", "gold3", "deeppink"), ncol = 1)
    col = c()
    for (i in 1:nrow(k.s)) {
        col = c(col, tutticolors[k.s[i, 2], ])
    }
    dirout = paste(getwd(), "/PLS_DA_", scaling, "/", sep = "")
    scor = paste(dirout, "ScorePlot_PLS_DA_", scaling, ".pdf", sep = "")
    pdf(scor)
    graphics::plot(score.x[, pcx], score.x[, pcy], col = col, pch = 19, xlab = c(xlab), ylab = c(ylab), 
        xlim = lim, ylim = lim, sub = paste("Cumulative Proportion of Variance Explained = ", cum, 
            "%", sep = ""), main = paste("PLS-DA Score Plot (", scaling, ")", sep = ""))
    axis(1, at = lim * 2, pos = c(0, 0), labels = FALSE, col = "grey", lwd = 0.7)
    axis(2, at = lim * 2, pos = c(0, 0), labels = FALSE, col = "grey", lwd = 0.7)
    
    
    if (Labels) {
        text(score.x[, pcx], score.x[, pcy], col = col, cex = 0.5, labels = rownames(score.x), pos = 1)
    } else {
        legend("topright", legend = levels(factor(slink[, 1])), bty = "n", pch = 19, col = seq_along(levels(factor(slink[, 
            1]))), text.col = seq_along(levels(factor(slink[, 1]))))
    }
    
    # requireNamespace(car) library(car)
    dataEllipse_sT(score.x[, pcx], score.x[, pcy], levels = c(0.95), add = TRUE, col = "black", lwd = 0.4, 
        plot.points = FALSE, center.cex = 0.2)
    dirout = paste(getwd(), "/PLS_DA_", scaling, "/", sep = "")
    # scor = paste(dirout, 'ScorePlot_PLS_DA_', scaling, '.pdf', sep='') dev.copy2pdf(file=scor)
    dev.off()
    Max.pc1 = 1.1 * (max(loading.x[, pcx]))
    Min.pc1 = 1.1 * (min(loading.x[, pcx]))
    Mpc1 = c(Min.pc1, Max.pc1)
    Max.pc2 = 1.1 * (max(loading.x[, pcy]))
    Min.pc2 = 1.1 * (min(loading.x[, pcy]))
    Mpc2 = c(Min.pc2, Max.pc2)
    load = paste(dirout, "W.cPlot_PLS_DA_", scaling, ".pdf", sep = "")
    # par(new = T)
    pdf(load)
    # dev.new()
    graphics::plot(loading.x[, pcx], loading.x[, pcy], xlim = Mpc1, ylim = Mpc2, xlab = paste("w*c values ", 
        pcx, sep = ""), ylab = paste("w*c values ", pcy, sep = ""), main = paste("PLS-DA Loading Plot (", 
        scaling, ")", sep = ""))
    text(loading.x[, pcx], loading.x[, pcy], labels = rownames(loading.x), cex = 0.7, pos = 1)
    axis(1, at = Mpc1 * 2, pos = c(0, 0), labels = FALSE, col = "grey", lwd = 0.7)
    axis(2, at = Mpc2 * 2, pos = c(0, 0), labels = FALSE, col = "grey", lwd = 0.7)
    # load = paste(dirout, 'W.cPlot_PLS_DA_', scaling, '.pdf', sep='') dev.copy2pdf(file=load)
    dev.off()
}