Add cpe calc to vignette and fix up BuildPathwayCrosstalkNetwork function

R/CytoscapeGraphingFunctions.R

@@ -279,7 +279,7 @@ cytoscape.graph.PCN.pathways <- function(PCN = pathway.crosstalk.network, net.na
   # edgeColors <- c("#9966FF", col2hex("green"))
   edgeColors <- c(gplots::col2hex(ggplot2::alpha("darkorchid1",0.33333)), gplots::col2hex(ggplot2::alpha("tomato", 0.25)), gplots::col2hex("green"))
   # edgeColors <- c("#9966FF", col2hex("green"))
-  edgeTypes <- c("PTM_cluster_evidence", "Protein_cluster_evidence", "pathway_Jaccard_similarity")
+  edgeTypes <- c("PTM_cluster_weights", "Protein_cluster_evidence", "pathway_Jaccard_similarity")
   RCy3::setEdgeColorMapping( 'interaction', edgeTypes, edgeColors, 'd', default.color="#FFFFFF")
   setEdgeWidths(ffactor = -1.2, log=TRUE) # Finally works!
   style.name <- paste("PCN style", (1+no.windows), sep=" ")

R/PathwayCrosstalkNetwork.R

@@ -59,24 +59,23 @@ ReadBioplanetFile <- function(bioplanet.file = "bioplanet.csv") {
 #'
 #' @param common.clusters The list of common clusters between all three distance metrics (Euclidean, Spearman, and SED). Can be made in MakeCorrelationNetwork
 #' @param bioplanet.file Either the path to a delimited Bioplanet file, or a named list of pathways where each list element is a character vector of gene symbols.
-#' @param createfile The path of where to create the edgelist file. Defaults to the working directory, if FALSE is provided, a file will not be created.
-#' @param PCN.edgelist.name Name of the PCN edgelist file to be created
 #' @return A list with these data structures at the given index: \enumerate{
 #' \item{Contains pathway source-target columns, along with the interaction type.}
 #' \item{Contains pathway source-target columns, with edge weights of their jaccard similarity and their Pathway-Pathway Evidence score.}
 #' \item{All pathways in the bioplanet database as a named list containing string vectors. Each vector is a pathway with strings associated with the genes in that pathway.}
+#' \item{The cluster-by-pathway CPE matrix used to compute PTM cluster weights.}
 #' }
 #' @export
 #'
 #' @examples
-#' Example_Output <- BuildPathwayCrosstalkNetwork(ex_common_clusters, ex_pathways_list, createfile = FALSE)
+#' Example_Output <- BuildPathwayCrosstalkNetwork(ex_common_clusters, ex_pathways_list)
 #' Example_Output[[1]][[3,]]
 #' Example_Output[[3]][[1:3]]
-BuildPathwayCrosstalkNetwork <- function(common.clusters, bioplanet.file = "bioplanet.csv", createfile = getwd(), PCN.edgelist.name = "PCN_file"){
+#' Example_Output[[4]][1:3, 1:3]
+BuildPathwayCrosstalkNetwork <- function(common.clusters, bioplanet.file = "bioplanet.csv"){
   message("Making PCN")
   start_time <- Sys.time()
   message(start_time)
-  if(is.character(createfile) && !dir.exists(createfile)) stop(paste("Could not find directory", createfile)) #If createfile is a path but an incorrect one
 
   # Accept either a file path or pre-built pathways.list
   if (is.character(bioplanet.file)) {
@@ -178,9 +177,13 @@ BuildPathwayCrosstalkNetwork <- function(common.clusters, bioplanet.file = "biop
   }
 
   ### Generate PCN network ###
-  # Get a vector of all the PTP weights for every pair of pathways using the CPE weights to filter. For a PTP weight to be non-NA, the PTP weight will be the sum of all clusters both pathways have nonzero CPEs in.
+  # Get a matrix of all the PTM Cluster Weights for every pair of pathways using the CPE
+  # weights. For a pathway-pathway PTM Cluster Weight to be non-NA, both pathways must
+  # have at least one non-NA CPE weight for the same cluster. Then, the PTM Cluster Weight
+  # is the sum of the CPE weights for all clusters that both pathways have a modified protein in,
+  # which is just the clusters that both pathways have non-NA CPE weights for.
 
-  PTPscore <- apply(PCNedgelist[,1:2], 1, function(x) {
+  PTMCW <- apply(PCNedgelist[,1:2], 1, function(x) {
     rows <- rowSums(!is.na(CPE.matrix[, x])) == 2
     if (any(rows)) {
       sum(CPE.matrix[rows, x], na.rm = TRUE)
@@ -189,61 +192,46 @@ BuildPathwayCrosstalkNetwork <- function(common.clusters, bioplanet.file = "biop
     }})
 
 
-  PTPscore[PTPscore== 0] <- NA # Safety check: Turn all 0s created in above line into NAs
+  PTMCW[PTMCW== 0] <- NA # Safety check: Turn all 0s created in above line into NAs
 
-  PCNedgelist <- cbind(PCNedgelist, PTPscore) #Bind all the columns together. Now Data structure is PATHWAY | PATHWAY | Jaccard | CPE
+  PCNedgelist <- cbind(PCNedgelist, PTMCW) #Bind all the columns together. Now Data structure is PATHWAY | PATHWAY | Jaccard | CPE
   PCNedgelist <- PCNedgelist[rowSums(is.na(PCNedgelist)) != 2, ] #Remove all rows that only have NA values for the jaccard and CPE values
   PCNedgelist <- as.data.frame(PCNedgelist)
-  names(PCNedgelist) <- c("source", "target", "pathway_Jaccard_similarity", "PTM_cluster_evidence")
-  # Sort by the highest PTM cluster evidence
-  # PCNedgelist <- PCNedgelist[order(PCNedgelist$PTM_cluster_evidence, decreasing = TRUE),]
+  names(PCNedgelist) <- c("source", "target", "pathway_Jaccard_similarity", "PTM_cluster_weights")
+
   PCNedgelist$pathway_Jaccard_similarity <- as.numeric(PCNedgelist$pathway_Jaccard_similarity)
-  PCNedgelist$PTM_cluster_evidence <- as.numeric(PCNedgelist$PTM_cluster_evidence)
+  PCNedgelist$PTM_cluster_weights <- as.numeric(PCNedgelist$PTM_cluster_weights)
   # Convert NA to 0
   PCNedgelist[is.na(PCNedgelist)] <- 0
   # Subset data frames
   zero_jaccard <- PCNedgelist[PCNedgelist$pathway_Jaccard_similarity == 0, ]
   nonzero_jaccard <- PCNedgelist[PCNedgelist$pathway_Jaccard_similarity > 0, ]
 
-  # Sort both by PTM_cluster_evidence in decreasing order
-  zero_jaccard <- zero_jaccard[order(zero_jaccard$PTM_cluster_evidence, decreasing=TRUE), ]
-  nonzero_jaccard <- nonzero_jaccard[order(nonzero_jaccard$PTM_cluster_evidence, decreasing=TRUE), ]
+  # Sort both by PTM_cluster_weights in decreasing order
+  zero_jaccard <- zero_jaccard[order(zero_jaccard$PTM_cluster_weights, decreasing=TRUE), ]
+  nonzero_jaccard <- nonzero_jaccard[order(nonzero_jaccard$PTM_cluster_weights, decreasing=TRUE), ]
 
   # Combine results: zero-jaccard block on top, then nonzero-jaccard block
   PCNedgelist <- rbind(zero_jaccard, nonzero_jaccard)
 
   # For Cytoscape graphing
   #Remove all rows that only have NA values for  CPE values
-  bioplanetCPEedges <- PCNedgelist[!is.na(PCNedgelist[,"PTM_cluster_evidence"]), c("source", "target", "PTM_cluster_evidence")]
+  bioplanetCPEedges <- PCNedgelist[!is.na(PCNedgelist[,"PTM_cluster_weights"]), c("source", "target", "PTM_cluster_weights")]
   # For Cytoscape it's useful to have both types of edges for plotting in different colors
 
   # Assign interaction, required for Cytoscape
-  bioplanetCPEedges$interaction <- "PTM_cluster_evidence"
+  bioplanetCPEedges$interaction <- "PTM_cluster_weights"
   # Create pathway crosstalk network with individual cluster and bioplanet edges
   jaccard.net <- bioplanetjaccardedges
   names(jaccard.net) <- c("source", "target", "Weight", "interaction")
   CPE.net <- bioplanetCPEedges
   names(CPE.net) <- c("source", "target", "Weight", "interaction")
   pathway.crosstalk.network <- rbind(CPE.net, jaccard.net)
 
-  ### Save edgefile for cytoscape plotting ###
-
-  if(is.character(createfile)){ #Don't need to check if directory exists since was done above
-    saved.dir <- getwd()
-    setwd(createfile)
-    filename <- paste(PCN.edgelist.name, ".csv", sep="") #Name of the file created with .csv appended
-    utils::write.csv(pathway.crosstalk.network, file = filename, row.names = FALSE) #Save to files for cytoscape...
-
-    cat(filename, "made in directory:", getwd()) #Tell the user where their files got put
-    setwd(saved.dir)
-  }
   end_time <- Sys.time()
   message(end_time)
   #calculate difference between start and end time
   total_time <- end_time - start_time
   message(noquote(paste("Total time: ", total_time, sep="")))
-  return(list(pathway.crosstalk.network, PCNedgelist, pathways.list))
-}
-
-
-
+  return(list(pathway.crosstalk.network, PCNedgelist, pathways.list, CPE.matrix))
+}

R/data_documentation.R

@@ -302,7 +302,7 @@
 #'  \item{source}{The name of a pathway interacting}
 #'  \item{target}{The name of the pathway the source interacts with}
 #'  \item{pathway_Jaccard_similarity}{A weight}
-#'  \item{PTM_cluster_evidence}{A weight}
+#'  \item{PTM_cluster_weights}{A weight}
 #' }
 #' @source "Produced by BuildPathwayCrosstalkNetwork"
 #' @examples

data-raw/processed-brca.R

@@ -48,6 +48,5 @@ usethis::use_data(BRCA_genemania.edges, overwrite = TRUE)
 bioplanet.file <- system.file("extdata", "bioplanet_pathway_June2025.csv", package = "PTMsToPathways")
 # clusters are already made in brca_clusterlist_data, so we can just use those
 common.clusters <- brca_clusterlist_data[[1]]
-BRCA_PCN.data <- BuildPathwayCrosstalkNetwork(common.clusters, bioplanet.file,
-                                         createfile = FALSE)
+BRCA_PCN.data <- BuildPathwayCrosstalkNetwork(common.clusters, bioplanet.file)
 usethis::use_data(BRCA_PCN.data, overwrite = TRUE)

tests/testthat/test_buildpathwaycrosstalknetwork.R

@@ -2,14 +2,14 @@
 
 test_that("BuildPathwayCrosstalkNetwork() gives right answer", {
 
-  output <- suppressMessages(BuildPathwayCrosstalkNetwork(ex_common_clusters, ex_pathways_list, createfile = FALSE))
+  output <- suppressMessages(BuildPathwayCrosstalkNetwork(ex_common_clusters, ex_pathways_list))
 
   pathway_crosstalk_network <- output[[1]]
 
   row_7 <- as.list(pathway_crosstalk_network[7,])
   row_25 <- as.list(pathway_crosstalk_network[25,])
 
-  exp_row_7 <- list("Lipid and lipoprotein metabolism", "Vitamin B12 metabolism", 0.3, "PTM_cluster_evidence")
+  exp_row_7 <- list("Lipid and lipoprotein metabolism", "Vitamin B12 metabolism", 0.3, "PTM_cluster_weights")
   exp_row_25 <- list("ERBB signaling pathway", "Validated nuclear estrogen receptor alpha network", 0.025974025974026, "pathway_Jaccard_similarity")
 
   expect_setequal(row_7, exp_row_7)
@@ -41,6 +41,14 @@ test_that("BuildPathwayCrosstalkNetwork() gives right answer", {
   expect_true("CRP" %in% pathways_list$'Selenium pathway')
   expect_true("FAM120B" %in% pathways_list$'RXR/VDR pathway')
   expect_false("ALB" %in% pathways_list$'Mitochondrial fatty acid beta-oxidation')
+
+  cpe_matrix <- output[[4]]
+
+  expect_true(is.matrix(cpe_matrix))
+  expect_equal(nrow(cpe_matrix), length(ex_common_clusters))
+  expect_equal(ncol(cpe_matrix), length(ex_pathways_list))
+  expect_equal(rownames(cpe_matrix), names(ex_common_clusters))
+  expect_equal(colnames(cpe_matrix), names(ex_pathways_list))
 
 
 })

vignettes/BRCANetworks.Rmd

@@ -319,8 +319,7 @@ following:
 
 ```{r eval = FALSE}
 bioplanet.file <- system.file("extdata", "bioplanet_pathway_June2025.csv", package = "PTMsToPathways")
-PCN.data <- BuildPathwayCrosstalkNetwork(common.clusters, bioplanet.file,
-                                         createfile = FALSE)
+PCN.data <- BuildPathwayCrosstalkNetwork(common.clusters, bioplanet.file)
 ```
 
 Or load in pre-computed results from within the PTMsToPathways package:

vignettes/CreatingNetworks.Rmd

@@ -412,14 +412,25 @@ The function [`BuildPathwayCrosstalkNetwork`](references/BuildPathwayCrosstalkNe
 takes the list of pathways (or a path to the pathway file) and the list of clusters generated in step 1 and returns
 the relationships between pathways in two forms: a dataframe of edges that can be uploaded to Cytoscape (the first element of the returned list)
 and a dataframe with columns for the two edge weights (the second element of the returned list). The third element of the returned list is the list of pathways
-(in case the user provided a path to the pathway file rather than a list of pathways). The function can be run as follows:
+(in case the user provided a path to the pathway file rather than a list of pathways). The fourth element is the cluster-by-pathway
+CPE matrix used for the PTM cluster weight calculations. The function can be run as follows:
 ```{r eval = TRUE, echo = TRUE}
-PCN.data <- BuildPathwayCrosstalkNetwork(common.clusters, bioplanet.file,
-                                         createfile = FALSE)
+PCN.data <- BuildPathwayCrosstalkNetwork(common.clusters, bioplanet.file)
 pathway.crosstalk.network <- PCN.data[[1]]
 PCNedgelist <- PCN.data[[2]]
 pathways.list <- PCN.data[[3]]
+CPE.matrix <- PCN.data[[4]]
 ```
+
+If you want to save the Cytoscape-ready edgefile, save the returned
+`pathway.crosstalk.network` dataframe explicitly:
+
+```{r eval = FALSE, echo = TRUE}
+utils::write.csv(pathway.crosstalk.network,
+                 file = "PCN_file.csv",
+                 row.names = FALSE)
+```
+
 The required columns for Cytoscape are:
 ```{r eval = TRUE, echo = TRUE}
 names(pathway.crosstalk.network)
@@ -448,7 +459,7 @@ dat <- PCNedgelist[PCNedgelist$source == "Axon guidance" & PCNedgelist$target ==
 knitr::kable(dat, align = 'l', digits = 2)
 ```
 
-## Metric calculations
+## Pathway-pathway weight calculations
 
 ### Jaccard Similarity
  The Jaccard similarity is defined as the size of
@@ -458,16 +469,16 @@ of unique genes that are in either pathway).
 
 ```{r eval = TRUE, echo = TRUE}
 size_union <- length(union(pathways.list[["Axon guidance"]], pathways.list[["ERBB signaling pathway"]]))
-size_union
 size_intersection <- length(intersect(pathways.list[["Axon guidance"]], pathways.list[["ERBB signaling pathway"]]))
-size_intersection
 jaccard_similarity <- size_intersection / size_union
 jaccard_similarity
 ```
 
 Notice that the found value is the same as the `pathway_Jaccard_similarity` value in the `PCNedgelist` dataframe.
 
-### Cluster-Pathway Evidence (CPE) Score
+### PTM Cluster Weight 
+
+P2P also calculates a PTM cluster based weight between pathways (`pathway_cluster_weight`). 
 
 For a given pathway $j$ ($pw_j$) and a given cluster $i$ ($cl_i$), we define the Cluster-Pathway Evidence (CPE),
 $$
@@ -508,28 +519,36 @@ So CPE for the first cluster and the `Axon guidance` pathway is:
 ```{r eval = TRUE, echo = TRUE}
 1/(1*8) + 1/(2*8)
 ```
+We can see that this is the same value as the CPE score for the first cluster and the `Axon guidance` pathway in the `CPE.matrix`:
+```{r eval = TRUE, echo = TRUE}
+CPE.matrix[1, "Axon guidance"]
+```
 
-Overall, two pathways are considered to be related if they both have positive CPE with the same cluster.
-In that case, the CPE score for the edge between those two pathways is the sum of the CPE scores for all
-clusters that contain 
+Overall, two pathways are considered to be related if they both have positive CPE with at least one cluster.
+In that case, the PTM cluster weight for the two pathways is the sum of the CPE
+scores for all clusters that have positive CPE with both pathways. So if a
+cluster has positive CPE with both pathways, it contributes to the PTM cluster
+weight between those two pathways.
 
-For example, here are the number of clusters that contain PTMs on one or more proteins in the `Axon guidance` pathway:
-```{r eval = TRUE, echo = TRUE}
-ag_clusts <- sapply(common.clusters, function(x) any(sapply(pathways.list[["Axon guidance"]], grepl,x=x)))
-sum(ag_clusts)
+For the `Axon guidance` and `ERBB signaling pathway` pair, we can see all of those clusters and their CPE scores with both pathways in the `CPE.matrix`:
+```{r eval = FALSE, echo = TRUE}
+cols <- c("Axon guidance", "ERBB signaling pathway")
+rows <- rowSums(!is.na(CPE.matrix[, cols])) == 2
+CPE.matrix[rows, cols]
 ```
-And here are the number of clusters that contain PTMs on one or more proteins in the `ERBB signaling pathway`:
-```{r eval = TRUE, echo = TRUE}
-erbb_clusts <- sapply(common.clusters, function(x) any(sapply(pathways.list[["ERBB signaling pathway"]], grepl,x=x)))
-sum(erbb_clusts)
+```{r eval = TRUE, echo = FALSE}
+cols <- c("Axon guidance", "ERBB signaling pathway")
+rows <- rowSums(!is.na(CPE.matrix[, cols])) == 2
+dat <- CPE.matrix[rows, cols]
+knitr::kable(dat, align = 'l', digits = 2)
 ```
-And overall, here are the number of clustesr that contain PTMs on at least one protein from both the `Axon guidance` pathway and the `ERBB signaling pathway`:
-```{r eval = TRUE, echo = TRUE}
-sum(ag_clusts & erbb_clusts)
+And the total is the same as the `pathway_cluster_weight` for this pair above.
+
+```{r eval= TRUE, echo = TRUE}
+sum(CPE.matrix[rows, cols])
 ```
-The CPE scores from these clusters to the two pathways would be summed to get the final CPE score for the edge between the two pathways.
 
-For more detail on the CPE calculation, see ["Ross et al.,
+For more detail on the PTM cluster weight calculation, see ["Ross et al.,
 2023"](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1010690).
 
 # Saving Data