R/comEigenGene.R

In YC3/mirNet: microRNA netowrk analyses

#' @title Computing eigen-genes or eigen-microRNAs from detected modules in co-expression network.
#'
#' @description \code{comEigenGene} computes eigen genes/microRNAs from each module detected in co-expression network, and then computes a P-value from the correlation of eigem gene/microRNA expression and phenotypes.
#'
#' @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 pheno A vector of sample phenotypes. Sample phenotype in a scientific research could be categorical (treatment/control), or continuous (age). If you have multiple phenotypes, use a list with names denoting the corresponding phenotype.
#' @param moduleList A list, entries are vectors of genes/microRNAs in detected modules.
#' @param type A string, denoting the type of summarization of expression of genes/microRNAs, 'egein' for principal component analysis, 'average' for an average of gene/microRNA expression, default to 'eigen'.
#' @param test A string, 'categorical' if sample phenotype if categorical, 'continuous' if sample phenotype is continuous, default to 'categorical'.
#' @param pc Numeric, the number of principal components you want to show in the figure, default to 1.
#'
#' @return A list with P-values as well as detected eigen genes/microRNAs from each module.
#'
#' @importFrom stats kruskal.test
#' @importFrom stats lm
#'
#' @export comEigenGene
#'
#' @examples
#' comEigenGene(data.norm, pheno.v, modules.l, type = 'eigen', test = 'categorical')



comEigenGene <- function(data, pheno, moduleList, type = 'eigen', test = 'categorical', pc = 1){

    L <- length(moduleList)
    eigenGene.l <- vector(length = L, mode = 'list')
    names(eigenGene.l) <- paste0('mod', 1:L)

    p.v <- vector(length = L)
    names(p.v) <- paste0('mod', 1:L)
    for(i in 1:L){

        if(length(moduleList[[i]]) > 1){

            exp <- data[moduleList[[i]], ]
            exp <- exp - apply(exp, 1, mean)

            if(type == 'eigen'){
                # SVD
                svd.o <- svd(exp)
                exp.v <- svd.o$v[, pc]
                eigenGene.l[[i]] <- exp.v
            }
            if(type == 'average'){
                exp.v <- colMeans(exp)
            }

            # association between eigen-miRNA expression and phenotypes
            if(test == 'categorical'){
                p.v[i] <- kruskal.test(exp.v ~ as.factor(pheno))$p.value
            }
            if(test == 'continuous'){
                lm.o <- lm(exp.v ~ pheno)
                p.v[i] <- summary(lm.o)$coeff[2, 4]
            }
        }
        else{
            p.v[i] <- NA
        }

     }
    list(p = p.v, eigen = eigenGene.l)
}