R/negative_cor.R

In AllenTiTaiWang/anamiR: An integrated analysis package of miRNA and mRNA expression data

#' Find the correlation coefficient between each gene and miRNA.
#'
#' This function will calculate the correlation coefficient between each gene
#' and miRNA from differential expressed data, which are produced by
#' \link{differExp_discrete} or \link{differExp_continuous}. After filtering
#' the positive and higher than cutoff value of correlation, this function
#' would return a matrix with seven columns, including miRNA, gene, correlation
#' coefficients and Fold change, P-adjust value for miRNA and gene. Each
#' row represents one potential miRNA-target gene interaction.
#'
#' @return matrix format with each row indicating one potential miRNA-target
#'    gene interaction and seven columns are miRNA, gene, correlation
#'    coefficient and Fold change, P-adjust value for miRNA and gene.
#'
#' @seealso \code{\link[stats]{cor}} for calculation of correlation.
#'
#' @param mrna_data differential expressed data in matrix format, with
#'    sample name in columns and gene symbol in rows, which is generated by
#'    \link{differExp_discrete} or \link{differExp_continuous}.
#' @param mirna_data differential expressed data in matrix format, with
#'    sample name in columns and miRNA in rows, which is generated by
#'    \link{differExp_discrete} or \link{differExp_continuous}, miRNA should
#'    be miRBase 21 version now.
#' @param method methods for calculating correlation coefficient, including
#'    "pearson", "spearman", "kendall". Default is "pearson". From function
#'    \code{\link[stats]{cor}}
#' @param cut.off an numeric value indicating a threshold of correlation
#'    coefficient for every potential miRNA-genes interactions. Default is -0.5,
#'    however, if no interaction pass the threshold, this function would add
#'    0.2 value in threshold until at least one interaction passed the threshold.
#'
#'
#' @examples
#' ## Use the internal dataset
#' data("mirna", package = "anamiR", envir = environment())
#' data("pheno.mirna", package = "anamiR", envir = environment())
#' data("mrna", package = "anamiR", envir = environment())
#' data("pheno.mrna", package = "anamiR", envir = environment())
#'
#' ## SummarizedExperiment class
#' require(SummarizedExperiment)
#' mirna_se <- SummarizedExperiment(
#'  assays = SimpleList(counts=mirna),
#'  colData = pheno.mirna)
#'
#' ## SummarizedExperiment class
#' require(SummarizedExperiment)
#' mrna_se <- SummarizedExperiment(
#'  assays = SimpleList(counts=mrna),
#'  colData = pheno.mrna)
#'
#' ## Finding differential miRNA from miRNA expression data with t.test
#' mirna_d <- differExp_discrete(
#'    se = mirna_se,
#'    class = "ER",
#'    method = "t.test"
#' )
#'
#' ## Finding differential mRNA from mRNA expression data with t.test
#' mrna_d <- differExp_discrete(
#'    se = mrna_se,
#'    class = "ER",
#'    method = "t.test"
#' )
#'
#' ## Correlation
#' cor <- negative_cor(mrna_data = mrna_d, mirna_data = mirna_d,
#'      method = "pearson"
#' )
#'
#' @import stats
#' @export
negative_cor <- function(
  mrna_data,
  mirna_data,
  method = c("pearson",
             "kendall",
             "spearman"),
  cut.off = -0.5
) {

  method <- match.arg(method)

  common_column <- intersect(colnames(mrna_data),
                             colnames(mirna_data))

  mrna_data <- as.data.frame(mrna_data[, common_column])
  mirna_data <- as.data.frame(mirna_data[, common_column])

  cal_cor <- function(data_1, data_2, cor_cut) {
    n <- 1
    corr <- list()
    for (i in seq_len(nrow(data_1))) {
      for (j in seq_len(nrow(data_2))) {
        mrna <- as.numeric(data_1[i, 1:(ncol(mrna_data) - 5)])
        mirna <- as.numeric(data_2[j, 1:(ncol(mirna_data) - 5)])
        tmp <- stats::cor(mrna, mirna, method = method)
        #print(tmp)
        if (tmp < cor_cut) {
          corr[[n]] <- row.names(data_2)[j]
          corr[[n]][2] <- row.names(data_1)[i]
          corr[[n]][3] <- tmp
          corr[[n]][4] <- data_2[["log-ratio"]][j]
          corr[[n]][5] <- data_2[["P-adjust"]][j]
          corr[[n]][6] <- data_2[["mean_case"]][j]
          corr[[n]][7] <- data_2[["mean_control"]][j]
          corr[[n]][8] <- data_1[["log-ratio"]][i]
          corr[[n]][9] <- data_1[["P-adjust"]][i]
          corr[[n]][10] <- data_1[["mean_case"]][i]
          corr[[n]][11] <- data_1[["mean_control"]][i]
          n <- n + 1
        }
      }
    }
    return(corr)
  }

  corr <- cal_cor(mrna_data, mirna_data, cut.off)
  corr <- do.call(rbind, corr)

  if (is.null(corr)) {
    cut.off <- cut.off + 0.2
    corr <- cal_cor(mrna_data, mirna_data, cut.off)
    corr <- do.call(rbind, corr)
  }
  last_column <- paste("Correlation(", cut.off, ")")
  colnames(corr) <- c("miRNA", "Gene", last_column,
                      "logratio(miRNA)", "P-adjust(miRNA)",
                      "mean_case(miRNA)", "mean_control(miRNA)",
                      "logratio(gene)", "P-adjust(gene)",
                      "mean_case(gene)", "mean_control(gene)")
  return(corr)
}