diff --git a/NAMESPACE b/NAMESPACE index 58ba02f6..693d1ba8 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,45 +1,11 @@ # Generated by roxygen2: do not edit by hand -export("%w/o%") export(BuildClusterFilteredNetwork) -export(BuildPathwayCrosstalkNetwork) export(EvaluateClusters) export(GetGeneMANIA.edges) export(GetKinsub.edges) export(GetSTRINGdb.edges) -export(GraphCfn) export(MakeClusterList) export(MakeCorrelationNetwork) export(MakeDBInput) -export(NodeEdgeKey) -export(ReadBioplanetFile) export(StandardizeGeneSymbols) -export(connectNodes.all) -export(cytoscape.graph.PCN.pathways) -export(filter.edges.0) -export(filter.edges.1) -export(filter.edges.between) -export(fix.excel) -export(get.co.clustered.ptms) -export(getCyEdgeNames) -export(graph.ptm.by.cluster) -export(harmonize_cfs) -export(make.cytoscape.node.file) -export(make.gene.data.from.ptmtable) -export(make.genepep.edges) -export(merge2cols) -export(mergeEdges) -export(name.peptide) -export(outersect) -export(ptms_to_cfn) -export(remove.autophos) -export(setCorrEdgeAppearance) -export(setEdgeWidths) -export(setNodeColorToRatios) -export(setNodeColorToRowz) -export(setNodeMapping) -export(setNodeSizeColorIndependently) -export(strip.cy.goo) -importFrom(gplots,heatmap.2) -importFrom(grDevices,colorRampPalette) -importFrom(grDevices,hcl) diff --git a/R/BuildClusterFilteredNetwork.R b/R/BuildClusterFilteredNetwork.R index b832b449..0b384c7a 100644 --- a/R/BuildClusterFilteredNetwork.R +++ b/R/BuildClusterFilteredNetwork.R @@ -20,43 +20,42 @@ #' @examples #' Example_Output <- BuildClusterFilteredNetwork(ex.gene.cccn, ex.stringdb.edges, ex.gm.edges) #' utils::head(Example_Output) - BuildClusterFilteredNetwork <- function(gene.cccn.edges, stringdb.edges = NULL, genemania.edges = NULL, kinsub.edges = NULL, db.filepaths = c(NULL)) { - - # Combine PPIs from different databases - # First Normalize Weights - if (!is.null(stringdb.edges)) { - stringdb.edges$Weight <- 100 * stringdb.edges$Weight / max(stringdb.edges$Weight, na.rm = TRUE) } # this returns a range of 0 to 100 - if (!is.null(genemania.edges)) { - genemania.edges$Weight <- 100 * genemania.edges$Weight / max(genemania.edges$Weight, na.rm = TRUE) } # this returns a range of 0 to 100 - - # Note: if additional PPI data is desired, duplicate the normalization of Weights above and add the additional edge file(s) here (and in the function header) - # Done see below - - # Combine gathered PPI edges into one data frame - combined.PPIs <- rbind(stringdb.edges, genemania.edges, kinsub.edges) - - if(!is.null(db.filepaths)){ - for(path in db.filepaths){ - db.edges <- utils::read.table(path) - db.edges$Weight <- 100 * db.edges$Weight / max(db.edges$Weight, na.rm = TRUE) # return a range of 0 to 100 - combined.PPIs <- rbind(combined.PPIs, db.edges) # add it to combined.PPIs - + # Combine PPIs from different databases + # First Normalize Weights + if (!is.null(stringdb.edges)) { + stringdb.edges$Weight <- 100 * stringdb.edges$Weight / max(stringdb.edges$Weight, na.rm = TRUE) + } # this returns a range of 0 to 100 + if (!is.null(genemania.edges)) { + genemania.edges$Weight <- 100 * genemania.edges$Weight / max(genemania.edges$Weight, na.rm = TRUE) + } # this returns a range of 0 to 100 + + # Note: if additional PPI data is desired, duplicate the normalization of Weights above and add the additional edge file(s) here (and in the function header) + # Done see below + + # Combine gathered PPI edges into one data frame + combined.PPIs <- rbind(stringdb.edges, genemania.edges, kinsub.edges) + + if (!is.null(db.filepaths)) { + for (path in db.filepaths) { + db.edges <- utils::read.table(path) + db.edges$Weight <- 100 * db.edges$Weight / max(db.edges$Weight, na.rm = TRUE) # return a range of 0 to 100 + combined.PPIs <- rbind(combined.PPIs, db.edges) # add it to combined.PPIs + } } - } - cfn1 <- merge(gene.cccn.edges[,c("source", "target")], combined.PPIs, by=c("source", "target")) + cfn1 <- merge(gene.cccn.edges[, c("source", "target")], combined.PPIs, by = c("source", "target")) - # Undirected edges may be reversed in their order so merge the other way around. - reversed <- combined.PPIs - reversed <- reversed[ , c("target", "source", setdiff(names(gene.cccn.edges), c("source", "target")))] - colnames(reversed)[seq_len(2)] <- c("source", "target") # Rename for merge compatibility - cfn2 <- merge(reversed[,c("source", "target")], combined.PPIs, by=c("source", "target")) + # Undirected edges may be reversed in their order so merge the other way around. + reversed <- combined.PPIs + reversed <- reversed[, c("target", "source", setdiff(names(gene.cccn.edges), c("source", "target")))] + colnames(reversed)[seq_len(2)] <- c("source", "target") # Rename for merge compatibility + cfn2 <- merge(reversed[, c("source", "target")], combined.PPIs, by = c("source", "target")) - # Combine both (removing redundant rows if needed) - cfn <- rbind(cfn1, cfn2) - cfn <- unique(cfn) + # Combine both (removing redundant rows if needed) + cfn <- rbind(cfn1, cfn2) + cfn <- unique(cfn) - return(list(combined.PPIs, cfn)) + return(list(combined.PPIs, cfn)) } -#______________________________________________________________________________________________________________________________ +# ______________________________________________________________________________________________________________________________ diff --git a/R/EvaluateClusters.R b/R/EvaluateClusters.R index f89d9ec9..2697e12d 100644 --- a/R/EvaluateClusters.R +++ b/R/EvaluateClusters.R @@ -113,220 +113,232 @@ #' #' # Unordered cohort data (e.g. BRCA): disable slope filter #' eval_brca <- EvaluateClusters( -#' brca_clusters, brca_tbl, -#' data.type = "ratio", -#' use.slope = FALSE, -#' index.mode = "density" +#' brca_clusters, brca_tbl, +#' data.type = "ratio", +#' use.slope = FALSE, +#' index.mode = "density" #' ) #' #' # LINCS data with ratio columns to strip, size-based index: #' eval_lincs <- EvaluateClusters( -#' lincs_clusters, lincs_tbl, -#' data.type = "ratio", -#' index.mode = "size", -#' use.slope = TRUE, -#' ratio.col.pattern = "to" +#' lincs_clusters, lincs_tbl, +#' data.type = "ratio", +#' index.mode = "size", +#' use.slope = TRUE, +#' ratio.col.pattern = "to" #' ) EvaluateClusters <- function( - clusterlist, - tbl.sc, - data.type = c("ratio", "intensity", "count"), - index.mode = c("density", "size"), - use.slope = TRUE, - ratio.col.pattern = NULL, - verbose = TRUE + clusterlist, + tbl.sc, + data.type = c("ratio", "intensity", "count"), + index.mode = c("density", "size"), + use.slope = TRUE, + ratio.col.pattern = NULL, + verbose = TRUE ) { + # ---- argument validation ------------------------------------------------ + data.type <- match.arg(data.type) + index.mode <- match.arg(index.mode) - # ---- argument validation ------------------------------------------------ - data.type <- match.arg(data.type) - index.mode <- match.arg(index.mode) - - if (!is.list(clusterlist) || length(clusterlist) == 0L) { - stop("'clusterlist' must be a non-empty list.") - } - if (!is.data.frame(tbl.sc) && !is.matrix(tbl.sc)) { - stop("'tbl.sc' must be a data frame or matrix.") - } - if (is.null(rownames(tbl.sc))) { - stop("'tbl.sc' must have row names identifying genes / PTMs.") - } - - # ---- helper: number of non-missing values in a vector ------------------ - .filled <- function(x) sum(!is.na(x)) - .nmissing <- function(x) sum( is.na(x)) - - # ---- helper: extract gene names from one cluster element --------------- - # Accepts either a character vector of names or a data frame with a - # PTMnames column (MakeClusterList output format). - .gene.names <- function(cl.el) { - if (is.data.frame(cl.el)) { - if ("PTMnames" %in% names(cl.el)) return(cl.el[["PTMnames"]]) - if ("Gene.Name" %in% names(cl.el)) return(cl.el[["Gene.Name"]]) - stop("Cluster data frame must contain a 'PTMnames' or 'Gene.Name' column.") + if (!is.list(clusterlist) || length(clusterlist) == 0L) { + stop("'clusterlist' must be a non-empty list.") } - if (is.character(cl.el)) return(cl.el) - stop("Each element of 'clusterlist' must be a character vector or data frame.") - } - - # ---- helper: linear slope of a numeric vector against its index -------- - # Returns NA if fewer than 2 non-missing values are available. - .get.slope <- function(x) { - idx <- which(!is.na(x)) - if (length(idx) < 2L) return(NA_real_) - stats::coef(stats::lm(x[idx] ~ idx))[2L] - } - - # ---- initialise results data frame ------------------------------------- - n.clusters <- length(clusterlist) - evaluation <- data.frame( - Group = seq_len(n.clusters), - no.genes = NA_integer_, - culled.by.slope = NA_integer_, - percent.singlesamplegenes = NA_real_, - no.samples = NA_integer_, - percent.singlegenesamples = NA_real_, - total.signal = NA_real_, - percent.NA = NA_real_, - intensity = NA_real_, - Index = NA_real_, - stringsAsFactors = FALSE - ) - - # ---- main loop ---------------------------------------------------------- - for (i in seq_len(n.clusters)) { - - if (verbose) message("Starting Group ", i) - - gene.names <- .gene.names(clusterlist[[i]]) - evaluation$no.genes[i] <- length(gene.names) - - # Extract the cluster sub-table. - # When the cluster has a single gene, data.frame() is needed to prevent - # drop to a vector. - if (length(gene.names) == 1L) { - at <- data.frame(tbl.sc[gene.names, , drop = FALSE]) - } else { - at <- data.frame(tbl.sc[rownames(tbl.sc) %in% gene.names, , drop = FALSE]) + if (!is.data.frame(tbl.sc) && !is.matrix(tbl.sc)) { + stop("'tbl.sc' must be a data frame or matrix.") } - - # Optionally remove ratio descriptor columns (e.g. LINCS "H1_to_control"). - if (!is.null(ratio.col.pattern)) { - ratio.cols <- grep(ratio.col.pattern, names(at)) - if (length(ratio.cols) > 0L) at <- at[, -ratio.cols, drop = FALSE] + if (is.null(rownames(tbl.sc))) { + stop("'tbl.sc' must have row names identifying genes / PTMs.") } - # Apply absolute value for ratio data so that signal magnitude is - # assessed without sign cancellation. - if (data.type == "ratio") at <- abs(at) - - # Restrict to samples that have at least one non-missing value. - has.data <- which(vapply(at, .filled, integer(1L)) > 0L) - n.has.data <- length(has.data) - - if (n.has.data > 1L) { - acol <- names(at)[has.data] - } else { - acol <- if (n.has.data == 1L) names(at)[has.data] else character(0L) + # ---- helper: number of non-missing values in a vector ------------------ + .filled <- function(x) sum(!is.na(x)) + .nmissing <- function(x) sum(is.na(x)) + + # ---- helper: extract gene names from one cluster element --------------- + # Accepts either a character vector of names or a data frame with a + # PTMnames column (MakeClusterList output format). + .gene.names <- function(cl.el) { + if (is.data.frame(cl.el)) { + if ("PTMnames" %in% names(cl.el)) { + return(cl.el[["PTMnames"]]) + } + if ("Gene.Name" %in% names(cl.el)) { + return(cl.el[["Gene.Name"]]) + } + stop("Cluster data frame must contain a 'PTMnames' or 'Gene.Name' column.") + } + if (is.character(cl.el)) { + return(cl.el) + } + stop("Each element of 'clusterlist' must be a character vector or data frame.") } - evaluation$no.samples[i] <- length(acol) - # Restrict columns to those with data before further calculations. - if (length(acol) > 0L) { - at <- at[, acol, drop = FALSE] + # ---- helper: linear slope of a numeric vector against its index -------- + # Returns NA if fewer than 2 non-missing values are available. + .get.slope <- function(x) { + idx <- which(!is.na(x)) + if (length(idx) < 2L) { + return(NA_real_) + } + stats::coef(stats::lm(x[idx] ~ idx))[2L] } - evaluation$total.signal[i] <- sum(at, na.rm = TRUE) - - # ---- degenerate case: single sample or single gene ------------------ - # Diagnostics cannot be computed meaningfully; mark as 100 % single. - if (length(acol) <= 1L || length(gene.names) == 1L) { - evaluation$culled.by.slope[i] <- length(gene.names) - evaluation$percent.NA[i] <- 0 - evaluation$percent.singlesamplegenes[i] <- 100 - evaluation$percent.singlegenesamples[i] <- 100 - + # ---- initialise results data frame ------------------------------------- + n.clusters <- length(clusterlist) + evaluation <- data.frame( + Group = seq_len(n.clusters), + no.genes = NA_integer_, + culled.by.slope = NA_integer_, + percent.singlesamplegenes = NA_real_, + no.samples = NA_integer_, + percent.singlegenesamples = NA_real_, + total.signal = NA_real_, + percent.NA = NA_real_, + intensity = NA_real_, + Index = NA_real_, + stringsAsFactors = FALSE + ) + + # ---- main loop ---------------------------------------------------------- + for (i in seq_len(n.clusters)) { + if (verbose) message("Starting Group ", i) + + gene.names <- .gene.names(clusterlist[[i]]) + evaluation$no.genes[i] <- length(gene.names) + + # Extract the cluster sub-table. + # When the cluster has a single gene, data.frame() is needed to prevent + # drop to a vector. + if (length(gene.names) == 1L) { + at <- data.frame(tbl.sc[gene.names, , drop = FALSE]) + } else { + at <- data.frame(tbl.sc[rownames(tbl.sc) %in% gene.names, , drop = FALSE]) + } + + # Optionally remove ratio descriptor columns (e.g. LINCS "H1_to_control"). + if (!is.null(ratio.col.pattern)) { + ratio.cols <- grep(ratio.col.pattern, names(at)) + if (length(ratio.cols) > 0L) at <- at[, -ratio.cols, drop = FALSE] + } + + # Apply absolute value for ratio data so that signal magnitude is + # assessed without sign cancellation. + if (data.type == "ratio") at <- abs(at) + + # Restrict to samples that have at least one non-missing value. + has.data <- which(vapply(at, .filled, integer(1L)) > 0L) + n.has.data <- length(has.data) + + if (n.has.data > 1L) { + acol <- names(at)[has.data] + } else { + acol <- if (n.has.data == 1L) names(at)[has.data] else character(0L) + } + evaluation$no.samples[i] <- length(acol) + + # Restrict columns to those with data before further calculations. + if (length(acol) > 0L) { + at <- at[, acol, drop = FALSE] + } + + evaluation$total.signal[i] <- sum(at, na.rm = TRUE) + + # ---- degenerate case: single sample or single gene ------------------ + # Diagnostics cannot be computed meaningfully; mark as 100 % single. + if (length(acol) <= 1L || length(gene.names) == 1L) { + evaluation$culled.by.slope[i] <- length(gene.names) + evaluation$percent.NA[i] <- 0 + evaluation$percent.singlesamplegenes[i] <- 100 + evaluation$percent.singlegenesamples[i] <- 100 + } else { + # ---- full diagnostic block ---------------------------------------- + + evaluation$percent.NA[i] <- + 100 * sum(vapply(at, .nmissing, integer(1L))) / + (nrow(at) * ncol(at)) + + # Genes present in only one sample. + single.sample.genes <- at[vapply(seq_len(nrow(at)), function(r) { + .filled(as.numeric(at[r, ])) == 1L + }, logical(1L)), , drop = FALSE] + evaluation$percent.singlesamplegenes[i] <- + 100 * nrow(single.sample.genes) / nrow(at) + + # Samples that contain only one gene. + single.gene.samples <- + sum(vapply(at, .filled, integer(1L)) == 1L) + evaluation$percent.singlegenesamples[i] <- + 100 * single.gene.samples / ncol(at) + + # ---- slope filter ------------------------------------------------- + if (use.slope) { + # Sort columns by decreasing total signal, then rows by decreasing + # total signal, to create a monotonically organised matrix for slope + # assessment. + col.order <- order(-vapply( + at, function(x) sum(x, na.rm = TRUE), + numeric(1L) + )) + row.order <- order(-vapply(seq_len(nrow(at)), function(r) { + sum(as.numeric(at[r, ]), na.rm = TRUE) + }, numeric(1L))) + cluster.mo <- at[row.order, col.order, drop = FALSE] + + slopes <- apply(cluster.mo, 1, .get.slope) + bad.slope <- names(which(is.na(slopes) | slopes > 0)) + + if (verbose && length(bad.slope) > 0L) { + message(" ", length(bad.slope), " gene(s) culled by slope filter.") + } + + evaluation$culled.by.slope[i] <- length(bad.slope) + } else { + evaluation$culled.by.slope[i] <- 0L + } + } # end full diagnostic block + } # end cluster loop + + # ---- compute index for all clusters ------------------------------------ + + cleargenes <- evaluation$no.genes - evaluation$culled.by.slope + realsamples <- evaluation$no.samples - + (evaluation$no.samples * evaluation$percent.singlegenesamples / 100) + intensity <- evaluation$total.signal - + (evaluation$total.signal * evaluation$percent.NA / 100) + + evaluation$intensity <- intensity + + if (index.mode == "density") { + # Reward data density per gene: dividing by no.genes penalises large + # diffuse clusters relative to compact, information-rich ones. + evaluation$Index <- + (1 + realsamples) * (1 + cleargenes) / + ((1 + evaluation$percent.NA) * evaluation$no.genes) } else { - # ---- full diagnostic block ---------------------------------------- - - evaluation$percent.NA[i] <- - 100 * sum(vapply(at, .nmissing, integer(1L))) / - (nrow(at) * ncol(at)) - - # Genes present in only one sample. - single.sample.genes <- at[vapply(seq_len(nrow(at)), function(r) - .filled(as.numeric(at[r, ])) == 1L, logical(1L)), , drop = FALSE] - evaluation$percent.singlesamplegenes[i] <- - 100 * nrow(single.sample.genes) / nrow(at) - - # Samples that contain only one gene. - single.gene.samples <- - sum(vapply(at, .filled, integer(1L)) == 1L) - evaluation$percent.singlegenesamples[i] <- - 100 * single.gene.samples / ncol(at) - - # ---- slope filter ------------------------------------------------- - if (use.slope) { - # Sort columns by decreasing total signal, then rows by decreasing - # total signal, to create a monotonically organised matrix for slope - # assessment. - col.order <- order(-vapply(at, function(x) sum(x, na.rm = TRUE), - numeric(1L))) - row.order <- order(-vapply(seq_len(nrow(at)), function(r) - sum(as.numeric(at[r, ]), na.rm = TRUE), numeric(1L))) - cluster.mo <- at[row.order, col.order, drop = FALSE] - - slopes <- apply(cluster.mo, 1, .get.slope) - bad.slope <- names(which(is.na(slopes) | slopes > 0)) - - if (verbose && length(bad.slope) > 0L) - message(" ", length(bad.slope), " gene(s) culled by slope filter.") - - evaluation$culled.by.slope[i] <- length(bad.slope) - } else { - evaluation$culled.by.slope[i] <- 0L - } - } # end full diagnostic block - - } # end cluster loop - - # ---- compute index for all clusters ------------------------------------ - - cleargenes <- evaluation$no.genes - evaluation$culled.by.slope - realsamples <- evaluation$no.samples - - (evaluation$no.samples * evaluation$percent.singlegenesamples / 100) - intensity <- evaluation$total.signal - - (evaluation$total.signal * evaluation$percent.NA / 100) - - evaluation$intensity <- intensity - - if (index.mode == "density") { - # Reward data density per gene: dividing by no.genes penalises large - # diffuse clusters relative to compact, information-rich ones. - evaluation$Index <- - (1 + realsamples) * (1 + cleargenes) / - ((1 + evaluation$percent.NA) * evaluation$no.genes) - } else { - # "size": reward large, signal-dense clusters. - evaluation$Index <- - intensity * (1 + realsamples) * (1 + cleargenes) / - (1 + evaluation$percent.NA) - } + # "size": reward large, signal-dense clusters. + evaluation$Index <- + intensity * (1 + realsamples) * (1 + cleargenes) / + (1 + evaluation$percent.NA) + } - # Drop the intensity column when using the size index (it is already - # embedded in the index formula and clutters the output). - output.cols <- c("Group", "no.genes", "culled.by.slope", - "percent.singlesamplegenes", "no.samples", - "percent.singlegenesamples", "total.signal", - "percent.NA", "Index") - if (index.mode == "density") { - output.cols <- append(output.cols, "intensity", - after = which(output.cols == "percent.NA")) - } + # Drop the intensity column when using the size index (it is already + # embedded in the index formula and clutters the output). + output.cols <- c( + "Group", "no.genes", "culled.by.slope", + "percent.singlesamplegenes", "no.samples", + "percent.singlegenesamples", "total.signal", + "percent.NA", "Index" + ) + if (index.mode == "density") { + output.cols <- append(output.cols, "intensity", + after = which(output.cols == "percent.NA") + ) + } - eval.sort <- evaluation[ - order(-evaluation$Index, evaluation$percent.NA), - output.cols - ] + eval.sort <- evaluation[ + order(-evaluation$Index, evaluation$percent.NA), + output.cols + ] - return(eval.sort) + return(eval.sort) } diff --git a/R/GatherPPIData.R b/R/GatherPPIData.R index 47625169..24725147 100644 --- a/R/GatherPPIData.R +++ b/R/GatherPPIData.R @@ -27,45 +27,44 @@ #' # sym.map <- StandardizeGeneSymbols(c("EPRS", "QARS", "DDR1", "DDR2")) #' # unique(sym.map$standard_symbol) StandardizeGeneSymbols <- function(genes, - species = 9606, - string.version = "12.0", - keep.unmapped = TRUE) { - - if (!requireNamespace("STRINGdb", quietly = TRUE)) { - stop("Please install STRINGdb: BiocManager::install('STRINGdb')") - } - - if (!is.character(genes) || length(genes) == 0) { - stop("`genes` must be a non-empty character vector.") - } - - genes <- unique(as.character(genes)) - input.df <- data.frame(Gene.Names = genes, stringsAsFactors = FALSE) - - string.db <- STRINGdb::STRINGdb$new( - version = string.version, - species = species, - score_threshold = 0, - network_type = "full", - link_data = "full", - input_directory = "" - ) - - mapped <- string.db$map(input.df, "Gene.Names", removeUnmappedRows = FALSE) - - out <- data.frame( - input_symbol = mapped$Gene.Names, - STRING_id = if ("STRING_id" %in% colnames(mapped)) mapped$STRING_id else NA_character_, - standard_symbol = if ("preferred_name" %in% colnames(mapped)) mapped$preferred_name else NA_character_, - mapped = !is.na(if ("STRING_id" %in% colnames(mapped)) mapped$STRING_id else NA_character_), - stringsAsFactors = FALSE - ) - - if (keep.unmapped) { - out$standard_symbol[is.na(out$standard_symbol)] <- out$input_symbol[is.na(out$standard_symbol)] - } - - out + species = 9606, + string.version = "12.0", + keep.unmapped = TRUE) { + if (!requireNamespace("STRINGdb", quietly = TRUE)) { + stop("Please install STRINGdb: BiocManager::install('STRINGdb')") + } + + if (!is.character(genes) || length(genes) == 0) { + stop("`genes` must be a non-empty character vector.") + } + + genes <- unique(as.character(genes)) + input.df <- data.frame(Gene.Names = genes, stringsAsFactors = FALSE) + + string.db <- STRINGdb::STRINGdb$new( + version = string.version, + species = species, + score_threshold = 0, + network_type = "full", + link_data = "full", + input_directory = "" + ) + + mapped <- string.db$map(input.df, "Gene.Names", removeUnmappedRows = FALSE) + + out <- data.frame( + input_symbol = mapped$Gene.Names, + STRING_id = if ("STRING_id" %in% colnames(mapped)) mapped$STRING_id else NA_character_, + standard_symbol = if ("preferred_name" %in% colnames(mapped)) mapped$preferred_name else NA_character_, + mapped = !is.na(if ("STRING_id" %in% colnames(mapped)) mapped$STRING_id else NA_character_), + stringsAsFactors = FALSE + ) + + if (keep.unmapped) { + out$standard_symbol[is.na(out$standard_symbol)] <- out$input_symbol[is.na(out$standard_symbol)] + } + + out } @@ -87,26 +86,26 @@ StandardizeGeneSymbols <- function(genes, #' #' @keywords internal .map_nodes_with_symbol_map <- function(gene.cccn.nodes, symbol.map = NULL) { - nodes <- unique(as.character(gene.cccn.nodes)) - - if (is.null(symbol.map)) { - return(nodes) - } - - required.map.cols <- c("input_symbol", "standard_symbol") - missing.map.cols <- setdiff(required.map.cols, colnames(symbol.map)) - if (length(missing.map.cols) > 0) { - stop( - "`symbol.map` is missing required columns: ", - paste(missing.map.cols, collapse = ", ") - ) - } + nodes <- unique(as.character(gene.cccn.nodes)) + + if (is.null(symbol.map)) { + return(nodes) + } - idx <- match(nodes, symbol.map$input_symbol) - mapped.nodes <- symbol.map$standard_symbol[idx] - mapped.nodes[is.na(mapped.nodes)] <- nodes[is.na(mapped.nodes)] + required.map.cols <- c("input_symbol", "standard_symbol") + missing.map.cols <- setdiff(required.map.cols, colnames(symbol.map)) + if (length(missing.map.cols) > 0) { + stop( + "`symbol.map` is missing required columns: ", + paste(missing.map.cols, collapse = ", ") + ) + } - unique(mapped.nodes) + idx <- match(nodes, symbol.map$input_symbol) + mapped.nodes <- symbol.map$standard_symbol[idx] + mapped.nodes[is.na(mapped.nodes)] <- nodes[is.na(mapped.nodes)] + + unique(mapped.nodes) } #' Make Database Input File #' @@ -120,10 +119,10 @@ StandardizeGeneSymbols <- function(genes, #' @export #' #' @examples -#' #MakeDBInput(ex.nodenames) -#' cat(ex.nodenames[[1]], sep = '\n') +#' # MakeDBInput(ex.nodenames) +#' cat(ex.nodenames[[1]], sep = "\n") MakeDBInput <- function(gene.cccn.nodes, file.path.name = "db_nodes.txt") { - utils::write.table(unique(c(gene.cccn.nodes[[1]], gene.cccn.nodes[[2]])), file = file.path.name, row.names = FALSE, col.names = FALSE, quote = FALSE) + utils::write.table(unique(c(gene.cccn.nodes[[1]], gene.cccn.nodes[[2]])), file = file.path.name, row.names = FALSE, col.names = FALSE, quote = FALSE) } # Pulls nodenames from the gene.cccn # @@ -174,238 +173,245 @@ MakeDBInput <- function(gene.cccn.nodes, file.path.name = "db_nodes.txt") { #' # symbol.map = sym.map #' # ) GetSTRINGdb.edges <- function(gene.cccn.edges, - gene.cccn.nodes, - local = FALSE, - string.local.path = "string_hs_hugo_full.tsv", - combined.score.threshold = 400, - include.transferred = TRUE, - symbol.map = NULL) { - - if (!is.character(gene.cccn.nodes) || length(gene.cccn.nodes) == 0) { - stop("`gene.cccn.nodes` must be a non-empty character vector.") - } - - if (!is.numeric(combined.score.threshold) || - length(combined.score.threshold) != 1 || - is.na(combined.score.threshold) || - combined.score.threshold < 0 || - combined.score.threshold > 1000) { - stop("`combined.score.threshold` must be a single number between 0 and 1000.") - } - - # -------------------------------------------------------------------------- - # Standardize nodes if a symbol map is supplied - # -------------------------------------------------------------------------- - nodes <- .map_nodes_with_symbol_map(gene.cccn.nodes, symbol.map) - - - # -------------------------------------------------------------------------- - # Helper for empty returns - # -------------------------------------------------------------------------- - .empty_edges <- function() { - data.frame( - source = character(0), - target = character(0), - interaction = character(0), - Weight = numeric(0), - stringsAsFactors = FALSE - ) - } - - # -------------------------------------------------------------------------- - # Local mode: fully offline - # -------------------------------------------------------------------------- - if (local) { - if (!file.exists(string.local.path)) { - stop("Local STRING file not found: ", string.local.path, - "\nGenerate it from protein.links.full.v12.0.txt.gz first.") + gene.cccn.nodes, + local = FALSE, + string.local.path = "string_hs_hugo_full.tsv", + combined.score.threshold = 400, + include.transferred = TRUE, + symbol.map = NULL) { + if (!is.character(gene.cccn.nodes) || length(gene.cccn.nodes) == 0) { + stop("`gene.cccn.nodes` must be a non-empty character vector.") } - message("Reading local STRING file: ", string.local.path) - dt <- utils::read.delim(string.local.path, stringsAsFactors = FALSE) + if (!is.numeric(combined.score.threshold) || + length(combined.score.threshold) != 1 || + is.na(combined.score.threshold) || + combined.score.threshold < 0 || + combined.score.threshold > 1000) { + stop("`combined.score.threshold` must be a single number between 0 and 1000.") + } + + # -------------------------------------------------------------------------- + # Standardize nodes if a symbol map is supplied + # -------------------------------------------------------------------------- + nodes <- .map_nodes_with_symbol_map(gene.cccn.nodes, symbol.map) + + + # -------------------------------------------------------------------------- + # Helper for empty returns + # -------------------------------------------------------------------------- + .empty_edges <- function() { + data.frame( + source = character(0), + target = character(0), + interaction = character(0), + Weight = numeric(0), + stringsAsFactors = FALSE + ) + } + + # -------------------------------------------------------------------------- + # Local mode: fully offline + # -------------------------------------------------------------------------- + if (local) { + if (!file.exists(string.local.path)) { + stop( + "Local STRING file not found: ", string.local.path, + "\nGenerate it from protein.links.full.v12.0.txt.gz first." + ) + } + + message("Reading local STRING file: ", string.local.path) + dt <- utils::read.delim(string.local.path, stringsAsFactors = FALSE) + + required.cols <- c( + "source", "target", + "experiments", "database", + "combined_score" + ) + + if (include.transferred) { + required.cols <- c( + required.cols, + "experiments_transferred", + "database_transferred" + ) + } + + missing.cols <- setdiff(required.cols, colnames(dt)) + if (length(missing.cols) > 0) { + stop( + "Local STRING file is missing required columns: ", + paste(missing.cols, collapse = ", ") + ) + } + + dt <- dt[ + dt$combined_score >= combined.score.threshold & + dt$source %in% nodes & + dt$target %in% nodes, + ] + + if (nrow(dt) == 0) { + warning("No local STRING edges passed the filters.") + return(.empty_edges()) + } + + if (include.transferred) { + keep <- dt$experiments > 0 | + dt$experiments_transferred > 0 | + dt$database > 0 | + dt$database_transferred > 0 + } else { + keep <- dt$experiments > 0 | dt$database > 0 + } + + dt <- dt[keep, ] + + if (nrow(dt) == 0) { + warning("No local STRING edges with selected evidence types passed the filters.") + return(.empty_edges()) + } + + dt$edgeType <- "STRINGdb" + + if (include.transferred) { + dt[dt$database_transferred > 0, "edgeType"] <- "database_transferred" + dt[dt$database > 0, "edgeType"] <- "database" + dt[dt$experiments_transferred > 0, "edgeType"] <- "experiments_transferred" + dt[dt$experiments > 0, "edgeType"] <- "experiments" + dt$Weight <- rowSums( + dt[, c( + "experiments", "experiments_transferred", + "database", "database_transferred" + ), drop = FALSE] + ) + } else { + dt[dt$database > 0, "edgeType"] <- "database" + dt[dt$experiments > 0, "edgeType"] <- "experiments" + dt$Weight <- rowSums( + dt[, c("experiments", "database"), drop = FALSE] + ) + } + # After computing dt$edgeType and dt$Weight + + # Canonicalize undirected edges: sort endpoints within each row + dt$u <- pmin(dt$source, dt$target) + dt$v <- pmax(dt$source, dt$target) + + stringdb.edges <- unique(dt[, c("u", "v", "edgeType", "Weight")]) + colnames(stringdb.edges) <- c("source", "target", "interaction", "Weight") + stringdb.edges$source <- as.character(stringdb.edges$source) + stringdb.edges$target <- as.character(stringdb.edges$target) + rownames(stringdb.edges) <- NULL + return(stringdb.edges) + } + + # -------------------------------------------------------------------------- + # Live mode + # -------------------------------------------------------------------------- + if (!requireNamespace("STRINGdb", quietly = TRUE)) { + stop("Please install STRINGdb: BiocManager::install('STRINGdb')") + } - required.cols <- c( - "source", "target", - "experiments", "database", - "combined_score" + nodenames <- data.frame(Gene.Names = nodes, stringsAsFactors = FALSE) + + string.db <- STRINGdb::STRINGdb$new( + version = "12.0", + species = 9606, + score_threshold = 0, + network_type = "full", + link_data = "full", + input_directory = "" ) - if (include.transferred) { - required.cols <- c( - required.cols, - "experiments_transferred", - "database_transferred" - ) + message("Querying STRINGdb for interactions between ", length(nodes), " genes...") + string.proteins <- string.db$get_proteins() + + message("Mapping genes to STRING IDs...") + mapped.genes <- string.db$map(nodenames, "Gene.Names", removeUnmappedRows = TRUE) + + if (nrow(mapped.genes) == 0) { + warning("No input genes could be mapped to STRING.") + return(.empty_edges()) } - missing.cols <- setdiff(required.cols, colnames(dt)) - if (length(missing.cols) > 0) { - stop("Local STRING file is missing required columns: ", - paste(missing.cols, collapse = ", ")) + message("Retrieving interactions for mapped genes...") + interactions <- string.db$get_interactions(mapped.genes$STRING_id) + + if (nrow(interactions) == 0) { + warning("STRINGdb returned no interactions.") + return(.empty_edges()) + } + + message("Formatting...") + interactions$Gene.1 <- vapply(interactions$from, function(x) { + string.proteins[match(x, string.proteins$protein_external_id), "preferred_name"] + }, FUN.VALUE = character(1)) + interactions$Gene.2 <- vapply(interactions$to, function(x) { + string.proteins[match(x, string.proteins$protein_external_id), "preferred_name"] + }, FUN.VALUE = character(1)) + + if (!"combined_score" %in% colnames(interactions)) { + stop("Live STRINGdb interactions did not include a `combined_score` column.") } - dt <- dt[ - dt$combined_score >= combined.score.threshold & - dt$source %in% nodes & - dt$target %in% nodes, + interactions <- interactions[ + interactions$combined_score >= combined.score.threshold & + interactions$Gene.1 %in% nodes & + interactions$Gene.2 %in% nodes, ] - if (nrow(dt) == 0) { - warning("No local STRING edges passed the filters.") - return(.empty_edges()) + if (nrow(interactions) == 0) { + warning("No live STRING edges passed the filters.") + return(.empty_edges()) } if (include.transferred) { - keep <- dt$experiments > 0 | - dt$experiments_transferred > 0 | - dt$database > 0 | - dt$database_transferred > 0 + str.e <- interactions[interactions$experiments > 0, ] + str.et <- interactions[interactions$experiments_transferred > 0, ] + str.d <- interactions[interactions$database > 0, ] + str.dt <- interactions[interactions$database_transferred > 0, ] + combined_interactions <- unique(rbind(str.e, str.et, str.d, str.dt)) } else { - keep <- dt$experiments > 0 | dt$database > 0 + str.e <- interactions[interactions$experiments > 0, ] + str.d <- interactions[interactions$database > 0, ] + combined_interactions <- unique(rbind(str.e, str.d)) } - dt <- dt[keep, ] - - if (nrow(dt) == 0) { - warning("No local STRING edges with selected evidence types passed the filters.") - return(.empty_edges()) + if (nrow(combined_interactions) == 0) { + warning("No live STRING edges with selected evidence types passed the filters.") + return(.empty_edges()) } - dt$edgeType <- "STRINGdb" + combined_interactions$edgeType <- "STRINGdb" if (include.transferred) { - dt[dt$database_transferred > 0, "edgeType"] <- "database_transferred" - dt[dt$database > 0, "edgeType"] <- "database" - dt[dt$experiments_transferred > 0, "edgeType"] <- "experiments_transferred" - dt[dt$experiments > 0, "edgeType"] <- "experiments" - dt$Weight <- rowSums( - dt[, c("experiments", "experiments_transferred", - "database", "database_transferred"), drop = FALSE] - ) + combined_interactions[combined_interactions$database_transferred > 0, "edgeType"] <- "database_transferred" + combined_interactions[combined_interactions$database > 0, "edgeType"] <- "database" + combined_interactions[combined_interactions$experiments_transferred > 0, "edgeType"] <- "experiments_transferred" + combined_interactions[combined_interactions$experiments > 0, "edgeType"] <- "experiments" + + combined_interactions$Weight <- rowSums( + combined_interactions[, c( + "experiments", "experiments_transferred", + "database", "database_transferred" + ), drop = FALSE] + ) } else { - dt[dt$database > 0, "edgeType"] <- "database" - dt[dt$experiments > 0, "edgeType"] <- "experiments" - dt$Weight <- rowSums( - dt[, c("experiments", "database"), drop = FALSE] - ) - } - # After computing dt$edgeType and dt$Weight + combined_interactions[combined_interactions$database > 0, "edgeType"] <- "database" + combined_interactions[combined_interactions$experiments > 0, "edgeType"] <- "experiments" - # Canonicalize undirected edges: sort endpoints within each row - dt$u <- pmin(dt$source, dt$target) - dt$v <- pmax(dt$source, dt$target) + combined_interactions$Weight <- rowSums( + combined_interactions[, c("experiments", "database"), drop = FALSE] + ) + } - stringdb.edges <- unique(dt[, c("u", "v", "edgeType", "Weight")]) + stringdb.edges <- unique(combined_interactions[, c("Gene.1", "Gene.2", "edgeType", "Weight")]) colnames(stringdb.edges) <- c("source", "target", "interaction", "Weight") - stringdb.edges$source <- as.character(stringdb.edges$source) - stringdb.edges$target <- as.character(stringdb.edges$target) rownames(stringdb.edges) <- NULL - return(stringdb.edges) - } - - # -------------------------------------------------------------------------- - # Live mode - # -------------------------------------------------------------------------- - if (!requireNamespace("STRINGdb", quietly = TRUE)) { - stop("Please install STRINGdb: BiocManager::install('STRINGdb')") - } - - nodenames <- data.frame(Gene.Names = nodes, stringsAsFactors = FALSE) - - string.db <- STRINGdb::STRINGdb$new( - version = "12.0", - species = 9606, - score_threshold = 0, - network_type = "full", - link_data = "full", - input_directory = "" - ) - - message("Querying STRINGdb for interactions between ", length(nodes), " genes...") - string.proteins <- string.db$get_proteins() - - message("Mapping genes to STRING IDs...") - mapped.genes <- string.db$map(nodenames, "Gene.Names", removeUnmappedRows = TRUE) - - if (nrow(mapped.genes) == 0) { - warning("No input genes could be mapped to STRING.") - return(.empty_edges()) - } - - message("Retrieving interactions for mapped genes...") - interactions <- string.db$get_interactions(mapped.genes$STRING_id) - - if (nrow(interactions) == 0) { - warning("STRINGdb returned no interactions.") - return(.empty_edges()) - } - - message("Formatting...") - interactions$Gene.1 <- vapply(interactions$from, function(x) { - string.proteins[match(x, string.proteins$protein_external_id), "preferred_name"] - }, FUN.VALUE=character(1)) - interactions$Gene.2 <- vapply(interactions$to, function(x) { - string.proteins[match(x, string.proteins$protein_external_id), "preferred_name"] - }, FUN.VALUE=character(1)) - - if (!"combined_score" %in% colnames(interactions)) { - stop("Live STRINGdb interactions did not include a `combined_score` column.") - } - - interactions <- interactions[ - interactions$combined_score >= combined.score.threshold & - interactions$Gene.1 %in% nodes & - interactions$Gene.2 %in% nodes, - ] - - if (nrow(interactions) == 0) { - warning("No live STRING edges passed the filters.") - return(.empty_edges()) - } - - if (include.transferred) { - str.e <- interactions[interactions$experiments > 0, ] - str.et <- interactions[interactions$experiments_transferred > 0, ] - str.d <- interactions[interactions$database > 0, ] - str.dt <- interactions[interactions$database_transferred > 0, ] - combined_interactions <- unique(rbind(str.e, str.et, str.d, str.dt)) - } else { - str.e <- interactions[interactions$experiments > 0, ] - str.d <- interactions[interactions$database > 0, ] - combined_interactions <- unique(rbind(str.e, str.d)) - } - - if (nrow(combined_interactions) == 0) { - warning("No live STRING edges with selected evidence types passed the filters.") - return(.empty_edges()) - } - - combined_interactions$edgeType <- "STRINGdb" - - if (include.transferred) { - combined_interactions[combined_interactions$database_transferred > 0, "edgeType"] <- "database_transferred" - combined_interactions[combined_interactions$database > 0, "edgeType"] <- "database" - combined_interactions[combined_interactions$experiments_transferred > 0, "edgeType"] <- "experiments_transferred" - combined_interactions[combined_interactions$experiments > 0, "edgeType"] <- "experiments" - - combined_interactions$Weight <- rowSums( - combined_interactions[, c("experiments", "experiments_transferred", - "database", "database_transferred"), drop = FALSE] - ) - } else { - combined_interactions[combined_interactions$database > 0, "edgeType"] <- "database" - combined_interactions[combined_interactions$experiments > 0, "edgeType"] <- "experiments" - - combined_interactions$Weight <- rowSums( - combined_interactions[, c("experiments", "database"), drop = FALSE] - ) - } - - stringdb.edges <- unique(combined_interactions[, c("Gene.1", "Gene.2", "edgeType", "Weight")]) - colnames(stringdb.edges) <- c("source", "target", "interaction", "Weight") - rownames(stringdb.edges) <- NULL - stringdb.edges + stringdb.edges } @@ -440,85 +446,88 @@ GetSTRINGdb.edges <- function(gene.cccn.edges, #' # symbol.map = sym.map #' # ) GetGeneMANIA.edges <- function(gm.results.path, - gene.cccn.nodes, - local = FALSE, - genemania.local.path = "hs_interactions_hugo.tsv", - gm.interaction.types = c("Pathway", "Physical Interactions"), - symbol.map = NULL) { - - nodes <- .map_nodes_with_symbol_map(gene.cccn.nodes, symbol.map) - - if (local) { - if (!file.exists(genemania.local.path)) { - stop("Local GeneMANIA file not found: ", genemania.local.path, - "\nGenerate it with scripts/genemania_hs_download.r or set local = FALSE.") - } - - message("Reading local GeneMANIA file: ", genemania.local.path) - dt <- utils::read.delim(genemania.local.path, stringsAsFactors = FALSE) - - required.cols <- c("Gene1", "Gene2", "Group", "Weight") - missing.cols <- setdiff(required.cols, colnames(dt)) - if (length(missing.cols) > 0) { - stop("Local GeneMANIA file is missing required columns: ", - paste(missing.cols, collapse = ", ")) - } - - dt <- dt[dt$Group %in% gm.interaction.types, ] - dt <- dt[dt$Gene1 %in% nodes & dt$Gene2 %in% nodes, ] - - if (nrow(dt) == 0) { - warning("No local GeneMANIA edges passed the filters.") - return(data.frame( - source = character(0), - target = character(0), - interaction = character(0), - Weight = numeric(0), - stringsAsFactors = FALSE - )) + gene.cccn.nodes, + local = FALSE, + genemania.local.path = "hs_interactions_hugo.tsv", + gm.interaction.types = c("Pathway", "Physical Interactions"), + symbol.map = NULL) { + nodes <- .map_nodes_with_symbol_map(gene.cccn.nodes, symbol.map) + + if (local) { + if (!file.exists(genemania.local.path)) { + stop( + "Local GeneMANIA file not found: ", genemania.local.path, + "\nGenerate it with scripts/genemania_hs_download.r or set local = FALSE." + ) + } + + message("Reading local GeneMANIA file: ", genemania.local.path) + dt <- utils::read.delim(genemania.local.path, stringsAsFactors = FALSE) + + required.cols <- c("Gene1", "Gene2", "Group", "Weight") + missing.cols <- setdiff(required.cols, colnames(dt)) + if (length(missing.cols) > 0) { + stop( + "Local GeneMANIA file is missing required columns: ", + paste(missing.cols, collapse = ", ") + ) + } + + dt <- dt[dt$Group %in% gm.interaction.types, ] + dt <- dt[dt$Gene1 %in% nodes & dt$Gene2 %in% nodes, ] + + if (nrow(dt) == 0) { + warning("No local GeneMANIA edges passed the filters.") + return(data.frame( + source = character(0), + target = character(0), + interaction = character(0), + Weight = numeric(0), + stringsAsFactors = FALSE + )) + } + + dt <- unique(dt[, c("Gene1", "Gene2", "Group", "Weight")]) + colnames(dt) <- c("source", "target", "interaction", "Weight") + rownames(dt) <- NULL + return(dt) } - dt <- unique(dt[, c("Gene1", "Gene2", "Group", "Weight")]) - colnames(dt) <- c("source", "target", "interaction", "Weight") - rownames(dt) <- NULL - return(dt) - } + # ---- Live mode: original Cytoscape export behaviour ---- - # ---- Live mode: original Cytoscape export behaviour ---- + all_lines <- readLines(gm.results.path) - all_lines <- readLines(gm.results.path) + start_line <- grep("Weight\\tType", all_lines) + end_line <- grep("^Gene\\s+GO ids", all_lines) + if (length(end_line) == 0) end_line <- length(all_lines) + 1 - start_line <- grep("Weight\\tType", all_lines) - end_line <- grep("^Gene\\s+GO ids", all_lines) - if (length(end_line) == 0) end_line <- length(all_lines) + 1 + network_lines <- all_lines[start_line[1]:(end_line[1] - 1)] - network_lines <- all_lines[start_line[1]:(end_line[1] - 1)] - - edgetable <- read.table( - text = network_lines, - header = TRUE, - stringsAsFactors = FALSE, - sep = "\t", - comment.char = "#", - na.strings = "", - quote = "", - fill = TRUE - ) + edgetable <- read.table( + text = network_lines, + header = TRUE, + stringsAsFactors = FALSE, + sep = "\t", + comment.char = "#", + na.strings = "", + quote = "", + fill = TRUE + ) - keeper <- edgetable$Type %in% gm.interaction.types - edgetable <- edgetable[keeper, ] - edgetable <- edgetable[, c("Gene.1", "Gene.2", "Type", "Weight")] - colnames(edgetable) <- c("source", "target", "interaction", "Weight") + keeper <- edgetable$Type %in% gm.interaction.types + edgetable <- edgetable[keeper, ] + edgetable <- edgetable[, c("Gene.1", "Gene.2", "Type", "Weight")] + colnames(edgetable) <- c("source", "target", "interaction", "Weight") - keep <- edgetable$source %in% nodes & edgetable$target %in% nodes - genemania.edges <- unique(edgetable[keep, ]) + keep <- edgetable$source %in% nodes & edgetable$target %in% nodes + genemania.edges <- unique(edgetable[keep, ]) - if (nrow(genemania.edges) == 0) { - warning("No GeneMANIA edges passed the filters.") - } + if (nrow(genemania.edges) == 0) { + warning("No GeneMANIA edges passed the filters.") + } - rownames(genemania.edges) <- NULL - return(genemania.edges) + rownames(genemania.edges) <- NULL + return(genemania.edges) } #' Get Kinase-Substrate Edges @@ -544,52 +553,51 @@ GetGeneMANIA.edges <- function(gm.results.path, #' # symbol.map = sym.map #' # ) GetKinsub.edges <- function(kinasesubstrate.filename = "Kinase_Substrate_Dataset.txt", - gene.cccn.nodes, - symbol.map = NULL) { - - nodes <- toupper(.map_nodes_with_symbol_map(gene.cccn.nodes, symbol.map)) - - kinasesubstrateraw <- read.table( - kinasesubstrate.filename, - header = TRUE, - skip = 3, - stringsAsFactors = FALSE, - sep = "\t", - na.strings = "", - fill = TRUE - ) - - kinsub <- kinasesubstrateraw - - if (any(is.na(kinsub$GENE))) { - kinsub[which(is.na(kinsub$GENE)), "GENE"] <- - as.character(kinsub[which(is.na(kinsub$GENE)), "KINASE"]) - } - - if (any(is.na(kinsub$SUB_GENE))) { - kinsub[which(is.na(kinsub$SUB_GENE)), "SUB_GENE"] <- - as.character(kinsub[which(is.na(kinsub$SUB_GENE)), "SUBSTRATE"]) - } - - kinase <- toupper(kinsub$GENE) - substrate <- toupper(kinsub$SUB_GENE) - - kinsub <- unique(data.frame( - source = kinase, - target = substrate, - stringsAsFactors = FALSE - )) - - kinsub.edges <- kinsub[kinsub$source %in% nodes & kinsub$target %in% nodes, ] - kinsub.edges$interaction <- "pp" - kinsub.edges$Weight <- 1 - - kinsub.edges <- remove.autophos(kinsub.edges) - - if (nrow(kinsub.edges) == 0) { - warning("No kinase-substrate edges passed the filters.") - } - - rownames(kinsub.edges) <- NULL - return(kinsub.edges) + gene.cccn.nodes, + symbol.map = NULL) { + nodes <- toupper(.map_nodes_with_symbol_map(gene.cccn.nodes, symbol.map)) + + kinasesubstrateraw <- read.table( + kinasesubstrate.filename, + header = TRUE, + skip = 3, + stringsAsFactors = FALSE, + sep = "\t", + na.strings = "", + fill = TRUE + ) + + kinsub <- kinasesubstrateraw + + if (any(is.na(kinsub$GENE))) { + kinsub[which(is.na(kinsub$GENE)), "GENE"] <- + as.character(kinsub[which(is.na(kinsub$GENE)), "KINASE"]) + } + + if (any(is.na(kinsub$SUB_GENE))) { + kinsub[which(is.na(kinsub$SUB_GENE)), "SUB_GENE"] <- + as.character(kinsub[which(is.na(kinsub$SUB_GENE)), "SUBSTRATE"]) + } + + kinase <- toupper(kinsub$GENE) + substrate <- toupper(kinsub$SUB_GENE) + + kinsub <- unique(data.frame( + source = kinase, + target = substrate, + stringsAsFactors = FALSE + )) + + kinsub.edges <- kinsub[kinsub$source %in% nodes & kinsub$target %in% nodes, ] + kinsub.edges$interaction <- "pp" + kinsub.edges$Weight <- 1 + + kinsub.edges <- remove.autophos(kinsub.edges) + + if (nrow(kinsub.edges) == 0) { + warning("No kinase-substrate edges passed the filters.") + } + + rownames(kinsub.edges) <- NULL + return(kinsub.edges) } diff --git a/R/MakeClusterList.R b/R/MakeClusterList.R index fbe8ebb6..9a425876 100644 --- a/R/MakeClusterList.R +++ b/R/MakeClusterList.R @@ -23,173 +23,177 @@ #' @export #' #' @examples -#' Example_Output <- MakeClusterList(ex_tiny_ptm_table) #Run function -#' Example_Output[[1]][1:3] #Display data -#' #Do we want to have one for adj.consensus? Doesn't seem like it'd be very helpful to view. -#' utils::head(Example_Output[[3]][, c(1,2,3,4,5)]) #Display data -MakeClusterList <- function(ptmtable, keeplength = 2, toolong = 3.5, tsne_perplexity = 15, tsne_max_iter = 5000){ - start_time <- Sys.time() - message("Starting correlation calculations and t-SNE.") - message("This may take a few minutes or hours for large data sets.") - - # Add if statement here to make sure functions are formatted correctly # - # Ensure ptmtable is a data frame with numeric values # - PTMnames <- rownames(ptmtable) - if (!all(vapply(ptmtable, is.numeric, FUN.VALUE = logical(1) ))) { - stop("All columns in 'ptmtable' must be numeric.") - } - - # SPEARMAN CALCULATION - - # Calculate Spearman correlation # - ptm.correlation.matrix <- stats::cor(t(ptmtable), use = "pairwise.complete.obs", method = "spearman") - # Note: this is the slowest step. We found use = "pairwise.complete.obs", - # method = "spearman" to perform the best according to evaluations with data - # with missing values, but it takes longer. - - # Replace diagonal with NA # - diag(ptm.correlation.matrix) <- NA - - # Calculate dissimilarity # - sp.diss.matrix <- 1 - abs(ptm.correlation.matrix) - - # Handle any remaining NA values by setting them to the maximum dissimilarity # - max.dist.sp <- max(sp.diss.matrix, na.rm = TRUE) - - #Set NA values to 100 * the max distance - sp.diss.matrix[is.na(sp.diss.matrix)] <- 100 * max.dist.sp - message("Spearman correlation calculation complete after ", round(Sys.time() - start_time, 2), " ", units(Sys.time() - start_time), " total.") - - # Run t-SNE # - tsne.results <- Rtsne::Rtsne(sp.diss.matrix, dims = 3, perplexity = tsne_perplexity, theta = 0.25, max_iter = tsne_max_iter, check_duplicates = FALSE, pca = FALSE) - spearman.coords <- tsne.results$Y - message("Spearman t-SNE calculation complete after ", round(Sys.time() - start_time, 2), " ", units(Sys.time() - start_time), " total.") - - #EUCLIDEAN CALCULATION - - # Convert the dataframe to a distance matrix using Euclidean distance # - ptmtable.dist <- as.matrix(stats::dist(ptmtable, method = "euclidean")) - - # Compute the maximum distance in the matrix, excluding NA values # - max.dist <- max(ptmtable.dist, na.rm = TRUE) - - # Replace NA values in the distance matrix with 100 times the maximum distance # - ptmtable.dist[is.na(ptmtable.dist)] <- 100 * max.dist - - # Normalize the distance matrix by scaling it to a range from 0 to 100. This becomes the distance matrix for euclidian distance which we will run Rtsne on# - eu.dist.calc <- 100 * ptmtable.dist / max(ptmtable.dist, na.rm = TRUE) - eu.dist.calc <- as.matrix(eu.dist.calc) #Fix eu.dist.calc RQ - message("Euclidean distance calculation complete after ", round(Sys.time() - start_time, 2), " ", units(Sys.time() - start_time), " total.") - - # Run t-SNE # - eu.ptms.tsne.list <- Rtsne::Rtsne(eu.dist.calc, dims = 3, perplexity = tsne_perplexity, theta = 0.25, max_iter = tsne_max_iter, check_duplicates = FALSE, pca = FALSE) - euclidean.coords <- eu.ptms.tsne.list$Y - message("Euclidean t-SNE calculation complete after ", round(Sys.time() - start_time, 2), " ", units(Sys.time() - start_time), " total.") - - #COMBINED CALCULATION - - #fix spearman thing; so do the exact same thing but no absolute value - sp.diss.calc <- 1 - ptm.correlation.matrix # range goes from (-1 to 1) to (0 to 2) - max.diss.sp <- max(sp.diss.calc, na.rm = TRUE) # find the max value (around 2) - sp.diss.calc[is.na(sp.diss.calc)] <- 50 * max.diss.sp # Set every NA value to 50 * the max dissimilarity - sp.diss.calc <- as.matrix(sp.diss.calc) # turn into a matrix - - #find average - combined.distance <- (sp.diss.calc + eu.dist.calc) / 2 - message("Combined distance calculation complete after ", round(Sys.time() - start_time, 2), " ", units(Sys.time() - start_time), " total.") - # Perform t-SNE on the combined distances # - tsne.result <- Rtsne::Rtsne(combined.distance, dims = 3, perplexity = tsne_perplexity, theta = 0.25, max_iter = tsne_max_iter, check_duplicates = FALSE, pca = FALSE) - sed.coords <- tsne.result$Y - message("SED t-SNE calculation complete after ", round(Sys.time() - start_time, 2), " ", units(Sys.time() - start_time), " total.") - - - #Nested function to analyze result - group <- NULL #Gets rid of check note - clustercreate <- function(tsne_coords, distance_name){ - - # Find clusters that are connected - # this returns a vector with cluster assignment for each PTM - result.disc2 <- vegan::distconnected(stats::dist(tsne_coords), toolong = toolong, trace = FALSE) - - # Create a plot of the clusters, with convex hulll around each cluster - # note that tsne_coords is a matrix with 3 columns (the t-SNE dimensions), but only first two will get plotted - vegan::ordiplot(tsne_coords, display = c(), main=paste("t-SNE Cluster Plot", distance_name)) - vegan::ordihull(tsne_coords, result.disc2, col="red", lwd=2) - - # make a dataframe with one column for PTMs and another for group assignment - result.span.df <- data.frame(PTMnames = PTMnames, group= result.disc2) - - # Convert data frame into a list of clusters (check doesn't like group but it's a column name) - result.span.list <- plyr::dlply(result.span.df, plyr::.(group)) # GROUP LIST ! - - end_time <- Sys.time() - total_time <- end_time - start_time - message("Clustering for ", distance_name, " complete after ", round(total_time, 2), " ", units(total_time), " total.") - return(result.span.list) - } #END of nested function - - #Create all tsne coords data sturcture and give it names - all.tsne.coords <- list(euclidean.coords, spearman.coords, sed.coords) - names(all.tsne.coords) <- c("Euclidean", "Spearman", "SED") - clusters.list <- mapply(clustercreate, all.tsne.coords, names(all.tsne.coords)) - - FindCommonClusters <- function(clusters.list, keeplength=3) { # >>>> NEW method - # For each distance metric's clusters: - # 1. Create a square matrix of all PTMs. - # 2. For each cluster, set all PTM–PTM pairs in the cluster to 1 (indicating co-membership). - # 3. The final matrix for a method has 1 for PTM pairs co-clustered in that method; 0 otherwise. - - co_membership_matrix <- function(clusters, all_ptms) { - # 1. square matrix of all PTMs - mat <- matrix(0, nrow = length(PTMnames), ncol = length(PTMnames), - dimnames = list(PTMnames, PTMnames)) - # 2. For each cluster, set all PTM–PTM pairs in the cluster to 1 (indicating co-membership). - for (cluster in clusters) { - ptms <- cluster$PTMnames - if (length(ptms) > 1) { - mat[ptms, ptms] <- 1 - } - } - diag(mat) <- 0 - return(mat) +#' Example_Output <- MakeClusterList(ex_tiny_ptm_table) # Run function +#' Example_Output[[1]][1:3] # Display data +#' # Do we want to have one for adj.consensus? Doesn't seem like it'd be very helpful to view. +#' utils::head(Example_Output[[3]][, c(1, 2, 3, 4, 5)]) # Display data +MakeClusterList <- function(ptmtable, keeplength = 2, toolong = 3.5, tsne_perplexity = 15, tsne_max_iter = 5000) { + start_time <- Sys.time() + message("Starting correlation calculations and t-SNE.") + message("This may take a few minutes or hours for large data sets.") + + # Add if statement here to make sure functions are formatted correctly # + # Ensure ptmtable is a data frame with numeric values # + PTMnames <- rownames(ptmtable) + if (!all(vapply(ptmtable, is.numeric, FUN.VALUE = logical(1)))) { + stop("All columns in 'ptmtable' must be numeric.") } - adjacency_matrices <- purrr::map(clusters.list, co_membership_matrix, all_ptms=all_ptms) - - # Step 2: Sum the Co-Membership Matrices Across Methods - adj.sum <- Reduce("+", adjacency_matrices) # values: 0 (never), 1, 2, 3 (co-clustered in all 3 methods) - - #Step 3: Build a Consensus Network (for Co-Clustering in All 3 Embeddings) - # Edges: Only keep edges where `adj.sum == 3` (meaning the PTM pair is in the same cluster in all three methods). - adj.consensus.matrix <- (adj.sum == 3) * 1 # This is a neat R trick to convert a logical matrix (`TRUE`/`FALSE`) to a numeric matrix (`1`/`0`) via multiplication. - g <- igraph::graph_from_adjacency_matrix(adj.consensus.matrix, mode="undirected", diag=FALSE) - - # Step 4: Extract Cliques (Consensus Clusters): since clusters partitioned the graph, every clique is a connected component - components <- igraph::components(g) - # filter out components smaller than keeplength parameter and create list of co-clusters - # from existing single-metric clusters - keep_ids <- which(components$csize >= keeplength) - clusters_in_all_three <- lapply(keep_ids, - function(id) {igraph::V(g)$name[components$membership == id]} - ) - if (length(clusters_in_all_three) > 0) { - names(clusters_in_all_three) <- paste0("ConsensusCluster", seq_along(clusters_in_all_three)) - } else { - warning("No cliques remain after filtering with keeplength = ", keeplength) - } + # SPEARMAN CALCULATION + + # Calculate Spearman correlation # + ptm.correlation.matrix <- stats::cor(t(ptmtable), use = "pairwise.complete.obs", method = "spearman") + # Note: this is the slowest step. We found use = "pairwise.complete.obs", + # method = "spearman" to perform the best according to evaluations with data + # with missing values, but it takes longer. + + # Replace diagonal with NA # + diag(ptm.correlation.matrix) <- NA + + # Calculate dissimilarity # + sp.diss.matrix <- 1 - abs(ptm.correlation.matrix) + + # Handle any remaining NA values by setting them to the maximum dissimilarity # + max.dist.sp <- max(sp.diss.matrix, na.rm = TRUE) + + # Set NA values to 100 * the max distance + sp.diss.matrix[is.na(sp.diss.matrix)] <- 100 * max.dist.sp + message("Spearman correlation calculation complete after ", round(Sys.time() - start_time, 2), " ", units(Sys.time() - start_time), " total.") + + # Run t-SNE # + tsne.results <- Rtsne::Rtsne(sp.diss.matrix, dims = 3, perplexity = tsne_perplexity, theta = 0.25, max_iter = tsne_max_iter, check_duplicates = FALSE, pca = FALSE) + spearman.coords <- tsne.results$Y + message("Spearman t-SNE calculation complete after ", round(Sys.time() - start_time, 2), " ", units(Sys.time() - start_time), " total.") + + # EUCLIDEAN CALCULATION + + # Convert the dataframe to a distance matrix using Euclidean distance # + ptmtable.dist <- as.matrix(stats::dist(ptmtable, method = "euclidean")) + + # Compute the maximum distance in the matrix, excluding NA values # + max.dist <- max(ptmtable.dist, na.rm = TRUE) + + # Replace NA values in the distance matrix with 100 times the maximum distance # + ptmtable.dist[is.na(ptmtable.dist)] <- 100 * max.dist + + # Normalize the distance matrix by scaling it to a range from 0 to 100. This becomes the distance matrix for euclidian distance which we will run Rtsne on# + eu.dist.calc <- 100 * ptmtable.dist / max(ptmtable.dist, na.rm = TRUE) + eu.dist.calc <- as.matrix(eu.dist.calc) # Fix eu.dist.calc RQ + message("Euclidean distance calculation complete after ", round(Sys.time() - start_time, 2), " ", units(Sys.time() - start_time), " total.") + + # Run t-SNE # + eu.ptms.tsne.list <- Rtsne::Rtsne(eu.dist.calc, dims = 3, perplexity = tsne_perplexity, theta = 0.25, max_iter = tsne_max_iter, check_duplicates = FALSE, pca = FALSE) + euclidean.coords <- eu.ptms.tsne.list$Y + message("Euclidean t-SNE calculation complete after ", round(Sys.time() - start_time, 2), " ", units(Sys.time() - start_time), " total.") + + # COMBINED CALCULATION + + # fix spearman thing; so do the exact same thing but no absolute value + sp.diss.calc <- 1 - ptm.correlation.matrix # range goes from (-1 to 1) to (0 to 2) + max.diss.sp <- max(sp.diss.calc, na.rm = TRUE) # find the max value (around 2) + sp.diss.calc[is.na(sp.diss.calc)] <- 50 * max.diss.sp # Set every NA value to 50 * the max dissimilarity + sp.diss.calc <- as.matrix(sp.diss.calc) # turn into a matrix + + # find average + combined.distance <- (sp.diss.calc + eu.dist.calc) / 2 + message("Combined distance calculation complete after ", round(Sys.time() - start_time, 2), " ", units(Sys.time() - start_time), " total.") + # Perform t-SNE on the combined distances # + tsne.result <- Rtsne::Rtsne(combined.distance, dims = 3, perplexity = tsne_perplexity, theta = 0.25, max_iter = tsne_max_iter, check_duplicates = FALSE, pca = FALSE) + sed.coords <- tsne.result$Y + message("SED t-SNE calculation complete after ", round(Sys.time() - start_time, 2), " ", units(Sys.time() - start_time), " total.") + + + # Nested function to analyze result + group <- NULL # Gets rid of check note + clustercreate <- function(tsne_coords, distance_name) { + # Find clusters that are connected + # this returns a vector with cluster assignment for each PTM + result.disc2 <- vegan::distconnected(stats::dist(tsne_coords), toolong = toolong, trace = FALSE) + + # Create a plot of the clusters, with convex hulll around each cluster + # note that tsne_coords is a matrix with 3 columns (the t-SNE dimensions), but only first two will get plotted + vegan::ordiplot(tsne_coords, display = c(), main = paste("t-SNE Cluster Plot", distance_name)) + vegan::ordihull(tsne_coords, result.disc2, col = "red", lwd = 2) + + # make a dataframe with one column for PTMs and another for group assignment + result.span.df <- data.frame(PTMnames = PTMnames, group = result.disc2) + + # Convert data frame into a list of clusters (check doesn't like group but it's a column name) + result.span.list <- plyr::dlply(result.span.df, plyr::.(group)) # GROUP LIST ! + + end_time <- Sys.time() + total_time <- end_time - start_time + message("Clustering for ", distance_name, " complete after ", round(total_time, 2), " ", units(total_time), " total.") + return(result.span.list) + } # END of nested function + + # Create all tsne coords data sturcture and give it names + all.tsne.coords <- list(euclidean.coords, spearman.coords, sed.coords) + names(all.tsne.coords) <- c("Euclidean", "Spearman", "SED") + clusters.list <- mapply(clustercreate, all.tsne.coords, names(all.tsne.coords)) + + FindCommonClusters <- function(clusters.list, keeplength = 3) { # >>>> NEW method + # For each distance metric's clusters: + # 1. Create a square matrix of all PTMs. + # 2. For each cluster, set all PTM–PTM pairs in the cluster to 1 (indicating co-membership). + # 3. The final matrix for a method has 1 for PTM pairs co-clustered in that method; 0 otherwise. + + co_membership_matrix <- function(clusters, all_ptms) { + # 1. square matrix of all PTMs + mat <- matrix(0, + nrow = length(PTMnames), ncol = length(PTMnames), + dimnames = list(PTMnames, PTMnames) + ) + # 2. For each cluster, set all PTM–PTM pairs in the cluster to 1 (indicating co-membership). + for (cluster in clusters) { + ptms <- cluster$PTMnames + if (length(ptms) > 1) { + mat[ptms, ptms] <- 1 + } + } + diag(mat) <- 0 + return(mat) + } - end_time <- Sys.time() - message("Consensus clustering complete after ", round(end_time - start_time, 2), " ", units(end_time - start_time), " total.") - return(list(adj.consensus.matrix, clusters_in_all_three)) - } + adjacency_matrices <- purrr::map(clusters.list, co_membership_matrix, all_ptms = all_ptms) + + # Step 2: Sum the Co-Membership Matrices Across Methods + adj.sum <- Reduce("+", adjacency_matrices) # values: 0 (never), 1, 2, 3 (co-clustered in all 3 methods) + + # Step 3: Build a Consensus Network (for Co-Clustering in All 3 Embeddings) + # Edges: Only keep edges where `adj.sum == 3` (meaning the PTM pair is in the same cluster in all three methods). + adj.consensus.matrix <- (adj.sum == 3) * 1 # This is a neat R trick to convert a logical matrix (`TRUE`/`FALSE`) to a numeric matrix (`1`/`0`) via multiplication. + g <- igraph::graph_from_adjacency_matrix(adj.consensus.matrix, mode = "undirected", diag = FALSE) + + # Step 4: Extract Cliques (Consensus Clusters): since clusters partitioned the graph, every clique is a connected component + components <- igraph::components(g) + # filter out components smaller than keeplength parameter and create list of co-clusters + # from existing single-metric clusters + keep_ids <- which(components$csize >= keeplength) + clusters_in_all_three <- lapply( + keep_ids, + function(id) { + igraph::V(g)$name[components$membership == id] + } + ) + if (length(clusters_in_all_three) > 0) { + names(clusters_in_all_three) <- paste0("ConsensusCluster", seq_along(clusters_in_all_three)) + } else { + warning("No cliques remain after filtering with keeplength = ", keeplength) + } + + end_time <- Sys.time() + message("Consensus clustering complete after ", round(end_time - start_time, 2), " ", units(end_time - start_time), " total.") + return(list(adj.consensus.matrix, clusters_in_all_three)) + } - # Find common clusters from existing single-metric clusters - clusters.common.list <- FindCommonClusters(clusters.list, keeplength) - adj.consensus.matrix <- clusters.common.list[[1]] - common.clusters <- clusters.common.list[[2]] + # Find common clusters from existing single-metric clusters + clusters.common.list <- FindCommonClusters(clusters.list, keeplength) + adj.consensus.matrix <- clusters.common.list[[1]] + common.clusters <- clusters.common.list[[2]] - message("MakeClusterList complete after ", round(Sys.time() - start_time, 2), " ", units(Sys.time() - start_time), " total.") + message("MakeClusterList complete after ", round(Sys.time() - start_time, 2), " ", units(Sys.time() - start_time), " total.") - return(list(common.clusters, adj.consensus.matrix, ptm.correlation.matrix)) + return(list(common.clusters, adj.consensus.matrix, ptm.correlation.matrix)) } diff --git a/R/MakeCorrelationNetwork.R b/R/MakeCorrelationNetwork.R index b9fb93fd..8760f33e 100644 --- a/R/MakeCorrelationNetwork.R +++ b/R/MakeCorrelationNetwork.R @@ -19,134 +19,136 @@ #' Example_Output[[1]][1:5, 1:5] #' Example_Output[[2]][1:5, 1:5] #' -MakeCorrelationNetwork <- function(adj.consensus.matrix, ptm.correlation.matrix){ - # Two nested functions for creating the PTM and gene CCCN, respectively - # Use the consensus adjacency matrix to filter PTM correlations, then create a graph and edge files for PTMs and genes - - start_time <- Sys.time() - - MakePTMCCCN <- function(adj.consensus.matrix, ptm.correlation.matrix) { - message("Making PTM CCCN") - # Use only PTM pairs that co-clustered in all three methods (adj.consensus.matrix == 1) - ptm.cccn <- ptm.correlation.matrix[sort(rownames(ptm.correlation.matrix)), sort(colnames(ptm.correlation.matrix))] - ptm.cccn.mask <- adj.consensus.matrix[sort(rownames(adj.consensus.matrix)), sort(colnames(adj.consensus.matrix))] - # Find active ptms before setting 0 to NA - active_flags <- (rowSums(adj.consensus.matrix) > 0) | (colSums(adj.consensus.matrix) > 0) - active_ptms <- rownames(adj.consensus.matrix)[active_flags] - ptm.cccn.mask[ptm.cccn.mask == 0] <- NA # Set 0 to NA for masking - - # Apply mask to Spearman correlations - ptm.cccn[is.na(ptm.cccn.mask)] <- NA - - # Subset the adjacency matrix to include only active PTMs - ptm.cccn.active <- ptm.cccn[active_ptms, active_ptms, drop=FALSE] - - # Remove self-loops - if(any(!is.na(diag(ptm.cccn)))) {diag(ptm.cccn) <- NA} - ptm.cccn <- ptm.cccn.active - - # Make edglist file: - # Replace NA with 0 for igraph compatibility - ptm.cccn0 <- ptm.cccn - ptm.cccn0[is.na(ptm.cccn0)] <- 0 - - # Create igraph object from correlation matrix - ptm.cccn.g <- igraph::graph_from_adjacency_matrix(ptm.cccn0, mode = "lower", diag = FALSE, weighted = TRUE) - - # Extract edge list as data.frame - ptm.cccn.edges <- as.data.frame(igraph::as_edgelist(ptm.cccn.g)) - names(ptm.cccn.edges) <- c("source", "target") - ptm.cccn.edges$Weight <- igraph::edge_attr(ptm.cccn.g, "weight") - - # Classify interaction by correlation value - ptm.cccn.edges$interaction <- "correlation" - ptm.cccn.edges$interaction[ptm.cccn.edges$Weight <= -0.5] <- "negative correlation" - ptm.cccn.edges$interaction[ptm.cccn.edges$Weight >= 0.5] <- "positive correlation" - - message("PTM CCCN complete after ", round(Sys.time() - start_time, 2), " ", units(Sys.time() - start_time), " total.") - return(list(ptm.cccn0, ptm.cccn.g, ptm.cccn.edges)) - } - # - ptm.cccn.list <- MakePTMCCCN(adj.consensus.matrix, ptm.correlation.matrix) - ptm.cccn <- ptm.cccn.list[[1]] - ptm.cccn.g <- ptm.cccn.list[[2]] - ptm.cccn.edges <- ptm.cccn.list[[3]] - # Build the Gene CCCN - # Start from igraph object saved from MakePTMCCCN() - # Double ddply Summing: By grouping and summing in both directions, you ensure the aggregation is performed for both genes in each pair, producing a correctly shaped and labeled gene–gene matrix. - MakeGeneCCCN <- function(ptm.cccn){ - message("Making Gene CCCN") - # ptm.cccn was returned above - gene.cccn <- data.frame(ptm.cccn, row.names = rownames(ptm.cccn), check.rows=TRUE, check.names=FALSE, fix.empty.names = FALSE) - # Check: identical(rownames(gene.cccn), colnames(gene.cccn)) # TRUE - gene.cccn$Gene.Name <- vapply(rownames(gene.cccn), function (x) unlist(strsplit(x, " ", fixed=TRUE))[1], FUN.VALUE = character(1)) - # Use only upper triangle so correlations are not duplicated during the next step - gene.cccn[lower.tri(gene.cccn)] <- NA - - # Sum correlations in one dimension, then the other dimension - gene.cccn2 <- dplyr::summarise( - dplyr::group_by(gene.cccn, .data$Gene.Name), - dplyr::across( - tidyselect::where(is.numeric), - ~sum(.x, na.rm = TRUE) - ) - ) - gene.cccn2 <- as.data.frame(gene.cccn2) - rownames(gene.cccn2) <- gene.cccn2$Gene.Name - gene.cccn2 <- gene.cccn2[, 2:ncol(gene.cccn2)] - # Transform to do the other dimension - gene.cccn2 <- data.frame(t(gene.cccn2)) - gene.cccn2$Gene <- vapply(rownames(gene.cccn2), function (x) unlist(strsplit(x, " ", fixed=TRUE))[1], FUN.VALUE = character(1)) - # Now sum the other dimension - gene.cccn3 <- dplyr::summarise( - dplyr::group_by(gene.cccn2, .data$Gene), - dplyr::across( - tidyselect::where(is.numeric), - ~sum(.x, na.rm = TRUE) - ) - ) - # R likes to put dots in column names, which is a problem for ambiguous gene names and gene names with hyphens - # so just work around the problem (once satisfied that the gene names actually match). - names(gene.cccn3)[2:ncol(gene.cccn3)] <- gene.cccn3$Gene - gene.cccn3 <- as.data.frame(gene.cccn3) - rownames(gene.cccn3) <- gene.cccn3$Gene - gene.cccn.matrix <- as.matrix(gene.cccn3[,2:ncol(gene.cccn3)]) - # Replace 0 with NA in the correlation matrix - gene.cccn.matrix[gene.cccn.matrix==0] <- NA - # Remove self-loops by setting diagonal to NA; diag() requires a matrix - diag(gene.cccn.matrix) <- NA - gene.cccn0 <- gene.cccn.matrix - gene.cccn0[is.na(gene.cccn0)] <- 0 # igraph doesn't like NAs - # For Graphing and making edge lists - gene.cccn.g <- igraph::graph_from_adjacency_matrix(gene.cccn0, mode = "lower", diag = FALSE, weighted = TRUE) - # Extract edge list as data.frame - gene.cccn.edges <- as.data.frame(igraph::as_edgelist(gene.cccn.g)) - names(gene.cccn.edges) <- c("source", "target") - gene.cccn.edges$Weight <- igraph::edge_attr(gene.cccn.g, "weight") - - # Classify interaction by correlation value - gene.cccn.edges$interaction <- "correlation" - gene.cccn.edges$interaction[gene.cccn.edges$Weight <= -0.5] <- "negative correlation" - gene.cccn.edges$interaction[gene.cccn.edges$Weight >= 0.5] <- "positive correlation" - - message("Gene CCCN complete after ", round(Sys.time() - start_time, 2), " ", units(Sys.time() - start_time), " total.") - return(list(gene.cccn.g, gene.cccn.edges, gene.cccn.matrix)) - } - - gene.cccn.list <- MakeGeneCCCN(ptm.cccn) - gene.cccn.g <- gene.cccn.list[[1]] - gene.cccn.edges <- gene.cccn.list[[2]] - gene.cccn <- gene.cccn.list[[3]] - - # Make a list of nodes for gathering PPI data - gene.cccn.nodes <- unique(c(gene.cccn.edges$source, gene.cccn.edges$target)) - if(length(gene.cccn.nodes) == 0) stop("No genes found (gene.cccn.nodes is empty)") - - ### Return Final Data Structure ### - return (list(ptm.cccn.edges, gene.cccn.edges, gene.cccn.nodes)) - - ### Graphing ### - # graph <- igraph::graph_from_adjacency_matrix(gene.cccn, mode = "lower", diag = FALSE, weighted = "Weight") - # plot(graph) - # Note: Cytoscape is recommended for graphing networks +MakeCorrelationNetwork <- function(adj.consensus.matrix, ptm.correlation.matrix) { + # Two nested functions for creating the PTM and gene CCCN, respectively + # Use the consensus adjacency matrix to filter PTM correlations, then create a graph and edge files for PTMs and genes + + start_time <- Sys.time() + + MakePTMCCCN <- function(adj.consensus.matrix, ptm.correlation.matrix) { + message("Making PTM CCCN") + # Use only PTM pairs that co-clustered in all three methods (adj.consensus.matrix == 1) + ptm.cccn <- ptm.correlation.matrix[sort(rownames(ptm.correlation.matrix)), sort(colnames(ptm.correlation.matrix))] + ptm.cccn.mask <- adj.consensus.matrix[sort(rownames(adj.consensus.matrix)), sort(colnames(adj.consensus.matrix))] + # Find active ptms before setting 0 to NA + active_flags <- (rowSums(adj.consensus.matrix) > 0) | (colSums(adj.consensus.matrix) > 0) + active_ptms <- rownames(adj.consensus.matrix)[active_flags] + ptm.cccn.mask[ptm.cccn.mask == 0] <- NA # Set 0 to NA for masking + + # Apply mask to Spearman correlations + ptm.cccn[is.na(ptm.cccn.mask)] <- NA + + # Subset the adjacency matrix to include only active PTMs + ptm.cccn.active <- ptm.cccn[active_ptms, active_ptms, drop = FALSE] + + # Remove self-loops + if (any(!is.na(diag(ptm.cccn)))) { + diag(ptm.cccn) <- NA + } + ptm.cccn <- ptm.cccn.active + + # Make edglist file: + # Replace NA with 0 for igraph compatibility + ptm.cccn0 <- ptm.cccn + ptm.cccn0[is.na(ptm.cccn0)] <- 0 + + # Create igraph object from correlation matrix + ptm.cccn.g <- igraph::graph_from_adjacency_matrix(ptm.cccn0, mode = "lower", diag = FALSE, weighted = TRUE) + + # Extract edge list as data.frame + ptm.cccn.edges <- as.data.frame(igraph::as_edgelist(ptm.cccn.g)) + names(ptm.cccn.edges) <- c("source", "target") + ptm.cccn.edges$Weight <- igraph::edge_attr(ptm.cccn.g, "weight") + + # Classify interaction by correlation value + ptm.cccn.edges$interaction <- "correlation" + ptm.cccn.edges$interaction[ptm.cccn.edges$Weight <= -0.5] <- "negative correlation" + ptm.cccn.edges$interaction[ptm.cccn.edges$Weight >= 0.5] <- "positive correlation" + + message("PTM CCCN complete after ", round(Sys.time() - start_time, 2), " ", units(Sys.time() - start_time), " total.") + return(list(ptm.cccn0, ptm.cccn.g, ptm.cccn.edges)) + } + # + ptm.cccn.list <- MakePTMCCCN(adj.consensus.matrix, ptm.correlation.matrix) + ptm.cccn <- ptm.cccn.list[[1]] + ptm.cccn.g <- ptm.cccn.list[[2]] + ptm.cccn.edges <- ptm.cccn.list[[3]] + # Build the Gene CCCN + # Start from igraph object saved from MakePTMCCCN() + # Double ddply Summing: By grouping and summing in both directions, you ensure the aggregation is performed for both genes in each pair, producing a correctly shaped and labeled gene–gene matrix. + MakeGeneCCCN <- function(ptm.cccn) { + message("Making Gene CCCN") + # ptm.cccn was returned above + gene.cccn <- data.frame(ptm.cccn, row.names = rownames(ptm.cccn), check.rows = TRUE, check.names = FALSE, fix.empty.names = FALSE) + # Check: identical(rownames(gene.cccn), colnames(gene.cccn)) # TRUE + gene.cccn$Gene.Name <- vapply(rownames(gene.cccn), function(x) unlist(strsplit(x, " ", fixed = TRUE))[1], FUN.VALUE = character(1)) + # Use only upper triangle so correlations are not duplicated during the next step + gene.cccn[lower.tri(gene.cccn)] <- NA + + # Sum correlations in one dimension, then the other dimension + gene.cccn2 <- dplyr::summarise( + dplyr::group_by(gene.cccn, .data$Gene.Name), + dplyr::across( + tidyselect::where(is.numeric), + ~ sum(.x, na.rm = TRUE) + ) + ) + gene.cccn2 <- as.data.frame(gene.cccn2) + rownames(gene.cccn2) <- gene.cccn2$Gene.Name + gene.cccn2 <- gene.cccn2[, 2:ncol(gene.cccn2)] + # Transform to do the other dimension + gene.cccn2 <- data.frame(t(gene.cccn2)) + gene.cccn2$Gene <- vapply(rownames(gene.cccn2), function(x) unlist(strsplit(x, " ", fixed = TRUE))[1], FUN.VALUE = character(1)) + # Now sum the other dimension + gene.cccn3 <- dplyr::summarise( + dplyr::group_by(gene.cccn2, .data$Gene), + dplyr::across( + tidyselect::where(is.numeric), + ~ sum(.x, na.rm = TRUE) + ) + ) + # R likes to put dots in column names, which is a problem for ambiguous gene names and gene names with hyphens + # so just work around the problem (once satisfied that the gene names actually match). + names(gene.cccn3)[2:ncol(gene.cccn3)] <- gene.cccn3$Gene + gene.cccn3 <- as.data.frame(gene.cccn3) + rownames(gene.cccn3) <- gene.cccn3$Gene + gene.cccn.matrix <- as.matrix(gene.cccn3[, 2:ncol(gene.cccn3)]) + # Replace 0 with NA in the correlation matrix + gene.cccn.matrix[gene.cccn.matrix == 0] <- NA + # Remove self-loops by setting diagonal to NA; diag() requires a matrix + diag(gene.cccn.matrix) <- NA + gene.cccn0 <- gene.cccn.matrix + gene.cccn0[is.na(gene.cccn0)] <- 0 # igraph doesn't like NAs + # For Graphing and making edge lists + gene.cccn.g <- igraph::graph_from_adjacency_matrix(gene.cccn0, mode = "lower", diag = FALSE, weighted = TRUE) + # Extract edge list as data.frame + gene.cccn.edges <- as.data.frame(igraph::as_edgelist(gene.cccn.g)) + names(gene.cccn.edges) <- c("source", "target") + gene.cccn.edges$Weight <- igraph::edge_attr(gene.cccn.g, "weight") + + # Classify interaction by correlation value + gene.cccn.edges$interaction <- "correlation" + gene.cccn.edges$interaction[gene.cccn.edges$Weight <= -0.5] <- "negative correlation" + gene.cccn.edges$interaction[gene.cccn.edges$Weight >= 0.5] <- "positive correlation" + + message("Gene CCCN complete after ", round(Sys.time() - start_time, 2), " ", units(Sys.time() - start_time), " total.") + return(list(gene.cccn.g, gene.cccn.edges, gene.cccn.matrix)) + } + + gene.cccn.list <- MakeGeneCCCN(ptm.cccn) + gene.cccn.g <- gene.cccn.list[[1]] + gene.cccn.edges <- gene.cccn.list[[2]] + gene.cccn <- gene.cccn.list[[3]] + + # Make a list of nodes for gathering PPI data + gene.cccn.nodes <- unique(c(gene.cccn.edges$source, gene.cccn.edges$target)) + if (length(gene.cccn.nodes) == 0) stop("No genes found (gene.cccn.nodes is empty)") + + ### Return Final Data Structure ### + return(list(ptm.cccn.edges, gene.cccn.edges, gene.cccn.nodes)) + + ### Graphing ### + # graph <- igraph::graph_from_adjacency_matrix(gene.cccn, mode = "lower", diag = FALSE, weighted = "Weight") + # plot(graph) + # Note: Cytoscape is recommended for graphing networks } diff --git a/R/data_documentation.R b/R/data_documentation.R index 66b011cb..9b881aef 100644 --- a/R/data_documentation.R +++ b/R/data_documentation.R @@ -123,11 +123,11 @@ #' PTM CCCN Edgelist -#' +#' #' PTM Co-Cluster Correlation Network as an Edgelist -#' +#' #' @keywords internal -#' +#' #' @format Data Frame of 5710 rows and 4 columns #' \describe{ #' \item{source}{A PTM} @@ -141,11 +141,11 @@ "ex_ptm_cccn_edges" #' Gene CCCN Edgelist -#' +#' #' Gene Co-Cluster Correlation Network as an Edgelist -#' +#' #' @keywords internal -#' +#' #' @format Data Frame of 2279 rows and 4 columns #' \describe{ #' \item{source}{A Gene} @@ -159,11 +159,11 @@ "ex_gene_cccn_edges" #' Gene list (nodes) -#' +#' #' List of Genes gotten from PTMs -#' +#' #' @keywords internal -#' +#' #' @format Character Vector of length 390 #' \describe{ #' \item{String}{Gene} @@ -196,11 +196,11 @@ #' Genemania Edges -#' +#' #' Edges from gene_cccn_nodes that are also found within the Genemania database -#' +#' #' @keywords internal -#' +#' #' @format Dataframe of 5 rows and 4 columns #' \describe{ #' \item{source}{First gene interacting} @@ -213,16 +213,15 @@ #' head(ex_genemania_edges) "ex_genemania_edges" -# DO KINSUB EDGES - +# DO KINSUB EDGES #' Combined PPIs -#' +#' #' Combines all Protein to Protein interaction data and combines their weights into one score from a range of 0 to 1 -#' +#' #' @keywords internal -#' +#' #' @format Dataframe of 8849 rows and 4 columns #' \describe{ #' \item{source}{First gene interacting} @@ -236,11 +235,11 @@ "ex_combined_ppi" #' Cfn -#' +#' #' The Cluster Filtered Network is a filtered version of the Combined PPIs but only for proteins found to interact via t-SNE clusters -#' +#' #' @keywords internal -#' +#' #' @format Dataframe of 196 rows and 4 columns #' \describe{ #' \item{source}{First gene interacting} @@ -256,11 +255,11 @@ #' Bioplanet -#' +#' #' A database file from the bioplanet website -#' +#' #' @keywords internal -#' +#' #' @format Dataframe of 196 rows and 4 columns #' \describe{ #' \item{PATHWAY_ID}{The pathway a gene belongs to} @@ -274,11 +273,11 @@ "ex_bioplanet" #' Pathway Crosstalk Network -#' +#' #' Specific genes from the bioplanet file that we found in Common Clusters and their respective cell signaling pathway relationships -#' +#' #' @keywords internal -#' +#' #' @format Dataframe of 37 rows and 4 columns #' \describe{ #' \item{source}{The name of a pathway interacting} @@ -292,11 +291,11 @@ "ex_pathway_crosstalk_network" #' PCN Edge List -#' +#' #' The Pathway Crosstalk Network but with direct values for both weight types -#' +#' #' @keywords internal -#' +#' #' @format Dataframe of 19 rows and 4 columns #' \describe{ #' \item{source}{The name of a pathway interacting} @@ -310,11 +309,11 @@ "ex_PCNedgelist" #' Pathways list -#' +#' #' A list of pathway names and the genes that make them up -#' +#' #' @keywords internal -#' +#' #' @format list of 9 character vectors #' \describe{ #' \item{Pathway name}{Genes in the pathway} @@ -325,18 +324,18 @@ "ex_pathways_list" #' BRCA Cluster List Data -#' +#' #' Output of MakeClusterList on the BRCA data -#' +#' #' @keywords internal -#' +#' #' @format list containing the 3 outputs of MakeClusterList #' \describe{ #' \item{Common Clusters}{231 lists of PTMs found to cluster, see Common Clusters} #' \item{Adjacency Matrix}{4237 x 4237, see Adjacency Consensus Matrix} #' \item{PTM correlation matrix}{4237 x 4237, see Correlation Matrix} -#' } -#' +#' } +#' #' @source "MakeClusterList on the BRCA data" #' @examples #' head("brca_clusterlist_data[[1]]") @@ -345,29 +344,29 @@ "brca_clusterlist_data" # BRCA CCCN Data -#' +#' #' Output of MakeCorrelationNetwork on the BRCA data -#' +#' #' @keywords internal "brca_CCCN_data" #' BRCA STRINGdb Edges -#' +#' #' Edges from the STRINGdb database for the genes in the BRCA dataset -#' +#' #' @keywords internal "BRCA_stringdb.edges" #' BRCA Genemania Edges -#' +#' #' Edges from the Genemania database for the genes in the BRCA dataset -#' +#' #' @keywords internal "BRCA_genemania.edges" #' BRCA PCN Data -#' +#' #' Output of BuildPathwayCrosstalkNetwork on the BRCA data -#' +#' #' @keywords internal -"BRCA_PCN.data" \ No newline at end of file +"BRCA_PCN.data" diff --git a/man/filter.edges.0.Rd b/man/filter.edges.0.Rd index 4f961d03..deeb68ba 100644 --- a/man/filter.edges.0.Rd +++ b/man/filter.edges.0.Rd @@ -7,6 +7,10 @@ filter.edges.0(nodenames, edge.file) } \arguments{ +\item{nodenames}{} + +\item{edge.file}{} + \item{cfn}{A version of ppi.network with only the edges that exist in cccn.matrix and have non-zero weights} \item{ptmtable}{A dataset for post-translational modifications. Formatted with numbered rows, and the first column containing PTM names. The rest of the column names should be drugs. Values are numeric values that represent how much the PTM has reacted to the drug.} diff --git a/man/make.genepep.edges.Rd b/man/make.genepep.edges.Rd index f0bb0c3b..b2df0f20 100644 --- a/man/make.genepep.edges.Rd +++ b/man/make.genepep.edges.Rd @@ -6,6 +6,9 @@ \usage{ make.genepep.edges(peptide.edgefile) } +\arguments{ +\item{peptide.edgefile}{} +} \description{ Title } diff --git a/man/remove.autophos.Rd b/man/remove.autophos.Rd index 7d8c3550..d2980675 100644 --- a/man/remove.autophos.Rd +++ b/man/remove.autophos.Rd @@ -6,6 +6,9 @@ \usage{ remove.autophos(edgefile) } +\arguments{ +\item{edgefile}{} +} \description{ Title } diff --git a/vignettes/CreatingNetworks.Rmd b/vignettes/CreatingNetworks.Rmd index 5e333101..40748bb6 100644 --- a/vignettes/CreatingNetworks.Rmd +++ b/vignettes/CreatingNetworks.Rmd @@ -259,11 +259,7 @@ also be downloaded and added to the PPI network. All three external databases have different interfaces for downloading data, so we show how to retrieve data from each of them below. -For this tutorial, we query STRINGdb and GeneMANIA directly. We also -provide the option (using the switch local = TRUE) to retrive edges from -human PPIs from these sources in pre-assembled files; the [BRCANetworks -vignette](BCRANetworks.html) demonstrates how to get the STRING-db and -GeneMANIA edges from the static downloaded networks. +For this tutorial, we query STRINGdb and GeneMANIA directly. We also provide the option (using the switch local = TRUE) to retrive edges from human PPIs from these sources in pre-assembled files; the [BRCANetworks vignette](BCRANetworks.html) demonstrates how to get the STRING-db and GeneMANIA edges from the static downloaded networks. ### 1. STRINGdb @@ -392,20 +388,24 @@ cfn.merged <- mergeEdges(cfn) # Step 5: Pathway Crosstalk Network -The final step is the creation of the Pathway Crosstalk Network (PCN), which -creates a set of pathway-pathway edges that have two weights: a Jaccard similarity -and a Cluster-Pathway Evidence score. +The final step is the creation of the Pathway Crosstalk Network (PCN). This step requires input of an external database from [NCATS BioPlanet](https://tripod.nih.gov/bioplanet/download/pathway.csv) that contains groups of genes (proteins) involved in various cellular -processes known as pathways. P2P provides a function [`ReadBioplanetFile`](references/ReadBioplanetFile.html) -that reads in the BioPlanet file and converts it into a list of pathways. The first pathways is displayed below: +processes known as pathways. `BuildPathwayCrosstalkNetwork` turns this +data file into a list of pathways and converts those pathways into a +list of pathway-pathway edges, each of which is assigned a Jaccard +similarity and a Cluster-Pathway Evidence score based on the common +clusters found in the gene co-cluster correlation network. Info about +the Cluster-Pathway Evidence score can be found +[here](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1010690) +For graphing in Cytoscape, the Cluster-Pathway Evidence and Jaccard +similarity edges are listed separately in the edgelist called +pathway.crosstalk.network. ```{r eval = TRUE, echo = TRUE} bioplanet.file <- system.file("extdata", "pathway.csv", package = "PTMsToPathways") -pathways.list <- ReadBioplanetFile(bioplanet.file) -pathways.list[[5]] ``` The function [`BuildPathwayCrosstalkNetwork`](references/BuildPathwayCrosstalkNetwork.html) @@ -442,7 +442,7 @@ And we can see some of the pathway crosstalk network edges below: pathway.crosstalk.network[1:5,] ``` -```{r eval = TRUE, echo = FALSE} +```{r eval = TRUE, echo = TRUE} dat <- pathway.crosstalk.network[1:5,] knitr::kable(dat, align = 'l', digits = 2) ``` diff --git a/vignettes/CytoscapeGraphing.Rmd b/vignettes/CytoscapeGraphing.Rmd index 75a379fc..0f07a5b8 100644 --- a/vignettes/CytoscapeGraphing.Rmd +++ b/vignettes/CytoscapeGraphing.Rmd @@ -40,6 +40,10 @@ knitr::opts_chunk$set( ``` +## + +## Graphing Networks in Cytoscape + ### Top-down approach starting with pathways It is possible to graph the entire PCN, CFN, and CCCNs in their entirety, though @@ -53,30 +57,26 @@ PTMs of interest in the R data objects, then make smaller, more interpretable graphs in Cytoscape using RCy3. For example, we will find names of all pathways in Bioplanet that contain EGFR. -To get started, load the PTMsToPathways package and read the Bioplanet file using the -built-in P2P function [ReadBioplanetFile()](reference/ReadBioplanetFile.html). +The data object `pathways.list` is a list, where the name of the list element is +the name of a Bioplanet pathway and each element is a character vector of the +genes in that pathway. Then we want to find interactions between the pathway +"Transmembrane transport of small molecules" and those pathways. The utility +functions used below are in CytoscapeGraphing.R ```{r eval = TRUE, echo = FALSE} library(PTMsToPathways) -pathways.list <- ReadBioplanetFile(system.file("extdata", "bioplanet_pathway_June2025.csv", package = "PTMsToPathways")) +funckey <- function_key ``` -The data object `pathways.list` is a list, where the name of the list element is -the name of a Bioplanet pathway and each element is a character vector of the -genes in that pathway. - ```{r eval = TRUE} -egfr_pathways <- names(pathways.list)[sapply(1:length(pathways.list), function(x) - {"EGFR" %in% pathways.list[[x]]})] +egfr_pathways <- names(ex_pathways_list)[sapply(1:length(ex_pathways_list), function(x) + {"EGFR" %in% ex_pathways_list[[x]]})] ``` We expect 83 pathways that contain EGFR, so let's check: ```{r eval = TRUE} -head(egfr_pathways) length(egfr_pathways) ``` -Then we want to find interactions between the pathway -"Transmembrane transport of small molecules" and those pathways. ```{r eval = TRUE} egfr_transporter.pcn <- filter.edges.between( @@ -160,7 +160,7 @@ sp1 <- connectNodes.all(c("FYN", 'MET'), ig.graph=NULL, # To include co-clustered PTMs in the network an extra step is necessary: sp1_plus <- get.co.clustered.ptms(sp1) -sp1_plus.nodes <- make.cytoscape.node.file(sp1_plus, function_key, ptmtable, +sp1_plus.nodes <- make.cytoscape.node.file(sp1_plus, funckey, ptmtable, include.gene.data = TRUE, include.coclustered.PTMs = TRUE)