diff --git a/Rpkg/.Rbuildignore b/Rpkg/.Rbuildignore index 7afa938..7e5c00c 100644 --- a/Rpkg/.Rbuildignore +++ b/Rpkg/.Rbuildignore @@ -22,6 +22,7 @@ ^\.gitignore$ ^\.gitattributes$ ^\.gitmodules$ +^coherence_internal\.md$ # Python package and related files ^python/ @@ -58,6 +59,8 @@ ^scratch/ ^tmp/ ^temp/ +^analysis/.*\.md$ +^analysis/.*\.html$ # Build artifacts ^.*\.tar\.gz$ diff --git a/Rpkg/analysis/R/hClust/RawModuleSubModuleDotPlotOfGOTerms-1.pdf b/Rpkg/analysis/R/hClust/RawModuleSubModuleDotPlotOfGOTerms-1.pdf deleted file mode 100644 index cb87c91..0000000 Binary files a/Rpkg/analysis/R/hClust/RawModuleSubModuleDotPlotOfGOTerms-1.pdf and /dev/null differ diff --git a/Rpkg/analysis/R/hClust/RawModuleSubModuleDotPlotOfGOTerms-1.png b/Rpkg/analysis/R/hClust/RawModuleSubModuleDotPlotOfGOTerms-1.png deleted file mode 100644 index e2eed1f..0000000 Binary files a/Rpkg/analysis/R/hClust/RawModuleSubModuleDotPlotOfGOTerms-1.png and /dev/null differ diff --git a/Rpkg/analysis/R/hClust/RawZHeatmapAll-1.pdf b/Rpkg/analysis/R/hClust/RawZHeatmapAll-1.pdf deleted file mode 100644 index d7182c6..0000000 Binary files a/Rpkg/analysis/R/hClust/RawZHeatmapAll-1.pdf and /dev/null differ diff --git a/Rpkg/analysis/R/hClust/RawZHeatmapAll-1.png b/Rpkg/analysis/R/hClust/RawZHeatmapAll-1.png deleted file mode 100644 index 4048ab8..0000000 Binary files a/Rpkg/analysis/R/hClust/RawZHeatmapAll-1.png and /dev/null differ diff --git a/Rpkg/analysis/R/hClust/RawZHeatmapSAINTContamModules-1.pdf b/Rpkg/analysis/R/hClust/RawZHeatmapSAINTContamModules-1.pdf deleted file mode 100644 index 8fe0186..0000000 Binary files a/Rpkg/analysis/R/hClust/RawZHeatmapSAINTContamModules-1.pdf and /dev/null differ diff --git a/Rpkg/analysis/R/hClust/RawZHeatmapSAINTContamModules-1.png b/Rpkg/analysis/R/hClust/RawZHeatmapSAINTContamModules-1.png deleted file mode 100644 index a585ae4..0000000 Binary files a/Rpkg/analysis/R/hClust/RawZHeatmapSAINTContamModules-1.png and /dev/null differ diff --git a/Rpkg/analysis/R/hClust/RawZHeatmapSAINTHighHighModules-1.pdf b/Rpkg/analysis/R/hClust/RawZHeatmapSAINTHighHighModules-1.pdf deleted file mode 100644 index f5d5b0e..0000000 Binary files a/Rpkg/analysis/R/hClust/RawZHeatmapSAINTHighHighModules-1.pdf and /dev/null differ diff --git a/Rpkg/analysis/R/hClust/RawZHeatmapSAINTHighHighModules-1.png b/Rpkg/analysis/R/hClust/RawZHeatmapSAINTHighHighModules-1.png deleted file mode 100644 index fab30fe..0000000 Binary files a/Rpkg/analysis/R/hClust/RawZHeatmapSAINTHighHighModules-1.png and /dev/null differ diff --git a/Rpkg/analysis/R/hClust/RawZHeatmapSAINTLowLowModules-1.pdf b/Rpkg/analysis/R/hClust/RawZHeatmapSAINTLowLowModules-1.pdf deleted file mode 100644 index 205ce1d..0000000 Binary files a/Rpkg/analysis/R/hClust/RawZHeatmapSAINTLowLowModules-1.pdf and /dev/null differ diff --git a/Rpkg/analysis/R/hClust/RawZHeatmapSAINTLowLowModules-1.png b/Rpkg/analysis/R/hClust/RawZHeatmapSAINTLowLowModules-1.png deleted file mode 100644 index 82a8d35..0000000 Binary files a/Rpkg/analysis/R/hClust/RawZHeatmapSAINTLowLowModules-1.png and /dev/null differ diff --git a/Rpkg/analysis/R/hClust/RawZHeatmapSAINTSignalModule-1.pdf b/Rpkg/analysis/R/hClust/RawZHeatmapSAINTSignalModule-1.pdf deleted file mode 100644 index 52e1e1e..0000000 Binary files a/Rpkg/analysis/R/hClust/RawZHeatmapSAINTSignalModule-1.pdf and /dev/null differ diff --git a/Rpkg/analysis/R/hClust/RawZHeatmapSAINTSignalModule-1.png b/Rpkg/analysis/R/hClust/RawZHeatmapSAINTSignalModule-1.png deleted file mode 100644 index 001e58a..0000000 Binary files a/Rpkg/analysis/R/hClust/RawZHeatmapSAINTSignalModule-1.png and /dev/null differ diff --git a/Rpkg/analysis/R/hClust/RawZHeatmapTMEMModule-1.pdf b/Rpkg/analysis/R/hClust/RawZHeatmapTMEMModule-1.pdf deleted file mode 100644 index 2bbe5c7..0000000 Binary files a/Rpkg/analysis/R/hClust/RawZHeatmapTMEMModule-1.pdf and /dev/null differ diff --git a/Rpkg/analysis/R/hClust/RawZHeatmapTMEMModule-1.png b/Rpkg/analysis/R/hClust/RawZHeatmapTMEMModule-1.png deleted file mode 100644 index 5ecbd39..0000000 Binary files a/Rpkg/analysis/R/hClust/RawZHeatmapTMEMModule-1.png and /dev/null differ diff --git a/Rpkg/vignettes/.gitignore b/Rpkg/vignettes/.gitignore index c833a2c..d3dc0e4 100644 --- a/Rpkg/vignettes/.gitignore +++ b/Rpkg/vignettes/.gitignore @@ -3,3 +3,4 @@ .RData .Ruserdata inst/doc +coherence_internal.md diff --git a/Rpkg/vignettes/HierarchicalClustering.Rmd b/Rpkg/vignettes/HierarchicalClustering.Rmd index d5b832c..b3dc71f 100644 --- a/Rpkg/vignettes/HierarchicalClustering.Rmd +++ b/Rpkg/vignettes/HierarchicalClustering.Rmd @@ -1,7 +1,8 @@ --- -title: "IMPS Hierarchical Clustering" +title: "Hierarchical Clustering vignette" author: "Erik Larsen" date: "`r Sys.Date()`" +always_allow_html: true output: html_document: toc: true @@ -75,7 +76,6 @@ library(visNetwork) options(stringsAsFactors = FALSE) ``` - ## Load the Pre-Computed Data ```{r Data Load} @@ -85,7 +85,6 @@ load("../../analysis/R/precompute_objects.rda") ```{r Second Data Load, echo = FALSE} rm(IPMS_mat_log) - ``` @@ -189,24 +188,24 @@ if (identical(fmt, "html")) { } ``` In order to assess this method and also to compare/validate the hierarchical SAINT -model, describing modules/clusters based on the proteins in them is useful. +model, describing modules/clusters based on the proteins in them is useful Similarly localized proteins could be expected to have similar abundances to -each other as well as the target (**TMEM184B**). +each other as well as the target (**TMEM184B**) Similar functioning proteins could also be expected to have similar abundances -to each other as well as the target. +to each other as well as the target Describing these modules/clusters attained by `hierarchical clustering` via gene ontology gene set enrichment analysis is, therefore, useful and standard -practice. +practice Furthermore, we may be able to assemble a `protein-protein interaction network` -from the results. +from the results I outline this below-- the methods and results for acquiring these modules, their constituent proteins, and how we make sense of them-- visualized as -`dotplots`, `heatmaps`, and a `network`. I include searchable tables printed last. +`dotplots`, `heatmaps`, and a `network`. I include searchable tables printed last ## Extract Modules @@ -657,13 +656,19 @@ global_enrichment_results_df |> strip.background = ggplot2::element_rect(color = 'black', fill = 'white')) + ggplot2::guides(#shape = ggplot2::guide_legend(override.aes = list(size = 4, shape = c(19, 17))), - color = ggplot2::guide_legend(override.aes = list(size = 4, shape = 19, color = c('navy', 'darkgoldenrod3', 'firebrick'))), - alpha = ggplot2::guide_legend(override.aes = list(size = 5, shape = 19)), + color = ggplot2::guide_legend(override.aes = list(size = 4, + shape = 19, + color = c('navy', + 'darkgoldenrod3', + 'firebrick'))), + alpha = ggplot2::guide_legend(override.aes = list(size = 5, + shape = 19)), size = ggplot2::guide_legend(override.aes = list(shape = 19)) # size = 'none' ) + - ggplot2::scale_color_manual(values = c('navy', 'darkgoldenrod3', 'firebrick')) + - ggplot2::scale_fill_manual(values = c('navy', 'darkgoldenrod3', 'firebrick')) + ggplot2::scale_fill_manual(values = c('navy', + 'darkgoldenrod3', + 'firebrick')) ``` ## Module GO (similarity) Network @@ -817,7 +822,9 @@ if (identical(fmt, "html")) { ```{r RawZHeatmapAll, render = T} pheatmap::pheatmap( IPMS_z_ordered, - color = grDevices::colorRampPalette(c('navy', 'white', 'darkgoldenrod'))((200)), + color = grDevices::colorRampPalette(c('navy', + 'white', + 'darkgoldenrod'))((200)), cluster_rows = hc, cluster_cols = F, show_rownames = F, @@ -852,7 +859,9 @@ pheatmap::pheatmap( dplyr::pull(Protein) ), ], - color = grDevices::colorRampPalette(c('navy', 'white', 'darkgoldenrod'))((200)), + color = grDevices::colorRampPalette(c('navy', + 'white', + 'darkgoldenrod'))((200)), cluster_rows = T, cluster_cols = F, show_rownames = T, @@ -901,7 +910,9 @@ pheatmap::pheatmap( grepl(classification, pattern = 'High SAINT, Low')) |> dplyr::pull(Protein)),], - color = grDevices::colorRampPalette(c('navy', 'white', 'darkgoldenrod'))((200)), + color = grDevices::colorRampPalette(c('navy', + 'white', + 'darkgoldenrod'))((200)), cluster_rows = T, cluster_cols = F, show_rownames = F, @@ -937,7 +948,9 @@ pheatmap::pheatmap( grepl(classification, pattern = 'High SAINT, High')) |> dplyr::pull(Protein)),], - color = grDevices::colorRampPalette(c('navy', 'white', 'darkgoldenrod'))((200)), + color = grDevices::colorRampPalette(c('navy', + 'white', + 'darkgoldenrod'))((200)), cluster_rows = T, cluster_cols = F, show_rownames = F, @@ -955,7 +968,10 @@ pheatmap::pheatmap( ] ), - annotation_colors = list(Module = c('darkorange', 'darkviolet', 'lightpink', 'skyblue')), + annotation_colors = list(Module = c('darkorange', + 'darkviolet', + 'lightpink', + 'skyblue')), angle_col = 90, use_raster = F, border_color = NA, @@ -973,7 +989,9 @@ pheatmap::pheatmap( grepl(classification, pattern = 'Low SAINT, Low')) |> dplyr::pull(Protein)),], - color = grDevices::colorRampPalette(c('navy', 'white', 'darkgoldenrod'))((200)), + color = grDevices::colorRampPalette(c('navy', + 'white', + 'darkgoldenrod'))((200)), cluster_rows = T, cluster_cols = F, show_rownames = F, @@ -1016,7 +1034,9 @@ pheatmap::pheatmap( grepl(classification, pattern = 'Low SAINT, High')) |> dplyr::pull(Protein)),], - color = grDevices::colorRampPalette(c('navy', 'white', 'darkgoldenrod'))((200)), + color = grDevices::colorRampPalette(c('navy', + 'white', + 'darkgoldenrod'))((200)), cluster_rows = T, cluster_cols = F, show_rownames = F, @@ -1058,7 +1078,7 @@ if (identical(fmt, "html")) { } else if (identical(fmt, "markdown")) { cat('
Results Tables\n') } else { - cat('#Results Tables\n') + cat('# Results Tables\n') } ``` @@ -1081,15 +1101,18 @@ DT::datatable( score_mass), ~paste0("module_", .x)) |> dplyr::filter(!is.na(GO_Term)) |> dplyr::filter(!is.na(GO_Category)) |> - dplyr::filter((grepl(module, pattern = '12') & - hSAINT_score > 0.5) | - (grepl(classification, - pattern = 'High SAINT, Low') & - hSAINT_score > 0.5) - ) - ) + dplyr::filter( (grepl(module, + pattern = '12') & + hSAINT_score > 0.5) | + (grepl(classification, + pattern = 'High SAINT, Low') & + hSAINT_score > 0.5) + ) +) ``` +This (complete) table exists as [hClust_protein_level_GO_GSEA.rda](../../analysis/R/hClust/hClust_protein_level_GO_GSEA.rda) in [Endo/analysis/R/hClust/](../../analysis/R/hClust) + ## Clustering Results Table by Module ```{r Print Module Level Results Table} @@ -1107,7 +1130,7 @@ DT::datatable( ) ``` - +This table exists as [hClust_module_level_GO_GSEA.rda](../../analysis/R/hClust/hClust_module_level_GO_GSEA.rda) in [Endo/analysis/R/hClust/](../../analysis/R/hClust) ```{r, results='asis', echo=FALSE} fmt <- knitr::pandoc_to() diff --git a/analysis/HierarchicalClustering.html b/analysis/HierarchicalClustering.html new file mode 100644 index 0000000..a2314f7 --- /dev/null +++ b/analysis/HierarchicalClustering.html @@ -0,0 +1,8366 @@ + + + + + + + + + + + + + + + +Hierarchical Clustering vignette + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + + + + + + +
+

Overview

+

This markdown was developed to formally document an ad-hoc cluster +analysis for TMEM184B immunoprecipitation tandem mass-spec experiments +for Dr. Martha Bhattacharya’s lab at the +University of Arizona for the endolysosomal paper under +revision (2025-02-28) for the Journal of Cell Science

+

This vignette documents the Hierarchical Clustering of +this IPMS data

+
+

Attach Packages

+
# data/wrangling
+library(tidyverse)
+library(Endo)
+library(DT)
+
+# clustering
+library(pheatmap)
+library(dynamicTreeCut)
+
+# GSEA/viz
+library(clusterProfiler)
+library(org.Hs.eg.db)
+library(ReactomePA)
+library(AnnotationDbi)
+
+# progress bar
+library(progressr)
+
+# parallelization
+library(future)
+library(future.apply)
+
+# network viz
+library(visNetwork)
+
+options(stringsAsFactors = FALSE)
+
+
+

Load the Pre-Computed Data

+
# Load precomputed objects from preprocessing
+load("../../analysis/R/precompute_objects.rda")
+
+
+
+

Hierarchical Clustering

+
+

Raw Counts

+
d <- stats::as.dist(1 - stats::cor(t(IPMS_mat), use = 'pairwise.complete.obs'))
+hc <- stats::hclust(d, method = "average")
+
sizes <- c(20, 30, 50, 80, 100)
+results <- list()
+
+for (m in sizes) {
+  cl <- dynamicTreeCut::cutreeDynamic(
+    dendro = hc,
+    distM = as.matrix(d),
+    deepSplit = 1,
+    pamRespectsDendro = TRUE,
+    minClusterSize = m
+  )
+  results[[as.character(m)]] <- cl
+  cat("minClusterSize =", m, "→", length(table(cl)), "modules\n")
+}
+
##  ..cutHeight not given, setting it to 0.976  ===>  99% of the (truncated) height range in dendro.
+##  ..done.
+## minClusterSize = 20 → 50 modules
+##  ..cutHeight not given, setting it to 0.976  ===>  99% of the (truncated) height range in dendro.
+##  ..done.
+## minClusterSize = 30 → 42 modules
+##  ..cutHeight not given, setting it to 0.976  ===>  99% of the (truncated) height range in dendro.
+##  ..done.
+## minClusterSize = 50 → 28 modules
+##  ..cutHeight not given, setting it to 0.976  ===>  99% of the (truncated) height range in dendro.
+##  ..done.
+## minClusterSize = 80 → 21 modules
+##  ..cutHeight not given, setting it to 0.976  ===>  99% of the (truncated) height range in dendro.
+##  ..done.
+## minClusterSize = 100 → 17 modules
+
rm(list = c('cl', 'm', 'sizes', 'results'))
+
    +

  • optimal minimum cluster size is 50 +(elbow; diminishing returns after)

    • +

  • # modules = 28

    • +
+
clusters_opt <- dynamicTreeCut::cutreeDynamic(hc,
+                                              distM = as.matrix(d),
+                                              deepSplit = 1,
+                                              pamRespectsDendro = T,
+                                              minClusterSize = 50)
+
##  ..cutHeight not given, setting it to 0.976  ===>  99% of the (truncated) height range in dendro.
+##  ..done.
+
names(clusters_opt) <- rownames(IPMS_mat)
+
module_ids <- sort(unique(clusters_opt))
+module_colors <- stats::setNames(
+  grDevices::rainbow(length(module_ids), s = 0.6, v = 0.9),
+  module_ids
+)
+
+annotation_colors <- list(Module = module_colors)
+

Sparse data is too difficult to discern from a global view

+

-> Z transform the raw counts after clustering to +plot

+
    +

  • the cluster object preserves the ordering of the proteins +(dendrogram)

    +
      +
    • Z-score therefore only changes the values of the columns and rows to +better visualize the sparse data
      • +
    • +
+
IPMS_z <- t(scale(t(IPMS_mat)))
+colnames(IPMS_z) <- colnames(IPMS_mat)
+
+condition <- sub('_.*', '', colnames(IPMS_mat))
+
+col_order <- unlist(
+  lapply(
+    split(seq_along(condition), condition), function(idx) {
+      submat <- IPMS_z[, idx, drop = FALSE]
+      hc_sub <- stats::hclust(stats::dist(t(submat)), method = 'ward.D2')
+      idx[hc_sub$order]
+    }))
+
+IPMS_z_ordered <- IPMS_z[, col_order]
+
+rm(list = c('condition', 'col_order'))
+
+
+
+

Module Descriptions

+

In order to assess this method and also to compare/validate the +hierarchical SAINT model, describing modules/clusters based on the +proteins in them is useful

+

Similarly localized proteins could be expected to have similar +abundances to each other as well as the target +(TMEM184B)

+

Similar functioning proteins could also be expected to have similar +abundances to each other as well as the target

+

Describing these modules/clusters attained by +hierarchical clustering via gene ontology gene set +enrichment analysis is, therefore, useful and standard practice

+

Furthermore, we may be able to assemble a +protein-protein interaction network from the results

+

I outline this below– the methods and results for acquiring these +modules, their constituent proteins, and how we make sense of them– +visualized as dotplots, heatmaps, and a +network. I include searchable tables printed last

+
+

Extract Modules

+

Extract gene names from the hierarchically-clustered modules and +convert to a dataframe

+
modules_list_idx <- lapply(sort(unique(clusters_opt)), function(m) {
+  which(clusters_opt == m)
+})
+names(modules_list_idx) <- sort(unique(clusters_opt))
+
+modules_list <- lapply(sort(unique(clusters_opt)), function(m) {
+  rownames(IPMS_mat)[clusters_opt == m]
+})
+names(modules_list) <- sort(unique(clusters_opt)) |> as.character()
+
+modules_df <- modules_list %>%
+  purrr::map_df(., .f = as.data.frame, .id = 'module') |> 
+  dplyr::rename(Protein = 2) |> 
+  dplyr::relocate(Protein, .before = 1)
+
+
+

Join Modules to hSAINT

+

Include score summaries (mean and median posterior probability +scores, ranks, and literature or CRAPome contamination rate)

+
module_scores <- modules_df |>  
+  dplyr::inner_join(
+    hSAINT
+  ) |> 
+  dplyr::arrange(desc(score)) |> 
+  dplyr::group_by(module) |> 
+  dplyr::summarize(mean_score = mean(score),
+                   median_score = stats::median(score),
+                   score_mass = sum(score),
+                   n = dplyr::n(),
+                   lit_flag_perc = (sum(lit_flag)/n)*100,
+                   CRAPome_flag_perc = (sum(CRAPome_flag)/n)*100) |> 
+  dplyr::arrange(desc(mean_score)) |> 
+  dplyr::mutate(mean_rank = dplyr::row_number()) |> 
+  dplyr::arrange(desc(median_score)) |> 
+  dplyr::mutate(median_rank = dplyr::row_number()) |> 
+  dplyr::arrange(desc(score_mass)) |> 
+  dplyr::mutate(total_rank = dplyr::row_number()) |> 
+  dplyr::arrange(mean_rank)
+

Join the module scores to the modules (per protein) dataframe and the +hSAINT results for classification

+
protein_module_df <- modules_df |> 
+  dplyr::inner_join(
+    hSAINT
+  ) |> 
+  dplyr::full_join(
+    module_scores |> 
+      dplyr::ungroup()
+  ) |> 
+  dplyr::mutate(
+    SAINT_classification = dplyr::case_when(mean_rank <= 5 ~ 'High SAINT',
+                                            mean_rank > 5 ~ 'Low SAINT'
+                                            ),
+    
+    contams_classification = dplyr::case_when(lit_flag_perc > 5 ~'High Contaminants',
+                                              TRUE ~ 'Low Contaminants')) |> 
+  dplyr::mutate(classification = paste0(SAINT_classification,
+                                        ", ",
+                                        contams_classification)) |> 
+  dplyr::select(-SAINT_classification,
+                -contams_classification) |> 
+  dplyr::relocate(classification, .after = module)
+
+
+

Module GO GSEA

+
EntrezIDs <- AnnotationDbi::select(
+  org.Hs.eg.db,
+  keys = modules_df |> dplyr::distinct(Protein) |> dplyr::pull(),
+  columns = 'ENTREZID',
+  keytype = 'SYMBOL'
+) |> 
+  dplyr::filter(!is.na(2))
+
+EntrezMap <- stats::setNames(EntrezIDs$ENTREZID, EntrezIDs$SYMBOL)
+

Create enrichment functions

+
background_genes <- EntrezIDs |> # use the EntrezIDs mapped from the modules_df
+  dplyr::pull(ENTREZID)
+
+get_GO <- function(genes, ontology, background_genes) {
+  clusterProfiler::enrichGO(
+    gene = genes,
+    universe = background_genes,
+    OrgDb = org.Hs.eg.db,
+    keyType = 'ENTREZID',
+    ont = ontology,
+    pAdjustMethod = 'BH',
+    qvalueCutoff = 0.05,
+    readable = T
+  )
+}
+
+get_kegg <- function(genes, EntrezMap) {
+  entrez <- EntrezMap[genes]
+  entrez <- entrez[!is.na(entrez)]
+  if (length(entrez)==0) {
+    return(NULL)
+  }
+  clusterProfiler::enrichKEGG(
+    gene = entrez,
+    organism = 'hsa',
+    pAdjustMethod = 'BH',
+    qvalueCutoff = 0.05
+  )
+}
+
+get_reactome <- function(genes, EntrezMap) {
+  entrez <- EntrezMap[genes]
+  entrez <- entrez[!is.na(entrez)]
+  if (length(entrez)==0) {
+    return(NULL)
+  }
+  ReactomePA::enrichPathway(
+    gene = entrez,
+    organism = 'human',
+    pAdjustMethod = 'BH',
+    qvalueCutoff = 0.05,
+    readable = T
+  )
+}
+

Run the functions on the dataset in parallel (multiple CPUs)

+
future::plan(future::multisession, workers = 6)
+progressr::handlers(progressr::handler_txtprogressbar())
+
+module_ids <- sort(unique(modules_df$module))
+
+progressr::with_progress({
+  p <- progressr::progressor(steps = length(module_ids))
+
+  enrichment_results <- future.apply::future_lapply(
+    module_ids,
+    function(mod) {
+
+      # tick the progress bar
+      p(sprintf("Module %s", mod))
+
+      genes <- modules_df |>
+        dplyr::filter(module == mod) |>
+        dplyr::inner_join(EntrezIDs |> dplyr::rename(Protein = 1)) |> 
+        dplyr::distinct(ENTREZID) |>
+        dplyr::pull(ENTREZID)
+
+      list(
+        module = mod,
+        go_bp  = get_GO(genes, ontology = "BP", universe = background_genes),
+        go_cc  = get_GO(genes, ontology = "CC", universe = background_genes),
+        go_mf  = get_GO(genes, ontology = "MF", universe = background_genes)
+        # add KEGG/Reactome later
+      )
+    }
+  )
+})
+
+names(enrichment_results) <- module_ids
+rm(p)
+

Gather enrichment results

+
global_enrichment_results_df <- lapply(enrichment_results, function(x) {
+  x %>%
+    purrr::map_df(., .f = as.data.frame, .id = 'GO_Category') |> 
+    dplyr::rename(module = 2, GO_Term = Description, GO_ID = ID) |> 
+    tidyr::fill(module, .direction = 'down') |> 
+    dplyr::filter(!is.na(GO_ID)) |> 
+    dplyr::mutate(GO_Category = gsub(GO_Category,
+                                     pattern = 'go_',
+                                     replacement = '') %>% 
+                    stringr::str_to_upper(.)) |> 
+    tibble::remove_rownames() |> 
+    dplyr::distinct()
+  }
+) %>%
+  purrr::map_df(., .f = as.data.frame)
+

The “High SAINT score” modules were among the modules +with significant enrichment of all the proteins detected in the +experiments

+

However, many modules did not have significant enrichment relative to +the entire set of proteins in the experiments, including +TMEM’s. Using hierarchical SAINT score as a +explanatory tool, this is likely because the modules without GO +enrichment contained a large mix of proteins scored as “noise” (low +posterior probability of being a true interactor) and signal (high +posterior probability of being a true interactor)

+

In order to better make sense of the clustering– describing the +“un-enriched” modules– “sub-module” GO GSEA was +conducted

+
+
+

SubModule GO GSEA

+

Sub-module enrichment analysis on modules without any enrichment (to +describe themes of proteins within these modules)

+
"%notin%" <- Negate("%in%")
+
+# index the modules without significant enrichment to describe sub-modules
+tier2_modules <- c(seq(1:28) |> as.character())[
+  which(c(seq(1:28) |> as.character()) %notin% c(global_enrichment_results_df |> 
+          dplyr::distinct(module) |> 
+          dplyr::pull()) == T)]
+
# function to gather sub-modules for each module that did not return enrichment results
+submodules_from_hclust <- function(
+  hc_object,
+  modules_df,
+  n_sub_per_module = 3L
+) {
+  # number of existing modules
+  n_modules <- length(unique(modules_df$module))
+
+  # total number of clusters in the finer cut
+  K_big <- n_modules * n_sub_per_module
+
+  # finer cut of the same tree
+  sub_labels_all <- stats::cutree(hc_object, k = K_big)
+
+  # assemble into a tibble
+  sub_df <- tibble::tibble(
+    Protein   = names(sub_labels_all),
+    sub_label = sub_labels_all
+  )
+
+  # attach existing module labels
+  sub_df <- dplyr::inner_join(
+    sub_df,
+    modules_df,
+    by = "Protein"
+  )
+
+  # define submodules as (module, sub_label) pairs
+  submodules_df <- sub_df %>%
+    dplyr::group_by(module, sub_label) %>%
+    dplyr::mutate(
+      submodule = paste0("sub", dplyr::cur_group_id())
+    ) %>%
+    dplyr::ungroup() %>%
+    dplyr::select(Protein, module, submodule)
+
+  return(submodules_df)
+}
+
# get the sub-modules
+submodules_df <- submodules_from_hclust(
+  hc_object  = hc,
+  modules_df = dplyr::filter(modules_df, module %in% tier2_modules),
+  n_sub_per_module = 3L
+)
+
# perform GO GSEA on the submodules
+
+submodule_enrichment_results <- submodules_df %>%
+  dplyr::inner_join(EntrezIDs |> dplyr::rename(Protein = 1)) |> 
+  dplyr::group_by(module, submodule) %>%
+  dplyr::group_modify(~ {
+
+    # submodule genes
+    genes <- .x$ENTREZID
+    
+    background_genes <- unique(submodules_df$ENTREZID[submodules_df$module == .y$module])
+
+    tibble::tibble(
+      BP = list(get_GO(genes, ontology = "BP", background_genes = background_genes)),
+      CC = list(get_GO(genes, ontology = "CC", background_genes = background_genes)),
+      MF = list(get_GO(genes, ontology = "MF", background_genes = background_genes))
+    )
+  }) %>%
+  dplyr::ungroup()
+rm(background_genes)
+
# wrangle the sub-module GO GSEA results
+submodule_enrichment_results_df <- submodule_enrichment_results |> 
+    tidyr::pivot_longer(
+      cols      = c('BP', 'CC', 'MF'),
+      names_to  = "GO_Category",
+      values_to = "values"
+    ) %>%
+    dplyr::mutate(
+      GO_tbl = purrr::map(values, ~ as.data.frame(.x))
+    ) |>
+    tidyr::unnest(GO_tbl) |>
+    dplyr::rename(
+      GO_ID   = ID,
+      GO_Term = Description
+    ) |>
+    dplyr::select(
+      module,
+      submodule,
+      GO_Category,
+      GO_Term,
+      GO_ID,
+      dplyr::everything(),
+      -values
+    )
+
modules_with_submodule_enrichment <- tier2_modules
+modules_without_submodule_enrichment <- generics::setdiff(modules_with_submodule_enrichment,
+                                                          submodule_enrichment_results_df |> 
+                                                            dplyr::distinct(module) |> 
+                                                            dplyr::pull())
+modules_with_submodule_enrichment <- submodule_enrichment_results_df |> 
+  dplyr::distinct(module) |> 
+  dplyr::pull()
+rm(tier2_modules)
+
+
+

Construct the Module GO (similarity) Network

+
go_network <- global_enrichment_results_df |> 
+  dplyr::full_join(submodule_enrichment_results_df) |> 
+  dplyr::full_join(module_scores) |> 
+  dplyr::mutate(Tier = dplyr::case_when(is.na(submodule) ~ 'Tier1',TRUE~'Tier2'),
+                node_to = GO_Term,
+                node_from = paste0(
+                  
+                  "Module", module,
+                  ifelse(is.na(submodule),
+                         '',
+                         paste0('_SubModule',
+                                submodule %>% 
+                                  stringr::str_extract(., pattern = '\\d+')))),
+                weight = -log10(p.adjust)) |> 
+  dplyr::filter(weight > -log10(0.05))
+
+g <- igraph::graph_from_data_frame(
+  d = go_network[, c("node_from", "node_to",
+                     "weight", "geneID", "p.adjust", "FoldEnrichment")],
+  directed = TRUE,
+  vertices =
+    dplyr::bind_rows(
+      go_network |> dplyr::transmute(name = node_from, score_mass, submodule),
+      go_network |> dplyr::transmute(name = node_to,
+                                     score_mass = NA_real_,
+                                     submodule = NA_character_)
+    ) |>
+    dplyr::distinct(name, .keep_all = TRUE)
+)
+
+
+
+

Plots

+
+

Module Dot Plot

+
global_enrichment_results_df |> 
+  dplyr::full_join(submodule_enrichment_results_df) |> 
+  dplyr::full_join(module_scores) |> 
+  dplyr::filter(!is.na(GO_Term)) |> 
+  dplyr::filter(!is.na(GO_Category)) |> 
+  dplyr::mutate(
+    Tier = dplyr::case_when(
+      is.na(submodule) ~ paste0('Mod. ', module),
+      TRUE ~ paste0('Mod. ', module,
+                    ':\nsubmod. ',
+                    gsub(submodule, pattern = 'sub', replacement = ''))
+      ),
+    is_enriched = dplyr::case_when(is.na(submodule) ~ 'Enriched',
+                                   TRUE ~ 'Not Enriched')
+    ) |> 
+  dplyr::group_by(module) |>  
+  dplyr::slice_head(n = 3) |>
+  dplyr::arrange(p.adjust) |>
+  ggplot2::ggplot(
+    ggplot2::aes(
+      x = -log10(p.adjust),
+      y = forcats::fct_rev(forcats::fct_infreq(GO_Term, mean_score)),
+      color = GO_Category,
+      size = FoldEnrichment,
+      alpha = mean_score,
+      shape = is_enriched,
+    )
+  ) +
+  ggplot2::geom_vline(xintercept = -log10(0.05),
+                      linetype = 'dashed',
+                      color = 'firebrick',
+                      size = 1,
+                      alpha = 0.5) +
+  ggplot2::geom_point() +
+  ggplot2::facet_grid(Tier ~ ., scales = 'free_y') +
+  ggplot2::labs(title = 'GO Enrichment Across\nModules and SubModules\nin TMEM184B IPMS\nHierarchical Clustering',
+                x = '-log10(Adj.P-val)',
+                y = 'GO Term',
+                alpha = 'Module\nMean\nHierarchical\nSAINT\nScore',
+                size = 'FoldEnrichment',
+                color = 'GO Category',
+                shape = 'IsEnrichedModule') +
+  ggplot2::theme_bw() +
+  ggplot2::theme(title = ggplot2::element_text(face = 'bold', size = 12),
+                 axis.text = ggplot2::element_text(face = 'bold', size = 10),
+                 strip.text = ggplot2::element_text(face = 'bold', size = 7), 
+                 strip.background = ggplot2::element_rect(color = 'black',
+                                                          fill = 'white')) +
+  ggplot2::guides(#shape = ggplot2::guide_legend(override.aes = list(size = 4, shape = c(19, 17))),
+                  color = ggplot2::guide_legend(override.aes = list(size = 4,
+                                                                    shape = 19,
+                                                                    color = c('navy',
+                                                                              'darkgoldenrod3',
+                                                                              'firebrick'))),
+                  alpha = ggplot2::guide_legend(override.aes = list(size = 5,
+                                                                    shape = 19)),
+                  size = ggplot2::guide_legend(override.aes = list(shape = 19))
+                  # size = 'none'
+                  ) +
+  ggplot2::scale_fill_manual(values = c('navy',
+                                        'darkgoldenrod3',
+                                        'firebrick'))
+

+
+
+

Module GO (similarity) Network

+
coords <- igraph::layout_with_fr(g)
+coords <- coords*200
+
+visNetwork::visNetwork(
+  nodes = data.frame(
+    id    = igraph::V(g)$name,
+    label = igraph::V(g)$name,
+    size = scales::rescale(
+      (\(x){ x[is.na(x)] <- 0; x })(igraph::V(g)$score_mass),
+      to = c(20, 60)
+    ),
+    x = coords[,1],
+    y = coords[,2],
+    group = base::interaction(
+                        {
+                          x <- igraph::V(g)$module
+                          if (length(x) == 0) rep("", igraph::vcount(g)) else x
+                        },
+                        {
+                          x <- igraph::V(g)$submodule
+                          if (length(x) == 0) {
+                            rep("", igraph::vcount(g))
+                          } else {
+                            ifelse(nzchar(x), x, igraph::V(g)$module)
+                          }
+                        },
+                        drop = TRUE,
+                        lex.order = TRUE
+                      ),
+    title = sapply(
+      igraph::V(g),
+      function(v) {
+        incoming_edges <- igraph::incident(g, v, mode = "all")
+        incoming_edges <- incoming_edges[
+          igraph::ends(g, incoming_edges)[,2] == igraph::V(g)[v]$name
+        ]
+        raw_gene_strings <- igraph::E(g)[incoming_edges]$geneID
+        raw_gene_strings <- unlist(raw_gene_strings, use.names = FALSE)
+
+        genes <- raw_gene_strings |>
+          as.character() |>
+          (\(x) x[!is.na(x)])() |>
+          strsplit("\\\\") |>
+          unlist() |>
+          (\(x) x[nzchar(x)])() |>
+          unique()
+
+        term <- igraph::V(g)[v]$name
+        saint_value <- igraph::V(g)[v]$score_mass
+        if (is.na(saint_value)) saint_value <- 0
+        is_module <- nzchar(igraph::V(g)[v]$submodule)
+
+        if (is_module) {
+          paste0(
+            "<div style='background-color:#ffffff;color:#1f4e79;",
+            "padding:6px 10px;border-radius:6px;",
+            "border:1px solid #1f4e79;font-size:12px;'>",
+            "<b>", term, "</b><br>",
+            "<b>SAINT mass:</b> ", signif(saint_value, 3), "<br>",
+            if (length(genes) > 0)
+              paste0("<u>Proteins:</u><br>", paste(genes, collapse = "<br>"))
+            else "",
+            "</div>"
+          )
+        } else {
+          paste0(
+            "<div style='background-color:#ffffff;color:#1f4e79;",
+            "padding:6px 10px;border-radius:6px;",
+            "border:1px solid #1f4e79;font-size:12px;'>",
+            "<b>", term, "</b><br>",
+            if (length(genes) > 0)
+              paste0("<u>Proteins:</u><br>", paste(genes, collapse = "<br>"))
+            else "",
+            "</div>"
+          )
+        }
+      }
+    )
+  ),
+  edges = data.frame(
+    from  = igraph::as_data_frame(g, "edges")$from,
+    to    = igraph::as_data_frame(g, "edges")$to,
+    value = igraph::E(g)$weight,
+    title = vapply(
+      seq_along(igraph::E(g)),
+      function(i) {
+        paste0(
+          "<div style='background-color:#ffffff;color:#1f4e79;",
+          "padding:6px 10px;border-radius:6px;",
+          "border:1px solid #1f4e79;font-size:12px;'>",
+          "<b>p.adjust:</b> ", signif(igraph::E(g)$p.adjust[i], 3), "<br>",
+          "<b>weight:</b> ", signif(igraph::E(g)$weight[i], 3),
+          "</div>"
+        )
+      },
+      FUN.VALUE = character(1)
+    )
+  )
+) |>
+  # visNetwork::visNodes(color = NULL) |>
+  visNetwork::visEdges(arrows = "to") |>
+  visNetwork::visNodes(
+    scaling = list(min = 5, max = 40),
+    shadow  = TRUE
+  ) |>
+  visNetwork::visLayout(randomSeed = 1) |>
+  visNetwork::visPhysics(enabled = FALSE) |>
+  visNetwork::visInteraction(
+    dragView = TRUE,
+    zoomView = TRUE
+  ) |>
+  visNetwork::visOptions(nodesIdSelection = TRUE)
+
+ +
rm(coords)
+

Modules with enriched GO Terms are blue

+

Here, only the enriched sub-modules were included: the modules with +that did not reach significant enrichment broadly, but did reach +significant enrichment when that module was broken up. These +Enriched sub-modules are +non-blue

+

Modules 2, 9, 15, +and 26 are of interest, having terms in line with what +might be expected from a TMEM experiment, and having robust +Hierarchical SAINT scores (hover over Module/SubModule +nodes)

+

TMEM184B, itself, was in Module 12, +which did NOT return any enriched GO Terms, likely due to its +membrane protein nature (and putative localization, where it may +progress through the endomembrane system and collide with a wide variety +of proteins). This is common among TM protein IPMS

+
+
+

Heatmaps

+
+

Raw Counts Heatmap

+
pheatmap::pheatmap(
+  IPMS_z_ordered,
+  color = grDevices::colorRampPalette(c('navy',
+                                        'white',
+                                        'darkgoldenrod'))((200)),
+  cluster_rows = hc,
+  cluster_cols = F,
+  show_rownames = F,
+  annotation_row = data.frame(Module = clusters_opt),
+  annotation_colors = annotation_colors,
+  angle_col = 90,
+  use_raster = F,
+  border_color = NA,
+  main = 'Hierarchically-Clustered\nTMEM184B IPMS TSCs, Z-scored'
+)
+

+
+
+

Raw Counts Heatmap (TMEM Module)

+
pheatmap::pheatmap(
+  IPMS_z_ordered[rownames(IPMS_z_ordered) %in% 
+                   
+                   c(protein_module_df |> 
+                       dplyr::full_join(global_enrichment_results_df) |> 
+                       dplyr::full_join(submodule_enrichment_results_df) |> 
+                       dplyr::select(-c(GO_Category:Count)) |> 
+                       dplyr::distinct() |> 
+                       dplyr::filter(
+                         module == protein_module_df |> 
+                           dplyr::full_join(global_enrichment_results_df) |> 
+                           dplyr::full_join(submodule_enrichment_results_df) |> 
+                           dplyr::select(-c(GO_Category:Count)) |> 
+                           dplyr::filter(Protein == 'TMEM184B') |> 
+                           dplyr::pull(module)
+                         ) |> 
+                       dplyr::pull(Protein)
+                     ),
+                 ],
+  color = grDevices::colorRampPalette(c('navy',
+                                        'white',
+                                        'darkgoldenrod'))((200)),
+  cluster_rows = T,
+  cluster_cols = F,
+  show_rownames = T,
+  fontsize_row = 5,
+  annotation_row = data.frame(
+    
+    Module12_SubModule = clusters_opt[
+      names(clusters_opt) %in% 
+        
+        c(protein_module_df |> 
+                       dplyr::full_join(global_enrichment_results_df) |> 
+                       dplyr::full_join(submodule_enrichment_results_df) |> 
+                       dplyr::select(-c(GO_Category:Count)) |> 
+                       dplyr::distinct() |> 
+                       dplyr::filter(
+                         module == protein_module_df |> 
+                           dplyr::full_join(global_enrichment_results_df) |> 
+                           dplyr::full_join(submodule_enrichment_results_df) |> 
+                           dplyr::select(-c(GO_Category:Count)) |> 
+                           dplyr::filter(Protein == 'TMEM184B') |> 
+                           dplyr::pull(module)
+                         ) |> dplyr::filter(!is.na(Protein)) |> 
+            dplyr::pull(Protein)
+                     )
+      ]
+    ),
+  annotation_colors = list(Module12_SubModule = 'darkorange2'),
+  angle_col = 90,
+  use_raster = F,
+  border_color = NA,
+  main = 'Hierarchically-Clustered TMEM184B IPMS TSCs\nZ-scored, TMEM Module'
+)
+

+

TMEM’s module (12) is not siginficantly +enriched above the proteins detected by the mass spectrometer.

+

Submodules within module 12 are not enriched

+
+
+

High hSAINT Score, Low Contamination Module Heatmap (Signal)

+
pheatmap::pheatmap(
+  IPMS_z_ordered[rownames(IPMS_z_ordered) %in% c(protein_module_df |> 
+                   dplyr::filter(Protein == 'TMEM184B' |
+                                   grepl(classification,
+                                         pattern = 'High SAINT, Low')) |> 
+                   dplyr::pull(Protein)),],
+  color = grDevices::colorRampPalette(c('navy',
+                                        'white',
+                                        'darkgoldenrod'))((200)),
+  cluster_rows = T,
+  cluster_cols = F,
+  show_rownames = F,
+  annotation_row = data.frame(
+    
+    Module = clusters_opt[
+      names(clusters_opt) %in% 
+        
+        c(protein_module_df |> 
+            dplyr::filter(Protein == 'TMEM184B' |
+                            grepl(classification,
+                                  pattern = 'High SAINT, Low')) |> 
+            dplyr::pull(Protein)
+          )
+      ]
+    ),
+  
+  annotation_colors = list(Module = c('firebrick', 'navy')),
+  angle_col = 90,
+  use_raster = F,
+  border_color = NA,
+  cutree_rows = 3,
+  main = 'Hierarchically-Clustered TMEM184B IPMS TSCs\nZ-scored, High SAINT Score, Low Contamination Cluster + TMEM'
+)
+

+
+
+

High hSAINT Score, High Contamination Modules Heatmap +(Sticky/Noise)

+
pheatmap::pheatmap(
+  IPMS_z_ordered[rownames(IPMS_z_ordered) %in% c(protein_module_df |> 
+                   dplyr::filter(Protein == 'TMEM184B' |
+                                   grepl(classification,
+                                         pattern = 'High SAINT, High')) |> 
+                   dplyr::pull(Protein)),],
+  color = grDevices::colorRampPalette(c('navy',
+                                        'white',
+                                        'darkgoldenrod'))((200)),
+  cluster_rows = T,
+  cluster_cols = F,
+  show_rownames = F,
+  annotation_row = data.frame(
+    
+    Module = clusters_opt[
+      names(clusters_opt) %in% 
+        
+        c(protein_module_df |> 
+            dplyr::filter(Protein == 'TMEM184B' |
+                            grepl(classification,
+                                  pattern = 'High SAINT, High')) |> 
+            dplyr::pull(Protein)
+          )
+      ]
+    ),
+  
+  annotation_colors = list(Module = c('darkorange',
+                                      'darkviolet',
+                                      'lightpink',
+                                      'skyblue')),
+  angle_col = 90,
+  use_raster = F,
+  border_color = NA,
+  cutree_rows = 2,
+  main = 'Hierarchically-Clustered TMEM184B IPMS TSCs\nZ-scored, High SAINT Score, High Contamination Clusters'
+)
+

+
+
+

Low hSAINT Score, Low Contamination Modules Heatmap (Low +Abundance)

+
pheatmap::pheatmap(
+  IPMS_z_ordered[rownames(IPMS_z_ordered) %in% c(protein_module_df |> 
+                   dplyr::filter(Protein == 'TMEM184B' |
+                                   grepl(classification,
+                                         pattern = 'Low SAINT, Low')) |> 
+                   dplyr::pull(Protein)),],
+  color = grDevices::colorRampPalette(c('navy',
+                                        'white',
+                                        'darkgoldenrod'))((200)),
+  cluster_rows = T,
+  cluster_cols = F,
+  show_rownames = F,
+  annotation_row = data.frame(
+    
+    Module = clusters_opt[
+      names(clusters_opt) %in% 
+        
+        c(protein_module_df |> 
+            dplyr::filter(Protein == 'TMEM184B' |
+                            grepl(classification,
+                                  pattern = 'Low SAINT, Low')) |> 
+            dplyr::pull(Protein)
+          )
+      ]
+    ),
+  
+  annotation_colors = list(Module = annotation_colors[[1]][
+    clusters_opt[
+      names(clusters_opt) %in% c(protein_module_df |> 
+                   dplyr::filter(Protein == 'TMEM184B' |
+                                   grepl(classification,
+                                         pattern = 'Low SAINT, Low')) |> 
+                   dplyr::pull(Protein))
+    ] |> as.numeric() |> unique()
+  ]),
+  angle_col = 90,
+  use_raster = F,
+  border_color = NA,
+  main = 'Hierarchically-Clustered TMEM184B IPMS TSCs\nZ-scored, Low SAINT Score, Low Contamination Cluster'
+)
+

+
+
+

Low hSAINT Score, High Contamination Modules Heatmap +(Background)

+
pheatmap::pheatmap(
+  IPMS_z_ordered[rownames(IPMS_z_ordered) %in% c(protein_module_df |> 
+                   dplyr::filter(Protein == 'TMEM184B' |
+                                   grepl(classification,
+                                         pattern = 'Low SAINT, High')) |> 
+                   dplyr::pull(Protein)),],
+  color = grDevices::colorRampPalette(c('navy',
+                                        'white',
+                                        'darkgoldenrod'))((200)),
+  cluster_rows = T,
+  cluster_cols = F,
+  show_rownames = F,
+  annotation_row = data.frame(
+    
+    Module = clusters_opt[
+      names(clusters_opt) %in% 
+        
+        c(protein_module_df |> 
+            dplyr::filter(Protein == 'TMEM184B' |
+                            grepl(classification,
+                                  pattern = 'Low SAINT, High')) |> 
+            dplyr::pull(Protein)
+          )
+      ]
+    ),
+  
+  annotation_colors = list(Module = annotation_colors[[1]][
+    clusters_opt[
+      names(clusters_opt) %in% c(protein_module_df |> 
+                   dplyr::filter(Protein == 'TMEM184B' |
+                                   grepl(classification,
+                                         pattern = 'Low SAINT, High')) |> 
+                   dplyr::pull(Protein))
+    ] |> as.numeric() |> unique()
+  ]),
+  angle_col = 90,
+  use_raster = F,
+  border_color = NA,
+  main = 'Hierarchically-Clustered TMEM184B IPMS TSCs\nZ-scored, Low SAINT Score, High Contamination Cluster'
+)
+

+
+
+
+
+

Results Tables

+
+

Clustering Results Table by Protein (Modules 9, 12)

+
DT::datatable(
+  protein_module_df |> 
+    dplyr::full_join(global_enrichment_results_df) |> 
+    dplyr::full_join(submodule_enrichment_results_df) |> 
+    dplyr::rename(hSAINT_score = score) |> 
+    dplyr::rename(module_size = n) %>%
+    dplyr::rename_with(.cols = c(dplyr::contains("mean_"),
+                                 dplyr::contains("median_"),
+                                 dplyr::contains("_perc"),
+                                 score_mass), ~paste0("module_", .x)) |> 
+    dplyr::filter(!is.na(GO_Term)) |> 
+    dplyr::filter(!is.na(GO_Category)) |> 
+    dplyr::filter( (grepl(module,
+                          pattern = '12') &
+                      hSAINT_score > 0.5) |
+                     (grepl(classification,
+                            pattern = 'High SAINT, Low') &
+                     hSAINT_score > 0.5)
+                   )
+)
+
+ +

This (complete) table exists as hClust_protein_level_GO_GSEA.rda +in Endo/analysis/R/hClust/

+
+
+

Clustering Results Table by Module

+
DT::datatable(
+  global_enrichment_results_df |> 
+  dplyr::full_join(submodule_enrichment_results_df) |> 
+  dplyr::full_join(module_scores) |> 
+  dplyr::rename(module_size = n) %>%
+  dplyr::rename_with(.cols = c(dplyr::contains("mean_"),
+                                 dplyr::contains("median_"),
+                                 dplyr::contains("_perc"),
+                                 score_mass), ~paste0("module_", .x)) |> 
+  dplyr::filter(!is.na(GO_Term)) |> 
+  dplyr::filter(!is.na(GO_Category))
+)
+
+ +

This table exists as hClust_module_level_GO_GSEA.rda +in Endo/analysis/R/hClust/

+
+
+
+

Summary

+

Module 9 has the highest volume of +hSAINT signal across modules, with some interesting +potential candidates

+

As noted above, it is not surprising Module 12, +TMEM’s module, was not enriched, containing a large amount +of proteins with low hSAINT scores, likely due to +TMEM’s nature as a transmembrane protein, exacerbated by +its localization to multiple organellar compartments

+

A more in-depth analysis of proteins across modules where there is +greater probability of interacting with TMEM is encouraged– +collecting groups of proteins, related by localization or function, is +useful for inferring data validity; however, data coherence +should not negate investigating proteins with high +hSAINT scores or are known/likely to interact with +TMEM or related processes

+

Collectively, it appears the hierarchical SAINT model +has done a fair job at delineating possible signal from noise

+
+ + + + +
+ + + + + + + + + + + + + + + diff --git a/analysis/HierarchicalClustering.md b/analysis/HierarchicalClustering.md new file mode 100644 index 0000000..1aa7907 --- /dev/null +++ b/analysis/HierarchicalClustering.md @@ -0,0 +1,1105 @@ +# Overview + +This markdown was developed to formally document an ad-hoc cluster +analysis for TMEM184B immunoprecipitation tandem mass-spec experiments +for **Dr. Martha Bhattacharya’s lab** at the **University of Arizona** +for the endolysosomal paper under revision (2025-02-28) for the Journal +of Cell Science + +This vignette documents the `Hierarchical Clustering` of this IPMS data + +## Attach Packages + +``` r +# data/wrangling +library(tidyverse) +library(Endo) +library(DT) + +# clustering +library(pheatmap) +library(dynamicTreeCut) + +# GSEA/viz +library(clusterProfiler) +library(org.Hs.eg.db) +library(ReactomePA) +library(AnnotationDbi) + +# progress bar +library(progressr) + +# parallelization +library(future) +library(future.apply) + +# network viz +library(visNetwork) + +options(stringsAsFactors = FALSE) +``` + +## Load the Pre-Computed Data + +``` r +# Load precomputed objects from preprocessing +load("../../analysis/R/precompute_objects.rda") +``` + +# Hierarchical Clustering + +## Raw Counts + +``` r +d <- stats::as.dist(1 - stats::cor(t(IPMS_mat), use = 'pairwise.complete.obs')) +hc <- stats::hclust(d, method = "average") +``` + +``` r +sizes <- c(20, 30, 50, 80, 100) +results <- list() + +for (m in sizes) { + cl <- dynamicTreeCut::cutreeDynamic( + dendro = hc, + distM = as.matrix(d), + deepSplit = 1, + pamRespectsDendro = TRUE, + minClusterSize = m + ) + results[[as.character(m)]] <- cl + cat("minClusterSize =", m, "→", length(table(cl)), "modules\n") +} +``` + + ## ..cutHeight not given, setting it to 0.976 ===> 99% of the (truncated) height range in dendro. + ## ..done. + ## minClusterSize = 20 → 50 modules + ## ..cutHeight not given, setting it to 0.976 ===> 99% of the (truncated) height range in dendro. + ## ..done. + ## minClusterSize = 30 → 42 modules + ## ..cutHeight not given, setting it to 0.976 ===> 99% of the (truncated) height range in dendro. + ## ..done. + ## minClusterSize = 50 → 28 modules + ## ..cutHeight not given, setting it to 0.976 ===> 99% of the (truncated) height range in dendro. + ## ..done. + ## minClusterSize = 80 → 21 modules + ## ..cutHeight not given, setting it to 0.976 ===> 99% of the (truncated) height range in dendro. + ## ..done. + ## minClusterSize = 100 → 17 modules + +``` r +rm(list = c('cl', 'm', 'sizes', 'results')) +``` + +- `optimal minimum cluster size` is **50** (elbow; diminishing returns + after) + +- `# modules` = **28** + +``` r +clusters_opt <- dynamicTreeCut::cutreeDynamic(hc, + distM = as.matrix(d), + deepSplit = 1, + pamRespectsDendro = T, + minClusterSize = 50) +``` + + ## ..cutHeight not given, setting it to 0.976 ===> 99% of the (truncated) height range in dendro. + ## ..done. + +``` r +names(clusters_opt) <- rownames(IPMS_mat) +``` + +``` r +module_ids <- sort(unique(clusters_opt)) +module_colors <- stats::setNames( + grDevices::rainbow(length(module_ids), s = 0.6, v = 0.9), + module_ids +) + +annotation_colors <- list(Module = module_colors) +``` + +Sparse data is too difficult to discern from a global view + +-\> Z transform the raw counts **after clustering** to plot + +- the cluster object preserves the ordering of the proteins (dendrogram) + + - Z-score therefore only changes the values of the columns and rows to + better visualize the sparse data + +``` r +IPMS_z <- t(scale(t(IPMS_mat))) +colnames(IPMS_z) <- colnames(IPMS_mat) + +condition <- sub('_.*', '', colnames(IPMS_mat)) + +col_order <- unlist( + lapply( + split(seq_along(condition), condition), function(idx) { + submat <- IPMS_z[, idx, drop = FALSE] + hc_sub <- stats::hclust(stats::dist(t(submat)), method = 'ward.D2') + idx[hc_sub$order] + })) + +IPMS_z_ordered <- IPMS_z[, col_order] + +rm(list = c('condition', 'col_order')) +``` + +# Module Descriptions + +In order to assess this method and also to compare/validate the +hierarchical SAINT model, describing modules/clusters based on the +proteins in them is useful + +Similarly localized proteins could be expected to have similar +abundances to each other as well as the target (**TMEM184B**) + +Similar functioning proteins could also be expected to have similar +abundances to each other as well as the target + +Describing these modules/clusters attained by `hierarchical clustering` +via gene ontology gene set enrichment analysis is, therefore, useful and +standard practice + +Furthermore, we may be able to assemble a +`protein-protein interaction network` from the results + +I outline this below– the methods and results for acquiring these +modules, their constituent proteins, and how we make sense of them– +visualized as `dotplots`, `heatmaps`, and a `network`. I include +searchable tables printed last + +## Extract Modules + +Extract gene names from the hierarchically-clustered modules and convert +to a dataframe + +``` r +modules_list_idx <- lapply(sort(unique(clusters_opt)), function(m) { + which(clusters_opt == m) +}) +names(modules_list_idx) <- sort(unique(clusters_opt)) + +modules_list <- lapply(sort(unique(clusters_opt)), function(m) { + rownames(IPMS_mat)[clusters_opt == m] +}) +names(modules_list) <- sort(unique(clusters_opt)) |> as.character() + +modules_df <- modules_list %>% + purrr::map_df(., .f = as.data.frame, .id = 'module') |> + dplyr::rename(Protein = 2) |> + dplyr::relocate(Protein, .before = 1) +``` + +## Join Modules to hSAINT + +Include score summaries (mean and median posterior probability scores, +ranks, and literature or CRAPome contamination rate) + +``` r +module_scores <- modules_df |> + dplyr::inner_join( + hSAINT + ) |> + dplyr::arrange(desc(score)) |> + dplyr::group_by(module) |> + dplyr::summarize(mean_score = mean(score), + median_score = stats::median(score), + score_mass = sum(score), + n = dplyr::n(), + lit_flag_perc = (sum(lit_flag)/n)*100, + CRAPome_flag_perc = (sum(CRAPome_flag)/n)*100) |> + dplyr::arrange(desc(mean_score)) |> + dplyr::mutate(mean_rank = dplyr::row_number()) |> + dplyr::arrange(desc(median_score)) |> + dplyr::mutate(median_rank = dplyr::row_number()) |> + dplyr::arrange(desc(score_mass)) |> + dplyr::mutate(total_rank = dplyr::row_number()) |> + dplyr::arrange(mean_rank) +``` + +Join the module scores to the modules (per protein) dataframe and the +`hSAINT` results for classification + +``` r +protein_module_df <- modules_df |> + dplyr::inner_join( + hSAINT + ) |> + dplyr::full_join( + module_scores |> + dplyr::ungroup() + ) |> + dplyr::mutate( + SAINT_classification = dplyr::case_when(mean_rank <= 5 ~ 'High SAINT', + mean_rank > 5 ~ 'Low SAINT' + ), + + contams_classification = dplyr::case_when(lit_flag_perc > 5 ~'High Contaminants', + TRUE ~ 'Low Contaminants')) |> + dplyr::mutate(classification = paste0(SAINT_classification, + ", ", + contams_classification)) |> + dplyr::select(-SAINT_classification, + -contams_classification) |> + dplyr::relocate(classification, .after = module) +``` + +## Module GO GSEA + +``` r +EntrezIDs <- AnnotationDbi::select( + org.Hs.eg.db, + keys = modules_df |> dplyr::distinct(Protein) |> dplyr::pull(), + columns = 'ENTREZID', + keytype = 'SYMBOL' +) |> + dplyr::filter(!is.na(2)) + +EntrezMap <- stats::setNames(EntrezIDs$ENTREZID, EntrezIDs$SYMBOL) +``` + +Create enrichment functions + +``` r +background_genes <- EntrezIDs |> # use the EntrezIDs mapped from the modules_df + dplyr::pull(ENTREZID) + +get_GO <- function(genes, ontology, background_genes) { + clusterProfiler::enrichGO( + gene = genes, + universe = background_genes, + OrgDb = org.Hs.eg.db, + keyType = 'ENTREZID', + ont = ontology, + pAdjustMethod = 'BH', + qvalueCutoff = 0.05, + readable = T + ) +} + +get_kegg <- function(genes, EntrezMap) { + entrez <- EntrezMap[genes] + entrez <- entrez[!is.na(entrez)] + if (length(entrez)==0) { + return(NULL) + } + clusterProfiler::enrichKEGG( + gene = entrez, + organism = 'hsa', + pAdjustMethod = 'BH', + qvalueCutoff = 0.05 + ) +} + +get_reactome <- function(genes, EntrezMap) { + entrez <- EntrezMap[genes] + entrez <- entrez[!is.na(entrez)] + if (length(entrez)==0) { + return(NULL) + } + ReactomePA::enrichPathway( + gene = entrez, + organism = 'human', + pAdjustMethod = 'BH', + qvalueCutoff = 0.05, + readable = T + ) +} +``` + +Run the functions on the dataset in parallel (multiple CPUs) + +``` r +future::plan(future::multisession, workers = 6) +progressr::handlers(progressr::handler_txtprogressbar()) + +module_ids <- sort(unique(modules_df$module)) + +progressr::with_progress({ + p <- progressr::progressor(steps = length(module_ids)) + + enrichment_results <- future.apply::future_lapply( + module_ids, + function(mod) { + + # tick the progress bar + p(sprintf("Module %s", mod)) + + genes <- modules_df |> + dplyr::filter(module == mod) |> + dplyr::inner_join(EntrezIDs |> dplyr::rename(Protein = 1)) |> + dplyr::distinct(ENTREZID) |> + dplyr::pull(ENTREZID) + + list( + module = mod, + go_bp = get_GO(genes, ontology = "BP", universe = background_genes), + go_cc = get_GO(genes, ontology = "CC", universe = background_genes), + go_mf = get_GO(genes, ontology = "MF", universe = background_genes) + # add KEGG/Reactome later + ) + } + ) +}) + +names(enrichment_results) <- module_ids +rm(p) +``` + +Gather enrichment results + +``` r +global_enrichment_results_df <- lapply(enrichment_results, function(x) { + x %>% + purrr::map_df(., .f = as.data.frame, .id = 'GO_Category') |> + dplyr::rename(module = 2, GO_Term = Description, GO_ID = ID) |> + tidyr::fill(module, .direction = 'down') |> + dplyr::filter(!is.na(GO_ID)) |> + dplyr::mutate(GO_Category = gsub(GO_Category, + pattern = 'go_', + replacement = '') %>% + stringr::str_to_upper(.)) |> + tibble::remove_rownames() |> + dplyr::distinct() + } +) %>% + purrr::map_df(., .f = as.data.frame) +``` + +The “`High SAINT score`” modules were among the modules with significant +enrichment of all the proteins detected in the experiments + +However, many modules did not have significant enrichment relative to +the entire set of proteins in the experiments, *including TMEM’s*. Using +`hierarchical SAINT score` as a explanatory tool, this is likely because +the modules without GO enrichment contained a large mix of proteins +scored as “noise” (low posterior probability of being a true interactor) +and signal (high posterior probability of being a true interactor) + +In order to better make sense of the clustering– describing the +**“un-enriched” modules**– “sub-module” GO GSEA was conducted + +## SubModule GO GSEA + +Sub-module enrichment analysis on modules without any enrichment (to +describe themes of proteins within these modules) + +``` r +"%notin%" <- Negate("%in%") + +# index the modules without significant enrichment to describe sub-modules +tier2_modules <- c(seq(1:28) |> as.character())[ + which(c(seq(1:28) |> as.character()) %notin% c(global_enrichment_results_df |> + dplyr::distinct(module) |> + dplyr::pull()) == T)] +``` + +``` r +# function to gather sub-modules for each module that did not return enrichment results +submodules_from_hclust <- function( + hc_object, + modules_df, + n_sub_per_module = 3L +) { + # number of existing modules + n_modules <- length(unique(modules_df$module)) + + # total number of clusters in the finer cut + K_big <- n_modules * n_sub_per_module + + # finer cut of the same tree + sub_labels_all <- stats::cutree(hc_object, k = K_big) + + # assemble into a tibble + sub_df <- tibble::tibble( + Protein = names(sub_labels_all), + sub_label = sub_labels_all + ) + + # attach existing module labels + sub_df <- dplyr::inner_join( + sub_df, + modules_df, + by = "Protein" + ) + + # define submodules as (module, sub_label) pairs + submodules_df <- sub_df %>% + dplyr::group_by(module, sub_label) %>% + dplyr::mutate( + submodule = paste0("sub", dplyr::cur_group_id()) + ) %>% + dplyr::ungroup() %>% + dplyr::select(Protein, module, submodule) + + return(submodules_df) +} +``` + +``` r +# get the sub-modules +submodules_df <- submodules_from_hclust( + hc_object = hc, + modules_df = dplyr::filter(modules_df, module %in% tier2_modules), + n_sub_per_module = 3L +) +``` + +``` r +# perform GO GSEA on the submodules + +submodule_enrichment_results <- submodules_df %>% + dplyr::inner_join(EntrezIDs |> dplyr::rename(Protein = 1)) |> + dplyr::group_by(module, submodule) %>% + dplyr::group_modify(~ { + + # submodule genes + genes <- .x$ENTREZID + + background_genes <- unique(submodules_df$ENTREZID[submodules_df$module == .y$module]) + + tibble::tibble( + BP = list(get_GO(genes, ontology = "BP", background_genes = background_genes)), + CC = list(get_GO(genes, ontology = "CC", background_genes = background_genes)), + MF = list(get_GO(genes, ontology = "MF", background_genes = background_genes)) + ) + }) %>% + dplyr::ungroup() +rm(background_genes) +``` + +``` r +# wrangle the sub-module GO GSEA results +submodule_enrichment_results_df <- submodule_enrichment_results |> + tidyr::pivot_longer( + cols = c('BP', 'CC', 'MF'), + names_to = "GO_Category", + values_to = "values" + ) %>% + dplyr::mutate( + GO_tbl = purrr::map(values, ~ as.data.frame(.x)) + ) |> + tidyr::unnest(GO_tbl) |> + dplyr::rename( + GO_ID = ID, + GO_Term = Description + ) |> + dplyr::select( + module, + submodule, + GO_Category, + GO_Term, + GO_ID, + dplyr::everything(), + -values + ) +``` + +``` r +modules_with_submodule_enrichment <- tier2_modules +modules_without_submodule_enrichment <- generics::setdiff(modules_with_submodule_enrichment, + submodule_enrichment_results_df |> + dplyr::distinct(module) |> + dplyr::pull()) +modules_with_submodule_enrichment <- submodule_enrichment_results_df |> + dplyr::distinct(module) |> + dplyr::pull() +rm(tier2_modules) +``` + +## Construct the Module GO (similarity) Network + +``` r +go_network <- global_enrichment_results_df |> + dplyr::full_join(submodule_enrichment_results_df) |> + dplyr::full_join(module_scores) |> + dplyr::mutate(Tier = dplyr::case_when(is.na(submodule) ~ 'Tier1',TRUE~'Tier2'), + node_to = GO_Term, + node_from = paste0( + + "Module", module, + ifelse(is.na(submodule), + '', + paste0('_SubModule', + submodule %>% + stringr::str_extract(., pattern = '\\d+')))), + weight = -log10(p.adjust)) |> + dplyr::filter(weight > -log10(0.05)) + +g <- igraph::graph_from_data_frame( + d = go_network[, c("node_from", "node_to", + "weight", "geneID", "p.adjust", "FoldEnrichment")], + directed = TRUE, + vertices = + dplyr::bind_rows( + go_network |> dplyr::transmute(name = node_from, score_mass, submodule), + go_network |> dplyr::transmute(name = node_to, + score_mass = NA_real_, + submodule = NA_character_) + ) |> + dplyr::distinct(name, .keep_all = TRUE) +) +``` + +# Plots + +## Module Dot Plot + +``` r +global_enrichment_results_df |> + dplyr::full_join(submodule_enrichment_results_df) |> + dplyr::full_join(module_scores) |> + dplyr::filter(!is.na(GO_Term)) |> + dplyr::filter(!is.na(GO_Category)) |> + dplyr::mutate( + Tier = dplyr::case_when( + is.na(submodule) ~ paste0('Mod. ', module), + TRUE ~ paste0('Mod. ', module, + ':\nsubmod. ', + gsub(submodule, pattern = 'sub', replacement = '')) + ), + is_enriched = dplyr::case_when(is.na(submodule) ~ 'Enriched', + TRUE ~ 'Not Enriched') + ) |> + dplyr::group_by(module) |> + dplyr::slice_head(n = 3) |> + dplyr::arrange(p.adjust) |> + ggplot2::ggplot( + ggplot2::aes( + x = -log10(p.adjust), + y = forcats::fct_rev(forcats::fct_infreq(GO_Term, mean_score)), + color = GO_Category, + size = FoldEnrichment, + alpha = mean_score, + shape = is_enriched, + ) + ) + + ggplot2::geom_vline(xintercept = -log10(0.05), + linetype = 'dashed', + color = 'firebrick', + size = 1, + alpha = 0.5) + + ggplot2::geom_point() + + ggplot2::facet_grid(Tier ~ ., scales = 'free_y') + + ggplot2::labs(title = 'GO Enrichment Across\nModules and SubModules\nin TMEM184B IPMS\nHierarchical Clustering', + x = '-log10(Adj.P-val)', + y = 'GO Term', + alpha = 'Module\nMean\nHierarchical\nSAINT\nScore', + size = 'FoldEnrichment', + color = 'GO Category', + shape = 'IsEnrichedModule') + + ggplot2::theme_bw() + + ggplot2::theme(title = ggplot2::element_text(face = 'bold', size = 12), + axis.text = ggplot2::element_text(face = 'bold', size = 10), + strip.text = ggplot2::element_text(face = 'bold', size = 7), + strip.background = ggplot2::element_rect(color = 'black', + fill = 'white')) + + ggplot2::guides(#shape = ggplot2::guide_legend(override.aes = list(size = 4, shape = c(19, 17))), + color = ggplot2::guide_legend(override.aes = list(size = 4, + shape = 19, + color = c('navy', + 'darkgoldenrod3', + 'firebrick'))), + alpha = ggplot2::guide_legend(override.aes = list(size = 5, + shape = 19)), + size = ggplot2::guide_legend(override.aes = list(shape = 19)) + # size = 'none' + ) + + ggplot2::scale_fill_manual(values = c('navy', + 'darkgoldenrod3', + 'firebrick')) +``` + +![](../../analysis/R/hClust/RawModuleSubModuleDotPlotOfGOTerms-1.png) + +## Module GO (similarity) Network + +``` r +coords <- igraph::layout_with_fr(g) +coords <- coords*200 + +visNetwork::visNetwork( + nodes = data.frame( + id = igraph::V(g)$name, + label = igraph::V(g)$name, + size = scales::rescale( + (\(x){ x[is.na(x)] <- 0; x })(igraph::V(g)$score_mass), + to = c(20, 60) + ), + x = coords[,1], + y = coords[,2], + group = base::interaction( + { + x <- igraph::V(g)$module + if (length(x) == 0) rep("", igraph::vcount(g)) else x + }, + { + x <- igraph::V(g)$submodule + if (length(x) == 0) { + rep("", igraph::vcount(g)) + } else { + ifelse(nzchar(x), x, igraph::V(g)$module) + } + }, + drop = TRUE, + lex.order = TRUE + ), + title = sapply( + igraph::V(g), + function(v) { + incoming_edges <- igraph::incident(g, v, mode = "all") + incoming_edges <- incoming_edges[ + igraph::ends(g, incoming_edges)[,2] == igraph::V(g)[v]$name + ] + raw_gene_strings <- igraph::E(g)[incoming_edges]$geneID + raw_gene_strings <- unlist(raw_gene_strings, use.names = FALSE) + + genes <- raw_gene_strings |> + as.character() |> + (\(x) x[!is.na(x)])() |> + strsplit("\\\\") |> + unlist() |> + (\(x) x[nzchar(x)])() |> + unique() + + term <- igraph::V(g)[v]$name + saint_value <- igraph::V(g)[v]$score_mass + if (is.na(saint_value)) saint_value <- 0 + is_module <- nzchar(igraph::V(g)[v]$submodule) + + if (is_module) { + paste0( + "
", + "", term, "
", + "SAINT mass: ", signif(saint_value, 3), "
", + if (length(genes) > 0) + paste0("Proteins:
", paste(genes, collapse = "
")) + else "", + "
" + ) + } else { + paste0( + "
", + "", term, "
", + if (length(genes) > 0) + paste0("Proteins:
", paste(genes, collapse = "
")) + else "", + "
" + ) + } + } + ) + ), + edges = data.frame( + from = igraph::as_data_frame(g, "edges")$from, + to = igraph::as_data_frame(g, "edges")$to, + value = igraph::E(g)$weight, + title = vapply( + seq_along(igraph::E(g)), + function(i) { + paste0( + "
", + "p.adjust: ", signif(igraph::E(g)$p.adjust[i], 3), "
", + "weight: ", signif(igraph::E(g)$weight[i], 3), + "
" + ) + }, + FUN.VALUE = character(1) + ) + ) +) |> + # visNetwork::visNodes(color = NULL) |> + visNetwork::visEdges(arrows = "to") |> + visNetwork::visNodes( + scaling = list(min = 5, max = 40), + shadow = TRUE + ) |> + visNetwork::visLayout(randomSeed = 1) |> + visNetwork::visPhysics(enabled = FALSE) |> + visNetwork::visInteraction( + dragView = TRUE, + zoomView = TRUE + ) |> + visNetwork::visOptions(nodesIdSelection = TRUE) +``` + +
+ + +``` r +rm(coords) +``` + +Modules with `enriched` GO Terms are `blue` + +Here, only the enriched sub-modules were included: the modules with that +did not reach significant enrichment broadly, but did reach significant +enrichment when that module was broken up. These `Enriched` +**sub-modules** are **non-blue** + +Modules **2**, **9**, **15**, and **26** are of interest, having terms +in line with what might be expected from a `TMEM` experiment, and having +robust `Hierarchical SAINT` scores (hover over Module/SubModule nodes) + +`TMEM184B`, itself, was in Module **12**, which did *NOT* return any +enriched GO Terms, likely due to its membrane protein nature (and +putative localization, where it may progress through the endomembrane +system and collide with a wide variety of proteins). This is common +among TM protein IPMS + +## Heatmaps + +### Raw Counts Heatmap + +``` r +pheatmap::pheatmap( + IPMS_z_ordered, + color = grDevices::colorRampPalette(c('navy', + 'white', + 'darkgoldenrod'))((200)), + cluster_rows = hc, + cluster_cols = F, + show_rownames = F, + annotation_row = data.frame(Module = clusters_opt), + annotation_colors = annotation_colors, + angle_col = 90, + use_raster = F, + border_color = NA, + main = 'Hierarchically-Clustered\nTMEM184B IPMS TSCs, Z-scored' +) +``` + +![](../../analysis/R/hClust/RawZHeatmapAll-1.png) + +### Raw Counts Heatmap (TMEM Module) + +``` r +pheatmap::pheatmap( + IPMS_z_ordered[rownames(IPMS_z_ordered) %in% + + c(protein_module_df |> + dplyr::full_join(global_enrichment_results_df) |> + dplyr::full_join(submodule_enrichment_results_df) |> + dplyr::select(-c(GO_Category:Count)) |> + dplyr::distinct() |> + dplyr::filter( + module == protein_module_df |> + dplyr::full_join(global_enrichment_results_df) |> + dplyr::full_join(submodule_enrichment_results_df) |> + dplyr::select(-c(GO_Category:Count)) |> + dplyr::filter(Protein == 'TMEM184B') |> + dplyr::pull(module) + ) |> + dplyr::pull(Protein) + ), + ], + color = grDevices::colorRampPalette(c('navy', + 'white', + 'darkgoldenrod'))((200)), + cluster_rows = T, + cluster_cols = F, + show_rownames = T, + fontsize_row = 5, + annotation_row = data.frame( + + Module12_SubModule = clusters_opt[ + names(clusters_opt) %in% + + c(protein_module_df |> + dplyr::full_join(global_enrichment_results_df) |> + dplyr::full_join(submodule_enrichment_results_df) |> + dplyr::select(-c(GO_Category:Count)) |> + dplyr::distinct() |> + dplyr::filter( + module == protein_module_df |> + dplyr::full_join(global_enrichment_results_df) |> + dplyr::full_join(submodule_enrichment_results_df) |> + dplyr::select(-c(GO_Category:Count)) |> + dplyr::filter(Protein == 'TMEM184B') |> + dplyr::pull(module) + ) |> dplyr::filter(!is.na(Protein)) |> + dplyr::pull(Protein) + ) + ] + ), + annotation_colors = list(Module12_SubModule = 'darkorange2'), + angle_col = 90, + use_raster = F, + border_color = NA, + main = 'Hierarchically-Clustered TMEM184B IPMS TSCs\nZ-scored, TMEM Module' +) +``` + +![](../../analysis/R/hClust/RawZHeatmapTMEMModule-1.png) + +`TMEM`’s module (**12**) is not siginficantly enriched above the +proteins detected by the mass spectrometer. + +Submodules within **module 12** are not enriched + +### High hSAINT Score, Low Contamination Module Heatmap (Signal) + +``` r +pheatmap::pheatmap( + IPMS_z_ordered[rownames(IPMS_z_ordered) %in% c(protein_module_df |> + dplyr::filter(Protein == 'TMEM184B' | + grepl(classification, + pattern = 'High SAINT, Low')) |> + dplyr::pull(Protein)),], + color = grDevices::colorRampPalette(c('navy', + 'white', + 'darkgoldenrod'))((200)), + cluster_rows = T, + cluster_cols = F, + show_rownames = F, + annotation_row = data.frame( + + Module = clusters_opt[ + names(clusters_opt) %in% + + c(protein_module_df |> + dplyr::filter(Protein == 'TMEM184B' | + grepl(classification, + pattern = 'High SAINT, Low')) |> + dplyr::pull(Protein) + ) + ] + ), + + annotation_colors = list(Module = c('firebrick', 'navy')), + angle_col = 90, + use_raster = F, + border_color = NA, + cutree_rows = 3, + main = 'Hierarchically-Clustered TMEM184B IPMS TSCs\nZ-scored, High SAINT Score, Low Contamination Cluster + TMEM' +) +``` + +![](../../analysis/R/hClust/RawZHeatmapSAINTSignalModule-1.png) + +### High hSAINT Score, High Contamination Modules Heatmap (Sticky/Noise) + +``` r +pheatmap::pheatmap( + IPMS_z_ordered[rownames(IPMS_z_ordered) %in% c(protein_module_df |> + dplyr::filter(Protein == 'TMEM184B' | + grepl(classification, + pattern = 'High SAINT, High')) |> + dplyr::pull(Protein)),], + color = grDevices::colorRampPalette(c('navy', + 'white', + 'darkgoldenrod'))((200)), + cluster_rows = T, + cluster_cols = F, + show_rownames = F, + annotation_row = data.frame( + + Module = clusters_opt[ + names(clusters_opt) %in% + + c(protein_module_df |> + dplyr::filter(Protein == 'TMEM184B' | + grepl(classification, + pattern = 'High SAINT, High')) |> + dplyr::pull(Protein) + ) + ] + ), + + annotation_colors = list(Module = c('darkorange', + 'darkviolet', + 'lightpink', + 'skyblue')), + angle_col = 90, + use_raster = F, + border_color = NA, + cutree_rows = 2, + main = 'Hierarchically-Clustered TMEM184B IPMS TSCs\nZ-scored, High SAINT Score, High Contamination Clusters' +) +``` + +![](../../analysis/R/hClust/RawZHeatmapSAINTHighHighModules-1.png) + +### Low hSAINT Score, Low Contamination Modules Heatmap (Low Abundance) + +``` r +pheatmap::pheatmap( + IPMS_z_ordered[rownames(IPMS_z_ordered) %in% c(protein_module_df |> + dplyr::filter(Protein == 'TMEM184B' | + grepl(classification, + pattern = 'Low SAINT, Low')) |> + dplyr::pull(Protein)),], + color = grDevices::colorRampPalette(c('navy', + 'white', + 'darkgoldenrod'))((200)), + cluster_rows = T, + cluster_cols = F, + show_rownames = F, + annotation_row = data.frame( + + Module = clusters_opt[ + names(clusters_opt) %in% + + c(protein_module_df |> + dplyr::filter(Protein == 'TMEM184B' | + grepl(classification, + pattern = 'Low SAINT, Low')) |> + dplyr::pull(Protein) + ) + ] + ), + + annotation_colors = list(Module = annotation_colors[[1]][ + clusters_opt[ + names(clusters_opt) %in% c(protein_module_df |> + dplyr::filter(Protein == 'TMEM184B' | + grepl(classification, + pattern = 'Low SAINT, Low')) |> + dplyr::pull(Protein)) + ] |> as.numeric() |> unique() + ]), + angle_col = 90, + use_raster = F, + border_color = NA, + main = 'Hierarchically-Clustered TMEM184B IPMS TSCs\nZ-scored, Low SAINT Score, Low Contamination Cluster' +) +``` + +![](../../analysis/R/hClust/RawZHeatmapSAINTLowLowModules-1.png) + +### Low hSAINT Score, High Contamination Modules Heatmap (Background) + +``` r +pheatmap::pheatmap( + IPMS_z_ordered[rownames(IPMS_z_ordered) %in% c(protein_module_df |> + dplyr::filter(Protein == 'TMEM184B' | + grepl(classification, + pattern = 'Low SAINT, High')) |> + dplyr::pull(Protein)),], + color = grDevices::colorRampPalette(c('navy', + 'white', + 'darkgoldenrod'))((200)), + cluster_rows = T, + cluster_cols = F, + show_rownames = F, + annotation_row = data.frame( + + Module = clusters_opt[ + names(clusters_opt) %in% + + c(protein_module_df |> + dplyr::filter(Protein == 'TMEM184B' | + grepl(classification, + pattern = 'Low SAINT, High')) |> + dplyr::pull(Protein) + ) + ] + ), + + annotation_colors = list(Module = annotation_colors[[1]][ + clusters_opt[ + names(clusters_opt) %in% c(protein_module_df |> + dplyr::filter(Protein == 'TMEM184B' | + grepl(classification, + pattern = 'Low SAINT, High')) |> + dplyr::pull(Protein)) + ] |> as.numeric() |> unique() + ]), + angle_col = 90, + use_raster = F, + border_color = NA, + main = 'Hierarchically-Clustered TMEM184B IPMS TSCs\nZ-scored, Low SAINT Score, High Contamination Cluster' +) +``` + +![](../../analysis/R/hClust/RawZHeatmapSAINTContamModules-1.png) + +# Results Tables + +## Clustering Results Table by Protein (Modules 9, 12) + +``` r +DT::datatable( + protein_module_df |> + dplyr::full_join(global_enrichment_results_df) |> + dplyr::full_join(submodule_enrichment_results_df) |> + dplyr::rename(hSAINT_score = score) |> + dplyr::rename(module_size = n) %>% + dplyr::rename_with(.cols = c(dplyr::contains("mean_"), + dplyr::contains("median_"), + dplyr::contains("_perc"), + score_mass), ~paste0("module_", .x)) |> + dplyr::filter(!is.na(GO_Term)) |> + dplyr::filter(!is.na(GO_Category)) |> + dplyr::filter( (grepl(module, + pattern = '12') & + hSAINT_score > 0.5) | + (grepl(classification, + pattern = 'High SAINT, Low') & + hSAINT_score > 0.5) + ) +) +``` + +
+ + +This (complete) table exists as +[hClust_protein_level_GO_GSEA.rda](../../analysis/R/hClust/hClust_protein_level_GO_GSEA.rda) +in [Endo/analysis/R/hClust/](../../analysis/R/hClust) + +## Clustering Results Table by Module + +``` r +DT::datatable( + global_enrichment_results_df |> + dplyr::full_join(submodule_enrichment_results_df) |> + dplyr::full_join(module_scores) |> + dplyr::rename(module_size = n) %>% + dplyr::rename_with(.cols = c(dplyr::contains("mean_"), + dplyr::contains("median_"), + dplyr::contains("_perc"), + score_mass), ~paste0("module_", .x)) |> + dplyr::filter(!is.na(GO_Term)) |> + dplyr::filter(!is.na(GO_Category)) +) +``` + +
+ + +This table exists as +[hClust_module_level_GO_GSEA.rda](../../analysis/R/hClust/hClust_module_level_GO_GSEA.rda) +in [Endo/analysis/R/hClust/](../../analysis/R/hClust) + +\#Summary + +**Module 9** has the highest volume of `hSAINT` signal across modules, +with some interesting potential candidates + +As noted above, it is not surprising **Module 12**, `TMEM`’s module, was +not enriched, containing a large amount of proteins with low `hSAINT` +scores, likely due to `TMEM`’s nature as a transmembrane protein, +exacerbated by its localization to multiple organellar compartments + +A more in-depth analysis of proteins across modules where there is +greater probability of interacting with `TMEM` is encouraged– collecting +groups of proteins, related by localization or function, is useful for +inferring data validity; however, data coherence **should not negate** +investigating proteins with high `hSAINT` scores or are known/likely to +interact with `TMEM` or related processes + +Collectively, it appears the `hierarchical SAINT` model has done a fair +job at delineating possible signal from noise diff --git a/analysis/R/hClust/RawModuleSubModuleDotPlotOfGOTerms-1.pdf b/analysis/R/hClust/RawModuleSubModuleDotPlotOfGOTerms-1.pdf index 1ce7ed4..706cd5e 100644 Binary files a/analysis/R/hClust/RawModuleSubModuleDotPlotOfGOTerms-1.pdf and b/analysis/R/hClust/RawModuleSubModuleDotPlotOfGOTerms-1.pdf differ diff --git a/analysis/R/hClust/RawModuleSubModuleDotPlotOfGOTerms-1.png b/analysis/R/hClust/RawModuleSubModuleDotPlotOfGOTerms-1.png index 5fe7afe..763a30f 100644 Binary files a/analysis/R/hClust/RawModuleSubModuleDotPlotOfGOTerms-1.png and b/analysis/R/hClust/RawModuleSubModuleDotPlotOfGOTerms-1.png differ diff --git a/analysis/R/hClust/RawZHeatmapAll-1.pdf b/analysis/R/hClust/RawZHeatmapAll-1.pdf index 37c5ac5..6d0fe09 100644 Binary files a/analysis/R/hClust/RawZHeatmapAll-1.pdf and b/analysis/R/hClust/RawZHeatmapAll-1.pdf differ diff --git a/analysis/R/hClust/RawZHeatmapSAINTContamModules-1.pdf b/analysis/R/hClust/RawZHeatmapSAINTContamModules-1.pdf index 587a5c9..ff8c9ab 100644 Binary files a/analysis/R/hClust/RawZHeatmapSAINTContamModules-1.pdf and b/analysis/R/hClust/RawZHeatmapSAINTContamModules-1.pdf differ diff --git a/analysis/R/hClust/RawZHeatmapSAINTHighHighModules-1.pdf b/analysis/R/hClust/RawZHeatmapSAINTHighHighModules-1.pdf index 51d2b1d..37bf8e1 100644 Binary files a/analysis/R/hClust/RawZHeatmapSAINTHighHighModules-1.pdf and b/analysis/R/hClust/RawZHeatmapSAINTHighHighModules-1.pdf differ diff --git a/analysis/R/hClust/RawZHeatmapSAINTLowLowModules-1.pdf b/analysis/R/hClust/RawZHeatmapSAINTLowLowModules-1.pdf index 8a0bc83..28d7180 100644 Binary files a/analysis/R/hClust/RawZHeatmapSAINTLowLowModules-1.pdf and b/analysis/R/hClust/RawZHeatmapSAINTLowLowModules-1.pdf differ diff --git a/analysis/R/hClust/RawZHeatmapSAINTSignalModule-1.pdf b/analysis/R/hClust/RawZHeatmapSAINTSignalModule-1.pdf index 86dc093..8bc5586 100644 Binary files a/analysis/R/hClust/RawZHeatmapSAINTSignalModule-1.pdf and b/analysis/R/hClust/RawZHeatmapSAINTSignalModule-1.pdf differ diff --git a/analysis/R/hClust/RawZHeatmapTMEMModule-1.pdf b/analysis/R/hClust/RawZHeatmapTMEMModule-1.pdf index 8fdde65..832f040 100644 Binary files a/analysis/R/hClust/RawZHeatmapTMEMModule-1.pdf and b/analysis/R/hClust/RawZHeatmapTMEMModule-1.pdf differ diff --git a/analysis/R/hClust/hClust_module_level_GO_GSEA.rda b/analysis/R/hClust/hClust_module_level_GO_GSEA.rda new file mode 100644 index 0000000..6a64e44 Binary files /dev/null and b/analysis/R/hClust/hClust_module_level_GO_GSEA.rda differ diff --git a/analysis/R/hClust/hClust_protein_level_GO_GSEA.rda b/analysis/R/hClust/hClust_protein_level_GO_GSEA.rda new file mode 100644 index 0000000..9c7c0e1 Binary files /dev/null and b/analysis/R/hClust/hClust_protein_level_GO_GSEA.rda differ