R/moduleCurve.R

In menggf/DEComplexDisease: A tool for differential expression analysis and DEGs based investigation to complex diseases by bi-clustering analysis

#' Display the patient-gene numbers during bi-clustering analysis
#'
#' Plot the curve of patient-gene number during bi-clustering analysis
#'
#' @docType methods
#'
#' @name module.curve
#' @param res.module a 'seed.module' or 'cluster.module' object returned by \code{\link{seed.module}} or \code{\link{cluster.module}}
#' @param mod the module to plot
#'
#' @author Guofeng Meng
#'
#' @references
#'
#' @details This function is used to display the patient and gene number during the bi-clustering analysis. It can be used for users to select the better gene or patient number for breakpoints.
#'
#' @return The plot for gene and patient number.
#'
#' @examples
#' module.curve(cluster.mod, 'M1')
#' @export

module.curve <- function(res.module, mod = names(res.module)[1]) {
    if (length(mod) > 1) {
        print("Warning: only the first elements of `mod` is used!")
        mod = mod[1]
    }
    if (!is(res.module, "seed.module") & !is(res.module, "cluster.module")) 
        stop("Error: reg.module: must be the output of 'seed.module' or 'cluster.module'!")
    mods = names(res.module)
    mods = mods[mods != "decd.specific" & mods != "decd.input" & mods != "decd.clustering"]
    if (!mod %in% mods) {
        stop("Error: mod: is not recognized")
    }
    
    pas = res.module[[mod]][["curve"]][["no.patient"]]
    ges = res.module[[mod]][["curve"]][["no.gene"]]
    sim = res.module[[mod]][["curve"]][["score"]]
    
    tag = sim != -1
    par(mar = c(4, 4, 2, 4) + 0.1)
    plot(pas, ges, pch = 16, xlab = "patients", ylab = "genes", col = "green")
    lines(pas, ges, col = "green", lwd = 2)
    points(pas[1], ges[1], col = "red", pch = "O")
    points(pas[length(pas)], ges[length(ges)], col = "brown", pch = "O")
    text(pas[length(pas)], ges[length(ges)], "max.patients", pos = 3, adj = 1)
    text(pas[1], ges[1], "max.genes", pos = 4)
    n = length(ges)
    wh1 = wh2 = wh3 = wh4 = -1
    if (n < 30) {
        wh1 = wh2 = wh3 = wh4 = round((1 + n)/2)
    } else {
        fit = lm(ges ~ poly(pas, 10, raw = TRUE))
        cf = fit$coefficients
        cf[is.na(cf)] = 0
        cf1 = .cof(cf)
        dt = .fv(pas, cf1)
        z = vapply(seq_len(n), function(m) sd(dt[seq_len(m)]), 0.1)
        wh1 = which.max(z)
        wh2 = which.max(ges * pas)
        wh3 = which.max(dt[seq_len(length(dt) - 20)])
        wh4 = which(sim == sim[tag][which.min(sim[tag])])
    }
    
    if (wh1 == wh2) {
        points(pas[wh1], ges[wh1], col = "red", pch = "O")
        text(pas[wh1], ges[wh1], "model", pos = 4)
    } else {
        points(pas[wh1], ges[wh1], col = "red", pch = "O")
        text(pas[wh1], ges[wh1], "slope.clustering", pos = 4)
        points(pas[wh2], ges[wh2], col = "red", pch = "O")
        text(pas[wh2], ges[wh2], "max.square", pos = 4)
        points(pas[wh3], ges[wh3], col = "red", pch = "O")
        text(pas[wh3], ges[wh3], "min.slope", pos = 2)
        points(pas[wh4], ges[wh4], col = "red", pch = "O")
        text(pas[wh4], ges[wh4], "min.similarity", pos = 2)
    }
    if (length(tag[tag]) != 0) {
        frac = (max(ges) - min(ges))/3
        frac2 = (max(sim[tag]) - min(sim[tag]))/3
        axis(side = 4, c(min(ges), min(ges) + frac, min(ges) + 2 * frac, max(ges)), 
            round(c(min(sim[tag]), min(sim[tag]) + frac2, min(sim[tag]) + 2 * frac2, 
                max(sim[tag])), digits = 2))
        par(new = TRUE)
        plot(pas[tag], sim[tag] * (max(ges) - min(ges))/(max(sim) - min(sim[tag])), 
            xlab = "", ylab = "", col = "brown", xaxt = "n", yaxt = "n")
        
        lines(pas[tag], sim[tag] * (max(ges) - min(ges))/(max(sim) - min(sim[tag])), 
            col = "brown", lty = 2)
        mtext("Similarity", side = 4, line = 3)
        legend("top", c("genes", "Similarity"), col = c("green", "brown"), lty = c(1, 
            2), lwd = 3)
    }
}