Ven hma dia phosphoproteomics

ven_hma_DIA-phosphoproteomics/01_Phosphosite_data_processing/01_Phosphosite_data_processing.R

@@ -0,0 +1,56 @@
+library(dplyr)
+library(MSnSet.utils)
+library(stringr)
+library(openxlsx)
+library(synapser)
+
+## 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('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))
+
+## 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", rowNames=FALSE)

ven_hma_DIA-phosphoproteomics/02_LimmaStats/02_LimmaStats.R

@@ -0,0 +1,41 @@
+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 
+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()
+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)

ven_hma_DIA-phosphoproteomics/03_KSEA_and_KSEAplot/03_KSEA_and_KSEAplot.R

@@ -0,0 +1,202 @@
+library(dplyr)
+library(ggplot2)
+library(openxlsx)
+library(synapser)
+
+## synapse login with .Renviron
+synLogin()
+
+## define folder ID and the target files
+folder_id <-"syn68733653"
+
+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")   
+RRef  <- df %>% filter(contrast == "Relapse-Refractory")  
+
+## subject NRRef 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))
+
+
+###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),
+               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)
+  )
+
+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
+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))
+
+
+###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),
+               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)
+  )
+
+ggsave('rel_ref_kins.png',height=8)

ven_hma_DIA-phosphoproteomics/04_GSEA_and_GSEAplot/04_GSEA_and_GSEAplot.R

@@ -0,0 +1,94 @@
+library(dplyr)
+library(ggplot2)
+library(openxlsx)
+library(msigdbr)
+library(fgsea)
+library(scales)
+library(synapser)
+
+## synapse login with .Renviron
+synLogin()
+
+## define folder ID and the target files
+folder_id <-"syn68733653"
+
+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")  
+
+# 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"
+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")
+  )