OGT_analysis

Packages

library(tidyverse)
library(ggplot2);theme_set(cowplot::theme_cowplot(font_size = 12) + theme(panel.grid.major = element_line(colour = "lightgrey", linewidth = 0.2), panel.grid.minor = element_line(colour = "lightgrey", linewidth = 0.2)))
library("reshape2")
library(ggrepel)
library(knitr)
library(RColorBrewer)
library(limma)
library(STRINGdb)
library(plotly)
library(shiny)

mutate <- dplyr::mutate
select <- dplyr::select
group_by <- dplyr::group_by
filter <- dplyr::filter
count <- dplyr::count

Graphics

options(
  ggplot2.discrete.colour = c("mediumvioletred", "cornflowerblue", "mediumseagreen", "goldenrod2", "lightsalmon2", "red3"),
  ggplot2.discrete.fill = c("mediumvioletred", "cornflowerblue", "mediumseagreen", "goldenrod2", "lightsalmon2", "red3")
)
tilted <- theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))
blank <- theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())

Metadata

TMT annotation: Channels 1-3: 0h 4-6: 2h 7-9: 4h 10-12: 24h 13-16: 48h 17-19: 7 days

# Object to store all dfs and matrices
res_ogt_glyco <- list()

res_ogt_glyco$TMT_info <- data.frame(
   sample = c(
    "sample-01", "sample-02", "sample-03",
    "sample-04", "sample-05", "sample-06",
    "sample-07", "sample-08", "sample-09",
    "sample-10", "sample-11", "sample-12",
    "sample-13", "sample-14", "sample-15", 
    "sample-16", "sample-17", "sample-18"
  ),
  time = factor(c(
    rep("0h", 3),
    rep("2h", 3),
    rep("4h", 3),
    rep("24h", 3),
    rep("48h", 3),
    rep("168h", 3)
  )),
  replicate = c(
    rep(c("rep1", "rep2", "rep3"), 6)
  )
)

Full proteome analysis

1. Load data

setwd("~/Documents/01_repos/help/OGT/")
res_ogt_glyco$FP <- read_tsv("protein.tsv")

# repair columnames
colnames(res_ogt_glyco$FP) <- gsub(" ", "\\.", colnames(res_ogt_glyco$FP))

# remove contaminants
res_ogt_glyco$FP <- res_ogt_glyco$FP %>% 
  filter(grepl("Homo sapiens", Organism) & !(grepl("contam", Protein))) %>% 
  filter(!(grepl("KRT", Gene)))

2. Normalisation

normalisation <- function(df) {
  # perform vsn normalisation on long df per protein and give back in right format
  m <- acast(df, Protein.ID ~ sample, value.var = "quant", fun.aggregate = function(x) sum(x, na.rm =T))
  m[m==0] <- NA
  m<-m[complete.cases(m),]
  m_norm <- limma::normalizeVSN(m)
  vsn::meanSdPlot(m_norm)
  df <- melt(m_norm, value.name = "quant_norm", varnames = c("Protein.ID", "sample")) %>% 
    inner_join(df, by = c("Protein.ID", "sample")) 
  
  return(tibble::tibble(df))
}

res_ogt_glyco$FP_normalised <- tt <- res_ogt_glyco$FP %>% 
  select(Protein.ID, Gene, matches(res_ogt_glyco$TMT_info$sample)) %>% 
  melt(variable.name = "sample", value.name = "quant") %>% 
  inner_join(res_ogt_glyco$TMT_info, by = c("sample")) %>% 
  drop_na() %>% 
  group_modify(~tibble(normalisation(.))) %>% 
  ungroup() %>% 
  mutate(time=factor(time, levels = c("0h", "2h", "4h", "24h", "48h", "168h")))

res_ogt_glyco$FP_normalised %>% 
  ggplot(aes(x = time, y = log2(quant), colour = replicate)) +
  geom_boxplot()

res_ogt_glyco$FP_normalised %>% 
  ggplot(aes(x = time, y = quant_norm, colour = replicate)) +
  geom_boxplot()

OGT and OGA

res_ogt_glyco$FP_normalised %>% 
  filter(Gene %in% c("OGA","OGT")) %>% 
  mutate(time = as.numeric(str_extract(time, "\\d+"))) %>% 
  filter(time != "168") %>% 
  ggplot(aes(x= time, y= quant_norm, colour = Gene, group = Gene)) +
  geom_point() +
  stat_summary(geom = "line", fun = mean, size = 1) +
  labs(y = "normalised TMT intensity", x = "OGT degron induction time [h]")

3. PCA

pca_plot <- function(data_matrix, x, y, pal) {
  m <- data_matrix 
  n_obs <- nrow(m)

  m <- m[matrixStats::rowVars(m, na.rm = T) %>%
    order(decreasing = T) %>%
    head((n_obs / 100) * 10), ]

  pca <- prcomp(t(m))
  
  data <- pca$x %>%
    as.data.frame() %>%
    rownames_to_column("sample") %>%
    select(sample, x = {{ x }}, y = {{ y }}) %>% 
     separate(sample, c("time",  "replicate"), sep = "_", remove = F)
  
  var <- (pca$sdev)^2/sum(pca$sdev^2)*100 
  vx <- as.numeric(str_extract(as.character({{ x }}) , "\\d"))
  vy <- as.numeric(str_extract(as.character({{ y }}) , "\\d"))

  p <- data %>%
     mutate(time=factor(time, levels = c("0h", "2h", "4h", "24h", "48h", "168h"))) %>% 
    ggplot(aes(x = x , y = y , colour = time, shape = replicate)) +
    geom_hline(yintercept = 0, colour = "lightgrey") +
    geom_vline(xintercept = 0, colour = "lightgrey") +
    geom_point(size = 3) +
    scale_colour_manual(values = pal) +
    cowplot::panel_border() +
    #lims(x = c(-10,30), y = c(-10,10)) +
    #guides(fill = "none", shape = "none") +
    labs(
      x = paste0(as.character({{ x }}), " (", round(var[vx], digits = 1), "%)"),
      y = paste0(as.character({{ y }}), " (", round(var[vy], digits = 1), "%)")
    )

  plot(p)
}

m <- res_ogt_glyco$FP_normalised %>%
  mutate(sample = paste( time, replicate, sep = "_")) %>%
  dcast(
    Protein.ID ~ sample,
    value.var = "quant_norm", fill = NA
  ) %>%
  drop_na() %>%
  column_to_rownames("Protein.ID") %>%
  as.matrix()

pca_plot(m, "PC1", "PC2", pal = c("#F2F0F7", "#DADAEB", "#BCBDDC", "#9E9AC8", "#756BB1", "#54278F"))

m <- res_ogt_glyco$FP_normalised %>%
  filter(time != "168h") %>% 
  mutate(sample = paste( time, replicate, sep = "_")) %>%
  dcast(
    Protein.ID ~ sample,
    value.var = "quant_norm", fill = NA
  ) %>%
  drop_na() %>%
  column_to_rownames("Protein.ID") %>%
  as.matrix()

pca_plot(m, "PC1", "PC2", pal = c("#F2F0F7", "#DADAEB", "#BCBDDC", "#9E9AC8", "#756BB1", "#54278F"))

4. Limma

res_ogt_glyco$eset_FP <- res_ogt_glyco$FP_normalised %>%
  ungroup() %>%
  mutate(ID = paste0(Protein.ID, "_", Gene),  sample = paste(time, replicate, sep = "_")) %>% 
  select(ID, sample, quant_norm) %>%
  dcast(ID ~ sample) %>%
  #select(ID, res_phospho_timecourse$TMT_info$sample) %>% 
  drop_na() %>%
  column_to_rownames("ID") %>%
  as.matrix()

t <- factor(str_extract(colnames(res_ogt_glyco$eset_FP), "\\d+h"), 
            levels = c("0h", "2h", "4h", "24h", "48h", "168h"))

res_ogt_glyco$design <- model.matrix(~ 1 + t)
res_ogt_glyco$design

   (Intercept) t2h t4h t24h t48h t168h
1            1   0   0    0    0     0
2            1   0   0    0    0     0
3            1   0   0    0    0     0
4            1   0   0    0    0     1
5            1   0   0    0    0     1
6            1   0   0    0    0     1
7            1   0   0    1    0     0
8            1   0   0    1    0     0
9            1   0   0    1    0     0
10           1   1   0    0    0     0
11           1   1   0    0    0     0
12           1   1   0    0    0     0
13           1   0   0    0    1     0
14           1   0   0    0    1     0
15           1   0   0    0    1     0
16           1   0   1    0    0     0
17           1   0   1    0    0     0
18           1   0   1    0    0     0
attr(,"assign")
[1] 0 1 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$t
[1] "contr.treatment"

fit <- lmFit(res_ogt_glyco$eset_FP, res_ogt_glyco$design)
fit_FP <- eBayes(fit)


res_ogt_glyco$limma_results <- lapply(colnames(fit_FP$coefficients), function(x) {
  limma::topTable(fit_FP, coef = x, number = Inf) %>%
    mutate(contrast = x, dataset = "FP") %>%
    rownames_to_column(var = "Protein.ID")
}) %>%
  bind_rows() %>%
  filter(contrast != "(Intercept)") %>%
  mutate(contrast = str_remove(contrast, "t")) %>% 
  mutate(contrast = factor(contrast, levels = c("0h", "2h", "4h", "24h", "48h", "168h"))) %>% 
  mutate(
    hit = ifelse((abs(logFC) > log2(1.5) & adj.P.Val < 0.05), "hit", "no hit"),
    direction = ifelse(logFC > 0 & hit == "hit", "up", "no change"),
    direction = ifelse(logFC < 0 & hit == "hit", "down", direction)
  ) %>% 
  separate(Protein.ID, into = c("Protein.ID", "Gene"), sep = "_")

Hits

ggplot(res_ogt_glyco$limma_results %>% filter(hit == "hit"), 
       aes(contrast, fill = direction)) +
  geom_bar(pos = "dodge") +
  scale_fill_manual(values = c("dodgerblue4", "indianred")) +
  theme(axis.text.x = element_text(angle = 45, hjust =1, vjust = 1)) +
  facet_wrap(~dataset, scales = "free") +
  labs(y = " # hits")

ggplot(res_ogt_glyco$limma_results %>% filter(hit == "hit" & contrast != "168h"), 
       aes(contrast, fill = direction)) +
  geom_bar(pos = "dodge") +
  scale_fill_manual(values = c("dodgerblue4", "indianred")) +
  labs(y = " # significantly regulated proteins", x = "OGT degron induction time")

DT::datatable(res_ogt_glyco$limma_results %>%  select(Gene, contrast, logFC, adj.P.Val),
  filter = "top",
)

res_ogt_glyco$limma_results %>% 
  group_by(contrast, direction) %>% 
  dplyr::count()

# A tibble: 15 × 3
# Groups:   contrast, direction [15]
   contrast direction     n
   <fct>    <chr>     <int>
 1 2h       down          3
 2 2h       no change  7763
 3 2h       up            2
 4 4h       down          2
 5 4h       no change  7761
 6 4h       up            5
 7 24h      down         91
 8 24h      no change  7650
 9 24h      up           27
10 48h      down        107
11 48h      no change  7607
12 48h      up           54
13 168h     down        639
14 168h     no change  6479
15 168h     up          650

#res_ogt_glyco$limma_results %>%  select(Gene, contrast, logFC, adj.P.Val, P.Value) %>%  openxlsx::write.xlsx("2025-02-27_OGT_FP_DEresults.xlsx")

Volcano

p <- ggplot(
  res_ogt_glyco$limma_results ,
  aes(x = logFC, y = -log10(P.Value), colour = direction, label = Gene)
) +
  geom_point(alpha = 0.5, stroke = 0) +
  scale_colour_manual(values = c("dodgerblue4", "lightgrey", "indianred")) +
  facet_wrap( ~ contrast, nrow = 1) +
  labs(x = "log2 fold-change", y = "-log10  p-value") +
  cowplot::panel_border()

plotly::ggplotly(p, tooltip = "label")

ggplot(
  res_ogt_glyco$limma_results %>%
    mutate(label = ifelse(Gene %in% c("OGT", "OGA"), Gene, "")) %>%
    mutate(strength = ifelse(Gene %in% c("OGT", "OGA"), "OGA|OGT", direction)),
  aes(x = logFC, y = -log10(P.Value), colour = direction, label = label)
) +
  geom_vline(xintercept = 0) +
  geom_vline(xintercept = c(0.6, -0.6), linetype = 3) +
  geom_point(aes(size = strength, alpha = strength)) +
  geom_text_repel(max.overlaps = Inf, size = 2, colour = "black", min.segment.length = 0.2, force_pull = 20) +
  scale_colour_manual(values = c("dodgerblue4", "lightgrey", "indianred")) +
  scale_size_manual(values = c(0.8,0.8, 2, 0.8)) +
   scale_alpha_manual(values = c(0.2,0.2, 1, 0.2)) +
  facet_wrap(~contrast, nrow = 1) +
  labs(x = "log2 fold-change", y = "-log10  p-value") +
  cowplot::theme_cowplot() +
  cowplot::panel_border() +
  guides(size = "none", alpha = "none")

OGT and OGA profiles

res_ogt_glyco$limma_results %>% 
  filter(Gene %in% c("OGA","OGT")) %>% 
  mutate(time = as.numeric(str_extract(contrast, "\\d+"))) %>% 
  filter(time != 168) %>% 
  bind_rows(data.frame(time = 0, logFC = 0, Gene = c("OGA", "OGT"))) %>% 
  ggplot(aes(x= time, y= logFC, colour = Gene, group = Gene)) +
  geom_hline(yintercept = 0) +
  geom_point() + 
  stat_summary(geom = "line", fun = mean, size = 1) +
  labs(y = "log2 fold-change against 0h", x = "OGT degron induction time [h]")

5. Clustering

hits <- res_ogt_glyco$limma_results %>% 
  filter(hit == "hit") 

cluster_matrix <- res_ogt_glyco$limma_results %>% 
  distinct(Protein.ID, logFC, contrast) %>%
  dcast(Protein.ID~contrast, value.var = "logFC") %>% 
  filter(Protein.ID %in% hits$Protein.ID) %>% 
  column_to_rownames("Protein.ID") %>% 
  as.matrix()

cluster_matrix  %>% 
  ComplexHeatmap::Heatmap(
  show_row_names = F,
  col = circlize::colorRamp2(c(-2, 0, 2), c("dodgerblue4", "white","indianred3")),
  name = "logFC"
  )

without timepoint 7 days

hits <- res_ogt_glyco$limma_results %>% 
  filter(hit == "hit" & contrast != "168h") 

cluster_matrix <- res_ogt_glyco$limma_results %>% 
  distinct(Protein.ID, logFC, contrast) %>%
  filter(contrast != "168h") %>% 
  dcast(Protein.ID~contrast, value.var = "logFC") %>% 
  filter(Protein.ID %in% hits$Protein.ID) %>% 
  column_to_rownames("Protein.ID") %>% 
  as.matrix()

cluster_matrix  %>% 
  ComplexHeatmap::Heatmap(
  show_row_names = F,
  cluster_columns = F,
  col = circlize::colorRamp2(c(-2, 0, 2), c("dodgerblue4", "white","indianred3")),
  name = "logFC"
  )

set.seed(1)
neural_clusters <- tt <- data.frame(
  cluster = cclust::cclust(
    cluster_matrix,
    dist = "euclidean", 
    method = "neuralgas", 
    centers = 4)$cluster,
  Protein.ID = rownames(cluster_matrix)
) %>%
  inner_join(
    res_ogt_glyco$limma_results
  ) %>%
  select(Protein.ID, Gene, cluster, contrast, logFC) %>% 
  mutate(time = as.numeric(str_extract(contrast, "\\d+"))) %>%
  bind_rows(data.frame(cluster = 1, contrast = "t0", time = (0), logFC = 0, Protein.ID = "dummy")) %>%
  complete(cluster, contrast, Protein.ID) %>%
  mutate(
    logFC = ifelse(contrast == "t0", 0, logFC),
    time = ifelse(contrast == "t0", 0, time),
    cluster = factor(cluster)
  ) %>%
  ungroup() %>%
  filter(Protein.ID != "dummy") %>% 
  drop_na()


neural_clusters %>%
  group_by(contrast, cluster, time) %>%
  summarise(logFC = median(logFC, na.rm = T)) %>%
  ggplot(
    aes(time, logFC, colour = cluster)
  ) +
  geom_point(size = 2) +
  geom_hline(yintercept = c(0.6, 0, -0.6), linetype = 3) +
  geom_hline(yintercept = c(0), linetype = 1) +
  geom_line(aes(group = cluster), size = 1) +
  scale_colour_manual(values = c("2" = "indianred", "1" = "darkred", "3" = "dodgerblue1", "4" = "dodgerblue4")) +
  labs(y = "median\nlog2 fold-change", x = "time [min]") +
  cowplot::theme_cowplot()

neural_clusters %>% 
  filter(time != 0) %>% 
  distinct(cluster, Protein.ID) %>% 
  ungroup() %>% 
  group_by(cluster) %>% 
  count

# A tibble: 4 × 2
# Groups:   cluster [4]
  cluster     n
  <fct>   <int>
1 1          53
2 2          16
3 3         126
4 4          13

6. Enrichment

hits <- res_ogt_glyco$limma_results %>% 
  filter(hit == "hit" ) 

#| cache: true
msigdf_decoupler <- msigdbr::msigdbr(species = "Homo sapiens") %>%
  filter(gs_subcat %in% c("CP:REACTOME", "GO:BP", "GO:CC", "GO:MF")) %>%
  transmute(source = str_replace(gs_name, "HALLMARK_", ""), mor = 1, target = gene_symbol) %>%
  distinct()

anno <- msigdf_decoupler %>% 
  filter(target %in% hits$Gene) %>% 
  group_by(source) %>% 
  summarise(members = paste(target, collapse="|"))

decoupler_results <- res_ogt_glyco$limma_results %>%
  distinct(Gene, contrast, logFC) %>%
  nest(data = c(Gene, logFC, contrast))%>%
  mutate(ER = map(data, function(df) {
    x = df %>%  acast(Gene ~ contrast, value.var = "logFC", fun.aggregate = mean)
    decoupleR::run_wmean(mat = x, network = msigdf_decoupler)
  }
)) %>%
  select(-data) %>%
  unnest(cols = c(ER)) %>%
  mutate(condition = factor(condition, levels = c("0h", "2h", "4h", "24h", "48h", "168h"))) %>% 
  inner_join(anno)

Visualise significant pathways of all modalities and timepoints together

int <- decoupler_results %>%
  filter(statistic == "norm_wmean" & p_value < 0.05) %>%
  group_by(condition) %>%
  slice_max(order_by = abs(score), n= 5)

p<- decoupler_results %>%
  filter(statistic == "norm_wmean" & source %in% int$source) %>%
  mutate(n_sig_proteins = str_count(members, "|")) %>% 
  ggplot(aes(y = reorder(source, score), x= score, colour = condition, label = n_sig_proteins)) +
  geom_vline(xintercept = 0, linetype = 1) +
  geom_vline(xintercept = c(-3,3), linetype =3) +
  geom_point() +
  scale_colour_manual(values = c("#F2F0F7", "#DADAEB", "#BCBDDC", "#9E9AC8", "#756BB1", "#54278F")) +
  labs(x="pathway enirchment score", y = "")
p

plotly::ggplotly(p, tooltip = "label")

int <- decoupler_results %>%
  filter(statistic == "norm_wmean" & p_value < 0.05 & abs(score) > 1.7) %>%
  filter(grepl("MITOCH|RIBOSOM", source))

p<- decoupler_results %>%
  filter(statistic == "norm_wmean" & source %in% int$source) %>%
  mutate(n_sig_proteins = str_count(members, "|")) %>% 
  ggplot(aes(y = reorder(source, score), x= score, colour = condition, label = n_sig_proteins)) +
  geom_vline(xintercept = 0, linetype = 1) +
  geom_vline(xintercept = c(-3,3), linetype =3) +
  geom_point() +
  scale_colour_manual(values = c("#F2F0F7", "#DADAEB", "#BCBDDC", "#9E9AC8", "#756BB1", "#54278F")) +
  labs(x="pathway enirchment score", y = "")
p

plotly::ggplotly(p, tooltip = "label")

DT::datatable(decoupler_results %>%   filter(statistic == "norm_wmean" & p_value < 0.05),
  filter = "top",
)

without 168h

decoupler_results <- res_ogt_glyco$limma_results %>%
  filter(contrast != "168h") %>% 
  distinct(Gene, contrast, logFC) %>%
  nest(data = c(Gene, logFC, contrast))%>%
  mutate(ER = map(data, function(df) {
    x = df %>%  acast(Gene ~ contrast, value.var = "logFC", fun.aggregate = mean)
    decoupleR::run_wmean(mat = x, network = msigdf_decoupler)
  }
)) %>%
  select(-data) %>%
  unnest(cols = c(ER)) %>%
  mutate(condition = factor(condition, levels = c("0h", "2h", "4h", "24h", "48h", "168h"))) %>% 
  inner_join(anno)

int <- decoupler_results %>%
  filter(statistic == "norm_wmean" & p_value < 0.05) %>%
  group_by(condition) %>%
  slice_max(order_by = abs(score), n= 5)

p<- decoupler_results %>%
  filter(statistic == "norm_wmean" & source %in% int$source) %>%
  mutate(n_sig_proteins = str_count(members, "|")) %>% 
  ggplot(aes(y = reorder(source, score), x= abs(score), fill = condition, label = n_sig_proteins)) +
  geom_vline(xintercept = 0, linetype = 1) +
  geom_vline(xintercept = c(3), linetype =3) +
  geom_point(size = 3, shape = 21, colour = "darkgrey") +
  #scale_size(range = c(2,4)) +
  scale_fill_manual(values = c("#F2F0F7", "#DADAEB", "#BCBDDC", "#9E9AC8", "#756BB1", "#54278F")) +
  labs(x="absolute pathway enirchment score", y = "")
p

Ribosome proteins

ribosome <- msigdf_decoupler %>%  filter(source == "GOCC_RIBOSOME")
mito <- msigdf_decoupler %>%  filter(source == "GOCC_MITOCHONDRION")

data1 <- res_ogt_glyco$limma_results %>%
  mutate(annotation = ifelse(Gene %in% ribosome$target, "ribosomal", "other")) %>%
  mutate(annotation = ifelse(Gene %in% mito$target, "mitochondrial", annotation)) %>%
  mutate(annotation = ifelse(Gene %in% c("OGT", "OGA"), "OGA|OGT", annotation)) %>%
  filter(annotation == "ribosomal") %>%
  filter(contrast != "168h")

data2 <- res_ogt_glyco$limma_results %>%
  mutate(annotation = ifelse(Gene %in% ribosome$target, "ribosomal", "other")) %>%
  mutate(annotation = ifelse(Gene %in% mito$target, "mitochondrial", annotation)) %>%
  mutate(annotation = ifelse(Gene %in% c("OGT", "OGA"), "OGA|OGT", annotation)) %>%
  filter(annotation == "other") %>%
  filter(contrast != "168h")

data3 <- res_ogt_glyco$limma_results %>%
  mutate(annotation = ifelse(Gene %in% ribosome$target, "ribosomal", "other")) %>%
  mutate(annotation = ifelse(Gene %in% mito$target, "mitochondrial", annotation)) %>%
  mutate(annotation = ifelse(Gene %in% c("OGT", "OGA"), "OGA|OGT", annotation)) %>%
  filter(annotation == "OGA|OGT") %>%
  filter(contrast != "168h")

data4 <- res_ogt_glyco$limma_results %>%
  mutate(annotation = ifelse(Gene %in% ribosome$target, "ribosomal", "other")) %>%
  mutate(annotation = ifelse(Gene %in% mito$target, "mitochondrial", annotation)) %>%
  mutate(annotation = ifelse(Gene %in% c("OGT", "OGA"), "OGA|OGT", annotation)) %>%
  filter(annotation == "mitochondrial") %>%
  filter(contrast != "168h")

 ggplot(data2, aes(x = logFC, y = -log10(P.Value), colour = annotation)) +
   geom_vline(xintercept = 0) +
   geom_vline(xintercept = c(0.6, -0.6), linetype = 3) +
   geom_point(data = data2, size = 0.5) +
   geom_point(data = data4, size = 0.8) +
   geom_point(data = data1, size = 0.8) +
   geom_point(data = data3, size = 2) +
   geom_label_repel(data = data3, aes(label = Gene), max.overlaps = Inf, size = 2, colour = "black", min.segment.length = 0.2, force_pull = 20) +
   scale_colour_manual(values = c("mitochondrial" = "orange", "ribosomal" = "cyan3", "other" = "grey", "OGA|OGT" = "black")) +
   # scale_colour_manual(values = c("dodgerblue4", "lightgrey", "indianred")) +
   facet_wrap(~contrast, nrow = 1) +
   labs(x = "log2 fold-change against 0h", y = "-log10  p-value") +
   cowplot::theme_cowplot() +
   cowplot::panel_border()

annotated_limma <- res_ogt_glyco$limma_results %>%
  mutate(annotation = ifelse(Gene %in% ribosome$target, "ribosomal", "other")) %>%
  mutate(annotation = ifelse(Gene %in% mito$target, "mitochondrial", annotation)) %>%
  mutate(annotation = ifelse(Gene %in% mito$target & Gene %in% ribosome$target, "both", annotation)) %>% 
  select(annotation, Gene, logFC, adj.P.Val, hit)

#annotated_limma %>% openxlsx::write.xlsx("annotated_DEA.xlsx")

ggplot(
  data2 %>%
    mutate(label = ifelse(Gene %in% c("OGT", "OGA"), Gene, "")) %>%
    mutate(strength = ifelse(Gene %in% c("OGT", "OGA"), "OGA|OGT", direction)),
  aes(x = logFC, y = -log10(P.Value), colour = annotation, label = label)
) +
  geom_vline(xintercept = 0) +
  geom_vline(xintercept = c(0.6, -0.6), linetype = 3) +
  geom_point(aes(size = strength, alpha = strength)) +
  geom_text_repel(max.overlaps = Inf, size = 2, colour = "black", min.segment.length = 0.2, force_pull = 20) +
   scale_colour_manual(values = c("ribosomal" = "red", "other" = "grey")) +
  scale_size_manual(values = c(0.8,0.8, 2, 0.8)) +
   scale_alpha_manual(values = c(0.2,0.2, 1, 0.2)) +
  facet_wrap(~contrast, nrow = 1) +
  labs(x = "log2 fold-change", y = "-log10  p-value") +
  cowplot::theme_cowplot() +
  cowplot::panel_border() +
  guides(size = "none", alpha = "none")

 res_ogt_glyco$limma_results %>%
     mutate(annotation = ifelse(Gene %in% ribosome$target, "ribosomal", "other")) %>%
   filter(contrast != "168h") %>% 
   filter(annotation == "ribosomal") %>% 
   acast(Gene ~ contrast, value.var = "logFC") %>% 
   ComplexHeatmap::Heatmap(show_row_names = F, col = circlize::colorRamp2(c(-0.5, 0, 0.5), c("dodgerblue4", "white","indianred3")),name = "log2FC")

mito <- msigdf_decoupler %>%  filter(source == "GOCC_MITOCHONDRION")
ribosome <- msigdf_decoupler %>%  filter(source == "GOCC_RIBOSOME") %>% 
  filter((target %in% mito$target))

data1 <- res_ogt_glyco$limma_results %>%
  mutate(annotation = ifelse(Gene %in% ribosome$target, "ribosomal+mito", "other")) %>%
  mutate(annotation = ifelse(Gene %in% c("OGT", "OGA"), "OGA|OGT", annotation)) %>%
  filter(annotation == "ribosomal+mito") %>%
  filter(contrast != "168h")

data2 <- res_ogt_glyco$limma_results %>%
  mutate(annotation = ifelse(Gene %in% ribosome$target, "ribosomal+mito", "other")) %>%
  mutate(annotation = ifelse(Gene %in% c("OGT", "OGA"), "OGA|OGT", annotation)) %>%
  filter(annotation == "other") %>%
  filter(contrast != "168h")

data3 <- res_ogt_glyco$limma_results %>%
  mutate(annotation = ifelse(Gene %in% ribosome$target, "ribosomal+mito", "other")) %>%
  mutate(annotation = ifelse(Gene %in% c("OGT", "OGA"), "OGA|OGT", annotation)) %>%
  filter(annotation == "OGA|OGT") %>%
  filter(contrast != "168h")



 ggplot(data2, aes(x = logFC, y = -log10(P.Value), colour = annotation)) +
   geom_vline(xintercept = 0) +
   geom_vline(xintercept = c(0.6, -0.6), linetype = 3) +
   geom_point(data = data2, size = 0.5) +
   geom_point(data = data1, size = 0.8) +
   geom_point(data = data3, size = 2) +
   geom_label_repel(data = data3, aes(label = Gene), max.overlaps = Inf, size = 2, colour = "black", min.segment.length = 0.2, force_pull = 20) +
   scale_colour_manual(values = c("mitochondrial" = "orange", "ribosomal+mito" = "red", "other" = "grey", "OGA|OGT" = "black")) +
   # scale_colour_manual(values = c("dodgerblue4", "lightgrey", "indianred")) +
   facet_wrap(~contrast, nrow = 1) +
   labs(x = "log2 fold-change against 0h", y = "-log10  p-value") +
   cowplot::theme_cowplot() +
   cowplot::panel_border()

O-glycoproteome

1. Load and format data

psm <- read_tsv("psm.tsv")

colnames(psm) <- gsub(" ", "\\.", colnames(psm))

psm_glyco <- psm %>%
  # remove N-glycosylation
  filter(!grepl("N.{1}T|N.{1}S", Peptide)) %>%
  filter(!(grepl("N.$", Peptide) & Next.AA %in% c("S", "T"))) %>%
  filter(!is.na(Total.Glycan.Composition)) %>%
  filter(Total.Glycan.Composition == "HexNAc(1) % 203.0794") %>% 
  select(Protein.ID, Gene, Protein.Description, Peptide, Start_AA = Protein.Start, End_AA = Protein.End, matches("sample")) %>% 
  melt(id.vars = c("Protein.ID", "Gene", "Protein.Description", "Peptide", "Start_AA", "End_AA"),
       variable.name = "sample", value.name = "intensity") %>% 
  inner_join(res_ogt_glyco$TMT_info) %>% 
  mutate(time=factor(time, levels = c("0h", "2h", "4h", "24h", "48h", "168h"))) %>% 
  drop_na()

psm_glyco %>% glimpse

Rows: 1,278
Columns: 10
$ Protein.ID          <chr> "Q9UI10", "P51610", "P51610", "Q9P2N6", "P49848", …
$ Gene                <chr> "EIF2B4", "HCFC1", "HCFC1", "KANSL3", "TAF6", "TAF…
$ Protein.Description <chr> "Translation initiation factor eIF-2B subunit delt…
$ Peptide             <chr> "KGEQGGPPPKASPSTAGETPSGVK", "SPITIITTK", "SPITIITT…
$ Start_AA            <dbl> 119, 794, 794, 857, 480, 480, 253, 771, 771, 564, …
$ End_AA              <dbl> 142, 802, 802, 866, 488, 488, 283, 793, 793, 584, …
$ sample              <chr> "sample-01", "sample-01", "sample-01", "sample-01"…
$ intensity           <dbl> 90266.141, 306371.750, 30600.857, 56383.891, 17486…
$ time                <fct> 0h, 0h, 0h, 0h, 0h, 0h, 0h, 0h, 0h, 0h, 0h, 0h, 0h…
$ replicate           <chr> "rep1", "rep1", "rep1", "rep1", "rep1", "rep1", "r…

psm_glyco %>%  distinct(Protein.ID, Gene) %>% write_tsv("plots/Oglyco_modproteins.tsv")

psm_glyco  %>% write_tsv("Oglyco_peptide_information.tsv")

set.seed(2)
psm_glyco %>% 
  filter(time != "168h") %>% 
  group_by(Gene, Start_AA, time, replicate) %>% 
  summarise(intensity = sum(intensity)) %>% 
  group_by(Gene, Start_AA, time) %>% 
  summarise(intensity = mean(intensity)) %>% 
  group_by(Gene, Start_AA) %>% 
  mutate(scaled_int = scale(log2(intensity))) %>% 
  drop_na() %>% 
  acast(Gene + Start_AA ~ time) %>% 
  #xscale() %>% 
  ComplexHeatmap::Heatmap(name = "scaled\nTMT intensity", cluster_columns = F, col = circlize::colorRamp2(c(-2, 0, 2), c("dodgerblue4", "white","indianred3")), rect_gp = grid::gpar(col = "white", lwd = 1), row_names_max_width = unit(20, "cm"))

psm_glyco %>% 
  filter(time != "168h") %>% 
  group_by(Gene, Start_AA, time, replicate) %>% 
  summarise(intensity = sum(intensity)) %>% 
  group_by(Gene, Start_AA, time) %>% 
  summarise(intensity = mean(intensity)) %>% 
  group_by(Gene, Start_AA) %>% 
  mutate(scaled_int = scale(log2(intensity))) %>% 
  drop_na() %>% 
  ggplot(aes(x=time, y=scaled_int)) +
  ggforce::geom_sina(colour = "grey") +
  geom_boxplot(fill=NA, width = 0.2) +
  ggpubr::stat_compare_means(method = "wilcox.test", ref.group = "0h", label = "p.format")

2. Visualise O-glycoproteins

glycoproteins <- unique(psm_glyco$Gene)

for (p in glycoproteins) {
  p = psm_glyco %>%
    filter(Gene == p) %>%
    ggplot(aes(x = time, y = log2(intensity), colour = Peptide, group = Peptide)) +
    geom_point() +
    geom_smooth(method = 'loess', formula = 'y ~ x', fill = "lightgrey") +
    labs(subtitle = p) +
    theme(legend.position = "bottom")
  plot(p)
}