Parameter values

params

## $rds150
## [1] "data/ipms_150_sce.rds"
## 
## $coords150
## [1] "data/visnetwork_coords_150.csv"
## 
## $rds500
## [1] "data/ipms_500_sce.rds"
## 
## $coords500
## [1] "data/visnetwork_coords_500.csv"
## 
## $rds150untag
## [1] "data/ipms_150_untag_sce.rds"
## 
## $interactors
## [1] "data/Gene_list_TFinteractors_Pombase20240226.txt"
## 
## $idmap
## [1] "data/id_mapping_table.txt"
## 
## $rad24_rad25_tbl
## [1] "data/FigS6_Rad24_Rad25_interactingTFs_list_motifs.txt"

Load required packages and helper functions

suppressPackageStartupMessages({
    library(visNetwork)
    library(SingleCellExperiment)
    library(dplyr)
    library(ggplot2)
    library(igraph)
    library(shiny)
    library(cowplot)
    library(ggplotify)
    library(ggraph)
    library(latex2exp)
    library(tidygraph)
    library(ggrepel)
    library(ComplexHeatmap)
    library(flextable)
    library(kableExtra)
    library(ggupset)
    library(tidyr)
    library(jsonlite)
})

source("params/plot_settings.R")
source("params/get_testres_function.R")
source("params/mapping_functions.R")
source("params/ipms_params.R")

Settings

w_trad <-  84 # heatmap body width (mm)
h_trad <- 178 # total heatmap height (mm)

fs <- 6  # small font size
fl <- 10 # large font size
ft <- 12 # title font size

Read data

Here we load SingleCellExperiment objects corresponding to the low- and high-salt IP-MS experiments, as well as an object corresponding to the low-salt experiments, where each condition was compared to the untagged controls rather than the large complement group as in the other two objects.

sce150 <- readRDS(params$rds150)
sce500 <- readRDS(params$rds500)
sce150untag <- readRDS(params$rds150untag)

We also load a table with mappings between different types of feature identifiers, and a table with known, annotated physical interactions (obtained from PomBase).

idmap <- read.delim(params$idmap)

interactors <- read.delim(params$interactors) |>
    dplyr::select(Systematic.ID, Gene.name, Physical.interactors) |>
    tidyr::separate_rows(Physical.interactors, sep = ",") |>
    dplyr::mutate(Gene.name = .capitalize(Gene.name),
                  Physical.interactors = .capitalize(Physical.interactors))

Helper function

## Overwrite the get_legend function from cowplot temporarily, 
## since the cowplot one doesn't work with ggplot2 3.5.0
## see https://github.com/wilkelab/cowplot/issues/202

get_legend <- function(plot, legend = NULL) {
    gt <- ggplotGrob(plot)
    
    pattern <- "guide-box"
    if (!is.null(legend)) {
        pattern <- paste0(pattern, "-", legend)
    }
    
    indices <- grep(pattern, gt$layout$name)
    
    not_empty <- !vapply(
        gt$grobs[indices], 
        inherits, what = "zeroGrob", 
        FUN.VALUE = logical(1)
    )
    indices <- indices[not_empty]
    
    if (length(indices) > 0) {
        return(gt$grobs[[indices[1]]])
    }
    return(NULL)
}

Calculate pairwise protein similarities

To define the similarity between two proteins, the following procedure is used: Let \(t_{ij}\) and \(padj_{ij}\) denote the t-statistic and adjusted p-value for protein \(i\) in comparison \(j\), respectively. Then define \[\tilde{t}_{ij} = \left\{\begin{array}{cl} 0 & \quad \textrm{if } t_{ij} <= 0 \textrm{ or } padj_{ij} > 0.1 \\ t_{ij} & \quad \textrm{otherwise} \end{array}\right.\]

The similarity between proteins \(i\) and \(k\) is then defined by \[sim_{ik} = \frac{\sum_{j=1}^N\tilde{t}_{ij}\tilde{t}_{kj}}{max\left(\sum_{j=1}^N\tilde{t}_{ij}, \sum_{j=1}^N\tilde{t}_{kj}\right)}\]

or, if written in another way: \[sim_{ik} = min\left(\frac{\sum_{j=1}^N\tilde{t}_{ij}\tilde{t}_{kj}}{\sum_{j=1}^N\tilde{t}_{ij}},\frac{\sum_{j=1}^N\tilde{t}_{ij}\tilde{t}_{kj}}{\sum_{j=1}^N\tilde{t}_{kj}}\right)\] it can be interpreted as the smallest of two weighted means.

The following helper function extracts t-statistics and adjusted p-values and calculates the truncated values mentioned above (\(\tilde{t}\)).

## Extract t-statistics and adjusted p-values and threshold
getTruncT <- function(sce, adjpThr = 0.1) {
    ## Extract adjusted p-values and log-fold changes
    tc <- .getTestCols(sce)
    
    adjpvals <- tc$adjp
    tstats <- tc$tstat
    colnames(adjpvals) <- .getProteinNameFromComparison(colnames(adjpvals))
    colnames(tstats) <- .getProteinNameFromComparison(colnames(tstats))
    adjpvals <- as.matrix(adjpvals[, grep("untagged|Untagged", colnames(adjpvals), invert = TRUE)])
    tstats <- as.matrix(tstats[, grep("untagged|Untagged", colnames(tstats), invert = TRUE)])
    stopifnot(colnames(adjpvals) == colnames(tstats),
              rownames(adjpvals) == rownames(tstats))

    ## Need to pass a relatively loose significance threshold in at least one IP
    ## (don't want things to be similar just because they are weakly but consistently
    ## changing in many IPs)
    tstats[adjpvals > adjpThr] <- 0
    tstats[tstats <= 0] <- 0
    tstats[is.na(tstats)] <- 0

    tstats
}

We also define a function to calculate the similarities between each pair of proteins.

## Helper function to get matrix of similarities
simsFromResults <- function(sce, adjpThr = 0.1) {
    
    tstats <- getTruncT(sce = sce, adjpThr = adjpThr)
    
    simil <- tstats %*% t(tstats)
    simil0 <- simil
    for (i in seq_len(nrow(simil))) {
        for (j in seq(i, ncol(simil))) {
            simil[i, j] <- simil0[i, j] / max(sum(tstats[i, ]), sum(tstats[j, ]))
            simil[j, i] <- simil0[j, i] / max(sum(tstats[i, ]), sum(tstats[j, ]))
        }
    }
    simil[is.na(simil)] <- 0
    
    ## Don't need self-connections
    diag(simil) <- 0

    stopifnot(rownames(simil) == colnames(simil))

    list(simil = simil, tstats = tstats, TFs = colnames(tstats))
}

To illustrate, let’s calculate the similarity for a specific example pair (Atf1 and Pcr1):

tmp <- getTruncT(sce = sce150, adjpThr = 0.1)[c("Atf1", "Pcr1"), ]
round(tmp[, colSums(tmp) != 0], 1)

##      Gaf1 Hsf1 Atf1 Pcr1
## Atf1  3.0  2.4 12.2  8.4
## Pcr1  2.3  2.5  7.2  6.9

\[sim_{\textrm{atf1,pcr1}} = \min\left(\frac{3\cdot2.3+2.4\cdot2.5+12.2\cdot7.2+8.4\cdot6.9}{2.3+2.5+7.2+6.9},\frac{3\cdot2.3+2.4\cdot2.5+12.2\cdot7.2+8.4\cdot6.9}{3+2.4+12.2+8.4}\right) = min\left(8.4,6.1\right) = 6.1\]

simsFromResults(sce150[c("Atf1", "Pcr1"), ], adjpThr = 0.1)$simil

##          Atf1     Pcr1
## Atf1 0.000000 6.114673
## Pcr1 6.114673 0.000000

Next, we calculate similarities for all pairs of proteins in each of the data sets.

simil150 <- simsFromResults(sce150, adjpThr = 0.1)
simil500 <- simsFromResults(sce500, adjpThr = 0.1)
maxSimil <- round(max(max(simil150$simil, na.rm = TRUE),
                      max(simil500$simil, na.rm = TRUE)), 1)

Generate networks

Given the matrix of similarities calculated above, we create a graph by first thresholding the similarities (setting all values below a given threshold to 0), and then considering the remaining (non-zero) values in the matrix as representing edges in the graph. Thus, two proteins in the graph are connected if their enrichment patterns across the IP experiments are similar. Only graph components with at least two proteins (i.e., no singletons) are displayed.

createNetwork <- function(similList, threshold, known_interactors, settings = "igraph") {
    ## Get similarity matrix and apply threshold
    simil <- similList$simil
    tfs <- similList$TFs
    simil[simil < threshold] <- 0
    simil[is.na(simil)] <- 0
    
    ## Create graph
    gr <- graph_from_adjacency_matrix(simil, weighted = TRUE)
    gr <- as.undirected(gr, mode = "collapse", edge.attr.comb = list(weight = "mean"))
    
    ## Remove singletons
    comps <- components(gr)
    too_small <- which(table(comps$membership) < 2)
    keep <- V(gr)[!(comps$membership %in% too_small)]
    gr <- induced_subgraph(gr, keep)
    
    ## Set attributes for vertices and edges
    ## Need col2hex here since visNetwork doesn't understand all R colors (e.g. gray70)
    V(gr)$color <- gplots::col2hex(bg_color)
    V(gr)$color[names(V(gr)) %in% tfs] <- main_colors[5]
    V(gr)$vertexclass <- "other"
    V(gr)$vertexclass[names(V(gr)) %in% tfs] <- "Bait"
    
    if (settings == "visNetwork") {
        V(gr)$shape <- "triangle"
        V(gr)$size <- 15
        V(gr)$label.cex = 2.08
    } else if (settings == "igraph") {
        V(gr)$shape <- "circle"
    } else {
        stop("Unknown 'settings' parameter")
    }
    
    V(gr)$shape[names(V(gr)) %in% tfs] <- "square"
    E(gr)$value <- E(gr)$weight
    E(gr)$width <- E(gr)$weight / max(simil) * 1.5
    
    ## Indicate annotated physical interactions
    eL <- as_edgelist(gr)
    for (i in seq_along(E(gr))) {
        v1 <- eL[i, 1]
        v2 <- eL[i, 2]
        tmp <- known_interactors |>
            dplyr::filter(((Systematic.ID == v1 | Gene.name == v1) & 
                               Physical.interactors == v2) | 
                              ((Systematic.ID == v2 | Gene.name == v2) & 
                                   Physical.interactors == v1))
        if (nrow(tmp) == 0) {
            ## Novel interaction
            ## One of the elements is a transcription factor -> color orange
            if (v1 %in% tfs || v2 %in% tfs) {
                E(gr)[i]$color <- main_colors[3]
                E(gr)[i]$edgeclass <- "Unannotated TF interaction"
            } else{
                E(gr)[i]$color <- gplots::col2hex(na_color)
                E(gr)[i]$edgeclass <- "Other"
            }
        } else {
            E(gr)[i]$color <- gplots::col2hex(na_color)
            E(gr)[i]$edgeclass <- "Other"
        }
    }
    E(gr)$edgeclass <- factor(E(gr)$edgeclass, 
                              levels = c("Unannotated TF interaction", 
                                         "Other"))
    gr
}

Next, we create networks for the low- and high-salt experiments.

gr150 <- createNetwork(similList = simil150, threshold = 6, 
                       known_interactors = interactors, settings = "igraph")
gr500 <- createNetwork(similList = simil500, threshold = 6, 
                       known_interactors = interactors, settings = "igraph")

We export networks (with a lower similarity threshold, to allow for interactive exploration) for use by TFexplorer.

gr150tfx <- createNetwork(similList = simil150, threshold = 3, 
                          known_interactors = interactors, settings = "igraph")
gr500tfx <- createNetwork(similList = simil500, threshold = 3, 
                          known_interactors = interactors, settings = "igraph")

## 150 mM NaCl
n150 <- as_tbl_graph(gr150tfx) |>
    activate(nodes) |>
    data.frame() |>
    mutate(idx = seq_along(name))
e150 <- as_tbl_graph(gr150tfx) |>
    activate(edges) |>
    data.frame() |>
    left_join(n150 |> select(name, idx) |> rename(from = idx, from_name = name), 
              by = "from") |>
    left_join(n150 |> select(name, idx) |> rename(to = idx, to_name = name),
              by = "to")
exp150 <- list(
    nodes = lapply(seq_len(nrow(n150)), function(i) {
        list(id = n150$name[i], 
             label = n150$name[i], 
             value = 1)
    }),
    edges = lapply(seq_len(nrow(e150)), function(i) {
        list(source = e150$from_name[i],
             target = e150$to_name[i],
             value = e150$value[i])
    })
)
exp150 <- toJSON(exp150,
                 pretty = TRUE, 
                 auto_unbox = TRUE)
write(exp150, file = "data/networkData_150mMSalt.json")

## 500 mM NaCl
n500 <- as_tbl_graph(gr500tfx) |>
    activate(nodes) |>
    data.frame() |>
    mutate(idx = seq_along(name))
e500 <- as_tbl_graph(gr500tfx) |>
    activate(edges) |>
    data.frame() |>
    left_join(n500 |> select(name, idx) |> rename(from = idx, from_name = name), 
              by = "from") |>
    left_join(n500 |> select(name, idx) |> rename(to = idx, to_name = name),
              by = "to")
exp500 <- list(
    nodes = lapply(seq_len(nrow(n500)), function(i) {
        list(id = n500$name[i], 
             label = n500$name[i], 
             value = 1)
    }),
    edges = lapply(seq_len(nrow(e500)), function(i) {
        list(source = e500$from_name[i],
             target = e500$to_name[i],
             value = e500$value[i])
    })
)
exp500 <- toJSON(exp500,
                 pretty = TRUE, 
                 auto_unbox = TRUE)
write(exp500, file = "data/networkData_500mMSalt.json")

14-3-3 protein interactions

We extract all IPs where either Rad24 or Rad25 are pulled down, using a relatively lenient adjusted p-value cutoff.

## Get all t-statistics and visualize
res150untag <- .getTestCols(sce150untag, adjp_cutoff = 0.01, logfc_cutoff = 1)
rad <- res150untag$tstat[c("Rad24", "Rad25"), ]
colnames(rad) <- .getProteinNameFromComparison(colnames(rad))
rad <- t(rad)
hmrad <- Heatmap(rad, 
                 cluster_rows = TRUE, 
                 cluster_columns = FALSE, 
                 col = circlize::colorRamp2(
                     breaks = seq(1.5, 13.0, length.out = 64),
                     colors = colorRampPalette(
                         binary_heatmap_colors[c("FALSE", "TRUE")])(64)),
                 border = TRUE, 
                 border_gp = gpar(lwd = 0.5),
                 name = "Mod.\nt-stat.",
                 na_col = binary_heatmap_colors["FALSE"],
                 column_names_gp = gpar(fontsize = fl),
                 row_names_gp = gpar(fontsize = fs),
                 show_row_names = TRUE, 
                 show_column_names = TRUE, 
                 use_raster = FALSE,
                 show_heatmap_legend = TRUE, 
                 width = unit(w_trad / 1.5, "mm"), 
                 heatmap_height = unit(1.4 * h_trad, "mm"),
                 column_title_gp = gpar(fontsize = ft),
                 heatmap_legend_param = list(
                     title_gp = gpar(fontsize = fl),
                     border = "gray10",
                     border_gp = gpar(lwd = 0.5)
                 ))
set.seed(42L)
hm_rad <- grid.grabExpr(
    hmrad <- draw(hmrad, merge_legend = TRUE,
                 heatmap_legend_side = "right",
                 annotation_legend_side = "right"),
    width = w_trad / 1.5, height = 1.4 * h_trad
)

## List of experiments where rad24 and/or rad25 were pulled down
comps <- colnames(res150untag$interactor)[colSums(
    res150untag$interactor[c("Rad24", "Rad25"), ], na.rm = TRUE) > 0]
explist <- .getUniProtIdFromComparison(comps, idmap = idmap)
cat(explist)

## Q10280 Q9URT4 P33520 P78926 Q01663 O13658 Q10328 O14283 O59741 Q09838 P41412 P78871 O14335 Q9P7B0 P36631 Q9C0Z1 Q9P7H9 Q10096 Q1MTM9 O94392 O60130 O59744 O74915 O59830

gnames <- .getProteinNameFromComparison(comps)
cat(gnames)

## Gaf1 Ams2 Res1 Map1 Pap1 Pho7 Phx1 Prr1 Prt1 Prz1 Res2 Rst2 Scr1 Stb3 Ste11 Toe3 Dal81 SPAC11D3.17 SPAC1327.01c SPAC2H10.01 Ntu1 Ntu2 SPCC757.04 SPCC965.10

## Read table
set_flextable_defaults(font.size = 8)
rad24_rad25_tbl <- read.delim(params$rad24_rad25_tbl)[, 1:2] |>
    setNames(c("TF", "R-x-x-[S/T]-x-P")) |>
    dplyr::mutate(TF = .capitalize(TF))
stopifnot(all(gnames %in% rad24_rad25_tbl$TF))
ftbl <- flextable(rad24_rad25_tbl) |>
  width(width =  c(0.5, 1.25)) |>
  set_table_properties(layout = "fixed")
ggtbl <- ggplot() +
    theme_void() +
    annotation_custom(grid::rasterGrob(as_raster(ftbl)),
                      xmin = -Inf, xmax = Inf, ymin = 0, ymax = 1)

data.frame(`Number of TFs interacting with at least one 14-3-3 protein` = 
               nrow(rad24_rad25_tbl),
           `Number of such TFs with at least one optimal 14-3-3 binding motif` = 
               sum(rad24_rad25_tbl$`R-x-x-[S/T]-x-P` != ""),
           check.names = FALSE) |>
    pivot_longer(cols = everything()) |>
    kbl(col.names = NULL) |>
    kable_styling()

Number of TFs interacting with at least one 14-3-3 protein	24
Number of such TFs with at least one optimal 14-3-3 binding motif	10

# UpSet plot
dfupset <- as.data.frame(t(res150untag$interactor[c("Rad24", "Rad25"), ])) |>
    tibble::rownames_to_column("TF") |>
    dplyr::mutate(TF = .getProteinNameFromComparison(TF))

ovL <- as.list(as.data.frame(dfupset[, c("Rad24", "Rad25")]))
ovL <- lapply(ovL, function(i) dfupset$TF[which(i)])

pd <- stack(ovL) |>
    group_by(values) |>
    summarise(ind = list(ind)) |>
    dplyr::left_join(rad24_rad25_tbl, by = join_by("values" == "TF")) |>
    dplyr::mutate(has_motif = !is.na(`R-x-x-[S/T]-x-P`) & `R-x-x-[S/T]-x-P` != "")

ggrad_upset <- ggplot(pd, aes(x = ind)) +
    geom_bar(aes(fill = has_motif)) +
    scale_fill_manual(values = c(`TRUE` = main_colors[3], `FALSE` = main_colors[5])) +
    scale_x_upset() +
    labs(x = element_blank(), y = "Number of IPs", fill = "Has motif: ") +
    theme_cowplot(12) +
    theme(axis.ticks.x = element_blank(),
          legend.position = "bottom") +
    theme_combmatrix(combmatrix.label.extra_spacing = 5,
                     combmatrix.label.make_space = TRUE,
                     combmatrix.label.width = unit(18, "mm")) +
    guides(fill = guide_legend(nrow = 1))

Pho7 volcano plot

md <- metadata(sce500)
volc <- md$testres$tests$Pho7_500_plate_vs_compl_Pho7_500_plate |>
    mutate(signif = adj.P.Val < adjpThreshold & logFC > log2fcThreshold)
(pho7volc <- ggplot(volc, aes(x = logFC, y = mlog10p)) + 
        geom_text_repel(data = volc |> filter(signif),
                        aes(label = pid), min.segment.length = 0,
                        max.overlaps = Inf) + 
        geom_point(aes(alpha = signif, color = signif), 
                   show.legend = c(alpha = FALSE, color = TRUE)) + 
        scale_color_manual(values = c(`TRUE` = main_colors[5], `FALSE` = na_color), 
                           labels = c(`TRUE` = paste0("adj. p < ", adjpThreshold, 
                                                      " and logFC > ", log2fcThreshold)),
                           na.value = "black", breaks = "TRUE") + 
        scale_alpha_manual(values = c(`TRUE` = 1, `FALSE` = 0.5)) + 
        geom_point(data = volc |> dplyr::filter(pid == "Pho7") |>
                       dplyr::mutate(type = "Bait"), 
                   mapping = aes(fill = type), color = main_colors[3], 
                   shape = 21, alpha = 1, size = 2) + 
        scale_fill_manual(values = c(Bait = main_colors[3])) + 
        labs(x = "log2(fold-change)", y = "-log10(p-value)", 
             title = "Pho7 vs complement (500 mM NaCl IP-MS)") + 
        theme_cowplot(ft) + 
        theme(legend.position = "bottom") + 
        guides(color = guide_legend(nrow = 1, byrow = TRUE, title = ""),
               fill = guide_legend(nrow = 1, byrow = TRUE, title = "")))

## Warning: Removed 828 rows containing missing values or values outside the scale
## range (`geom_point()`).

Figure 6

## Create small example similarity matrix and graph, for three selected proteins
poi <- c("Atf1", "Pcr1", "Arp8")
simil150example <- simsFromResults(sce150[poi, ])
gr150example <- createNetwork(similList = simil150example, threshold = 6,
                              known_interactors = interactors, 
                              settings = "igraph")

## General formula
## See https://uliniemann.com/blog/2022-02-21-math-annotations-in-ggplot2-with-latex2exp/
eqn1 <- data.frame(x = 0, y = 0, 
                   label = r"( $sim_{ik} = min\left(\frac{\sum_{j=1}^N\tilde{t}_{ij}\tilde{t}_{kj}}{\sum_{j=1}^N\tilde{t}_{ij}},\frac{\sum_{j=1}^N\tilde{t}_{ij}\tilde{t}_{kj}}{\sum_{j=1}^N\tilde{t}_{kj}}\right)$ )")

ggeqn1 <- ggplot(eqn1) +
    geom_text(aes(x = x, y = y, 
                  label = TeX(label, output = "character"), 
                  color = label), 
              parse = TRUE,
              key_glyph = draw_key_blank,
              size = 16/.pt) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_reverse(expand = c(0, 0)) + 
    scale_color_manual(values = "transparent", labels = TeX(eqn1$label,
                                                            italic = TRUE), 
                       name = "") + 
    guides(color = guide_legend(label.position = "top")) + 
    theme_void() + 
    theme(legend.text = element_text(size = 11))

## Specific example (atf1, pcr1)
eqn2 <- data.frame(x = 0, y = 0, 
                   label = r"( $sim_{Atf1,Pcr1} = min\left(\frac{3 \cdot 2.3 + 2.4 \cdot 2.5 + 12.2 \cdot 7.2 + 8.4 \cdot 6.9}{2.3 + 2.5 + 7.2 + 6.9},\frac{3 \cdot 2.3 + 2.4 \cdot 2.5 + 12.2 \cdot 7.2 + 8.4 \cdot 6.9}{3 + 2.4 + 12.2 + 8.4}\right) = min\left(8.4, 6.1\right) = 6.1$ )")
# eqn2 <- data.frame(x = 0, y = 0, 
#                    label = r"( $sim_{\textrm{Atf1,Pcr1}} = \min\left(\frac{3 \cdot 2.3 + 2.4 \cdot 2.5 + 12.2 \cdot 7.2 + 8.4 \cdot 6.9}{2.3 + 2.5 + 7.2 + 6.9},\frac{3 \cdot 2.3 + 2.4 \cdot 2.5 + 12.2 \cdot 7.2 + 8.4 \cdot 6.9}{3 + 2.4 + 12.2 + 8.4}\right) = min\left(8.4, 6.1\right) = 6.1$ )")

ggeqn2 <- ggplot(eqn2) +
    geom_text(aes(x = x, y = y, 
                  label = TeX(label, output = "character")), 
              parse = TRUE, 
              size = 11/.pt) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_reverse(expand = c(0, 0)) + 
    theme_void()

## Truncated t-statistics
tstats_vis <- data.frame(simil150example$tstats) |>
    tibble::rownames_to_column("Protein") |>
    tidyr::pivot_longer(names_to = "IP", values_to = "Value", -Protein) |>
    dplyr::mutate(Protein = factor(Protein, levels = poi))

## Reorder IPs for better visualization
enr <- unique(tstats_vis$IP[tstats_vis$Value > 0])
notenr <- setdiff(tstats_vis$IP, enr)
allips <- c(notenr, enr)
intlen <- round(length(notenr) / (length(enr) + 1))
idx <- seq(from = intlen, by = intlen, length.out = length(enr))
levs <- allips[order(c(seq_along(notenr), idx + 0.5))]
tstats_vis$IP <- factor(tstats_vis$IP, levels = levs)

## Plot truncated tstats across all IPs and label relevant points
labels <- tstats_vis |>
    dplyr::filter(Value > 0) |>
    dplyr::pull(IP) |>
    unique()

tplot <- ggplot(tstats_vis, aes(x = IP, y = Value, color = Protein, fill = Protein)) + 
        geom_col() + 
        geom_text(
            aes(label = ifelse(IP %in% labels & Value > 0, round(Value, digits = 1), "")),
            fontface ="plain", color = "black", size = 2.5, nudge_y = 2) +
        scale_color_manual(values = c(main_colors[1], main_colors[5], "black")) +
        scale_fill_manual(values = c(main_colors[1], main_colors[5], "black")) +
        scale_y_continuous(expand = c(0, 1)) + 
        facet_wrap(~ Protein, ncol = 1) +
        labs(x = "IP", y = "Thresholded mod. t-stat.") + 
        theme_cowplot(font_size = 10) + 
        theme(legend.position = "none", 
              axis.text.x = element_blank(), 
              axis.ticks.x = element_blank(), 
              strip.background = element_rect(fill = NA, colour = NA), 
              strip.text = element_text(hjust = 0, size = 10, face = "bold"))

## Similarity matrix
similsub <- simil150example$simil[poi, poi] |>
    as.data.frame() |>
    tibble::rownames_to_column("row") |>
    tidyr::pivot_longer(names_to = "column", values_to = "Similarity", -row) |>
    dplyr::mutate(Similarity = round(Similarity, digits = 1),
                  row = factor(row, levels = rev(poi)),
                  column = factor(column, levels = poi)) |>
    dplyr::mutate(simlab = ifelse(Similarity == 0, "", Similarity)) |>
    dplyr::mutate(Similarity = ifelse(Similarity == 0, NA, Similarity))
similplot <- ggplot(similsub, aes(x = column, y = row, fill = Similarity)) + 
    geom_tile() + 
    geom_text(aes(label = simlab, color = Similarity > 4), 
              fontface = "bold") +
    scale_color_manual(values = c(`TRUE` = "white", `FALSE` = "black")) + 
    scale_fill_gradient(low = "white", high = binary_heatmap_colors["TRUE"], 
                    na.value = na_color, 
                    limits = range(similsub$Similarity, na.rm = TRUE)) + 
    theme_cowplot() + 
    theme(axis.line = element_blank(), 
          axis.ticks = element_blank(),
          legend.position = "bottom",
          legend.title = element_text(size = 10, vjust = 0.82),
          legend.margin = margin(0, 0, 0, 0),
          legend.box.spacing = unit(0, "pt")) + 
    labs(x = "", y = "") + 
    guides(color = "none")

## Induced graph
ggr150example <- ggraph(as_tbl_graph(gr150example), 
                        layout = cbind(x = c(-300, 800), y = c(800, -300)), 
                        weights = weight) + 
    geom_edge_link(aes(color = edgeclass, width = weight)) +
    scale_edge_color_manual(values = c(`Unannotated TF interaction` = main_colors[3], 
                                       Other = gplots::col2hex(na_color)),
                            name = "") + 
    scale_edge_width(range = c(1.2, 1.9), guide = "none") + 
    geom_node_point(size = 4, aes(fill = vertexclass, shape = vertexclass), 
                    color = "black", stroke = 0.1) + 
    scale_shape_manual(values = c(other = 21, Bait = 22), name = "") + 
    scale_fill_manual(values = c(other = gplots::col2hex(bg_color), Bait = main_colors[5]),
                      name = "") + 
    geom_node_text(aes(label = name), size = 3.5, repel = TRUE, 
                   nudge_x = 0, nudge_y = -5) + 
    theme_graph(base_family = "sans", plot_margin = margin(30, 10, 40, 10)) + 
    guides(fill = guide_legend(override.aes = list(size = 4, alpha = 1))) + 
    coord_cartesian(clip = "off") 

## Put together
figs6c <- cowplot::plot_grid(
    cowplot::plot_grid(
        tplot, 
        get_legend(ggeqn1, "right"), 
        similplot + theme(legend.position = "bottom"), 
        ggr150example + theme(legend.position = "none"), 
        align = "v", axis = "t", 
        labels = c("C", "", "", ""), nrow = 1, vjust = -0.5,
        rel_widths = c(1.1, 0.9, 1, 0.5)
    ),
    ggeqn2, ncol = 1, labels = "", 
    rel_heights = c(3, 1)
)

Put together

fig6b <- ggdraw() + 
    draw_image("schematics/Fig6B_Pho7_protein_domains_cartoon_Arial_RGB.png")

coords150 <- read.csv(params$coords150)
coords500 <- read.csv(params$coords500)
coords150$y <- -coords150$y
coords500$y <- -coords500$y

.plotGraph <- function(grph, coords, nudge_y) {
    ggraph(as_tbl_graph(grph), 
           layout = coords, 
           weights = weight) + 
    geom_edge_link(aes(color = edgeclass, width = weight)) +
    scale_edge_color_manual(values = c(`Unannotated TF interaction` = main_colors[3], 
                                       Other = gplots::col2hex(na_color)),
                            breaks = "Unannotated TF interaction",
                            name = "") + 
    scale_edge_width(range = c(0.2, 0.9), guide = "none") + 
    geom_node_point(size = 2, aes(fill = vertexclass, shape = vertexclass), 
                    color = "black", stroke = 0.1) + 
    scale_shape_manual(values = c(other = 21, Bait = 22), name = "",
                       breaks = "Bait") + 
    scale_fill_manual(values = c(other = gplots::col2hex(bg_color), Bait = main_colors[5]),
                      name = "", breaks = "Bait") + 
    geom_node_text(aes(label = name), size = 2.5, repel = FALSE, 
                   nudge_x = 0, nudge_y = nudge_y) + 
    theme_graph(base_family = "sans", plot_margin = margin(10, 10, 10, 10)) + 
    guides(fill = guide_legend(override.aes = list(size = 4, alpha = 1)))
}
ggr150 <- .plotGraph(grph = gr150, coords = as.matrix(coords150[, c("x", "y")]),
                     nudge_y = -45)
ggr500 <- .plotGraph(grph = gr500, coords = as.matrix(coords500[, c("x", "y")]),
                     nudge_y = -35)

cowplot::plot_grid(
    cowplot::plot_grid(
        pho7volc, 
        cowplot::plot_grid(
            cowplot::plot_grid(NULL, fig6b, NULL, nrow = 1,
                               rel_widths = c(0.25, 1, 0.1)),  
            ggrad_upset,
            ncol = 1, labels = c("B", "C"),
            rel_heights = c(1, 1.75)), 
        nrow = 1, labels = c("A", "")
    ),
    NULL,
    cowplot::plot_grid(
        cowplot::plot_grid(
            ggr150 + theme(legend.position = "none") + 
                ggtitle("150 mM NaCl IP-MS") + 
                theme(plot.title = element_text(hjust = 0.075, size = ft)), 
            ggr500 + theme(legend.position = "none") + 
                ggtitle("500 mM NaCl IP-MS") +
                theme(plot.title = element_text(hjust = 0.075, size = ft)), 
            nrow = 1, labels = c("D", "E"), vjust = 1.8
        ),
        get_legend(ggr150 + theme(legend.position = "bottom"), "bottom"),
        nrow = 2, labels = c("", ""), rel_heights = c(1, 0.1)
    ),
    rel_heights = c(3.5, 0.3, 4.25), 
    ncol = 1
)

Supplementary figure