Datasets
Models
Diff of
/R/DIscBIO-generic-DEGanalysis.R
[000000]
..
[28e211]
Switch to side-by-side view
--- a +++ b/R/DIscBIO-generic-DEGanalysis.R @@ -0,0 +1,448 @@ +#' @title Determining differentially expressed genes (DEGs) between all +#' individual clusters. +#' @description This function defines DEGs between all individual clusters +#' generated by either K-means or model based clustering. +#' @param object \code{DISCBIO} class object. +#' @param Clustering Clustering has to be one of the following: +#' ["K-means","MB"]. Default is "K-means" +#' @param K A numeric value of the number of clusters. +#' @param fdr A numeric value of the false discovery rate. Default is 0.05. +#' @param name A string vector showing the name to be used to save the resulted +#' tables. +#' @param export A logical vector that allows writing the final gene list in +#' excel file. Default is TRUE. +#' @param quiet if `TRUE`, suppresses intermediate text output +#' @param plot if `TRUE`, plots are generated +#' @param filename_deg Name of the exported DEG table +#' @param filename_sigdeg Name of the exported sigDEG table +#' @importFrom graphics title +#' @importFrom utils write.csv capture.output +#' @param ... additional parameters to be passed to samr() +#' @return A list containing two tables. +setGeneric( + name = "DEGanalysis", + def = function(object, + K, + Clustering = "K-means", + fdr = 0.05, + name = "Name", + export = FALSE, + quiet = FALSE, + plot = TRUE, + filename_deg = "DEGsTable", + filename_sigdeg = "sigDEG", + ...) { + standardGeneric("DEGanalysis") + } +) + +#' @export +#' @rdname DEGanalysis +setMethod( + f = "DEGanalysis", + signature = "DISCBIO", + definition = function(object, + K, + Clustering = "K-means", + fdr = 0.05, + name = "Name", + export = FALSE, + quiet = FALSE, + plot = TRUE, + filename_deg, + filename_sigdeg, + ...) { + # Validation + if (!(Clustering %in% c("K-means", "MB"))) { + stop("Clustering has to be either K-means or MB") + } + if (Clustering == "K-means") { + Cluster_ID <- object@cpart + if (length(object@cpart) < 1) { + stop("run Clustexp before running DEGanalysisM") + } + } + if (Clustering == "MB") { + Cluster_ID <- object@MBclusters$clusterid + if (length(object@MBclusters$clusterid) < 1) { + stop("run ExprmclustMB before running DEGanalysisM") + } + } + + # Defining initial objects + gene_list <- object@FinalGeneList + gene_names <- rownames(object@expdata) + idx_genes <- is.element(gene_names, gene_list) + gene_names2 <- gene_names[idx_genes] + dataset <- object@expdata[gene_names2, ] + Nam <- colnames(dataset) + num <- 1:K + num1 <- paste("CL", num, sep = "") + + for (n in num) { + Nam <- ifelse((Cluster_ID == n), num1[n], Nam) + } + colnames(dataset) <- Nam + + if (!quiet) { + message("The dataset is ready for differential expression analysis") + } + num1 <- paste("CL", num, sep = "") + clustName <- paste("Cl", num, sep = "") + ClusterLength <- vector() + for (i in num) { + d1 <- clustName[i] + d2 <- dataset[, which(names(dataset) == num1[i])] + ClusterLength[i] <- length(d2) + assign(d1, d2) + } + ccdf <- data.frame(clustName, ClusterLength) + ccdff <- ccdf[order(ClusterLength), ] + clustName <- ccdff[, 1] + if (!quiet) { + print(clustName) + } + first <- vector() + second <- vector() + if (K < 2) { + stop("K has to be at least 2") + } else if (K == 2) { + first <- c(paste0(clustName[1])) + second <- c(paste0(clustName[2])) + } else if (K == 3) { + first <- c( + rep(paste0(clustName[1]), 2), rep(paste0(clustName[2]), 1) + ) + second <- c( + paste0(clustName[2]), + paste0(clustName[3]), + paste0(clustName[3]) + ) + } else if (K == 4) { + first <- c( + rep(paste0(clustName[1]), 3), + rep(paste0(clustName[2]), 1), + rep(paste0(clustName[4]), 1), + rep(paste0(clustName[3]), 1) + ) + second <- c( + paste0(clustName[2]), + paste0(clustName[3]), + paste0(clustName[4]), + paste0(clustName[3]), + paste0(clustName[2]), + paste0(clustName[4]) + ) + } else if (K == 5) { + first <- c( + rep(paste0(clustName[1]), 4), + rep(paste0(clustName[2]), 3), + rep(paste0(clustName[3]), 2), + rep(paste0(clustName[5]), 1) + ) + second <- c( + paste0(clustName[2]), + paste0(clustName[3]), + paste0(clustName[4]), + paste0(clustName[5]), + paste0(clustName[3]), + paste0(clustName[4]), + paste0(clustName[5]), + paste0(clustName[4]), + paste0(clustName[5]), + paste0(clustName[4]) + ) + } else if (K == 6) { + first <- c( + rep(paste0(clustName[1]), 3), + rep(paste0(clustName[5]), 1), + rep(paste0(clustName[1]), 1), + rep(paste0(clustName[2]), 2), + rep(paste0(clustName[5]), 1), + rep(paste0(clustName[2]), 1), + rep(paste0(clustName[3]), 1), + rep(paste0(clustName[5]), 1), + rep(paste0(clustName[3]), 1), + rep(paste0(clustName[5]), 1), + rep(paste0(clustName[4]), 1), + rep(paste0(clustName[5]), 1) + ) + second <- c( + paste0(clustName[2]), + paste0(clustName[3]), + paste0(clustName[4]), + paste0(clustName[1]), + paste0(clustName[6]), + paste0(clustName[3]), + paste0(clustName[4]), + paste0(clustName[2]), + paste0(clustName[6]), + paste0(clustName[4]), + paste0(clustName[3]), + paste0(clustName[6]), + paste0(clustName[4]), + paste0(clustName[6]), + paste0(clustName[6]) + ) + } + o <- 1:K + oo <- o[-length(o)] + com <- sum(oo) + if (!quiet) message("Number of comparisons: ", com * 2, "\n") + comNUM <- paste("comp", seq_len(com), sep = "") + DEGsTable <- data.frame() + DEGsE <- vector() + DEGsS <- vector() + for (i in seq_len(com)) { + FinalDEGsL <- data.frame() + FinalDEGsU <- data.frame() + FDRl <- vector() + FDRu <- vector() + + d1 <- comNUM[i] + d2 <- cbind(get(first[i]), get(second[i])) + assign(d1, d2) + len <- + c(length(get(first[i])[1, ]), length(get(second[i])[1, ])) + y <- c(rep(1:2, len)) + L <- as.matrix(get(comNUM[i])) + gname <- rownames(get(comNUM[i])) + x <- L + data <- list(x = x, y = y, geneid = gname) + if (quiet) { + invisible(capture.output( + samr.obj <- sammy( + data, + resp.type = "Two class unpaired", + assay.type = "seq", + testStatistic = "wilcoxon", + random.seed = 15, + ... + ) + )) + } else { + samr.obj <- sammy( + data, + resp.type = "Two class unpaired", + assay.type = "seq", + testStatistic = "wilcoxon", + random.seed = 15, + ... + ) + } + if (quiet) { + invisible(capture.output( + delta.table <- samr.compute.delta.table(samr.obj) + )) + } else { + delta.table <- samr.compute.delta.table(samr.obj) + } + wm <- which.min(delta.table[, 5]) + if (delta.table[wm, 5] <= fdr) { + w <- which(delta.table[, 5] <= fdr) + delta <- delta.table[w[1], 1] - 0.001 + + if (plot) { + samr.plot(samr.obj, delta) + title(paste0( + "DEGs in the ", second[i], " in ", first[i], " VS ", + second[i] + )) + } + siggenes.table <- samr.compute.siggenes.table( + samr.obj, delta, data, delta.table + ) + FDRl <- + as.numeric(siggenes.table$genes.lo[, 8]) / 100 + FDRu <- + as.numeric(siggenes.table$genes.up[, 8]) / 100 + + siggenes.table$genes.lo[, 8] <- FDRl + siggenes.table$genes.up[, 8] <- FDRu + + DEGsTable[i, 1] <- + paste0(first[i], " VS ", second[i]) + DEGsTable[i, 2] <- second[i] + DEGsTable[i, 3] <- length(FDRu) + DEGsTable[i, 4] <- paste0( + "Up-regulated-", name, second[i], "in", first[i], "VS", + second[i], ".csv" + ) + DEGsTable[i, 5] <- length(FDRl) + DEGsTable[i, 6] <- paste0( + "Low-regulated-", name, second[i], "in", first[i], "VS", + second[i], ".csv" + ) + DEGsTable[i + com, 1] <- + paste0(first[i], " VS ", second[i]) + DEGsTable[i + com, 2] <- first[i] + DEGsTable[i + com, 3] <- length(FDRu) + DEGsTable[i + com, 4] <- paste0( + "Low-regulated-", name, first[i], "in", first[i], "VS", + second[i], ".csv" + ) + DEGsTable[i + com, 5] <- length(FDRl) + DEGsTable[i + com, 6] <- paste0( + "Up-regulated-", name, first[i], "in", first[i], "VS", + second[i], ".csv" + ) + + if (length(FDRl) > 0) { + genes <- siggenes.table$genes.lo[, 3] + if (quiet) { + suppressMessages( + geneList <- + AnnotationDbi::select( + org.Hs.eg.db, + keys = keys(org.Hs.eg.db), + columns = c("SYMBOL", "ENSEMBL") + ) + ) + GL <- c(1, "MTRNR2", "ENSG00000210082") + geneList <- rbind(geneList, GL) + } else { + geneList <- + AnnotationDbi::select( + org.Hs.eg.db, + keys = keys(org.Hs.eg.db), + columns = c("SYMBOL", "ENSEMBL") + ) + GL <- c(1, "MTRNR2", "ENSG00000210082") + geneList <- rbind(geneList, GL) + } + FinalDEGsL <- cbind(genes, siggenes.table$genes.lo) + gene_list <- geneList[, 3] + idx_genes <- is.element(gene_list, genes) + genes2 <- geneList[idx_genes, ] + FinalDEGsL <- merge( + FinalDEGsL, + genes2, + by.x = "genes", + by.y = "ENSEMBL", + all.x = TRUE + ) + FinalDEGsL[, 3] <- FinalDEGsL[, 11] + FinalDEGsL <- FinalDEGsL[, c(-1, -10, -11)] + FinalDEGsL <- FinalDEGsL[order(FinalDEGsL[, 8]), ] + FinalDEGsL[is.na(FinalDEGsL[, 2]), c(2, 3)] <- + FinalDEGsL[is.na(FinalDEGsL[, 2]), 3] + if (!quiet) { + message( + "Low-regulated genes in the ", second[i], + " in ", first[i], " VS ", second[i], "\n" + ) + } + if (export) { + write.csv( + FinalDEGsL, + file = paste0( + "Low-regulated-", name, second[i], "in", + first[i], "VS", second[i], ".csv" + ) + ) + write.csv( + FinalDEGsL, + file = paste0( + "Up-regulated-", name, first[i], "in", first[i], + "VS", second[i], ".csv" + ) + ) + } + DEGsS <- c(DEGsS, FinalDEGsL[, 2]) + DEGsE <- c(DEGsE, as.character(FinalDEGsL[, 3])) + } + if (length(FDRu) > 0) { + genes <- siggenes.table$genes.up[, 3] + if (quiet) { + suppressMessages( + geneList <- + AnnotationDbi::select( + org.Hs.eg.db, + keys = keys(org.Hs.eg.db), + columns = c("SYMBOL", "ENSEMBL") + ) + ) + GL <- c(1, "MTRNR2", "ENSG00000210082") + GL1 <- c(1, "MTRNR1", "ENSG00000211459") + geneList <- rbind(geneList, GL, GL1) + } else { + geneList <- + AnnotationDbi::select( + org.Hs.eg.db, + keys = keys(org.Hs.eg.db), + columns = c("SYMBOL", "ENSEMBL") + ) + GL <- c(1, "MTRNR2", "ENSG00000210082") + GL1 <- c(1, "MTRNR1", "ENSG00000211459") + geneList <- rbind(geneList, GL, GL1) + } + FinalDEGsU <- cbind(genes, siggenes.table$genes.up) + gene_list <- geneList[, 3] + idx_genes <- is.element(gene_list, genes) + genes2 <- geneList[idx_genes, ] + FinalDEGsU <- merge( + FinalDEGsU, genes2, + by.x = "genes", + by.y = "ENSEMBL", all.x = TRUE + ) + FinalDEGsU[, 3] <- FinalDEGsU[, 11] + FinalDEGsU <- FinalDEGsU[, c(-1, -10, -11)] + FinalDEGsU <- FinalDEGsU[order(FinalDEGsU[, 8]), ] + FinalDEGsU[is.na(FinalDEGsU[, 2]), c(2, 3)] <- + FinalDEGsU[is.na(FinalDEGsU[, 2]), 3] + if (!quiet) { + message( + "Up-regulated genes in the ", second[i], " in ", + first[i], " VS ", second[i], "\n" + ) + } + if (export) { + write.csv( + FinalDEGsU, + file = paste0( + "Up-regulated-", name, second[i], "in", + first[i], "VS", second[i], ".csv" + ) + ) + write.csv( + FinalDEGsU, + file = paste0( + "Low-regulated-", name, first[i], "in", + first[i], "VS", second[i], ".csv" + ) + ) + } + DEGsS <- c(DEGsS, FinalDEGsU[, 2]) + DEGsE <- c(DEGsE, as.character(FinalDEGsU[, 3])) + } + } else { + DEGsTable[i, 1] <- paste0(first[i], " VS ", second[i]) + DEGsTable[i, 2] <- second[i] + DEGsTable[i, 3] <- NA + DEGsTable[i, 4] <- NA + DEGsTable[i, 5] <- NA + DEGsTable[i, 6] <- NA + + DEGsTable[i + com, 1] <- paste0(first[i], " VS ", second[i]) + DEGsTable[i + com, 2] <- first[i] + DEGsTable[i + com, 3] <- NA + DEGsTable[i + com, 4] <- NA + DEGsTable[i + com, 5] <- NA + DEGsTable[i + com, 6] <- NA + } + } + + if (!quiet & export) { + message("The results of DEGs are saved in your directory") + } + colnames(DEGsTable) <- c( + "Comparisons", "Target cluster", "Gene number", "File name", + "Gene number", "File name" + ) + if (export) write.csv(DEGsTable, file = paste0(filename_deg, ".csv")) + if (!quiet) print(DEGsTable) + sigDEG <- cbind(DEGsE, DEGsS) + if (export) write.csv(sigDEG, file = paste0(filename_sigdeg, ".csv")) + return(list(sigDEG, DEGsTable)) + } +)