From 9a8f30bcda55078c5b7fde463ad6e0c1e3c4adf6 Mon Sep 17 00:00:00 2001
From: LDAY26 <le.day@pnnl.gov>
Date: Wed, 29 Apr 2026 12:20:01 -0700
Subject: [PATCH 01/12] Phosphosite raw data cleaning and normalization

---
 .../01_Phosphosite_data_processing.R          | 49 +++++++++++++++++++
 1 file changed, 49 insertions(+)
 create mode 100644 ven_hma_DIA-phosphoproteomics/01_Phosphosite_data_processing/01_Phosphosite_data_processing.R

diff --git a/ven_hma_DIA-phosphoproteomics/01_Phosphosite_data_processing/01_Phosphosite_data_processing.R b/ven_hma_DIA-phosphoproteomics/01_Phosphosite_data_processing/01_Phosphosite_data_processing.R
new file mode 100644
index 0000000..1c78f30
--- /dev/null
+++ b/ven_hma_DIA-phosphoproteomics/01_Phosphosite_data_processing/01_Phosphosite_data_processing.R
@@ -0,0 +1,49 @@
+library(dplyr)
+library(stringr)
+library(openxlsx)
+
+## load phosphosite level data and metadata
+meta <- read.xlsx("PTRC_metadata_Exp28_removesamples.xlsx",check.names=FALSE)
+df <- read.xlsx("PTRC_EXP28_InSilico_DiaNN_phosphosites_90_removesamples.xlsx", check.names=FALSE)
+
+## make sure columns containing intensity values are in numeric
+df <- df%>% mutate(across('PTRC_Exp28_Phos_01':'PTRC_Exp28_Phos_41', as.numeric))
+
+## different uniprot protein accession IDs could be mapped to the same gene name. Create a new SITE2 column that is formatted as GeneName-Residue#, e.g.TADA2A-S6
+df <- df %>% mutate(SITE = str_c(Gene.Names, "-", Residue, Site))
+
+## subset site and sample columns
+df <- df[,9:45]
+## sum rows with same SITE ID
+df <- df %>% group_by(SITE) %>%
+  summarise(SITE=dplyr::first(SITE),
+            across(everything(), sum, na.rm=TRUE))
+
+## replace 0 with NA, followed by log2 transformation
+df[df == 0] <- NA
+df[,2:37] <- log(df[, 2:37], 2)
+# check data distribution prior to median centering
+boxplot(df[,2:37], cex.axis=1, las=2) 
+
+# median centering
+Zero_Center_Norm <- function(df) {
+  med_norm <- function (df) 
+  {
+    norm.coeff <- apply(df, 2, median, na.rm = TRUE)# collect median of each sample from specified dataframe
+    df1 <- sweep(df, 2, norm.coeff, "-") #subtract the median from each respective column in dataframe
+    avg_of_median <- mean(norm.coeff) # calculate average of medians of each sample in group
+    df1 <- df1 + avg_of_median #add average of averages back to each subtracted sample value
+    return(df1)
+  }
+  df <- med_norm(df)
+  return(df)
+}
+
+df[,2:37] <- Zero_Center_Norm(df[,2:37])
+
+# check data distribution after median centering
+boxplot(df[,2:37], cex.axis=1, las=2) 
+# remove sites with NAs across all samples.
+df <- df[rowSums(!is.na(df[ , 2:37])) > 0, ] 
+
+write.xlsx(df, "PTRC_EXP28_InSilico_Cleaned_PhosphositeData.xlsx")

From bc91d0ef5e8305b7c25dd032e683f8efb2465e1b Mon Sep 17 00:00:00 2001
From: LDAY26 <le.day@pnnl.gov>
Date: Wed, 29 Apr 2026 12:33:16 -0700
Subject: [PATCH 02/12] Use Limma to compare patients with different outcomes

---
 .../02_LimmaStats/02_LimmaStats.R             | 30 +++++++++++++++++++
 1 file changed, 30 insertions(+)
 create mode 100644 ven_hma_DIA-phosphoproteomics/02_LimmaStats/02_LimmaStats.R

diff --git a/ven_hma_DIA-phosphoproteomics/02_LimmaStats/02_LimmaStats.R b/ven_hma_DIA-phosphoproteomics/02_LimmaStats/02_LimmaStats.R
new file mode 100644
index 0000000..dbf2ac0
--- /dev/null
+++ b/ven_hma_DIA-phosphoproteomics/02_LimmaStats/02_LimmaStats.R
@@ -0,0 +1,30 @@
+library(dplyr)
+library(MSnSet.utils)
+library(openxlsx)
+
+## load log2 transformed, median centered phosphosite level data 
+df <- read.xlsx("PTRC_EXP28_InSilico_Cleaned_PhosphositeData.xlsx", check.names=FALSE)
+meta <- read.xlsx("PTRC_metadata_Exp28_removesamples.xlsx",check.names=FALSE)
+
+
+## build MSnSet
+exprs <- df %>% arrange(., SITE) %>% tibble::column_to_rownames(var="SITE") %>% as.matrix()
+fData <- data.frame(rownames(exprs)) %>% tibble::column_to_rownames(var="rownames.exprs.")
+exprs <- exprs[ , paste(meta$Sample, sep = "")]
+pData <- meta %>% mutate(SampleID = Sample) %>% tibble::column_to_rownames(var="Sample")
+
+dfSet <- MSnSet(exprs = exprs,
+                pData = pData, fData = fData)
+
+# set up contrasts to compare 1) responders VS nonresponders 2) no-relapse VS refractory, relapse VS refractory, relapse VS no-relapse
+contrasts1 <- c("ResponseGroupResponders-ResponseGroupNonResponders")
+contrasts2 <- c("subcohortNorelapse-subcohortRefractory", "subcohortRelapse-subcohortRefractory", "subcohortRelapse-subcohortNorelapse")
+
+tests1 <- limma_contrasts(eset = dfSet, model.str = "~ 0 + ResponseGroup", coef.str = "ResponseGroup",
+                          contrasts = contrasts1, trend = TRUE, robust = TRUE, plot = FALSE)
+
+tests2 <- limma_contrasts(eset = dfSet, model.str = "~ 0 + subcohort", coef.str = "subcohort",
+                          contrasts = contrasts2, trend = TRUE, robust = TRUE, plot = FALSE)
+
+write.xlsx(tests1, "PTRC_EXP28_Responder VS NonResponder.xlsx", rowNames=FALSE)
+write.xlsx(tests2, "PTRC_EXP28_Refractory VS relapse and no-relapse.xlsx", rowNames=FALSE)

From ab4d638d063fe957fb7d132d71192637ae897352 Mon Sep 17 00:00:00 2001
From: LDAY26 <le.day@pnnl.gov>
Date: Wed, 29 Apr 2026 13:04:35 -0700
Subject: [PATCH 03/12] Kinase Substrate Enrichment Analysis and associated
 plot

---
 .../03_KSEA_and_KSEAplot.R                    | 61 +++++++++++++++++++
 1 file changed, 61 insertions(+)
 create mode 100644 ven_hma_DIA-phosphoproteomics/03_KSEA_and_KSEAplot/03_KSEA_and_KSEAplot.R

diff --git a/ven_hma_DIA-phosphoproteomics/03_KSEA_and_KSEAplot/03_KSEA_and_KSEAplot.R b/ven_hma_DIA-phosphoproteomics/03_KSEA_and_KSEAplot/03_KSEA_and_KSEAplot.R
new file mode 100644
index 0000000..26d0752
--- /dev/null
+++ b/ven_hma_DIA-phosphoproteomics/03_KSEA_and_KSEAplot/03_KSEA_and_KSEAplot.R
@@ -0,0 +1,61 @@
+library(dplyr)
+library(ggplot2)
+library(openxlsx)
+
+
+# load Limma result
+df <- read.xlsx("PTRC_EXP28_Refractory VS relapse and no-relapse.xlsx", check.names=FALSE)
+## load kinase substrate database
+KSDB <- read.csv('PSP&NetworKIN_Kinase_Substrate_Dataset_July2016 1.csv', stringsAsFactors = FALSE)
+
+## parse out 3 comparisons from Limma output file
+NRRef <- df %>% filter(contrast == "Norelapse-Refractory")   
+RRef  <- df %>% filter(contrast == "Relapse-Refractory")  
+NRR   <- df %>% filter(contrast == "Relapse-Norelapse")  
+
+## subject NRRef, RRef, NRR to following codes for KSEA analysis and plotting.
+## prepare input for KSEA
+fold_change <- NRRef$logFC
+fold_change <- 2**fold_change
+
+PhosInp <- data.frame(Protein = "NULL", Gene = NRRef$feature, Peptide = "NULL", 
+                      Residue.Both = NRRef$feature, p = "NULL", FC = fold_change) %>%
+  dplyr::mutate(Residue.Both = sub("^.*-", "", Residue.Both)) %>%
+  dplyr::mutate(Gene = sub("^(.*)-[^-]*$", "\\1", Gene))
+
+ksea_res_full <- KSEAapp::KSEA.Scores(KSDB, PhosInp, NetworKIN = TRUE, NetworKIN.cutoff = Inf) 
+
+ksea_res <- ksea_res_full %>%
+  dplyr::select(Kinase.Gene, m, p.value, FDR, z.score) %>%
+  dplyr::rename(kinase = Kinase.Gene, z_score = z.score, p_value = p.value,
+                adj_p_val = FDR, site_size = m)
+
+# filter out kinases enriched by >=5 phosphosites and have p-value <0.05
+kinase <- ksea_res %>% filter(site_size >=5 & p_value < 0.05)
+
+kinase <- kinase %>%
+  arrange(desc(z_score)) %>%
+  mutate(kinase=factor(kinase, levels=kinase))
+
+
+# plot KSEA results with lollipop plot
+ggplot(kinase, aes(x = z_score, y = kinase)) +
+  geom_segment(aes(x = 0, xend = z_score, y = kinase, yend = kinase),
+               size=1.5) +
+  # points, size by m, color by FDR
+  geom_point(aes(size = site_size, colour = p_value)) +
+  scale_size(range = c(7, 12), name = "phosphosite\nsubstrates") +
+  scale_colour_viridis_c(option = "plasma", direction = -1, name = "p-value") +
+  geom_vline(xintercept=0, color="black", size=1.5)+
+  labs(x = "z-score", y = NULL, title = "") +
+  theme_minimal(base_size = 12) +
+  theme(
+    panel.grid.major.y = element_line(color="grey90", linewidth=0.3),
+    panel.grid.minor = element_line(color="grey90", linewidth=0.3),
+    axis.title.x=element_text(size=30, face="bold"),
+    axis.text.y=element_text(size=27, face="bold", color="black"),
+    axis.text.x=element_text(size=30, face="bold", color="black"),
+    legend.title=element_text(size=22, face="bold"),
+    legend.text=element_text(size=20)
+  )
+

From d2a27b1089ee157b5d35ab55ae39b026eef5526c Mon Sep 17 00:00:00 2001
From: LDAY26 <le.day@pnnl.gov>
Date: Wed, 29 Apr 2026 13:07:01 -0700
Subject: [PATCH 04/12] Update 01_Phosphosite_data_processing.R

---
 .../01_Phosphosite_data_processing.R                         | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/ven_hma_DIA-phosphoproteomics/01_Phosphosite_data_processing/01_Phosphosite_data_processing.R b/ven_hma_DIA-phosphoproteomics/01_Phosphosite_data_processing/01_Phosphosite_data_processing.R
index 1c78f30..b4a2ddf 100644
--- a/ven_hma_DIA-phosphoproteomics/01_Phosphosite_data_processing/01_Phosphosite_data_processing.R
+++ b/ven_hma_DIA-phosphoproteomics/01_Phosphosite_data_processing/01_Phosphosite_data_processing.R
@@ -2,8 +2,7 @@ library(dplyr)
 library(stringr)
 library(openxlsx)
 
-## load phosphosite level data and metadata
-meta <- read.xlsx("PTRC_metadata_Exp28_removesamples.xlsx",check.names=FALSE)
+## load phosphosite level data 
 df <- read.xlsx("PTRC_EXP28_InSilico_DiaNN_phosphosites_90_removesamples.xlsx", check.names=FALSE)
 
 ## make sure columns containing intensity values are in numeric
@@ -46,4 +45,4 @@ boxplot(df[,2:37], cex.axis=1, las=2)
 # remove sites with NAs across all samples.
 df <- df[rowSums(!is.na(df[ , 2:37])) > 0, ] 
 
-write.xlsx(df, "PTRC_EXP28_InSilico_Cleaned_PhosphositeData.xlsx")
+write.xlsx(df, "PTRC_EXP28_InSilico_Cleaned_PhosphositeData.xlsx", rowNames=FALSE)

From 4b01f48bd2352e3c4614b782c3c39209f0e80afe Mon Sep 17 00:00:00 2001
From: LDAY26 <le.day@pnnl.gov>
Date: Wed, 29 Apr 2026 13:23:56 -0700
Subject: [PATCH 05/12] Gene Set Enrichment Analysis and associated plot

---
 .../04_GSEA_and_GSEAplot.R                    | 88 +++++++++++++++++++
 1 file changed, 88 insertions(+)
 create mode 100644 ven_hma_DIA-phosphoproteomics/04_GSEA_and_GSEAplot/04_GSEA_and_GSEAplot.R

diff --git a/ven_hma_DIA-phosphoproteomics/04_GSEA_and_GSEAplot/04_GSEA_and_GSEAplot.R b/ven_hma_DIA-phosphoproteomics/04_GSEA_and_GSEAplot/04_GSEA_and_GSEAplot.R
new file mode 100644
index 0000000..ce5290f
--- /dev/null
+++ b/ven_hma_DIA-phosphoproteomics/04_GSEA_and_GSEAplot/04_GSEA_and_GSEAplot.R
@@ -0,0 +1,88 @@
+library(dplyr)
+library(ggplot2)
+library(openxlsx)
+library(msigdbr)
+library(fgsea)
+library(scales)
+
+
+# load Limma result
+df <- read.xlsx("PTRC_EXP28_Refractory VS relapse and no-relapse.xlsx", check.names=FALSE)
+
+## parse out 3 comparisons from Limma output file
+NRRef <- df %>% filter(contrast == "Norelapse-Refractory")   
+RRef  <- df %>% filter(contrast == "Relapse-Refractory")  
+NRR   <- df %>% filter(contrast == "Relapse-Norelapse")  
+
+# suject NRRef, RRef, NRR to the following code for GSEA
+# prepare input for GSEA
+names(NRRef)[names(NRRef) == "logFC"] <- "log2FC"
+names(NRRef)[names(NRRef) == "P.Value"] <- "pvalue"
+NRRef$gene <- sub("-.*", "", NRRef$feature)
+
+gene_level_NRRef <- NRRef[,c(1,4,9)] %>%
+  filter(!is.na(log2FC), !is.na(pvalue), !is.na(gene)) %>%
+  group_by(gene) %>%
+  slice_max(order_by = abs(log2FC), n = 1, with_ties = FALSE) %>%
+  ungroup()
+
+# Build ranking statistic for GSEA: sign(log2FC) * -log10(p)
+gene_level_NRRef <- gene_level_NRRef %>%
+  mutate(rank_stat = sign(log2FC) * -log10(pvalue))
+
+# Create named numeric vector, sorted decreasing
+ranks_NRRef <- gene_level_NRRef$rank_stat
+names(ranks_NRRef) <- gene_level_NRRef$gene
+ranks_NRRef <- sort(ranks_NRRef, decreasing = TRUE)
+
+# Hallmark pathways (H collection), human
+m_df <- msigdbr(species = "Homo sapiens", collection = "H")
+
+# Convert to a list: names = pathway, each element = vector of genes
+pathways <- split(m_df$gene_symbol, m_df$gs_name)
+
+set.seed(123)
+
+fgsea_res_NRRef <- fgsea(
+  pathways = pathways,
+  stats    = ranks_NRRef,
+  minSize  = 10,
+  maxSize  = 500,
+  nperm    = 10000
+)
+
+# calculate GeneRatio for plotting purpose
+fgsea_res_NRRef <- fgsea_res_NRRef %>%
+  mutate(
+    GeneRatio = lengths(leadingEdge) / size,
+    Count     = lengths(leadingEdge)       
+  )
+
+# filter out pathways with p-value <0.05 
+pfgsea <- fgsea_res_NRRef %>% filter(pval <0.05) %>% arrange(desc(GeneRatio))
+# clean up pathway labels
+pfgsea$pathway <-sub("^HALLMARK_", "", pfgsea$pathway)
+pfgsea$pathway <-gsub("_", " ", pfgsea$pathway)
+pfgsea$pathway <- factor(pfgsea$pathway, levels=rev(pfgsea$pathway))
+
+ggplot(pfgsea, aes(x = GeneRatio, y = pathway)) +
+  geom_point(aes(size = -log10(pval), color = NES)) +
+  scale_y_discrete(labels = label_wrap(20))+
+  scale_color_gradient2(low="blue", mid="white", high="red", midpoint=0, name="NES") +
+  scale_size_continuous(name = "-log10(pval)", range=c(6,10)) +
+  labs(
+    x = "GeneRatio",
+    y = NULL,
+    title = ""
+  ) +
+  theme_bw() +
+  theme(
+    axis.text.y = element_text(size = 22, face="bold", color="black"),
+    plot.title  = element_text(hjust = 0.5),
+    axis.title.x=element_text(size=30, face="bold"),
+    axis.title.y=element_text(size=27, face="bold", color="black"),
+    axis.text.x=element_text(size=27, face="bold", color="black"),
+    legend.title=element_text(size=22, face="bold"),
+    legend.text=element_text(size=20),
+    plot.margin = unit(c(1, 1, 1, 1), "cm")
+  )

From 1019156c34e25d29fee7ebfe18c7f4eadfb3a2b7 Mon Sep 17 00:00:00 2001
From: LDAY26 <le.day@pnnl.gov>
Date: Wed, 29 Apr 2026 14:19:20 -0700
Subject: [PATCH 06/12] Analysis on Human:Ecoli phosphopeptide dilution
 experiment

---
 .../Dilution_experiment/Dilution_experiment.R | 134 ++++++++++++++++++
 1 file changed, 134 insertions(+)
 create mode 100644 ven_hma_DIA-phosphoproteomics/Dilution_experiment/Dilution_experiment.R

diff --git a/ven_hma_DIA-phosphoproteomics/Dilution_experiment/Dilution_experiment.R b/ven_hma_DIA-phosphoproteomics/Dilution_experiment/Dilution_experiment.R
new file mode 100644
index 0000000..ec356e7
--- /dev/null
+++ b/ven_hma_DIA-phosphoproteomics/Dilution_experiment/Dilution_experiment.R
@@ -0,0 +1,134 @@
+library(ggplot2)
+library(dplyr)
+library(tidyr)
+library(openxlsx)
+library(stringr)
+library(ggside)
+
+# load dilution experiment phosphosite level data generated through in silico, combined, or hybrid searches
+df <- read.xlsx("PTRC_insilico_site_human.xlsx", check.names=FALSE)
+
+## make sure columns containing intensity values are in numeric
+df[,8:22] <-lapply(df[,8:22], as.numeric) 
+
+## different uniprot protein accession IDs could be mapped to the same gene name. Create a new SITE2 column that is formatted as GeneName-Residue#, e.g.TADA2A-S6
+df <- df %>% mutate(SITE = str_c(Gene.Names, "-", Residue, Site))
+
+## subset site and sample columns
+df <- df[,7:22]
+## sum rows with same SITE ID
+df <- df %>% group_by(SITE) %>%
+  summarise(SITE=dplyr::first(SITE),
+            across(everything(), sum, na.rm=TRUE))
+
+## log2 transform
+df[df == 0] <- NA
+df[,2:16] <- log(df[, 2:16], 2)  
+## check data distribution
+boxplot(df[,2:16], cex.axis=1, las=2)
+
+Zero_Center_Norm <- function(df) {
+  med_norm <- function (df) 
+  {
+    norm.coeff <- apply(df, 2, median, na.rm = TRUE)# collect median of each sample from specified dataframe
+    df1 <- sweep(df, 2, norm.coeff, "-") #subtract the median from each respective column in dataframe
+    avg_of_median <- mean(norm.coeff) # calculate average of medians of each sample in group
+    df1 <- df1 + avg_of_median #add average of averages back to each subtracted sample value
+    return(df1)
+  }
+  df <- med_norm(df)
+  return(df)
+}
+
+df[,2:16] <- Zero_Center_Norm(df[,2:16])
+
+## check data distribution after median centering
+boxplot(df[,2:16], cex.axis=1, las=2)
+
+## subset 3 replicates from each dilution point. Focus on sites identified across all 3 replicates. Calculate the mean across replicates per site. 
+df_D1 <- df[,c(1,2,7,12)]
+df_D1$CountD1 <- apply(df_D1[,2:4], 1, function(x) sum(is.na(x)))
+df_D1 <- filter(df_D1, CountD1 == 0)
+df_D1$D1 <- rowMeans(df_D1[,2:4])
+
+
+df_D2 <- df[,c(1,3,8,13)]
+df_D2$CountD2 <- apply(df_D2[,2:4], 1, function(x) sum(is.na(x)))
+df_D2 <- filter(df_D2, CountD2 ==0)
+df_D2$D2 <- rowMeans(df_D2[,2:4])
+
+df_D3 <- df[,c(1,4,9,14)]
+df_D3$CountD3 <- apply(df_D3[,2:4], 1, function(x) sum(is.na(x)))
+df_D3 <- filter(df_D3, CountD3 == 0)
+df_D3$D3 <- rowMeans(df_D3[,2:4])
+
+df_D4 <- df[,c(1,5,10,15)]
+df_D4$CountD4 <- apply(df_D4[,2:4], 1, function(x) sum(is.na(x)))
+df_D4 <- filter(df_D4, CountD4 == 0)
+df_D4$D4 <- rowMeans(df_D4[,2:4])
+
+df_D5 <- df[,c(1,6,11,16)]
+df_D5$CountD5 <- apply(df_D5[,2:4], 1, function(x) sum(is.na(x)))
+df_D5 <- filter(df_D5, CountD5 == 0)
+df_D5$D5 <- rowMeans(df_D5[,2:4])
+
+## merge sites identified across all dilution points.
+M2 <- merge(df_D1[,c(1,6)], df_D2[,c(1,6)], by="SITE")
+M3 <- merge(M2, df_D3[,c(1,6)], by ="SITE")
+M4 <- merge(M3, df_D4[,c(1,6)], by ="SITE")
+M5 <- merge(M4, df_D5[,c(1,6)], by ="SITE")
+
+## use middle dilution point 50:50 as reference. Calculate log2FC between each dilution and the 50:50 reference.
+M5$logFC_D1 <- M5$D1-M5$D3
+M5$logFC_D2 <- M5$D2-M5$D3
+M5$logFC_D3 <- M5$D3-M5$D3
+M5$logFC_D4 <- M5$D4-M5$D3
+M5$logFC_D5 <- M5$D5-M5$D3
+
+# quick overview of log2FC
+boxplot(M5[,7:11], cex.axis=1, las=2)
+
+# pivot data longer for ggplot 
+M5 <- M5[,c(1,7:11)]
+dfp <- M5 %>% pivot_longer(cols= 2:6,
+                           names_to="Dilution",
+                           values_to="intensity")
+## plot log2FC 
+ggplot(dfp, aes(x = Dilution, y = intensity)) +
+  geom_hline(yintercept=c(1,0.585,0,-1,-3.322), color=c("red","khaki4", "springgreen3", "deepskyblue","magenta"), linetype="dashed", linewidth=1)+
+  geom_boxplot(
+    aes(color=Dilution),
+    fill=NA,
+    position=position_dodge(width=0.6),
+    outlier.shape=NA,
+    width=0.5,
+    lwd=1)+
+  geom_point(
+    aes(color=Dilution),
+    position=position_jitterdodge(
+      jitter.width = 0.3, 
+      dodge.width = 0.6
+    ), 
+    alpha = 0.03, size = 0.3) +
+  scale_x_discrete(
+    name = "Dilution",
+    labels= c("logFC_D1"= "100:0",
+              "logFC_D2"= "75:25",
+              "logFC_D3"="50:50",
+              "logFC_D4"="25:75",
+              "logFC_D5"="5:95"))+
+  labs(y = "logFC",
+       color = "Dilution"
+  )+
+  theme(panel.grid.major = element_blank(),
+        panel.grid.minor = element_blank(),
+        panel.border     = element_blank(),
+        panel.background = element_blank(),
+        axis.line = element_line(color="black", size=1),
+        axis.title.x=element_text(size=30, face="bold"),
+        axis.title.y=element_text(size=30, face="bold"),
+        axis.text.y=element_text(size=27, face="bold", color="black"),
+        axis.text.x=element_text(size=30, face="bold", color="black"),
+        legend.title=element_text(size=27, face="bold"),
+        legend.text=element_text(size=20))
+

From 35ebf6dd693377ebe8b4638fccf2300924e7ef95 Mon Sep 17 00:00:00 2001
From: LDAY26 <le.day@pnnl.gov>
Date: Tue, 9 Jun 2026 13:25:28 -0700
Subject: [PATCH 07/12] Update 03_KSEA_and_KSEAplot.R

---
 .../03_KSEA_and_KSEAplot.R                    | 100 ++++++++++++++++--
 1 file changed, 93 insertions(+), 7 deletions(-)

diff --git a/ven_hma_DIA-phosphoproteomics/03_KSEA_and_KSEAplot/03_KSEA_and_KSEAplot.R b/ven_hma_DIA-phosphoproteomics/03_KSEA_and_KSEAplot/03_KSEA_and_KSEAplot.R
index 26d0752..5ef723b 100644
--- a/ven_hma_DIA-phosphoproteomics/03_KSEA_and_KSEAplot/03_KSEA_and_KSEAplot.R
+++ b/ven_hma_DIA-phosphoproteomics/03_KSEA_and_KSEAplot/03_KSEA_and_KSEAplot.R
@@ -1,19 +1,60 @@
 library(dplyr)
 library(ggplot2)
 library(openxlsx)
+library(synapser)
 
+## synapse login with .Renviron
+synLogin()
 
-# load Limma result
-df <- read.xlsx("PTRC_EXP28_Refractory VS relapse and no-relapse.xlsx", check.names=FALSE)
-## load kinase substrate database
-KSDB <- read.csv('PSP&NetworKIN_Kinase_Substrate_Dataset_July2016 1.csv', stringsAsFactors = FALSE)
+## define folder ID and the target files
+folder_id <-"syn68733653"
+file1_name <- "PTRC_EXP28_Phospho_Stats_Results.xlsx"
+file2_name <- "PTRC_EXP28_KSEA_Dataset_July2016 1.csv"
 
-## parse out 3 comparisons from Limma output file
+## find target files in this folder
+children <- synGetChildren(folder_id)
+children_list <- as.list(children)
+
+file1_id <- NULL
+file2_id <- NULL
+
+for (child in children_list) {
+  if (child$name == file1_name) {
+    file1_id <- child$id
+  } else if (child$name == file2_name) {
+    file2_id <- child$id
+  }
+}
+
+## download and load files
+if (!is.null(file1_id)) {
+  cat("Downloading", file1_name, "...\n")
+  file1_entity <- synGet(file1_id)
+  
+  # Load using openxlsx function
+  df <- read.xlsx(file1_entity$path)
+  cat("Successfully loaded 'df_phospho_stats'!\n")
+} else {
+  cat("Error: Could not find", file1_name, "in the folder.\n")
+}
+
+if (!is.null(file2_id)) {
+  cat("Downloading", file2_name, "...\n")
+  file2_entity <- synGet(file2_id)
+  
+  # Load into data frame
+  KSDB <- read.csv(file2_entity$path)
+  cat("Successfully loaded 'KSDB'!\n")
+} else {
+  cat("Error: Could not find", file2_name, "in the folder.\n")
+}
+
+
+## parse out 2 comparisons from Limma output file
 NRRef <- df %>% filter(contrast == "Norelapse-Refractory")   
 RRef  <- df %>% filter(contrast == "Relapse-Refractory")  
-NRR   <- df %>% filter(contrast == "Relapse-Norelapse")  
 
-## subject NRRef, RRef, NRR to following codes for KSEA analysis and plotting.
+## subject NRRef to following codes for KSEA analysis and plotting.
 ## prepare input for KSEA
 fold_change <- NRRef$logFC
 fold_change <- 2**fold_change
@@ -59,3 +100,48 @@ ggplot(kinase, aes(x = z_score, y = kinase)) +
     legend.text=element_text(size=20)
   )
 
+## subject RRef to following codes for KSEA analysis and plotting.
+## prepare input for KSEA
+fold_change <- RRef$logFC
+fold_change <- 2**fold_change
+
+PhosInp <- data.frame(Protein = "NULL", Gene = RRef$feature, Peptide = "NULL", 
+                      Residue.Both = RRef$feature, p = "NULL", FC = fold_change) %>%
+  dplyr::mutate(Residue.Both = sub("^.*-", "", Residue.Both)) %>%
+  dplyr::mutate(Gene = sub("^(.*)-[^-]*$", "\\1", Gene))
+
+ksea_res_full <- KSEAapp::KSEA.Scores(KSDB, PhosInp, NetworKIN = TRUE, NetworKIN.cutoff = Inf) 
+
+ksea_res <- ksea_res_full %>%
+  dplyr::select(Kinase.Gene, m, p.value, FDR, z.score) %>%
+  dplyr::rename(kinase = Kinase.Gene, z_score = z.score, p_value = p.value,
+                adj_p_val = FDR, site_size = m)
+
+# filter out kinases enriched by >=5 phosphosites and have p-value <0.05
+kinase <- ksea_res %>% filter(site_size >=5 & p_value < 0.05)
+
+kinase <- kinase %>%
+  arrange(desc(z_score)) %>%
+  mutate(kinase=factor(kinase, levels=kinase))
+
+
+# plot KSEA results with lollipop plot
+ggplot(kinase, aes(x = z_score, y = kinase)) +
+  geom_segment(aes(x = 0, xend = z_score, y = kinase, yend = kinase),
+               size=1.5) +
+  # points, size by m, color by FDR
+  geom_point(aes(size = site_size, colour = p_value)) +
+  scale_size(range = c(7, 12), name = "phosphosite\nsubstrates") +
+  scale_colour_viridis_c(option = "plasma", direction = -1, name = "p-value") +
+  geom_vline(xintercept=0, color="black", size=1.5)+
+  labs(x = "z-score", y = NULL, title = "") +
+  theme_minimal(base_size = 12) +
+  theme(
+    panel.grid.major.y = element_line(color="grey90", linewidth=0.3),
+    panel.grid.minor = element_line(color="grey90", linewidth=0.3),
+    axis.title.x=element_text(size=30, face="bold"),
+    axis.text.y=element_text(size=27, face="bold", color="black"),
+    axis.text.x=element_text(size=30, face="bold", color="black"),
+    legend.title=element_text(size=22, face="bold"),
+    legend.text=element_text(size=20)
+  )
\ No newline at end of file

From 7019974306e80d0611938bd6b6771e5ea8557503 Mon Sep 17 00:00:00 2001
From: LDAY26 <le.day@pnnl.gov>
Date: Tue, 9 Jun 2026 16:06:21 -0700
Subject: [PATCH 08/12] Update 01_Phosphosite_data_processing.R

---
 .../01_Phosphosite_data_processing.R               | 14 +++++++++++---
 1 file changed, 11 insertions(+), 3 deletions(-)

diff --git a/ven_hma_DIA-phosphoproteomics/01_Phosphosite_data_processing/01_Phosphosite_data_processing.R b/ven_hma_DIA-phosphoproteomics/01_Phosphosite_data_processing/01_Phosphosite_data_processing.R
index b4a2ddf..dd74aae 100644
--- a/ven_hma_DIA-phosphoproteomics/01_Phosphosite_data_processing/01_Phosphosite_data_processing.R
+++ b/ven_hma_DIA-phosphoproteomics/01_Phosphosite_data_processing/01_Phosphosite_data_processing.R
@@ -1,12 +1,20 @@
 library(dplyr)
+library(MSnSet.utils)
 library(stringr)
 library(openxlsx)
+library(synapser)
 
-## load phosphosite level data 
-df <- read.xlsx("PTRC_EXP28_InSilico_DiaNN_phosphosites_90_removesamples.xlsx", check.names=FALSE)
+## synapse login with .Renviron
+synLogin()
+
+file_id <- synFindEntityId("PTRC_EXP28_InSilico_DiaNN_phosphosites_90.xlsx", parent = "syn68733653")
+
+# Download and read the file
+file_entity <- synGet(file_id)
+df <- read.xlsx(file_entity$path)
 
 ## make sure columns containing intensity values are in numeric
-df <- df%>% mutate(across('PTRC_Exp28_Phos_01':'PTRC_Exp28_Phos_41', as.numeric))
+df <- df%>% mutate(across('Phos_01':'Phos_41', as.numeric))
 
 ## different uniprot protein accession IDs could be mapped to the same gene name. Create a new SITE2 column that is formatted as GeneName-Residue#, e.g.TADA2A-S6
 df <- df %>% mutate(SITE = str_c(Gene.Names, "-", Residue, Site))

From a8ae445fdfc79cfa15aeecdfb68267be5821e299 Mon Sep 17 00:00:00 2001
From: LDAY26 <le.day@pnnl.gov>
Date: Tue, 9 Jun 2026 16:18:53 -0700
Subject: [PATCH 09/12] Update 03_KSEA_and_KSEAplot.R

---
 .../03_KSEA_and_KSEAplot.R                    | 47 +++----------------
 1 file changed, 7 insertions(+), 40 deletions(-)

diff --git a/ven_hma_DIA-phosphoproteomics/03_KSEA_and_KSEAplot/03_KSEA_and_KSEAplot.R b/ven_hma_DIA-phosphoproteomics/03_KSEA_and_KSEAplot/03_KSEA_and_KSEAplot.R
index 5ef723b..c596bfc 100644
--- a/ven_hma_DIA-phosphoproteomics/03_KSEA_and_KSEAplot/03_KSEA_and_KSEAplot.R
+++ b/ven_hma_DIA-phosphoproteomics/03_KSEA_and_KSEAplot/03_KSEA_and_KSEAplot.R
@@ -8,47 +8,14 @@ synLogin()
 
 ## define folder ID and the target files
 folder_id <-"syn68733653"
-file1_name <- "PTRC_EXP28_Phospho_Stats_Results.xlsx"
-file2_name <- "PTRC_EXP28_KSEA_Dataset_July2016 1.csv"
-
-## find target files in this folder
-children <- synGetChildren(folder_id)
-children_list <- as.list(children)
-
-file1_id <- NULL
-file2_id <- NULL
-
-for (child in children_list) {
-  if (child$name == file1_name) {
-    file1_id <- child$id
-  } else if (child$name == file2_name) {
-    file2_id <- child$id
-  }
-}
-
-## download and load files
-if (!is.null(file1_id)) {
-  cat("Downloading", file1_name, "...\n")
-  file1_entity <- synGet(file1_id)
-  
-  # Load using openxlsx function
-  df <- read.xlsx(file1_entity$path)
-  cat("Successfully loaded 'df_phospho_stats'!\n")
-} else {
-  cat("Error: Could not find", file1_name, "in the folder.\n")
-}
-
-if (!is.null(file2_id)) {
-  cat("Downloading", file2_name, "...\n")
-  file2_entity <- synGet(file2_id)
-  
-  # Load into data frame
-  KSDB <- read.csv(file2_entity$path)
-  cat("Successfully loaded 'KSDB'!\n")
-} else {
-  cat("Error: Could not find", file2_name, "in the folder.\n")
-}
 
+file1_id <- synFindEntityId("PTRC_EXP28_Phospho_Stats_Results.xlsx", parent = folder_id)
+file1_entity <- synGet(file1_id)
+df <- read.xlsx(file1_entity$path)
+
+file2_id <- synFindEntityId("PTRC_EXP28_KSEA_Dataset_July2016 1.csv", parent = folder_id)
+file2_entity <- synGet(file2_id)
+KSDB <- read.csv(file2_entity$path)
 
 ## parse out 2 comparisons from Limma output file
 NRRef <- df %>% filter(contrast == "Norelapse-Refractory")   

From 33cdd962fd54811089b17575ee52ff683d855a26 Mon Sep 17 00:00:00 2001
From: LDAY26 <le.day@pnnl.gov>
Date: Tue, 9 Jun 2026 16:23:53 -0700
Subject: [PATCH 10/12] Update 02_LimmaStats.R

---
 .../02_LimmaStats/02_LimmaStats.R                 | 15 +++++++++++++--
 1 file changed, 13 insertions(+), 2 deletions(-)

diff --git a/ven_hma_DIA-phosphoproteomics/02_LimmaStats/02_LimmaStats.R b/ven_hma_DIA-phosphoproteomics/02_LimmaStats/02_LimmaStats.R
index dbf2ac0..ea8701c 100644
--- a/ven_hma_DIA-phosphoproteomics/02_LimmaStats/02_LimmaStats.R
+++ b/ven_hma_DIA-phosphoproteomics/02_LimmaStats/02_LimmaStats.R
@@ -1,11 +1,22 @@
 library(dplyr)
 library(MSnSet.utils)
 library(openxlsx)
+library(synapser)
+
+## synapse login with .Renviron
+synLogin()
+
+## define folder ID and the target files
+folder_id <-"syn68733653"
 
 ## load log2 transformed, median centered phosphosite level data 
-df <- read.xlsx("PTRC_EXP28_InSilico_Cleaned_PhosphositeData.xlsx", check.names=FALSE)
-meta <- read.xlsx("PTRC_metadata_Exp28_removesamples.xlsx",check.names=FALSE)
+file1_id <- synFindEntityId("PTRC_EXP28_InSilico_Cleaned_PhosphositeData.xlsx", parent = folder_id)
+file1_entity <- synGet(file1_id)
+df <- read.xlsx(file1_entity$path)
 
+file2_id <- synFindEntityId("PTRC_metadata_Exp28_removesamples.xlsx", parent = folder_id)
+file2_entity <- synGet(file2_id)
+meta <- read.xlsx(file2_entity$path)
 
 ## build MSnSet
 exprs <- df %>% arrange(., SITE) %>% tibble::column_to_rownames(var="SITE") %>% as.matrix()

From ed305aff96084ca3c9ae867b60d85d0bf152926c Mon Sep 17 00:00:00 2001
From: LDAY26 <le.day@pnnl.gov>
Date: Tue, 9 Jun 2026 16:29:13 -0700
Subject: [PATCH 11/12] Update 04_GSEA_and_GSEAplot.R

---
 .../04_GSEA_and_GSEAplot/04_GSEA_and_GSEAplot.R  | 16 +++++++++++-----
 1 file changed, 11 insertions(+), 5 deletions(-)

diff --git a/ven_hma_DIA-phosphoproteomics/04_GSEA_and_GSEAplot/04_GSEA_and_GSEAplot.R b/ven_hma_DIA-phosphoproteomics/04_GSEA_and_GSEAplot/04_GSEA_and_GSEAplot.R
index ce5290f..32caddc 100644
--- a/ven_hma_DIA-phosphoproteomics/04_GSEA_and_GSEAplot/04_GSEA_and_GSEAplot.R
+++ b/ven_hma_DIA-phosphoproteomics/04_GSEA_and_GSEAplot/04_GSEA_and_GSEAplot.R
@@ -4,17 +4,23 @@ library(openxlsx)
 library(msigdbr)
 library(fgsea)
 library(scales)
+library(synapser)
 
+## synapse login with .Renviron
+synLogin()
 
-# load Limma result
-df <- read.xlsx("PTRC_EXP28_Refractory VS relapse and no-relapse.xlsx", check.names=FALSE)
+## define folder ID and the target files
+folder_id <-"syn68733653"
 
-## parse out 3 comparisons from Limma output file
+file_id <- synFindEntityId("PTRC_EXP28_Phospho_Stats_Results.xlsx", parent = folder_id)
+file_entity <- synGet(file_id)
+df <- read.xlsx(file_entity$path)
+
+## parse out 2 comparisons from Limma output file
 NRRef <- df %>% filter(contrast == "Norelapse-Refractory")   
 RRef  <- df %>% filter(contrast == "Relapse-Refractory")  
-NRR   <- df %>% filter(contrast == "Relapse-Norelapse")  
 
-# suject NRRef, RRef, NRR to the following code for GSEA
+# suject NRRef, RRef to the following code for GSEA
 # prepare input for GSEA
 names(NRRef)[names(NRRef) == "logFC"] <- "log2FC"
 names(NRRef)[names(NRRef) == "P.Value"] <- "pvalue"

From f942673423251763faf497a274b32a430a1e9353 Mon Sep 17 00:00:00 2001
From: Sara JC Gosline <sara.gosline@pnnl.gov>
Date: Wed, 10 Jun 2026 10:08:33 -0700
Subject: [PATCH 12/12] added additional plotting ideas.

---
 .../03_KSEA_and_KSEAplot.R                    | 94 ++++++++++++++++++-
 1 file changed, 91 insertions(+), 3 deletions(-)

diff --git a/ven_hma_DIA-phosphoproteomics/03_KSEA_and_KSEAplot/03_KSEA_and_KSEAplot.R b/ven_hma_DIA-phosphoproteomics/03_KSEA_and_KSEAplot/03_KSEA_and_KSEAplot.R
index 5ef723b..6a8e1a6 100644
--- a/ven_hma_DIA-phosphoproteomics/03_KSEA_and_KSEAplot/03_KSEA_and_KSEAplot.R
+++ b/ven_hma_DIA-phosphoproteomics/03_KSEA_and_KSEAplot/03_KSEA_and_KSEAplot.R
@@ -64,21 +64,59 @@ PhosInp <- data.frame(Protein = "NULL", Gene = NRRef$feature, Peptide = "NULL",
   dplyr::mutate(Residue.Both = sub("^.*-", "", Residue.Both)) %>%
   dplyr::mutate(Gene = sub("^(.*)-[^-]*$", "\\1", Gene))
 
-ksea_res_full <- KSEAapp::KSEA.Scores(KSDB, PhosInp, NetworKIN = TRUE, NetworKIN.cutoff = Inf) 
+ksea_res_full <- KSEAapp::KSEA.Scores(KSDB, PhosInp, NetworKIN = TRUE, 
+                                      NetworKIN.cutoff = Inf) 
 
 ksea_res <- ksea_res_full %>%
   dplyr::select(Kinase.Gene, m, p.value, FDR, z.score) %>%
   dplyr::rename(kinase = Kinase.Gene, z_score = z.score, p_value = p.value,
                 adj_p_val = FDR, site_size = m)
 
+
 # filter out kinases enriched by >=5 phosphosites and have p-value <0.05
-kinase <- ksea_res %>% filter(site_size >=5 & p_value < 0.05)
+kinase <- ksea_res %>% filter(site_size >= 5 & p_value < 0.05)
 
 kinase <- kinase %>%
   arrange(desc(z_score)) %>%
   mutate(kinase=factor(kinase, levels=kinase))
 
 
+###ADDED BY SARA
+##now get complete scores
+ksea_comp <- KSEAapp::KSEA.Complete(KSDB, PhosInp, NetworKIN = TRUE, NetworKIN.cutoff = Inf, m.cutoff = 5, p.cutoff = 0.05) 
+
+links <- readr::read_csv('Kinase-Substrate Links.csv')  |>
+  dplyr::rename(kinase = 'Kinase.Gene') |>
+  right_join(kinase)
+
+links <- links |> 
+  rowwise() |>
+  mutate(site = paste0(c(`Substrate.Gene`, `Substrate.Mod`), collapse = '-'))
+
+
+links1 <- links |>
+  mutate(comparison = 'Norelapse-Refractory')
+
+##we can also look at one kinase of interest
+
+links |> subset(kinase == 'AURKA') |>
+  ggplot(aes(x=reorder(site, log2FC), y = log2FC, fill = p_value)) + geom_bar(stat='identity') +
+  coord_flip()
+ggsave('nr_ref_aurka.png',height=9)
+
+links$site[abs(links$log2FC) < 1.5] <- ""
+ 
+ggplot(links, aes(x = reorder(kinase,z_score), y = log2FC, col = log2FC)) + 
+  geom_boxplot(outliers=FALSE) + 
+  geom_jitter() + 
+  ggrepel::geom_label_repel(aes(label = site)) +
+  coord_flip()
+
+ggsave('nr_ref_subs.png', height=9)
+
+####end add
+
+
 # plot KSEA results with lollipop plot
 ggplot(kinase, aes(x = z_score, y = kinase)) +
   geom_segment(aes(x = 0, xend = z_score, y = kinase, yend = kinase),
@@ -100,6 +138,9 @@ ggplot(kinase, aes(x = z_score, y = kinase)) +
     legend.text=element_text(size=20)
   )
 
+ggsave('nr_ref_kins.png',height=9)
+
+
 ## subject RRef to following codes for KSEA analysis and plotting.
 ## prepare input for KSEA
 fold_change <- RRef$logFC
@@ -125,6 +166,51 @@ kinase <- kinase %>%
   mutate(kinase=factor(kinase, levels=kinase))
 
 
+###ADDED BY SARA
+##now get complete scores
+ksea_comp <- KSEAapp::KSEA.Complete(KSDB, PhosInp, NetworKIN = TRUE, NetworKIN.cutoff = Inf, m.cutoff = 5, p.cutoff = 0.05) 
+
+links <- readr::read_csv('Kinase-Substrate Links.csv')  |>
+  dplyr::rename(kinase = 'Kinase.Gene') |>
+  right_join(kinase)
+
+links <- links |> 
+  rowwise() |>
+  mutate(site = paste0(c(`Substrate.Gene`, `Substrate.Mod`), collapse = '-')) 
+
+links2 <- links |>
+  mutate(comparison = 'Relapse-Refractory')
+
+links |> 
+  subset(kinase == 'AURKB') |>
+  ggplot(aes(x = reorder(site, log2FC), y = log2FC, fill = p_value)) + 
+  geom_bar(stat = 'identity') +
+  coord_flip()
+
+ggsave('nr_ref_aurkb.png',height=9)
+
+rbind(links1, links2) |>
+  subset(kinase %in% c('AURKB','AURKA')) |>
+  ggplot(aes(x=reorder(site, log2FC), y = log2FC, col = kinase, shape = comparison)) +
+  geom_jitter() +
+  coord_flip()
+
+ggsave('aurk_test.png')
+
+links$site[abs(links$log2FC) < 1.5] <- ""
+
+ggplot(links, aes(x = reorder(kinase,z_score), y = log2FC, col = log2FC)) + 
+  geom_boxplot(outliers=FALSE) + 
+  geom_jitter() + 
+  ggrepel::geom_label_repel(aes(label = site)) +
+  coord_flip()
+
+##
+
+####end add
+ggsave('rel_ref_subs.png',height=9)
+
+
 # plot KSEA results with lollipop plot
 ggplot(kinase, aes(x = z_score, y = kinase)) +
   geom_segment(aes(x = 0, xend = z_score, y = kinase, yend = kinase),
@@ -144,4 +230,6 @@ ggplot(kinase, aes(x = z_score, y = kinase)) +
     axis.text.x=element_text(size=30, face="bold", color="black"),
     legend.title=element_text(size=22, face="bold"),
     legend.text=element_text(size=20)
-  )
\ No newline at end of file
+  )
+
+ggsave('rel_ref_kins.png',height=8)