R/plot_ExpSd_Sig.R

In YC3/mirNet: microRNA netowrk analyses

#' @title Plot the stadard deviation against average experssion of genes/microRNAs in the context of differential expression tests.
#'
#' @description \code{plot_ExpSd_Sig} plot the standard deviation of genes/microRNAs against their average expression level. With significantly up/down regulated genes/microRNAs highlighted.
#'
#' @param data A matrix, the normalized gene/microRNA expression dataset, should be a numeric matrix, with rows referring to genes/microRNAs and columns to samples.
#' @param stat The statistic table output from function \code{\link[mirNet]{DoLmFit}}, \code{\link[mirNet]{pairComp}} or \code{\link[mirNet]{DoWilcoxonTest}}.
#' @param pheno A vector of sample phenotypes. Sample phenotype in a scientific research could be treatment/control, normal/cancer or smoker/non-smoker. Different phenotypes should each be encoded as 0/1 when inputting to \code{DoTest}, for example, Normal-0; Cancer-1. Currently, only this function accepts two-class pheno data.
#' @param c_fdr The name of the column stores adjusted p-value.
#' @param fdr_cutoff The cutoff of adjusted p-value.
#' @param c_dir The name of the column stores a statistic that is indicative of the direction of alteration.
#' @param c_logFC The name of the column stores log2 fold change.
#' @param logFC_cutoff The cutoff of log2 fold change.
#' @param filename The name of the output figure.
#'
#' @return A scatter plot with significant genes/microRNAs highlighted.
#'
#' @seealso \code{\link[mirNet]{plot_ExpSd}}
#'
#' @export plot_ExpSd_Sig
#'
#' @examples
#' plot_ExpSd_Sig(data.m, stat, pheno, c_fdr = 'fdr', fdr_cutoff = 0.05, c_dir = 'logFC', c_logFC = 'logFC', logFC_cutoff = 1, filename = 'examples')




plot_ExpSd_Sig <- function(data, stat, pheno, c_fdr = 'fdr', fdr_cutoff = 0.05, c_dir = 'logFC', c_logFC = 'logFC', logFC_cutoff = 1, filename){
    # mRNA differential expression results
    up <- rownames(stat)[intersect(intersect(which(stat[, c_fdr] <= fdr_cutoff), which(stat[, c_dir] > 0)), which(abs(stat[, c_logFC]) > logFC_cutoff))]
    dn <- rownames(stat)[intersect(intersect(which(stat[, c_fdr] <= fdr_cutoff), which(stat[, c_dir] < 0)), which(abs(stat[, c_logFC]) > logFC_cutoff))]
    nonsig <- setdiff(rownames(data), c(up, dn))

    data.c <- data[, pheno %in% c('01', 1)]
    up.mean.c <- rowMeans(data.c)[up]
    dn.mean.c <- rowMeans(data.c)[dn]
    no.mean.c <- rowMeans(data.c)[nonsig]
    up.sd.c <- apply(data.c, 1, sd)[up]
    dn.sd.c <- apply(data.c, 1, sd)[dn]
    no.sd.c <- apply(data.c, 1, sd)[nonsig]
    df.c <- data.frame(AveExp = c(no.mean.c, up.mean.c, dn.mean.c), SD = c(no.sd.c, up.sd.c, dn.sd.c), change = rep(c('NONSIG', 'UP', 'DN'), times = c(length(no.mean.c), length(up.mean.c), length(dn.mean.c))))

    data.n <- data[, pheno %in% c('11', 0)]
    up.mean.n <- rowMeans(data.n)[up]
    dn.mean.n <- rowMeans(data.n)[dn]
    no.mean.n <- rowMeans(data.n)[nonsig]
    up.sd.n <- apply(data.n, 1, sd)[up]
    dn.sd.n <- apply(data.n, 1, sd)[dn]
    no.sd.n <- apply(data.n, 1, sd)[nonsig]
    df.n <- data.frame(AveExp = c(no.mean.n, up.mean.n, dn.mean.n), SD = c(no.sd.n, up.sd.n, dn.sd.n), change = rep(c('NONSIG', 'UP', 'DN'), times = c(length(no.mean.n), length(up.mean.n), length(dn.mean.n))))

    df <- rbind(df.c, df.n) %>% mutate(status = rep(c('Caner', 'Normal'), times = c(nrow(df.c), nrow(df.n))))
    df$change <- factor(df$change, levels = c(c('NONSIG', 'UP', 'DN')))


    pdf(paste0(filename, '.pdf'), width = 10, height = 4)
    p <- df %>% ggplot(aes(x = AveExp, y = SD, color = change)) + geom_point(alpha = 0.7) + facet_wrap(.~status, ncol = 2, scale = 'free') + scale_color_manual(values = c('grey', 'red', 'blue')) + theme_bw()
    plot(p)
    dev.off()
}