In nkurzaw/phosphoTPP: Functions and analysis workflow for phosphoTPP datasets
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.path = "inst/figures/README-"
)
reRunFitting <- FALSE
phosphoTPP
Analysis of the phosphoTPP dataset by Potel et al. and comparison to Huang et al.
Potel*, Kurzawa*, Becher*, et al. 2020. "Impact of Phosphorylation on
Thermal Stability of Proteins." BioRxiv, https://doi.org/10.1101/2020.01.14.903849
Introduction
This vignette illustrates the analysis of the phosphoTPP dataset by @potel_2020 and compares it to the one by @huang_2019.
Step-by-step walk through the analysis
Download of required data sets
The supplementary tables 41592_2019_499_MOESM3_ESM.xlsx and 41592_2019_499_MOESM4_ESM.xlsx by @huang_2019 can be downloaded from the journal website: https://www.nature.com/articles/s41592-019-0499-3#Sec26 .
The supplementary tables 56035_2_data_set_478471_pxb8gk.xlsx and 56035_2_data_set_478475_pxb8l1-tableS3.xlsx by @ochoa_2020 can also be downloaded from the respective journal website: https://www.nature.com/articles/s41587-019-0344-3#Sec30 .
The processed data by @potel_2020 and reprocessed data by @huang_2019 can be downloaded from Mendeley: https://data.mendeley.com/datasets/4vwzvxfcnd/4
The table 1-s2.0-S0092867418303854-mmc4.xlsx (used for supplementary figures comparing different TPP-TR experiments) can be downloaded from the respective journal website: https://www.sciencedirect.com/science/article/pii/S0092867418303854#app2
All tables should be saved in the dat folder of the package.
Setup
Install the phosphoTPP package from github
if (!requireNamespace("devtools", quietly = TRUE))
install.packages("devtools")
devtools::install_github("nkurzaw/phosphoTPP")
We then load the required libraries
library(phosphoTPP)
library(dplyr)
library(readr)
library(tidyr)
library(readxl)
library(GGally)
library(ggpmisc)
Set some plotting options
# set plotting sizes
anno_size <- 2
pval_size <- 2
rsq_size <- 2
# define plotting theme
theme_paper <-
theme(axis.title = element_text(size = 8),
axis.text = element_text(size = 6),
legend.text = element_text(size = 6),
title = element_text(size = 9))
And we set the thresholds for filtering for high-quality peptide identifications and quantifications and the filtering criteria for melting curves.
# set thresholds for quality filtering
s2i_thres <- 0.5
p2t_thres <- 4
mascot_score_thres <- 20
fdr_at_score_thres <- 0.01
r2_thres <- 0.8
plateau_thres <- 0.2
Walk through the analysis
Load the config table which tells us which TMT channels were used for which temperature treatment
cfg_tab <- readxl::read_xlsx(file.path("dat", "TPP-TR_config_phospho.xlsx"))
cfg_tab
We then create the temperature annotation data frames for the phosphopeptide and unmodified samples
temperature_anno <- cfg_tab %>%
dplyr::select(matches('[0-9]{3}')) %>%
gather(key, temperature) %>%
mutate(channel = paste("sig", key, sep = "_")) %>%
mutate(fc_channel = paste("rel_fc", key, sep = "_"))
temperature_anno_prot <- cfg_tab %>%
dplyr::select(matches('[0-9]{3}')) %>%
gather(key, temperature) %>%
mutate(channel = paste("signal_sum", key, sep = "_")) %>%
mutate(fc_channel = paste("rel_fc", key, sep = "_"))
We now identify the different raw input files
# read in lists of phospho and nbt protein and peptides tables per replicate
phospho_protein_files <- list.files("dat", pattern = 'phospho.+protein', full.names = TRUE)
phospho_protein_files
nbf_protein_files <- list.files("dat", pattern = 'NBF.+protein', full.names = TRUE)
nbf_protein_files
phospho_peptide_files <- list.files("dat", pattern = 'phospho.+peptide.+txt', full.names = TRUE)
phospho_peptide_files
phosho_mxq_peptide_files <- list.files("dat", pattern = 'pTPP.+txt', full.names = TRUE)
phosho_mxq_peptide_files
nbf_peptide_files <- list.files("dat", pattern = 'NBF.+peptide', full.names = TRUE)
nbf_peptide_files
Then, we read in the individual files into a list of replicate data for the respective data types ('phosphoproteins', unmodified proteins, phosphopeptides, unmodified peptides, etc.)
# read data files into lists of replicates
phospho_protein_tabs <- lapply(seq_len(length(phospho_protein_files)), function(file_i){
read_delim(phospho_protein_files[file_i], delim = "\t") %>%
mutate(replicate = file_i)
})
nbf_protein_tabs <- lapply(seq_len(length(nbf_protein_files)), function(file_i){
read_delim(nbf_protein_files[file_i], delim = "\t") %>%
mutate(replicate = file_i)
})
phospho_peptide_tabs <- lapply(seq_len(length(phospho_peptide_files)), function(file_i){
ibq_tab <- read_delim(phospho_peptide_files[file_i], delim = "\t") %>%
filter_at(vars(starts_with("sig_")), all_vars(!is.na(.)))
mxq_tab <- read_delim(phosho_mxq_peptide_files[file_i], delim = "\t")
combined_tab <- left_join(
ibq_tab %>% dplyr::select(protein_id, sequence, modifications,
matches("sig_"), s2i, p2t, msms_id),
mxq_tab %>% dplyr::select(protein_id_mxq = Proteins, sequence = Sequence, msms_id = `MS/MS scan number`,
mod_sequence = `Modified sequence`, phospho_prob = `Phospho (STY) Probabilities`,
phospho_score_diff = `Phospho (STY) Score Diffs`, phospho_site_STY = `Phospho (STY)`),
by = c("sequence", "msms_id")) %>%
filter(!is.na(protein_id_mxq)) %>%
group_by(sequence, msms_id) %>%
filter(!duplicated(mod_sequence)) %>%
ungroup %>%
filter(getProb(phospho_prob)) %>%
filter(p2t >= p2t_thres, s2i >= s2i_thres) %>%
mutate(replicate = file_i) %>%
mutate(rel_fc_126 = sig_126 / sig_126,
rel_fc_127L = sig_127L / sig_126,
rel_fc_127H = sig_127H / sig_126,
rel_fc_128L = sig_128L / sig_126,
rel_fc_128H = sig_128H / sig_126,
rel_fc_129L = sig_129L / sig_126,
rel_fc_129H = sig_129H / sig_126,
rel_fc_130L = sig_130L / sig_126,
rel_fc_130H = sig_130H / sig_126,
rel_fc_131L = sig_131L / sig_126)
return(combined_tab)
})
phospho_non_phospho_peptide_tabs <- lapply(seq_len(length(phospho_peptide_files)), function(file_i){
read_delim(phospho_peptide_files[file_i], delim = "\t") %>%
filter(is_decoy == 0, score > mascot_score_thres & is_unique == 1, rank == 1,
fdr_at_score < fdr_at_score_thres, !is.na(modifications),
grepl("TMT", modifications), !grepl("Phospho", modifications),
in_quantification_of_protein == 1, !grepl("##", protein_id)) %>%
filter_at(vars(starts_with("sig_")), all_vars(!is.na(.))) %>%
mutate(replicate = file_i) %>%
mutate(rel_fc_126 = sig_126 / sig_126,
rel_fc_127L = sig_127L / sig_126,
rel_fc_127H = sig_127H / sig_126,
rel_fc_128L = sig_128L / sig_126,
rel_fc_128H = sig_128H / sig_126,
rel_fc_129L = sig_129L / sig_126,
rel_fc_129H = sig_129H / sig_126,
rel_fc_130L = sig_130L / sig_126,
rel_fc_130H = sig_130H / sig_126,
rel_fc_131L = sig_131L / sig_126)
})
nbf_peptide_tabs <- lapply(seq_len(length(nbf_peptide_files)), function(file_i){
read_delim(nbf_peptide_files[file_i], delim = "\t") %>%
filter(is_decoy == 0, score > mascot_score_thres & is_unique == 1, rank == 1,
fdr_at_score < fdr_at_score_thres, !is.na(modifications),
grepl("TMT", modifications),
in_quantification_of_protein == 1, !grepl("##", protein_id)) %>%
filter_at(vars(starts_with("sig_")), all_vars(!is.na(.))) %>%
mutate(replicate = file_i) %>%
mutate(rel_fc_126 = sig_126 / sig_126,
rel_fc_127L = sig_127L / sig_126,
rel_fc_127H = sig_127H / sig_126,
rel_fc_128L = sig_128L / sig_126,
rel_fc_128H = sig_128H / sig_126,
rel_fc_129L = sig_129L / sig_126,
rel_fc_129H = sig_129H / sig_126,
rel_fc_130L = sig_130L / sig_126,
rel_fc_130H = sig_130H / sig_126,
rel_fc_131L = sig_131L / sig_126)
})
phospho_protein_tabs <- readRDS("prerun/phospho_protein_tabs.rds")
nbf_protein_tabs <- readRDS("prerun/nbf_protein_tabs.rds")
phospho_peptide_tabs <- readRDS("prerun/phospho_peptide_tabs.rds")
phospho_non_phospho_peptide_tabs <- readRDS("prerun/phospho_non_phospho_peptide_tabs.rds")
nbf_peptide_tabs <- readRDS("prerun/nbf_peptide_tabs.rds")
Next, we compute our first level normalization factors which normalize all peptides based on non-phosphorylated peptides found in matching replicates of phospho-enriched and non-bound fraction samples.
n_rep <- min(length(phospho_peptide_files),
length(nbf_peptide_files))
per_exp_norm_factors <- lapply(seq_len(n_rep), function(rep_i){
norm_in_df <- bind_rows(nbf_peptide_tabs[[rep_i]] %>%
mutate(dataset_id = paste("nbf", replicate, sep = "_")),
phospho_non_phospho_peptide_tabs[[rep_i]] %>%
mutate(dataset_id = paste("phospho", replicate, sep = "_")))
nbf_peptide_norm_factors <- getNormFactors4Data(
pep_tab = norm_in_df,
temperature_anno = temperature_anno,
filter_criteria = list("rel_fc_127H_low" = 0,
"rel_fc_129H_low" = 0.4,
"rel_fc_129H_high" = 0.6,
"rel_fc_130H_low" = 0,
"rel_fc_130H_high" = 0.3,
"rel_fc_131L_low" = 0,
"rel_fc_131L_high" = 0.2),
n_rep = 2)
median_norm_factors <- nbf_peptide_norm_factors
out_df <- tibble(
temperature = nbf_peptide_norm_factors[[1]]$temperature,
median_norm_factor = nbf_peptide_norm_factors[[1]]$median_rel_value /
nbf_peptide_norm_factors[[2]]$median_rel_value)
return(out_df)
})
per_exp_norm_factors
Then, we retrieve our second level of normalization factors which perform the normalization described by @savitski_2014 (originally done on protein level) which normalize peptide fold changes in the different temperature channels towards the melting curve of median values retrieved from high-quality melting peptides.
norm_in_df <-
bind_rows(nbf_peptide_tabs) %>%
mutate(dataset_id = paste("nbf", replicate, sep = "_"))#,
nbf_peptide_norm_factors <- getNormFactors4Data(
pep_tab = norm_in_df,
temperature_anno = temperature_anno,
filter_criteria = list("rel_fc_127H_low" = 1,
"rel_fc_129H_low" = 0.4,
"rel_fc_129H_high" = 0.6,
"rel_fc_130H_low" = 0,
"rel_fc_130H_high" = 0.3,
"rel_fc_131L_low" = 0,
"rel_fc_131L_high" = 0.2),
n_rep = n_rep)
nbf_peptide_norm_factors
Now, we apply both normalization factors and visualize the melting curves for our phosphopeptides and non-bound fraction peptides.
nbf_peptide_filtered <- bind_rows(lapply(seq_len(length(nbf_peptide_tabs)),
function(tab_i){
nbf_peptide_tabs[[tab_i]] %>%
dplyr::select(protein_id, sequence, modifications,
matches("rel_fc_"), replicate, score) %>%
gather(key, value, -protein_id, -sequence, -modifications,
-replicate, -score) %>%
group_by(protein_id, sequence, replicate, modifications, key) %>%
filter(score == max(score)) %>%
left_join(temperature_anno, by = c("key" = "fc_channel")) %>%
left_join(nbf_peptide_norm_factors[[tab_i]],
by = "temperature") %>%
mutate(rel_value = value * norm_factor) %>%
left_join(nbf_protein_tabs[[tab_i]] %>% dplyr::select(protein_id, gene_name),
by = c("protein_id"))
}))
phospho_peptide_tab_filtered <- bind_rows(lapply(seq_len(length(phospho_peptide_tabs)), function(tab_i){
phospho_peptide_tabs[[tab_i]] %>%
dplyr::select(-matches("sig_")) %>%
gather(key, value, matches("rel_fc_")) %>%
group_by(protein_id, sequence, mod_sequence, phospho_site_STY,
replicate, key) %>%
summarise(value = mean(value, na.rm = TRUE),
mean_s2i = mean(s2i, na.rm = TRUE),
mean_p2t = mean(p2t, na.rm = TRUE)) %>%
ungroup %>%
left_join(temperature_anno %>% dplyr::select(-key),
by = c("key" = "fc_channel")) %>%
left_join(per_exp_norm_factors[[tab_i]],
by = "temperature") %>%
left_join(nbf_peptide_norm_factors[[tab_i]],
by = "temperature") %>%
mutate(rel_value = value * median_norm_factor * norm_factor) %>%
left_join(phospho_protein_tabs[[tab_i]] %>%
dplyr::select(protein_id, gene_name),
by = c("protein_id")) %>%
filter(!grepl("##", protein_id))
}))
p_qc1 <- nbf_peptide_filtered %>%
ggplot(aes(temperature, rel_value)) +
geom_boxplot(aes(group = temperature)) +
facet_wrap(~replicate) +
coord_cartesian(ylim = c(0, 2.5))
p_qc2 <- phospho_peptide_tab_filtered %>%
ggplot(aes(temperature, rel_value)) +
geom_boxplot(aes(group = temperature)) +
facet_wrap(~replicate) +
coord_cartesian(ylim = c(0, 2.5))
cowplot::plot_grid(p_qc1 + ggtitle("nbf"),
p_qc2 + ggtitle("phospho"),
ncol = 2)
We then fit sigmoids to the melting profiles of the non-bound fractions peptides
# fit curves to nbf dataset
out_nbf_df <- nbf_peptide_filtered %>%
group_by(gene_name, replicate) %>%
do({
suppressWarnings(fitMeltcurveModelAndEval(df = .))
}) %>%
ungroup
out_nbf_df <- readRDS("prerun/out_nbf_df.rds")
We then filter for hiqh quality fit according to the criteria defined by @savitski_2014 and inspect reproducibility by visualization of melting point scatter plots between the different replicates
# filter for high quality fits
nbf_tm_rep_df <- out_nbf_df %>%
filter(rSq > r2_thres, plateau < plateau_thres) %>%
dplyr::select(gene_name, replicate, meltPoint) %>%
spread(replicate, meltPoint)
pm <- GGally::ggpairs(nbf_tm_rep_df, 2:6,
upper = "blank", diag = NULL,
lower = list(continuous = GGally::wrap(lowerFn)))
GGally::ggmatrix(list(pm[2, 1], NULL, NULL, NULL,
pm[3, 1], pm[3,2], NULL, NULL,
pm[4, 1], pm[4,2], pm[4,3], NULL,
pm[5, 1], pm[5,2], pm[5,3], pm[5,4]), 4, 4,
xAxisLabels = c("Tm rep1", "Tm rep2", "Tm rep3", "Tm rep4"),
yAxisLabels = c("Tm rep2", "Tm rep3","Tm rep4","Tm rep5"))
Next, we fit sigmoids to the melting profiles of the phosphopeptides
# fit phospho peptide melting curves
out_phospho_df <- phospho_peptide_tab_filtered %>%
dplyr::select(protein_id, gene_name, sequence,
mod_sequence, mean_s2i, mean_p2t,
phospho_site_STY, replicate,
rel_value, temperature) %>%
group_by(protein_id, gene_name, sequence,
mod_sequence, mean_s2i, mean_p2t,
phospho_site_STY, replicate) %>%
do({
suppressWarnings(fitMeltcurveModelAndEval(df = .))
}) %>%
ungroup
out_phospho_df <- readRDS("prerun/out_phospho_df.rds")
And again we inspect reproducibility by visualizing scatter plots of meling points estimated in the different replicates
phospho_tm_rep_df <- out_phospho_df %>%
filter(rSq > r2_thres, plateau < plateau_thres) %>%
dplyr::select(gene_name, sequence, mod_sequence, phospho_site_STY, replicate, meltPoint) %>%
spread(replicate, meltPoint)
pm_phospho <- ggpairs(phospho_tm_rep_df, 5:9,
upper = "blank", diag = NULL,
lower = list(continuous = wrap(lowerFn)))
ggmatrix(list(pm_phospho[2, 1], NULL, NULL, NULL,
pm_phospho[3, 1], pm_phospho[3,2], NULL, NULL,
pm_phospho[4, 1], pm_phospho[4,2], pm_phospho[4,3], NULL,
pm_phospho[5, 1], pm_phospho[5,2], pm_phospho[5,3], pm_phospho[5,4]), 4, 4,
xAxisLabels = c("Tm rep1", "Tm rep2", "Tm rep3", "Tm rep4"),
yAxisLabels = c("Tm rep2", "Tm rep3","Tm rep4","Tm rep5"))
Now we join our results of the melting curve fitting for phosphopeptides and unmodified proteins and apply the test for differential melting points established by @savitski_2014
nbf_tm_rep_fil <- out_nbf_df %>%
filter(rSq > r2_thres, plateau < plateau_thres) %>%
group_by(gene_name) %>%
filter(n() > 2) %>%
ungroup %>%
dplyr::select(gene_name, replicate, meltPoint) %>%
spread(replicate, meltPoint)
phospho_tm_rep_fil <- out_phospho_df %>%
filter(rSq > r2_thres, plateau < plateau_thres) %>%
group_by(gene_name, sequence, mod_sequence, phospho_site_STY) %>%
filter(n() > 2) %>%
ungroup() %>%
dplyr::select(gene_name, sequence, mod_sequence, phospho_site_STY, replicate, meltPoint) %>%
spread(replicate, meltPoint)
combo_tm_df <- left_join(phospho_tm_rep_fil, nbf_tm_rep_fil,
by = "gene_name") %>%
rowwise() %>%
mutate(id = paste(gene_name, mod_sequence, phospho_site_STY, sep = "_")) %>%
ungroup()
combo_mean_tm_significant_df <- combo_tm_df %>%
rowwise() %>%
mutate(mean_phospho = mean(c(`1.x`, `2.x`, `3.x`, `4.x`, `5.x`), na.rm = TRUE),
mean_nbf = mean(c(`1.y`, `2.y`, `3.y`, `4.y`, `5.y`), na.rm = TRUE)) %>%
ungroup() %>%
rowwise() %>%
mutate(delta_meltPoint_1 = `1.x` - `1.y`,
delta_meltPoint_2 = `2.x` - `2.y`,
delta_meltPoint_3 = `3.x` - `3.y`,
delta_meltPoint_4 = `4.x` - `4.y`,
delta_meltPoint_5 = `5.x` - `5.y`) %>%
ungroup() %>%
mutate(delta_meltPoint_1_scaled = scale(delta_meltPoint_1),
delta_meltPoint_2_scaled = scale(delta_meltPoint_2),
delta_meltPoint_3_scaled = scale(delta_meltPoint_3),
delta_meltPoint_4_scaled = scale(delta_meltPoint_4),
delta_meltPoint_5_scaled = scale(delta_meltPoint_5)) %>%
mutate(p_value_1 = 1 - pnorm(abs(delta_meltPoint_1_scaled)),
p_value_2 = 1 - pnorm(abs(delta_meltPoint_2_scaled)),
p_value_3 = 1 - pnorm(abs(delta_meltPoint_3_scaled)),
p_value_4 = 1 - pnorm(abs(delta_meltPoint_4_scaled)),
p_value_5 = 1 - pnorm(abs(delta_meltPoint_5_scaled))) %>%
ungroup() %>%
mutate(p_adj_1 = p.adjust(p_value_1, method = "BH"),
p_adj_2 = p.adjust(p_value_2, method = "BH"),
p_adj_3 = p.adjust(p_value_3, method = "BH"),
p_adj_4 = p.adjust(p_value_4, method = "BH"),
p_adj_5 = p.adjust(p_value_5, method = "BH")) %>%
rowwise() %>%
mutate(significant = countAtLeast(x = c(p_adj_1, p_adj_2, p_adj_3, p_adj_4, p_adj_5),
val_with_sign =
c(delta_meltPoint_1_scaled,
delta_meltPoint_2_scaled,
delta_meltPoint_3_scaled,
delta_meltPoint_4_scaled,
delta_meltPoint_5_scaled),
min_n = 2)) %>%
ungroup()
To maximize overlap with @huang_2019 , we here choose gene names overlapping their data in case there are ambiguous gene names
## load Huang et al. phospho dataset
huang_phospho_tab <- readxl::read_xlsx(file.path("dat", "41592_2019_499_MOESM4_ESM.xlsx"), skip = 2)
colnames(huang_phospho_tab)[15] <- "AveragePhosphoTm"
colnames(huang_phospho_tab)[21] <- "deltaTmPval"
colnames(huang_phospho_tab)[22] <- "HighQuality"
huang_phospho_hq <- filter(
huang_phospho_tab,
HighQuality == "Yes")
combo_mean_tm_fcs_choosen_gene_names <-
chooseCoherentGeneNames(in_df = combo_mean_tm_significant_df,
their_data = huang_phospho_tab)
Then, we annotate the phosphorylation sites found by mapping back the modified sequence to protein sequences obtained from Uniprot
uniprot_seq_df <- read_csv("dat/uniprot_seq_df.csv")
combo_mean_tm_pSite_anno <- left_join(
combo_mean_tm_fcs_choosen_gene_names,
uniprot_seq_df,
by = "gene_name") %>%
rowwise() %>%
mutate(Gene_pSite = getPhosphoSiteId(
gene_name = gene_name, string = mod_sequence, seq = seq)) %>%
ungroup %>%
mutate(significant_huang = Gene_pSite %in%
filter(huang_phospho_hq, as.numeric(deltaTmPval) < 0.05)$Gene_pSite)
p_1b <- ggplot(combo_mean_tm_pSite_anno, aes(mean_nbf, mean_phospho)) +
geom_point(shape = 16, alpha = 0.2, size = 0.5) +
geom_line(aes(x,y), data = tibble(x = 40:65, y = 40:65),
color = "darkgray", size = 0.25) +
geom_point(data = filter(combo_mean_tm_pSite_anno, significant),
color = "#b2182b",
shape = 16, size = 0.5) +
geom_point(data = filter(combo_mean_tm_pSite_anno, significant_huang),
color = "#1f78b4",
shape = 16, size = 0.5) +
geom_point(data = filter(combo_mean_tm_pSite_anno, significant_huang, significant),
color = "#ff7f00",
shape = 16, size = 0.5) +
coord_cartesian(xlim = c(40, 65), ylim = c(40, 65)) +
ggpmisc::stat_poly_eq(formula = y ~ x, parse = TRUE, size = rsq_size) +
labs(x = expression('T'[m]^unmodified* ' '* '('*~degree*C*')'),
y = expression('T'[m]^phospho* ' '* '('*~degree*C*')')) +
theme_bw() +
theme_paper
Now, we load the data on machine learning-derived functional scores for phosphosites by @ochoa_2020 and join them to our dataset
# annotate ochoa et al. functional scores
ochoa_et_al_anno <- readxl::read_xlsx(file.path("dat", "56035_2_data_set_478471_pxb8gk.xlsx"))
ochoa_et_al_mapping <- read_delim(file.path("dat", "ochoa_mapping.txt"), delim = "\t")
ochoa_et_al_func_anno <- readxl::read_xlsx(file.path("dat", "56035_2_data_set_478475_pxb8l1-tableS3.xlsx"))
ochoa_et_al_func_anno_psite <- left_join(
ochoa_et_al_func_anno,
ochoa_et_al_anno,
by = c("uniprot", "position")
)
ochoa_et_al_anno_fil <- left_join(
ochoa_et_al_func_anno_psite %>% dplyr::select(uniprot, position, residue, functional_score, is_disopred, isHotspot, isInterface),
ochoa_et_al_mapping %>% dplyr::select(uniprot = From, gene_name = To),
by = "uniprot") %>%
mutate(Gene_pSite = paste0(gene_name, "_", "p", residue, position))
combo_mean_tm_pSite_anno_functional_scores <- left_join(
combo_mean_tm_pSite_anno %>%
rowwise() %>% mutate(mean_delta_meltPoint = mean(c(delta_meltPoint_1, delta_meltPoint_2, delta_meltPoint_3,
delta_meltPoint_4, delta_meltPoint_5), na.rm = TRUE)) %>%
ungroup,
ochoa_et_al_anno_fil %>% na.omit %>% dplyr::select(Gene_pSite, functional_score, is_disopred, isHotspot, isInterface),
by = "Gene_pSite")
We join the functional scores also to the data by @huang_2019
huang_raw_df <- readxl::read_xlsx(file.path("dat", "41592_2019_499_MOESM3_ESM.xlsx"), skip = 2)
colnames(huang_raw_df)[14] <- "AverageBulkTm"
huang_bulk <- huang_raw_df %>%
dplyr::select(uniprot_id = `Uniprot Accession`,
average_tm_nbf = AverageBulkTm) %>%
mutate(average_tm_nbf = as.numeric(average_tm_nbf))
huang_phospho_hq_potel_anno <-
huang_phospho_hq %>%
dplyr::select(uniprot_id = `Uniprot Accession`,
Gene_pSite,
average_tm_phospho = AveragePhosphoTm,
p_value = deltaTmPval) %>%
mutate(significant_potel = Gene_pSite %in%
filter(combo_mean_tm_pSite_anno_functional_scores, significant)$Gene_pSite) %>%
mutate(average_tm_phospho = as.numeric(average_tm_phospho),
p_value = as.numeric(p_value)) %>%
left_join(huang_bulk, by = "uniprot_id")
p_1a <- ggplot(huang_phospho_hq_potel_anno, aes(average_tm_nbf, average_tm_phospho)) +
geom_point(shape = 16, alpha = 0.2, size = 0.5) +
geom_line(aes(x,y), data = tibble(x = 40:65, y = 40:65),
color = "darkgray", size = 0.25) +
geom_point(data = filter(huang_phospho_hq_potel_anno, p_value < 0.05),
color = "#1f78b4", size = 0.5,
shape = 16) +
geom_point(data = filter(huang_phospho_hq_potel_anno, significant_potel),
color = "#b2182b", size = 0.5,
shape = 16) +
geom_point(data = filter(huang_phospho_hq_potel_anno, p_value < 0.05, significant_potel),
color = "#ff7f00", size = 0.5,
shape = 16) +
coord_cartesian(xlim = c(40, 65), ylim = c(40, 65)) +
ggpmisc::stat_poly_eq(formula = y ~ x, parse = TRUE, size = rsq_size) +
labs(x = expression('T'[m]^unmodified*' '* '('*~degree*C*')'),
y = expression('T'[m]^phospho*' '* '('*~degree*C*')')) +
theme_bw() +
theme_paper
And now we visualize the distributions of functional scores for the significantly vs. non-significantly differently thermally shifted phosphosites for the results by @huang_2019 and @potel_2020
# ochoa functional score on Huang et al data
huang_phospho_hq_functional_score <- left_join(
huang_phospho_hq %>%
dplyr::select(Gene_pSite, average_deltaTm = AveragePhosphoTm,
p_value = deltaTmPval) %>%
mutate(average_deltaTm = as.numeric(average_deltaTm),
p_value = as.numeric(p_value)) %>%
mutate(significant = (p_value < 0.05)),
ochoa_et_al_anno_fil %>% na.omit %>% dplyr::select(Gene_pSite, functional_score),
by = "Gene_pSite")
# combine both datasets to save space for plot
func_score_combo_df <- bind_rows(
filter(combo_mean_tm_pSite_anno_functional_scores, !is.na(functional_score)) %>%
dplyr::select(Gene_pSite, significant, functional_score) %>%
mutate(group = "This study"),
filter(huang_phospho_hq_functional_score, !is.na(functional_score)) %>%
dplyr::select(Gene_pSite, significant, functional_score)%>%
mutate(group = "Huang et al.")
)
p_1c <- ggplot(func_score_combo_df , aes(significant, functional_score)) +
geom_violin(aes(fill = significant), alpha = 0.75, color = NA) +
geom_boxplot(outlier.shape = NA, width = 0.25, alpha = 0.5) +
ggsignif::geom_signif(comparisons = list(c("TRUE", "FALSE")),
textsize = anno_size, size = 0.25) +
geom_text(stat = "count", aes(label = ..count..), y = 0.01,
size = anno_size) +
scale_fill_manual(values = c("gray10", "gray70")) +
facet_wrap(~ group) +
labs(x = "found significantly de/stabilized", y = "functional score by Ochoa et al.") +
coord_cartesian(ylim = c(-0.1, 1.1)) +
theme_bw() +
theme(legend.position = "none",
strip.text = element_text(size = 8)) +
theme_paper
We now plot examples melting curves showing significantly altered melting profiles in the phosphopeptide data compared to the unmodified proteins
# examples from this study
p_1d <- plotDiffPhosphoMeltcurve(
g_name = "LYN", p_seq = "_VIEDNEpYTAR_",
nbf_df = nbf_peptide_filtered,
phospho_df = phospho_peptide_tab_filtered,
Gene_pSite = "LYN_pY397") +
theme_paper
p_1e <- plotDiffPhosphoMeltcurve(
g_name = "LMNA", p_seq = "_SGAQASSpTPLpSPTR_",
nbf_df = nbf_peptide_filtered,
phospho_df = phospho_peptide_tab_filtered,
Gene_pSite = "LMNA_pT19pS22") +
theme_paper
p_1f <- plotDiffPhosphoMeltcurve(
g_name = "CARHSP1", p_seq = "_GNVVPpSPLPTRR_",
nbf_df = nbf_peptide_filtered,
phospho_df = phospho_peptide_tab_filtered,
Gene_pSite = "CARHSP1_pS41") +
theme_paper
Finally, we combine all plots created above and plot them together
# assemble full figure
cowplot::plot_grid(
cowplot::plot_grid(p_1a, p_1b, p_1c,
labels = c("a", "b", "c"), ncol = 3, rel_widths = c(0.3, 0.3, 0.3)),
cowplot::plot_grid(p_1d, p_1e, p_1f, labels = c("d", "e", "f"),
ncol = 3, rel_widths = c(0.3, 0.3, 0.3)),
ncol = 1
)
Make output tables
supp1_tab <- phospho_peptide_tab_filtered %>%
dplyr::select(gene_name, protein_id,
sequence, mod_sequence,
phospho_site_STY, key,
replicate, temperature,
rel_value) %>%
unite(combo_var, c(key, replicate, temperature),
sep = "_") %>%
spread(combo_var, rel_value)
write_csv(supp1_tab, path = "phosphoTPP-supplementary_data1-supplementary_table_phospho_peptide_fold_changes.csv")
supp2_tab <- nbf_peptide_filtered %>%
ungroup %>%
dplyr::select(gene_name, protein_id,
sequence, modifications,
key, replicate, temperature,
rel_value) %>%
unite(combo_var, c(key, replicate, temperature),
sep = "_") %>%
group_by(gene_name, protein_id,
sequence, modifications,
combo_var) %>%
summarize(rel_value = mean(rel_value, na.rm = TRUE)) %>%
ungroup %>%
spread(combo_var, rel_value)
write_csv(supp2_tab, path = "phosphoTPP-supplementary_data2-supplementary_table_non-modified_peptide_fold_changes.csv")
phospho_median_fc_df <- chooseCoherentGeneNames(
phospho_peptide_tab_filtered,
combo_mean_tm_pSite_anno_functional_scores) %>%
group_by(gene_name,
mod_sequence,
key, temperature) %>%
summarize(median_rel_value = median(rel_value, na.rm = TRUE)) %>%
ungroup %>%
unite(combo_var, c(key, temperature),
sep = "_") %>%
spread(combo_var, median_rel_value) %>%
dplyr::select(gene_name, mod_sequence,
phospho_median_fc_37 = rel_fc_126_37,
phospho_median_fc_40.4 = rel_fc_127L_40.4,
phospho_median_fc_44 = rel_fc_127H_44,
phospho_median_fc_46.9 = rel_fc_128L_46.9,
phospho_median_fc_49.8 = rel_fc_128H_49.8,
phospho_median_fc_52.9 = rel_fc_129L_52.9,
phospho_median_fc_55.5 = rel_fc_129H_55.5,
phospho_median_fc_58.6 = rel_fc_130L_58.6,
phospho_median_fc_62 = rel_fc_130H_62,
phospho_median_fc_66.3 = rel_fc_131L_66.3)
nbf_median_fc_df <- chooseCoherentGeneNames(
nbf_peptide_filtered,
combo_mean_tm_pSite_anno_functional_scores) %>%
group_by(gene_name,
key, temperature) %>%
summarize(median_rel_value = median(rel_value, na.rm = TRUE)) %>%
ungroup %>%
unite(combo_var, c(key, temperature),
sep = "_") %>%
spread(combo_var, median_rel_value) %>%
dplyr::select(gene_name,
unmodified_median_fc_37 = rel_fc_126_37,
unmodified_median_fc_40.4 = rel_fc_127L_40.4,
unmodified_median_fc_44 = rel_fc_127H_44,
unmodified_median_fc_46.9 = rel_fc_128L_46.9,
unmodified_median_fc_49.8 = rel_fc_128H_49.8,
unmodified_median_fc_52.9 = rel_fc_129L_52.9,
unmodified_median_fc_55.5 = rel_fc_129H_55.5,
unmodified_median_fc_58.6 = rel_fc_130L_58.6,
unmodified_median_fc_62 = rel_fc_130H_62,
unmodified_median_fc_66.3 = rel_fc_131L_66.3)
supp3_tab <- combo_mean_tm_pSite_anno_functional_scores %>%
filter(!is.na(mean_nbf), !is.na(mean_phospho)) %>%
dplyr::select(gene_name, Gene_pSite, mod_sequence,
Tm_mean_phospho = mean_phospho,
Tm_mean_unmodified = mean_nbf,
mean_delta_meltPoint,
significant,
functional_score,
Tm_phospho_rep1 = `1.x`,
Tm_phospho_rep2 = `2.x`,
Tm_phospho_rep3 = `3.x`,
Tm_phospho_rep4 = `4.x`,
Tm_phospho_rep5 = `5.x`,
Tm_unmodified_rep1 = `1.y`,
Tm_unmodified_rep2 = `2.y`,
Tm_unmodified_rep3 = `3.y`,
Tm_unmodified_rep4 = `4.y`,
Tm_unmodified_rep5 = `5.y`,
delta_meltPoint_1,
delta_meltPoint_2,
delta_meltPoint_3,
delta_meltPoint_4,
delta_meltPoint_5,
p_adj_1, p_adj_2,
p_adj_3, p_adj_4,
p_adj_5) %>%
left_join(phospho_median_fc_df,
by = c("gene_name", "mod_sequence")) %>%
left_join(nbf_median_fc_df,
by = c("gene_name"))
write_csv(supp3_tab, path = "phosphoTPP-supplementary_data3-supplementary_table_phospho_non-modified_comparisons.csv")
Supplementary figures
\renewcommand{\figurename}{\textbf{Supplementary figure}}
lowerFn <- function(data, mapping, dot.size) {
p <- ggplot(data = data, mapping = mapping) +
geom_point(alpha = 0.15, size = 0.5) +
geom_line(aes(x,y), data = tibble(x = 40:70, y = 40:70),
color = "darkgray") +
ggpmisc::stat_poly_eq(formula = y ~ x, parse = TRUE, size = 2.5) +
theme_bw() +
theme(axis.text = element_text(size = 6))
p
}
lowerFnReg <- function(data, mapping, dot.size) {
p <- ggplot(data = data, mapping = mapping) +
geom_point(alpha = 0.15) +
geom_line(aes(x,y), data = tibble(x = 40:70, y = 40:70),
color = "darkgray") +
ggpmisc::stat_poly_eq(formula = y ~ x, parse = TRUE) +
theme_bw()
p
}
technical_rep_theme <- theme(
panel.background = element_rect(fill = "lightblue",
colour = "lightblue",
size = 0.5, linetype = "solid"))
huang_data_combo <- left_join(
huang_phospho_hq %>%
dplyr::select(uniprot_id = `Uniprot Accession`,
Gene_pSite,
HighQuality,
p_value = deltaTmPval,
average_tm_phospho = AveragePhosphoTm) %>%
filter(HighQuality == "Yes") %>%
mutate(average_tm_phospho = as.numeric(average_tm_phospho),
p_value = as.numeric(p_value)),
huang_bulk,
by = "uniprot_id")
ggplot(huang_data_combo, aes(average_tm_nbf, average_tm_phospho)) +
geom_point(shape = 16, alpha = 0.2) +
geom_line(aes(x,y), data = tibble(x = 42.5:65, y = 42.5:65),
color = "darkgray") +
coord_fixed(xlim = c(42.5, 65), ylim = c(42.5, 65)) +
ggpmisc::stat_poly_eq(formula = y ~ x, parse = TRUE) +
labs(x = expression('T'[m]^unmodified* ' ' * '('*~degree*C*')'),
y = expression('T'[m]^phospho* ' ' * '('*~degree*C*')')) +
theme_bw()
huang_unmod_mat_numeric <- data.frame(
apply(huang_raw_df[,2:13],
2, as.numeric))
huang_unmod_pm <- ggpairs(huang_unmod_mat_numeric,
upper = "blank", diag = NULL,
lower = list(continuous = wrap(lowerFn)))
ggmatrix(list(huang_unmod_pm[2, 1] + technical_rep_theme, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
huang_unmod_pm[3, 1], huang_unmod_pm[3,2], NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
huang_unmod_pm[4,1], huang_unmod_pm[4,2], huang_unmod_pm[4,3] + technical_rep_theme,
NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
huang_unmod_pm[5,1], huang_unmod_pm[5,2], huang_unmod_pm[5,3], huang_unmod_pm[5,4], NULL, NULL, NULL, NULL, NULL, NULL, NULL,
huang_unmod_pm[6,1], huang_unmod_pm[6,2], huang_unmod_pm[6,3], huang_unmod_pm[6,4],
huang_unmod_pm[6,5]+ technical_rep_theme, NULL, NULL, NULL, NULL, NULL, NULL,
huang_unmod_pm[7,1], huang_unmod_pm[7,2], huang_unmod_pm[7,3], huang_unmod_pm[7,4],
huang_unmod_pm[7,5], huang_unmod_pm[7,6], NULL, NULL, NULL, NULL, NULL,
huang_unmod_pm[8,1], huang_unmod_pm[8,2], huang_unmod_pm[8,3], huang_unmod_pm[8,4], huang_unmod_pm[8,5],
huang_unmod_pm[8,6], huang_unmod_pm[8,7]+ technical_rep_theme, NULL, NULL, NULL, NULL,
huang_unmod_pm[9,1], huang_unmod_pm[9,2], huang_unmod_pm[9,3], huang_unmod_pm[9,4], huang_unmod_pm[9,5],
huang_unmod_pm[9,6], huang_unmod_pm[9,7], huang_unmod_pm[9,8], NULL, NULL, NULL,
huang_unmod_pm[10,1], huang_unmod_pm[10,2], huang_unmod_pm[10,3], huang_unmod_pm[10,4], huang_unmod_pm[10,5],
huang_unmod_pm[10,6], huang_unmod_pm[10,7], huang_unmod_pm[10,8], huang_unmod_pm[10,9] + technical_rep_theme, NULL, NULL,
huang_unmod_pm[11,1], huang_unmod_pm[11,2], huang_unmod_pm[11,3], huang_unmod_pm[11,4], huang_unmod_pm[11,5],
huang_unmod_pm[11,6], huang_unmod_pm[11,7], huang_unmod_pm[11,8], huang_unmod_pm[11,9], huang_unmod_pm[11,10], NULL,
huang_unmod_pm[12,1], huang_unmod_pm[12,2], huang_unmod_pm[12,3], huang_unmod_pm[12,4], huang_unmod_pm[12,5],
huang_unmod_pm[12,6], huang_unmod_pm[12,7], huang_unmod_pm[12,8], huang_unmod_pm[12,9], huang_unmod_pm[12,10],
huang_unmod_pm[12,11] + technical_rep_theme),
11, 11,
xAxisLabels = c("Tm rep1_1", "Tm rep1_2", "Tm rep2_1", "Tm rep2_2", "Tm rep3_1", "Tm rep3_2",
"Tm rep4_1", "Tm rep4_2" ,"Tm rep5_1", "Tm rep5_2", "Tm rep6_1"),
yAxisLabels = c("Tm rep1_2", "Tm rep2_1", "Tm rep2_2", "Tm rep3_1", "Tm rep3_2", "Tm rep4_1",
"Tm rep4_2", "Tm rep5_1", "Tm rep5_2", "Tm rep6_1", "Tm rep6_2")) +
theme(strip.text = element_text(size = 6))
huang_mat_numeric <- data.frame(
apply(huang_phospho_hq[,5:14],
2, as.numeric))
huang_pm <- ggpairs(huang_mat_numeric,
upper = "blank", diag = NULL,
lower = list(continuous = wrap(lowerFn)))
ggmatrix(list(huang_pm[2, 1] + technical_rep_theme, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
huang_pm[3, 1], huang_pm[3,2], NULL, NULL, NULL, NULL, NULL, NULL, NULL,
huang_pm[4,1], huang_pm[4,2], huang_pm[4,3] + technical_rep_theme, NULL, NULL, NULL, NULL, NULL, NULL,
huang_pm[5,1], huang_pm[5,2], huang_pm[5,3], huang_pm[5,4], NULL, NULL, NULL, NULL, NULL,
huang_pm[6,1], huang_pm[6,2], huang_pm[6,3], huang_pm[6,4], huang_pm[6,5]+ technical_rep_theme, NULL, NULL, NULL, NULL,
huang_pm[7,1], huang_pm[7,2], huang_pm[7,3], huang_pm[7,4], huang_pm[7,5], huang_pm[7,6], NULL, NULL, NULL,
huang_pm[8,1], huang_pm[8,2], huang_pm[8,3], huang_pm[8,4], huang_pm[8,5],
huang_pm[8,6], huang_pm[8,7]+ technical_rep_theme, NULL, NULL,
huang_pm[9,1], huang_pm[9,2], huang_pm[9,3], huang_pm[9,4], huang_pm[9,5],
huang_pm[9,6], huang_pm[9,7], huang_pm[9,8], NULL,
huang_pm[10,1], huang_pm[10,2], huang_pm[10,3], huang_pm[10,4], huang_pm[10,5],
huang_pm[10,6], huang_pm[10,7], huang_pm[10,8], huang_pm[10,9] + technical_rep_theme),
9, 9,
xAxisLabels = c("Tm rep1_1", "Tm rep1_2", "Tm rep2_1", "Tm rep2_2", "Tm rep3_1", "Tm rep3_2",
"Tm rep4_1", "Tm rep4_2" ,"Tm rep5_1"),
yAxisLabels = c("Tm rep1_2", "Tm rep2_1", "Tm rep2_2", "Tm rep3_1", "Tm rep3_2", "Tm rep4_1",
"Tm rep4_2", "Tm rep5_1", "Tm rep5_2")) +
theme(strip.text = element_text(size = 6))
pm <- ggpairs(nbf_tm_rep_df, 2:6,
upper = "blank", diag = NULL,
lower = list(continuous = wrap(lowerFnReg)))
ggmatrix(list(pm[2, 1], NULL, NULL, NULL,
pm[3, 1], pm[3,2], NULL, NULL,
pm[4, 1], pm[4,2], pm[4,3], NULL,
pm[5, 1], pm[5,2], pm[5,3], pm[5,4]), 4, 4,
xAxisLabels = c("Tm rep1", "Tm rep2", "Tm rep3", "Tm rep4"),
yAxisLabels = c("Tm rep2", "Tm rep3","Tm rep4","Tm rep5"))
pm_phospho <- ggpairs(phospho_tm_rep_df, 5:9,
upper = "blank", diag = NULL,
lower = list(continuous = wrap(lowerFnReg)))
ggmatrix(list(pm_phospho[2, 1], NULL, NULL, NULL,
pm_phospho[3, 1], pm_phospho[3,2], NULL, NULL,
pm_phospho[4, 1], pm_phospho[4,2], pm_phospho[4,3], NULL,
pm_phospho[5, 1], pm_phospho[5,2], pm_phospho[5,3], pm_phospho[5,4]), 4, 4,
xAxisLabels = c("Tm rep1", "Tm rep2", "Tm rep3", "Tm rep4"),
yAxisLabels = c("Tm rep2", "Tm rep3","Tm rep4","Tm rep5"))
colnames(huang_phospho_hq)[19] <- "deltaTm"
delta_tm_df <- bind_rows(
huang_phospho_hq %>%
mutate(delta_tm = as.numeric(deltaTm)) %>%
dplyr::select(delta_tm) %>%
mutate(dataset = "Huang et al."),
combo_mean_tm_significant_df %>%
rowwise() %>%
mutate(mean_delta_meltPoint = mean(c(delta_meltPoint_1, delta_meltPoint_2, delta_meltPoint_3,
delta_meltPoint_4, delta_meltPoint_5), na.rm = TRUE)) %>%
ungroup %>%
dplyr::select(delta_tm = mean_delta_meltPoint) %>%
mutate(dataset = "This study")
)
ggplot(delta_tm_df, aes(delta_tm, fill = dataset)) +
geom_density(alpha = 0.25) +
xlab(expression(Delta*'T'[m])) +
theme_bw() +
theme(legend.position = "bottom")
mxq_unmod_df_huang <-
read_delim("dat/Huang_nonphospho_evidence.txt", delim = "\t")
mxq_unmod_df_potel <-
read_delim("dat/Potel_nonphospho_evidence.txt", delim = "\t")
mxq_phospho_df_huang <-
read_delim("dat/Huang_phospho_evidence.txt", delim = "\t")
mxq_phospho_df_potel <-
read_delim("dat/Potel_phospho_evidence.txt", delim = "\t")
mxq_combo_df <- bind_rows(
mxq_unmod_df_huang %>% dplyr::select(PIF) %>% mutate(dataset = "Huang et al. non-modified", PIF = as.numeric(PIF)),
mxq_phospho_df_huang %>% dplyr::select(PIF) %>% mutate(dataset = "Huang et al. enriched phosphopeptides", PIF = as.numeric(PIF)),
mxq_unmod_df_potel %>% dplyr::select(PIF) %>% mutate(dataset = "This study non-modified", PIF = as.numeric(PIF)),
mxq_phospho_df_potel %>% dplyr::select(PIF) %>% mutate(dataset = "This study enriched phosphopeptides", PIF = as.numeric(PIF))
)
ggplot(mxq_combo_df, aes(dataset, PIF)) +
geom_violin(alpha = 0.25, fill = "gray", color = "gray") +
geom_boxplot(width = 0.15, outlier.shape = NA) +
theme_bw() +
labs(x = "") +
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 90, hjust =1))
nbf_mean_tm_df <- nbf_tm_rep_fil %>%
ungroup %>%
rowwise() %>%
mutate(mean_Tm = mean(c(`1`, `2`, `3`, `4`, `5`), na.rm = TRUE)) %>%
ungroup %>%
dplyr::select(gene_name, mean_Tm)
huang_mapping <- read_delim("dat/huang_mapping.txt", delim = "\t")
huang_unmod_mean_tm_df <- huang_bulk %>%
left_join(huang_mapping %>% dplyr::select(uniprot_id = From, gene_name = To),
by = "uniprot_id") %>%
dplyr::select(gene_name, huang_mean_tm = average_tm_nbf) %>%
group_by(gene_name) %>%
summarize(huang_mean_tm = mean(huang_mean_tm, na.rm = TRUE)) %>%
ungroup %>%
filter(!is.nan(huang_mean_tm))
tan_hq_meltCurve_param_df <-
readRDS("dat/tan_hq_meltCurve_param_df.rds") %>%
dplyr::select(gene_name, tan_Tm = meltPoint)
becher_cell_hela <- read_xlsx("dat/1-s2.0-S0092867418303854-mmc4.xlsx",
sheet = 2) %>%
dplyr::select(gene_name, becher_cell_hela_mean_Tm = Tm.median) %>%
mutate(becher_cell_hela_mean_Tm = as.numeric(becher_cell_hela_mean_Tm))
combo_tm_df <-
tan_hq_meltCurve_param_df %>%
full_join(becher_cell_hela, by = "gene_name") %>%
full_join(nbf_mean_tm_df, by = "gene_name") %>%
full_join(huang_unmod_mean_tm_df, by = "gene_name")
pm_others <- ggpairs(combo_tm_df, 2:5,
upper = "blank", diag = NULL,
lower = list(continuous = wrap(lowerFnReg)))
ggmatrix(list(pm_others[2, 1], NULL, NULL,
pm_others[3, 1], pm_others[3,2], NULL,
pm_others[4, 1], pm_others[4,2], pm_others[4,3]),
3, 3,
xAxisLabels = c("Tan et al. (K562)", "Becher et al. (HeLa)",
"This study (HeLa)"),
yAxisLabels = c("Becher et al. (HeLa)",
"This study (HeLa)", "Huang et al. (HEK293T)")) +
theme(strip.text = element_text(size = 8))
Delta Tm comparison between the two studies
mean_deltaTm_df <- combo_mean_tm_pSite_anno %>%
rowwise() %>%
mutate(mean_delta_meltPoint = mean(
c(delta_meltPoint_1, delta_meltPoint_2, delta_meltPoint_3,
delta_meltPoint_4, delta_meltPoint_5), na.rm = TRUE)) %>%
ungroup %>%
dplyr::select(gene_name, Gene_pSite, mod_sequence,
phospho_site_STY, mean_delta_meltPoint,
significant)
delta_tm <- left_join(
mean_deltaTm_df %>%
dplyr::select(Gene_pSite, mean_delta_meltPoint,
significant),
huang_phospho_hq_potel_anno %>%
mutate(mean_delta_meltPoint_huang = average_tm_phospho -
average_tm_nbf) %>%
dplyr::select(Gene_pSite, mean_delta_meltPoint_huang, p_value),
by = "Gene_pSite"
) %>%
na.omit()
ggplot(delta_tm, aes(mean_delta_meltPoint_huang, mean_delta_meltPoint)) +
geom_point(alpha = 0.5, shape = 16) +
geom_point(data = filter(delta_tm, significant),
color = "#b2182b",
shape = 16) +
geom_point(data = filter(delta_tm, p_value < 0.05),
color = "#1f78b4",
shape = 16) +
geom_point(data = filter(delta_tm, p_value < 0.05, significant),
color = "#ff7f00",
shape = 16) +
geom_smooth(method = "lm", color = "darkgray", size = 1) +
ggrepel::geom_text_repel(aes(label = Gene_pSite),
data = filter(delta_tm, p_value < 0.05, significant)) +
ggpmisc::stat_poly_eq(formula = y ~ x, parse = TRUE) +
labs(x = expression(Delta*' T'[m]^Potel*' '* '('*~degree*C*')'),
y = expression(Delta* ' T'[m]^Huang*' '* '('*~degree*C*')')) +
theme_bw()
Blue dots correspond to phosphosites found significantly shifted by Huang et al., red dots were found significantly shifted by Potel et al., salmon coloured dots were found by both.
sessionInfo()
References
nkurzaw/phosphoTPP documentation built on Aug. 21, 2021, 5:12 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.path = "inst/figures/README-" ) reRunFitting <- FALSE
phosphoTPP
Analysis of the phosphoTPP dataset by Potel et al. and comparison to Huang et al.
Potel*, Kurzawa*, Becher*, et al. 2020. "Impact of Phosphorylation on Thermal Stability of Proteins." BioRxiv, https://doi.org/10.1101/2020.01.14.903849
Introduction
This vignette illustrates the analysis of the phosphoTPP dataset by @potel_2020 and compares it to the one by @huang_2019.
Step-by-step walk through the analysis
Download of required data sets
The supplementary tables 41592_2019_499_MOESM3_ESM.xlsx and 41592_2019_499_MOESM4_ESM.xlsx by @huang_2019 can be downloaded from the journal website: https://www.nature.com/articles/s41592-019-0499-3#Sec26 .
The supplementary tables 56035_2_data_set_478471_pxb8gk.xlsx and 56035_2_data_set_478475_pxb8l1-tableS3.xlsx by @ochoa_2020 can also be downloaded from the respective journal website: https://www.nature.com/articles/s41587-019-0344-3#Sec30 .
The processed data by @potel_2020 and reprocessed data by @huang_2019 can be downloaded from Mendeley: https://data.mendeley.com/datasets/4vwzvxfcnd/4
The table 1-s2.0-S0092867418303854-mmc4.xlsx (used for supplementary figures comparing different TPP-TR experiments) can be downloaded from the respective journal website: https://www.sciencedirect.com/science/article/pii/S0092867418303854#app2
All tables should be saved in the dat folder of the package.
Setup
Install the phosphoTPP package from github
if (!requireNamespace("devtools", quietly = TRUE)) install.packages("devtools") devtools::install_github("nkurzaw/phosphoTPP")
We then load the required libraries
library(phosphoTPP) library(dplyr) library(readr) library(tidyr) library(readxl) library(GGally) library(ggpmisc)
Set some plotting options
# set plotting sizes anno_size <- 2 pval_size <- 2 rsq_size <- 2 # define plotting theme theme_paper <- theme(axis.title = element_text(size = 8), axis.text = element_text(size = 6), legend.text = element_text(size = 6), title = element_text(size = 9))
And we set the thresholds for filtering for high-quality peptide identifications and quantifications and the filtering criteria for melting curves.
# set thresholds for quality filtering s2i_thres <- 0.5 p2t_thres <- 4 mascot_score_thres <- 20 fdr_at_score_thres <- 0.01 r2_thres <- 0.8 plateau_thres <- 0.2
Walk through the analysis
Load the config table which tells us which TMT channels were used for which temperature treatment
cfg_tab <- readxl::read_xlsx(file.path("dat", "TPP-TR_config_phospho.xlsx")) cfg_tab
We then create the temperature annotation data frames for the phosphopeptide and unmodified samples
temperature_anno <- cfg_tab %>% dplyr::select(matches('[0-9]{3}')) %>% gather(key, temperature) %>% mutate(channel = paste("sig", key, sep = "_")) %>% mutate(fc_channel = paste("rel_fc", key, sep = "_")) temperature_anno_prot <- cfg_tab %>% dplyr::select(matches('[0-9]{3}')) %>% gather(key, temperature) %>% mutate(channel = paste("signal_sum", key, sep = "_")) %>% mutate(fc_channel = paste("rel_fc", key, sep = "_"))
We now identify the different raw input files
# read in lists of phospho and nbt protein and peptides tables per replicate phospho_protein_files <- list.files("dat", pattern = 'phospho.+protein', full.names = TRUE) phospho_protein_files
nbf_protein_files <- list.files("dat", pattern = 'NBF.+protein', full.names = TRUE) nbf_protein_files
phospho_peptide_files <- list.files("dat", pattern = 'phospho.+peptide.+txt', full.names = TRUE) phospho_peptide_files
phosho_mxq_peptide_files <- list.files("dat", pattern = 'pTPP.+txt', full.names = TRUE) phosho_mxq_peptide_files
nbf_peptide_files <- list.files("dat", pattern = 'NBF.+peptide', full.names = TRUE) nbf_peptide_files
Then, we read in the individual files into a list of replicate data for the respective data types ('phosphoproteins', unmodified proteins, phosphopeptides, unmodified peptides, etc.)
# read data files into lists of replicates phospho_protein_tabs <- lapply(seq_len(length(phospho_protein_files)), function(file_i){ read_delim(phospho_protein_files[file_i], delim = "\t") %>% mutate(replicate = file_i) }) nbf_protein_tabs <- lapply(seq_len(length(nbf_protein_files)), function(file_i){ read_delim(nbf_protein_files[file_i], delim = "\t") %>% mutate(replicate = file_i) }) phospho_peptide_tabs <- lapply(seq_len(length(phospho_peptide_files)), function(file_i){ ibq_tab <- read_delim(phospho_peptide_files[file_i], delim = "\t") %>% filter_at(vars(starts_with("sig_")), all_vars(!is.na(.))) mxq_tab <- read_delim(phosho_mxq_peptide_files[file_i], delim = "\t") combined_tab <- left_join( ibq_tab %>% dplyr::select(protein_id, sequence, modifications, matches("sig_"), s2i, p2t, msms_id), mxq_tab %>% dplyr::select(protein_id_mxq = Proteins, sequence = Sequence, msms_id = `MS/MS scan number`, mod_sequence = `Modified sequence`, phospho_prob = `Phospho (STY) Probabilities`, phospho_score_diff = `Phospho (STY) Score Diffs`, phospho_site_STY = `Phospho (STY)`), by = c("sequence", "msms_id")) %>% filter(!is.na(protein_id_mxq)) %>% group_by(sequence, msms_id) %>% filter(!duplicated(mod_sequence)) %>% ungroup %>% filter(getProb(phospho_prob)) %>% filter(p2t >= p2t_thres, s2i >= s2i_thres) %>% mutate(replicate = file_i) %>% mutate(rel_fc_126 = sig_126 / sig_126, rel_fc_127L = sig_127L / sig_126, rel_fc_127H = sig_127H / sig_126, rel_fc_128L = sig_128L / sig_126, rel_fc_128H = sig_128H / sig_126, rel_fc_129L = sig_129L / sig_126, rel_fc_129H = sig_129H / sig_126, rel_fc_130L = sig_130L / sig_126, rel_fc_130H = sig_130H / sig_126, rel_fc_131L = sig_131L / sig_126) return(combined_tab) }) phospho_non_phospho_peptide_tabs <- lapply(seq_len(length(phospho_peptide_files)), function(file_i){ read_delim(phospho_peptide_files[file_i], delim = "\t") %>% filter(is_decoy == 0, score > mascot_score_thres & is_unique == 1, rank == 1, fdr_at_score < fdr_at_score_thres, !is.na(modifications), grepl("TMT", modifications), !grepl("Phospho", modifications), in_quantification_of_protein == 1, !grepl("##", protein_id)) %>% filter_at(vars(starts_with("sig_")), all_vars(!is.na(.))) %>% mutate(replicate = file_i) %>% mutate(rel_fc_126 = sig_126 / sig_126, rel_fc_127L = sig_127L / sig_126, rel_fc_127H = sig_127H / sig_126, rel_fc_128L = sig_128L / sig_126, rel_fc_128H = sig_128H / sig_126, rel_fc_129L = sig_129L / sig_126, rel_fc_129H = sig_129H / sig_126, rel_fc_130L = sig_130L / sig_126, rel_fc_130H = sig_130H / sig_126, rel_fc_131L = sig_131L / sig_126) }) nbf_peptide_tabs <- lapply(seq_len(length(nbf_peptide_files)), function(file_i){ read_delim(nbf_peptide_files[file_i], delim = "\t") %>% filter(is_decoy == 0, score > mascot_score_thres & is_unique == 1, rank == 1, fdr_at_score < fdr_at_score_thres, !is.na(modifications), grepl("TMT", modifications), in_quantification_of_protein == 1, !grepl("##", protein_id)) %>% filter_at(vars(starts_with("sig_")), all_vars(!is.na(.))) %>% mutate(replicate = file_i) %>% mutate(rel_fc_126 = sig_126 / sig_126, rel_fc_127L = sig_127L / sig_126, rel_fc_127H = sig_127H / sig_126, rel_fc_128L = sig_128L / sig_126, rel_fc_128H = sig_128H / sig_126, rel_fc_129L = sig_129L / sig_126, rel_fc_129H = sig_129H / sig_126, rel_fc_130L = sig_130L / sig_126, rel_fc_130H = sig_130H / sig_126, rel_fc_131L = sig_131L / sig_126) })
phospho_protein_tabs <- readRDS("prerun/phospho_protein_tabs.rds") nbf_protein_tabs <- readRDS("prerun/nbf_protein_tabs.rds") phospho_peptide_tabs <- readRDS("prerun/phospho_peptide_tabs.rds") phospho_non_phospho_peptide_tabs <- readRDS("prerun/phospho_non_phospho_peptide_tabs.rds") nbf_peptide_tabs <- readRDS("prerun/nbf_peptide_tabs.rds")
Next, we compute our first level normalization factors which normalize all peptides based on non-phosphorylated peptides found in matching replicates of phospho-enriched and non-bound fraction samples.
n_rep <- min(length(phospho_peptide_files), length(nbf_peptide_files)) per_exp_norm_factors <- lapply(seq_len(n_rep), function(rep_i){ norm_in_df <- bind_rows(nbf_peptide_tabs[[rep_i]] %>% mutate(dataset_id = paste("nbf", replicate, sep = "_")), phospho_non_phospho_peptide_tabs[[rep_i]] %>% mutate(dataset_id = paste("phospho", replicate, sep = "_"))) nbf_peptide_norm_factors <- getNormFactors4Data( pep_tab = norm_in_df, temperature_anno = temperature_anno, filter_criteria = list("rel_fc_127H_low" = 0, "rel_fc_129H_low" = 0.4, "rel_fc_129H_high" = 0.6, "rel_fc_130H_low" = 0, "rel_fc_130H_high" = 0.3, "rel_fc_131L_low" = 0, "rel_fc_131L_high" = 0.2), n_rep = 2) median_norm_factors <- nbf_peptide_norm_factors out_df <- tibble( temperature = nbf_peptide_norm_factors[[1]]$temperature, median_norm_factor = nbf_peptide_norm_factors[[1]]$median_rel_value / nbf_peptide_norm_factors[[2]]$median_rel_value) return(out_df) }) per_exp_norm_factors
Then, we retrieve our second level of normalization factors which perform the normalization described by @savitski_2014 (originally done on protein level) which normalize peptide fold changes in the different temperature channels towards the melting curve of median values retrieved from high-quality melting peptides.
norm_in_df <- bind_rows(nbf_peptide_tabs) %>% mutate(dataset_id = paste("nbf", replicate, sep = "_"))#, nbf_peptide_norm_factors <- getNormFactors4Data( pep_tab = norm_in_df, temperature_anno = temperature_anno, filter_criteria = list("rel_fc_127H_low" = 1, "rel_fc_129H_low" = 0.4, "rel_fc_129H_high" = 0.6, "rel_fc_130H_low" = 0, "rel_fc_130H_high" = 0.3, "rel_fc_131L_low" = 0, "rel_fc_131L_high" = 0.2), n_rep = n_rep) nbf_peptide_norm_factors
Now, we apply both normalization factors and visualize the melting curves for our phosphopeptides and non-bound fraction peptides.
nbf_peptide_filtered <- bind_rows(lapply(seq_len(length(nbf_peptide_tabs)), function(tab_i){ nbf_peptide_tabs[[tab_i]] %>% dplyr::select(protein_id, sequence, modifications, matches("rel_fc_"), replicate, score) %>% gather(key, value, -protein_id, -sequence, -modifications, -replicate, -score) %>% group_by(protein_id, sequence, replicate, modifications, key) %>% filter(score == max(score)) %>% left_join(temperature_anno, by = c("key" = "fc_channel")) %>% left_join(nbf_peptide_norm_factors[[tab_i]], by = "temperature") %>% mutate(rel_value = value * norm_factor) %>% left_join(nbf_protein_tabs[[tab_i]] %>% dplyr::select(protein_id, gene_name), by = c("protein_id")) })) phospho_peptide_tab_filtered <- bind_rows(lapply(seq_len(length(phospho_peptide_tabs)), function(tab_i){ phospho_peptide_tabs[[tab_i]] %>% dplyr::select(-matches("sig_")) %>% gather(key, value, matches("rel_fc_")) %>% group_by(protein_id, sequence, mod_sequence, phospho_site_STY, replicate, key) %>% summarise(value = mean(value, na.rm = TRUE), mean_s2i = mean(s2i, na.rm = TRUE), mean_p2t = mean(p2t, na.rm = TRUE)) %>% ungroup %>% left_join(temperature_anno %>% dplyr::select(-key), by = c("key" = "fc_channel")) %>% left_join(per_exp_norm_factors[[tab_i]], by = "temperature") %>% left_join(nbf_peptide_norm_factors[[tab_i]], by = "temperature") %>% mutate(rel_value = value * median_norm_factor * norm_factor) %>% left_join(phospho_protein_tabs[[tab_i]] %>% dplyr::select(protein_id, gene_name), by = c("protein_id")) %>% filter(!grepl("##", protein_id)) })) p_qc1 <- nbf_peptide_filtered %>% ggplot(aes(temperature, rel_value)) + geom_boxplot(aes(group = temperature)) + facet_wrap(~replicate) + coord_cartesian(ylim = c(0, 2.5)) p_qc2 <- phospho_peptide_tab_filtered %>% ggplot(aes(temperature, rel_value)) + geom_boxplot(aes(group = temperature)) + facet_wrap(~replicate) + coord_cartesian(ylim = c(0, 2.5)) cowplot::plot_grid(p_qc1 + ggtitle("nbf"), p_qc2 + ggtitle("phospho"), ncol = 2)
We then fit sigmoids to the melting profiles of the non-bound fractions peptides
# fit curves to nbf dataset out_nbf_df <- nbf_peptide_filtered %>% group_by(gene_name, replicate) %>% do({ suppressWarnings(fitMeltcurveModelAndEval(df = .)) }) %>% ungroup
out_nbf_df <- readRDS("prerun/out_nbf_df.rds")
We then filter for hiqh quality fit according to the criteria defined by @savitski_2014 and inspect reproducibility by visualization of melting point scatter plots between the different replicates
# filter for high quality fits nbf_tm_rep_df <- out_nbf_df %>% filter(rSq > r2_thres, plateau < plateau_thres) %>% dplyr::select(gene_name, replicate, meltPoint) %>% spread(replicate, meltPoint) pm <- GGally::ggpairs(nbf_tm_rep_df, 2:6, upper = "blank", diag = NULL, lower = list(continuous = GGally::wrap(lowerFn))) GGally::ggmatrix(list(pm[2, 1], NULL, NULL, NULL, pm[3, 1], pm[3,2], NULL, NULL, pm[4, 1], pm[4,2], pm[4,3], NULL, pm[5, 1], pm[5,2], pm[5,3], pm[5,4]), 4, 4, xAxisLabels = c("Tm rep1", "Tm rep2", "Tm rep3", "Tm rep4"), yAxisLabels = c("Tm rep2", "Tm rep3","Tm rep4","Tm rep5"))
Next, we fit sigmoids to the melting profiles of the phosphopeptides
# fit phospho peptide melting curves out_phospho_df <- phospho_peptide_tab_filtered %>% dplyr::select(protein_id, gene_name, sequence, mod_sequence, mean_s2i, mean_p2t, phospho_site_STY, replicate, rel_value, temperature) %>% group_by(protein_id, gene_name, sequence, mod_sequence, mean_s2i, mean_p2t, phospho_site_STY, replicate) %>% do({ suppressWarnings(fitMeltcurveModelAndEval(df = .)) }) %>% ungroup
out_phospho_df <- readRDS("prerun/out_phospho_df.rds")
And again we inspect reproducibility by visualizing scatter plots of meling points estimated in the different replicates
phospho_tm_rep_df <- out_phospho_df %>% filter(rSq > r2_thres, plateau < plateau_thres) %>% dplyr::select(gene_name, sequence, mod_sequence, phospho_site_STY, replicate, meltPoint) %>% spread(replicate, meltPoint) pm_phospho <- ggpairs(phospho_tm_rep_df, 5:9, upper = "blank", diag = NULL, lower = list(continuous = wrap(lowerFn))) ggmatrix(list(pm_phospho[2, 1], NULL, NULL, NULL, pm_phospho[3, 1], pm_phospho[3,2], NULL, NULL, pm_phospho[4, 1], pm_phospho[4,2], pm_phospho[4,3], NULL, pm_phospho[5, 1], pm_phospho[5,2], pm_phospho[5,3], pm_phospho[5,4]), 4, 4, xAxisLabels = c("Tm rep1", "Tm rep2", "Tm rep3", "Tm rep4"), yAxisLabels = c("Tm rep2", "Tm rep3","Tm rep4","Tm rep5"))
Now we join our results of the melting curve fitting for phosphopeptides and unmodified proteins and apply the test for differential melting points established by @savitski_2014
nbf_tm_rep_fil <- out_nbf_df %>% filter(rSq > r2_thres, plateau < plateau_thres) %>% group_by(gene_name) %>% filter(n() > 2) %>% ungroup %>% dplyr::select(gene_name, replicate, meltPoint) %>% spread(replicate, meltPoint) phospho_tm_rep_fil <- out_phospho_df %>% filter(rSq > r2_thres, plateau < plateau_thres) %>% group_by(gene_name, sequence, mod_sequence, phospho_site_STY) %>% filter(n() > 2) %>% ungroup() %>% dplyr::select(gene_name, sequence, mod_sequence, phospho_site_STY, replicate, meltPoint) %>% spread(replicate, meltPoint) combo_tm_df <- left_join(phospho_tm_rep_fil, nbf_tm_rep_fil, by = "gene_name") %>% rowwise() %>% mutate(id = paste(gene_name, mod_sequence, phospho_site_STY, sep = "_")) %>% ungroup() combo_mean_tm_significant_df <- combo_tm_df %>% rowwise() %>% mutate(mean_phospho = mean(c(`1.x`, `2.x`, `3.x`, `4.x`, `5.x`), na.rm = TRUE), mean_nbf = mean(c(`1.y`, `2.y`, `3.y`, `4.y`, `5.y`), na.rm = TRUE)) %>% ungroup() %>% rowwise() %>% mutate(delta_meltPoint_1 = `1.x` - `1.y`, delta_meltPoint_2 = `2.x` - `2.y`, delta_meltPoint_3 = `3.x` - `3.y`, delta_meltPoint_4 = `4.x` - `4.y`, delta_meltPoint_5 = `5.x` - `5.y`) %>% ungroup() %>% mutate(delta_meltPoint_1_scaled = scale(delta_meltPoint_1), delta_meltPoint_2_scaled = scale(delta_meltPoint_2), delta_meltPoint_3_scaled = scale(delta_meltPoint_3), delta_meltPoint_4_scaled = scale(delta_meltPoint_4), delta_meltPoint_5_scaled = scale(delta_meltPoint_5)) %>% mutate(p_value_1 = 1 - pnorm(abs(delta_meltPoint_1_scaled)), p_value_2 = 1 - pnorm(abs(delta_meltPoint_2_scaled)), p_value_3 = 1 - pnorm(abs(delta_meltPoint_3_scaled)), p_value_4 = 1 - pnorm(abs(delta_meltPoint_4_scaled)), p_value_5 = 1 - pnorm(abs(delta_meltPoint_5_scaled))) %>% ungroup() %>% mutate(p_adj_1 = p.adjust(p_value_1, method = "BH"), p_adj_2 = p.adjust(p_value_2, method = "BH"), p_adj_3 = p.adjust(p_value_3, method = "BH"), p_adj_4 = p.adjust(p_value_4, method = "BH"), p_adj_5 = p.adjust(p_value_5, method = "BH")) %>% rowwise() %>% mutate(significant = countAtLeast(x = c(p_adj_1, p_adj_2, p_adj_3, p_adj_4, p_adj_5), val_with_sign = c(delta_meltPoint_1_scaled, delta_meltPoint_2_scaled, delta_meltPoint_3_scaled, delta_meltPoint_4_scaled, delta_meltPoint_5_scaled), min_n = 2)) %>% ungroup()
To maximize overlap with @huang_2019 , we here choose gene names overlapping their data in case there are ambiguous gene names
## load Huang et al. phospho dataset huang_phospho_tab <- readxl::read_xlsx(file.path("dat", "41592_2019_499_MOESM4_ESM.xlsx"), skip = 2) colnames(huang_phospho_tab)[15] <- "AveragePhosphoTm" colnames(huang_phospho_tab)[21] <- "deltaTmPval" colnames(huang_phospho_tab)[22] <- "HighQuality" huang_phospho_hq <- filter( huang_phospho_tab, HighQuality == "Yes") combo_mean_tm_fcs_choosen_gene_names <- chooseCoherentGeneNames(in_df = combo_mean_tm_significant_df, their_data = huang_phospho_tab)
Then, we annotate the phosphorylation sites found by mapping back the modified sequence to protein sequences obtained from Uniprot
uniprot_seq_df <- read_csv("dat/uniprot_seq_df.csv") combo_mean_tm_pSite_anno <- left_join( combo_mean_tm_fcs_choosen_gene_names, uniprot_seq_df, by = "gene_name") %>% rowwise() %>% mutate(Gene_pSite = getPhosphoSiteId( gene_name = gene_name, string = mod_sequence, seq = seq)) %>% ungroup %>% mutate(significant_huang = Gene_pSite %in% filter(huang_phospho_hq, as.numeric(deltaTmPval) < 0.05)$Gene_pSite) p_1b <- ggplot(combo_mean_tm_pSite_anno, aes(mean_nbf, mean_phospho)) + geom_point(shape = 16, alpha = 0.2, size = 0.5) + geom_line(aes(x,y), data = tibble(x = 40:65, y = 40:65), color = "darkgray", size = 0.25) + geom_point(data = filter(combo_mean_tm_pSite_anno, significant), color = "#b2182b", shape = 16, size = 0.5) + geom_point(data = filter(combo_mean_tm_pSite_anno, significant_huang), color = "#1f78b4", shape = 16, size = 0.5) + geom_point(data = filter(combo_mean_tm_pSite_anno, significant_huang, significant), color = "#ff7f00", shape = 16, size = 0.5) + coord_cartesian(xlim = c(40, 65), ylim = c(40, 65)) + ggpmisc::stat_poly_eq(formula = y ~ x, parse = TRUE, size = rsq_size) + labs(x = expression('T'[m]^unmodified* ' '* '('*~degree*C*')'), y = expression('T'[m]^phospho* ' '* '('*~degree*C*')')) + theme_bw() + theme_paper
Now, we load the data on machine learning-derived functional scores for phosphosites by @ochoa_2020 and join them to our dataset
# annotate ochoa et al. functional scores ochoa_et_al_anno <- readxl::read_xlsx(file.path("dat", "56035_2_data_set_478471_pxb8gk.xlsx")) ochoa_et_al_mapping <- read_delim(file.path("dat", "ochoa_mapping.txt"), delim = "\t") ochoa_et_al_func_anno <- readxl::read_xlsx(file.path("dat", "56035_2_data_set_478475_pxb8l1-tableS3.xlsx")) ochoa_et_al_func_anno_psite <- left_join( ochoa_et_al_func_anno, ochoa_et_al_anno, by = c("uniprot", "position") ) ochoa_et_al_anno_fil <- left_join( ochoa_et_al_func_anno_psite %>% dplyr::select(uniprot, position, residue, functional_score, is_disopred, isHotspot, isInterface), ochoa_et_al_mapping %>% dplyr::select(uniprot = From, gene_name = To), by = "uniprot") %>% mutate(Gene_pSite = paste0(gene_name, "_", "p", residue, position)) combo_mean_tm_pSite_anno_functional_scores <- left_join( combo_mean_tm_pSite_anno %>% rowwise() %>% mutate(mean_delta_meltPoint = mean(c(delta_meltPoint_1, delta_meltPoint_2, delta_meltPoint_3, delta_meltPoint_4, delta_meltPoint_5), na.rm = TRUE)) %>% ungroup, ochoa_et_al_anno_fil %>% na.omit %>% dplyr::select(Gene_pSite, functional_score, is_disopred, isHotspot, isInterface), by = "Gene_pSite")
We join the functional scores also to the data by @huang_2019
huang_raw_df <- readxl::read_xlsx(file.path("dat", "41592_2019_499_MOESM3_ESM.xlsx"), skip = 2) colnames(huang_raw_df)[14] <- "AverageBulkTm" huang_bulk <- huang_raw_df %>% dplyr::select(uniprot_id = `Uniprot Accession`, average_tm_nbf = AverageBulkTm) %>% mutate(average_tm_nbf = as.numeric(average_tm_nbf)) huang_phospho_hq_potel_anno <- huang_phospho_hq %>% dplyr::select(uniprot_id = `Uniprot Accession`, Gene_pSite, average_tm_phospho = AveragePhosphoTm, p_value = deltaTmPval) %>% mutate(significant_potel = Gene_pSite %in% filter(combo_mean_tm_pSite_anno_functional_scores, significant)$Gene_pSite) %>% mutate(average_tm_phospho = as.numeric(average_tm_phospho), p_value = as.numeric(p_value)) %>% left_join(huang_bulk, by = "uniprot_id") p_1a <- ggplot(huang_phospho_hq_potel_anno, aes(average_tm_nbf, average_tm_phospho)) + geom_point(shape = 16, alpha = 0.2, size = 0.5) + geom_line(aes(x,y), data = tibble(x = 40:65, y = 40:65), color = "darkgray", size = 0.25) + geom_point(data = filter(huang_phospho_hq_potel_anno, p_value < 0.05), color = "#1f78b4", size = 0.5, shape = 16) + geom_point(data = filter(huang_phospho_hq_potel_anno, significant_potel), color = "#b2182b", size = 0.5, shape = 16) + geom_point(data = filter(huang_phospho_hq_potel_anno, p_value < 0.05, significant_potel), color = "#ff7f00", size = 0.5, shape = 16) + coord_cartesian(xlim = c(40, 65), ylim = c(40, 65)) + ggpmisc::stat_poly_eq(formula = y ~ x, parse = TRUE, size = rsq_size) + labs(x = expression('T'[m]^unmodified*' '* '('*~degree*C*')'), y = expression('T'[m]^phospho*' '* '('*~degree*C*')')) + theme_bw() + theme_paper
And now we visualize the distributions of functional scores for the significantly vs. non-significantly differently thermally shifted phosphosites for the results by @huang_2019 and @potel_2020
# ochoa functional score on Huang et al data huang_phospho_hq_functional_score <- left_join( huang_phospho_hq %>% dplyr::select(Gene_pSite, average_deltaTm = AveragePhosphoTm, p_value = deltaTmPval) %>% mutate(average_deltaTm = as.numeric(average_deltaTm), p_value = as.numeric(p_value)) %>% mutate(significant = (p_value < 0.05)), ochoa_et_al_anno_fil %>% na.omit %>% dplyr::select(Gene_pSite, functional_score), by = "Gene_pSite") # combine both datasets to save space for plot func_score_combo_df <- bind_rows( filter(combo_mean_tm_pSite_anno_functional_scores, !is.na(functional_score)) %>% dplyr::select(Gene_pSite, significant, functional_score) %>% mutate(group = "This study"), filter(huang_phospho_hq_functional_score, !is.na(functional_score)) %>% dplyr::select(Gene_pSite, significant, functional_score)%>% mutate(group = "Huang et al.") ) p_1c <- ggplot(func_score_combo_df , aes(significant, functional_score)) + geom_violin(aes(fill = significant), alpha = 0.75, color = NA) + geom_boxplot(outlier.shape = NA, width = 0.25, alpha = 0.5) + ggsignif::geom_signif(comparisons = list(c("TRUE", "FALSE")), textsize = anno_size, size = 0.25) + geom_text(stat = "count", aes(label = ..count..), y = 0.01, size = anno_size) + scale_fill_manual(values = c("gray10", "gray70")) + facet_wrap(~ group) + labs(x = "found significantly de/stabilized", y = "functional score by Ochoa et al.") + coord_cartesian(ylim = c(-0.1, 1.1)) + theme_bw() + theme(legend.position = "none", strip.text = element_text(size = 8)) + theme_paper
We now plot examples melting curves showing significantly altered melting profiles in the phosphopeptide data compared to the unmodified proteins
# examples from this study p_1d <- plotDiffPhosphoMeltcurve( g_name = "LYN", p_seq = "_VIEDNEpYTAR_", nbf_df = nbf_peptide_filtered, phospho_df = phospho_peptide_tab_filtered, Gene_pSite = "LYN_pY397") + theme_paper p_1e <- plotDiffPhosphoMeltcurve( g_name = "LMNA", p_seq = "_SGAQASSpTPLpSPTR_", nbf_df = nbf_peptide_filtered, phospho_df = phospho_peptide_tab_filtered, Gene_pSite = "LMNA_pT19pS22") + theme_paper p_1f <- plotDiffPhosphoMeltcurve( g_name = "CARHSP1", p_seq = "_GNVVPpSPLPTRR_", nbf_df = nbf_peptide_filtered, phospho_df = phospho_peptide_tab_filtered, Gene_pSite = "CARHSP1_pS41") + theme_paper
Finally, we combine all plots created above and plot them together
# assemble full figure cowplot::plot_grid( cowplot::plot_grid(p_1a, p_1b, p_1c, labels = c("a", "b", "c"), ncol = 3, rel_widths = c(0.3, 0.3, 0.3)), cowplot::plot_grid(p_1d, p_1e, p_1f, labels = c("d", "e", "f"), ncol = 3, rel_widths = c(0.3, 0.3, 0.3)), ncol = 1 )
Make output tables
supp1_tab <- phospho_peptide_tab_filtered %>% dplyr::select(gene_name, protein_id, sequence, mod_sequence, phospho_site_STY, key, replicate, temperature, rel_value) %>% unite(combo_var, c(key, replicate, temperature), sep = "_") %>% spread(combo_var, rel_value) write_csv(supp1_tab, path = "phosphoTPP-supplementary_data1-supplementary_table_phospho_peptide_fold_changes.csv") supp2_tab <- nbf_peptide_filtered %>% ungroup %>% dplyr::select(gene_name, protein_id, sequence, modifications, key, replicate, temperature, rel_value) %>% unite(combo_var, c(key, replicate, temperature), sep = "_") %>% group_by(gene_name, protein_id, sequence, modifications, combo_var) %>% summarize(rel_value = mean(rel_value, na.rm = TRUE)) %>% ungroup %>% spread(combo_var, rel_value) write_csv(supp2_tab, path = "phosphoTPP-supplementary_data2-supplementary_table_non-modified_peptide_fold_changes.csv") phospho_median_fc_df <- chooseCoherentGeneNames( phospho_peptide_tab_filtered, combo_mean_tm_pSite_anno_functional_scores) %>% group_by(gene_name, mod_sequence, key, temperature) %>% summarize(median_rel_value = median(rel_value, na.rm = TRUE)) %>% ungroup %>% unite(combo_var, c(key, temperature), sep = "_") %>% spread(combo_var, median_rel_value) %>% dplyr::select(gene_name, mod_sequence, phospho_median_fc_37 = rel_fc_126_37, phospho_median_fc_40.4 = rel_fc_127L_40.4, phospho_median_fc_44 = rel_fc_127H_44, phospho_median_fc_46.9 = rel_fc_128L_46.9, phospho_median_fc_49.8 = rel_fc_128H_49.8, phospho_median_fc_52.9 = rel_fc_129L_52.9, phospho_median_fc_55.5 = rel_fc_129H_55.5, phospho_median_fc_58.6 = rel_fc_130L_58.6, phospho_median_fc_62 = rel_fc_130H_62, phospho_median_fc_66.3 = rel_fc_131L_66.3) nbf_median_fc_df <- chooseCoherentGeneNames( nbf_peptide_filtered, combo_mean_tm_pSite_anno_functional_scores) %>% group_by(gene_name, key, temperature) %>% summarize(median_rel_value = median(rel_value, na.rm = TRUE)) %>% ungroup %>% unite(combo_var, c(key, temperature), sep = "_") %>% spread(combo_var, median_rel_value) %>% dplyr::select(gene_name, unmodified_median_fc_37 = rel_fc_126_37, unmodified_median_fc_40.4 = rel_fc_127L_40.4, unmodified_median_fc_44 = rel_fc_127H_44, unmodified_median_fc_46.9 = rel_fc_128L_46.9, unmodified_median_fc_49.8 = rel_fc_128H_49.8, unmodified_median_fc_52.9 = rel_fc_129L_52.9, unmodified_median_fc_55.5 = rel_fc_129H_55.5, unmodified_median_fc_58.6 = rel_fc_130L_58.6, unmodified_median_fc_62 = rel_fc_130H_62, unmodified_median_fc_66.3 = rel_fc_131L_66.3) supp3_tab <- combo_mean_tm_pSite_anno_functional_scores %>% filter(!is.na(mean_nbf), !is.na(mean_phospho)) %>% dplyr::select(gene_name, Gene_pSite, mod_sequence, Tm_mean_phospho = mean_phospho, Tm_mean_unmodified = mean_nbf, mean_delta_meltPoint, significant, functional_score, Tm_phospho_rep1 = `1.x`, Tm_phospho_rep2 = `2.x`, Tm_phospho_rep3 = `3.x`, Tm_phospho_rep4 = `4.x`, Tm_phospho_rep5 = `5.x`, Tm_unmodified_rep1 = `1.y`, Tm_unmodified_rep2 = `2.y`, Tm_unmodified_rep3 = `3.y`, Tm_unmodified_rep4 = `4.y`, Tm_unmodified_rep5 = `5.y`, delta_meltPoint_1, delta_meltPoint_2, delta_meltPoint_3, delta_meltPoint_4, delta_meltPoint_5, p_adj_1, p_adj_2, p_adj_3, p_adj_4, p_adj_5) %>% left_join(phospho_median_fc_df, by = c("gene_name", "mod_sequence")) %>% left_join(nbf_median_fc_df, by = c("gene_name")) write_csv(supp3_tab, path = "phosphoTPP-supplementary_data3-supplementary_table_phospho_non-modified_comparisons.csv")
Supplementary figures
\renewcommand{\figurename}{\textbf{Supplementary figure}}
lowerFn <- function(data, mapping, dot.size) { p <- ggplot(data = data, mapping = mapping) + geom_point(alpha = 0.15, size = 0.5) + geom_line(aes(x,y), data = tibble(x = 40:70, y = 40:70), color = "darkgray") + ggpmisc::stat_poly_eq(formula = y ~ x, parse = TRUE, size = 2.5) + theme_bw() + theme(axis.text = element_text(size = 6)) p } lowerFnReg <- function(data, mapping, dot.size) { p <- ggplot(data = data, mapping = mapping) + geom_point(alpha = 0.15) + geom_line(aes(x,y), data = tibble(x = 40:70, y = 40:70), color = "darkgray") + ggpmisc::stat_poly_eq(formula = y ~ x, parse = TRUE) + theme_bw() p } technical_rep_theme <- theme( panel.background = element_rect(fill = "lightblue", colour = "lightblue", size = 0.5, linetype = "solid"))
huang_data_combo <- left_join( huang_phospho_hq %>% dplyr::select(uniprot_id = `Uniprot Accession`, Gene_pSite, HighQuality, p_value = deltaTmPval, average_tm_phospho = AveragePhosphoTm) %>% filter(HighQuality == "Yes") %>% mutate(average_tm_phospho = as.numeric(average_tm_phospho), p_value = as.numeric(p_value)), huang_bulk, by = "uniprot_id") ggplot(huang_data_combo, aes(average_tm_nbf, average_tm_phospho)) + geom_point(shape = 16, alpha = 0.2) + geom_line(aes(x,y), data = tibble(x = 42.5:65, y = 42.5:65), color = "darkgray") + coord_fixed(xlim = c(42.5, 65), ylim = c(42.5, 65)) + ggpmisc::stat_poly_eq(formula = y ~ x, parse = TRUE) + labs(x = expression('T'[m]^unmodified* ' ' * '('*~degree*C*')'), y = expression('T'[m]^phospho* ' ' * '('*~degree*C*')')) + theme_bw()
huang_unmod_mat_numeric <- data.frame( apply(huang_raw_df[,2:13], 2, as.numeric)) huang_unmod_pm <- ggpairs(huang_unmod_mat_numeric, upper = "blank", diag = NULL, lower = list(continuous = wrap(lowerFn))) ggmatrix(list(huang_unmod_pm[2, 1] + technical_rep_theme, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, huang_unmod_pm[3, 1], huang_unmod_pm[3,2], NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, huang_unmod_pm[4,1], huang_unmod_pm[4,2], huang_unmod_pm[4,3] + technical_rep_theme, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, huang_unmod_pm[5,1], huang_unmod_pm[5,2], huang_unmod_pm[5,3], huang_unmod_pm[5,4], NULL, NULL, NULL, NULL, NULL, NULL, NULL, huang_unmod_pm[6,1], huang_unmod_pm[6,2], huang_unmod_pm[6,3], huang_unmod_pm[6,4], huang_unmod_pm[6,5]+ technical_rep_theme, NULL, NULL, NULL, NULL, NULL, NULL, huang_unmod_pm[7,1], huang_unmod_pm[7,2], huang_unmod_pm[7,3], huang_unmod_pm[7,4], huang_unmod_pm[7,5], huang_unmod_pm[7,6], NULL, NULL, NULL, NULL, NULL, huang_unmod_pm[8,1], huang_unmod_pm[8,2], huang_unmod_pm[8,3], huang_unmod_pm[8,4], huang_unmod_pm[8,5], huang_unmod_pm[8,6], huang_unmod_pm[8,7]+ technical_rep_theme, NULL, NULL, NULL, NULL, huang_unmod_pm[9,1], huang_unmod_pm[9,2], huang_unmod_pm[9,3], huang_unmod_pm[9,4], huang_unmod_pm[9,5], huang_unmod_pm[9,6], huang_unmod_pm[9,7], huang_unmod_pm[9,8], NULL, NULL, NULL, huang_unmod_pm[10,1], huang_unmod_pm[10,2], huang_unmod_pm[10,3], huang_unmod_pm[10,4], huang_unmod_pm[10,5], huang_unmod_pm[10,6], huang_unmod_pm[10,7], huang_unmod_pm[10,8], huang_unmod_pm[10,9] + technical_rep_theme, NULL, NULL, huang_unmod_pm[11,1], huang_unmod_pm[11,2], huang_unmod_pm[11,3], huang_unmod_pm[11,4], huang_unmod_pm[11,5], huang_unmod_pm[11,6], huang_unmod_pm[11,7], huang_unmod_pm[11,8], huang_unmod_pm[11,9], huang_unmod_pm[11,10], NULL, huang_unmod_pm[12,1], huang_unmod_pm[12,2], huang_unmod_pm[12,3], huang_unmod_pm[12,4], huang_unmod_pm[12,5], huang_unmod_pm[12,6], huang_unmod_pm[12,7], huang_unmod_pm[12,8], huang_unmod_pm[12,9], huang_unmod_pm[12,10], huang_unmod_pm[12,11] + technical_rep_theme), 11, 11, xAxisLabels = c("Tm rep1_1", "Tm rep1_2", "Tm rep2_1", "Tm rep2_2", "Tm rep3_1", "Tm rep3_2", "Tm rep4_1", "Tm rep4_2" ,"Tm rep5_1", "Tm rep5_2", "Tm rep6_1"), yAxisLabels = c("Tm rep1_2", "Tm rep2_1", "Tm rep2_2", "Tm rep3_1", "Tm rep3_2", "Tm rep4_1", "Tm rep4_2", "Tm rep5_1", "Tm rep5_2", "Tm rep6_1", "Tm rep6_2")) + theme(strip.text = element_text(size = 6))
huang_mat_numeric <- data.frame( apply(huang_phospho_hq[,5:14], 2, as.numeric)) huang_pm <- ggpairs(huang_mat_numeric, upper = "blank", diag = NULL, lower = list(continuous = wrap(lowerFn))) ggmatrix(list(huang_pm[2, 1] + technical_rep_theme, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, huang_pm[3, 1], huang_pm[3,2], NULL, NULL, NULL, NULL, NULL, NULL, NULL, huang_pm[4,1], huang_pm[4,2], huang_pm[4,3] + technical_rep_theme, NULL, NULL, NULL, NULL, NULL, NULL, huang_pm[5,1], huang_pm[5,2], huang_pm[5,3], huang_pm[5,4], NULL, NULL, NULL, NULL, NULL, huang_pm[6,1], huang_pm[6,2], huang_pm[6,3], huang_pm[6,4], huang_pm[6,5]+ technical_rep_theme, NULL, NULL, NULL, NULL, huang_pm[7,1], huang_pm[7,2], huang_pm[7,3], huang_pm[7,4], huang_pm[7,5], huang_pm[7,6], NULL, NULL, NULL, huang_pm[8,1], huang_pm[8,2], huang_pm[8,3], huang_pm[8,4], huang_pm[8,5], huang_pm[8,6], huang_pm[8,7]+ technical_rep_theme, NULL, NULL, huang_pm[9,1], huang_pm[9,2], huang_pm[9,3], huang_pm[9,4], huang_pm[9,5], huang_pm[9,6], huang_pm[9,7], huang_pm[9,8], NULL, huang_pm[10,1], huang_pm[10,2], huang_pm[10,3], huang_pm[10,4], huang_pm[10,5], huang_pm[10,6], huang_pm[10,7], huang_pm[10,8], huang_pm[10,9] + technical_rep_theme), 9, 9, xAxisLabels = c("Tm rep1_1", "Tm rep1_2", "Tm rep2_1", "Tm rep2_2", "Tm rep3_1", "Tm rep3_2", "Tm rep4_1", "Tm rep4_2" ,"Tm rep5_1"), yAxisLabels = c("Tm rep1_2", "Tm rep2_1", "Tm rep2_2", "Tm rep3_1", "Tm rep3_2", "Tm rep4_1", "Tm rep4_2", "Tm rep5_1", "Tm rep5_2")) + theme(strip.text = element_text(size = 6))
pm <- ggpairs(nbf_tm_rep_df, 2:6, upper = "blank", diag = NULL, lower = list(continuous = wrap(lowerFnReg))) ggmatrix(list(pm[2, 1], NULL, NULL, NULL, pm[3, 1], pm[3,2], NULL, NULL, pm[4, 1], pm[4,2], pm[4,3], NULL, pm[5, 1], pm[5,2], pm[5,3], pm[5,4]), 4, 4, xAxisLabels = c("Tm rep1", "Tm rep2", "Tm rep3", "Tm rep4"), yAxisLabels = c("Tm rep2", "Tm rep3","Tm rep4","Tm rep5"))
pm_phospho <- ggpairs(phospho_tm_rep_df, 5:9, upper = "blank", diag = NULL, lower = list(continuous = wrap(lowerFnReg))) ggmatrix(list(pm_phospho[2, 1], NULL, NULL, NULL, pm_phospho[3, 1], pm_phospho[3,2], NULL, NULL, pm_phospho[4, 1], pm_phospho[4,2], pm_phospho[4,3], NULL, pm_phospho[5, 1], pm_phospho[5,2], pm_phospho[5,3], pm_phospho[5,4]), 4, 4, xAxisLabels = c("Tm rep1", "Tm rep2", "Tm rep3", "Tm rep4"), yAxisLabels = c("Tm rep2", "Tm rep3","Tm rep4","Tm rep5"))
colnames(huang_phospho_hq)[19] <- "deltaTm" delta_tm_df <- bind_rows( huang_phospho_hq %>% mutate(delta_tm = as.numeric(deltaTm)) %>% dplyr::select(delta_tm) %>% mutate(dataset = "Huang et al."), combo_mean_tm_significant_df %>% rowwise() %>% mutate(mean_delta_meltPoint = mean(c(delta_meltPoint_1, delta_meltPoint_2, delta_meltPoint_3, delta_meltPoint_4, delta_meltPoint_5), na.rm = TRUE)) %>% ungroup %>% dplyr::select(delta_tm = mean_delta_meltPoint) %>% mutate(dataset = "This study") ) ggplot(delta_tm_df, aes(delta_tm, fill = dataset)) + geom_density(alpha = 0.25) + xlab(expression(Delta*'T'[m])) + theme_bw() + theme(legend.position = "bottom")
mxq_unmod_df_huang <- read_delim("dat/Huang_nonphospho_evidence.txt", delim = "\t") mxq_unmod_df_potel <- read_delim("dat/Potel_nonphospho_evidence.txt", delim = "\t") mxq_phospho_df_huang <- read_delim("dat/Huang_phospho_evidence.txt", delim = "\t") mxq_phospho_df_potel <- read_delim("dat/Potel_phospho_evidence.txt", delim = "\t") mxq_combo_df <- bind_rows( mxq_unmod_df_huang %>% dplyr::select(PIF) %>% mutate(dataset = "Huang et al. non-modified", PIF = as.numeric(PIF)), mxq_phospho_df_huang %>% dplyr::select(PIF) %>% mutate(dataset = "Huang et al. enriched phosphopeptides", PIF = as.numeric(PIF)), mxq_unmod_df_potel %>% dplyr::select(PIF) %>% mutate(dataset = "This study non-modified", PIF = as.numeric(PIF)), mxq_phospho_df_potel %>% dplyr::select(PIF) %>% mutate(dataset = "This study enriched phosphopeptides", PIF = as.numeric(PIF)) ) ggplot(mxq_combo_df, aes(dataset, PIF)) + geom_violin(alpha = 0.25, fill = "gray", color = "gray") + geom_boxplot(width = 0.15, outlier.shape = NA) + theme_bw() + labs(x = "") + theme(legend.position = "bottom", axis.text.x = element_text(angle = 90, hjust =1))
nbf_mean_tm_df <- nbf_tm_rep_fil %>% ungroup %>% rowwise() %>% mutate(mean_Tm = mean(c(`1`, `2`, `3`, `4`, `5`), na.rm = TRUE)) %>% ungroup %>% dplyr::select(gene_name, mean_Tm) huang_mapping <- read_delim("dat/huang_mapping.txt", delim = "\t") huang_unmod_mean_tm_df <- huang_bulk %>% left_join(huang_mapping %>% dplyr::select(uniprot_id = From, gene_name = To), by = "uniprot_id") %>% dplyr::select(gene_name, huang_mean_tm = average_tm_nbf) %>% group_by(gene_name) %>% summarize(huang_mean_tm = mean(huang_mean_tm, na.rm = TRUE)) %>% ungroup %>% filter(!is.nan(huang_mean_tm)) tan_hq_meltCurve_param_df <- readRDS("dat/tan_hq_meltCurve_param_df.rds") %>% dplyr::select(gene_name, tan_Tm = meltPoint) becher_cell_hela <- read_xlsx("dat/1-s2.0-S0092867418303854-mmc4.xlsx", sheet = 2) %>% dplyr::select(gene_name, becher_cell_hela_mean_Tm = Tm.median) %>% mutate(becher_cell_hela_mean_Tm = as.numeric(becher_cell_hela_mean_Tm)) combo_tm_df <- tan_hq_meltCurve_param_df %>% full_join(becher_cell_hela, by = "gene_name") %>% full_join(nbf_mean_tm_df, by = "gene_name") %>% full_join(huang_unmod_mean_tm_df, by = "gene_name") pm_others <- ggpairs(combo_tm_df, 2:5, upper = "blank", diag = NULL, lower = list(continuous = wrap(lowerFnReg))) ggmatrix(list(pm_others[2, 1], NULL, NULL, pm_others[3, 1], pm_others[3,2], NULL, pm_others[4, 1], pm_others[4,2], pm_others[4,3]), 3, 3, xAxisLabels = c("Tan et al. (K562)", "Becher et al. (HeLa)", "This study (HeLa)"), yAxisLabels = c("Becher et al. (HeLa)", "This study (HeLa)", "Huang et al. (HEK293T)")) + theme(strip.text = element_text(size = 8))
Delta Tm comparison between the two studies
mean_deltaTm_df <- combo_mean_tm_pSite_anno %>% rowwise() %>% mutate(mean_delta_meltPoint = mean( c(delta_meltPoint_1, delta_meltPoint_2, delta_meltPoint_3, delta_meltPoint_4, delta_meltPoint_5), na.rm = TRUE)) %>% ungroup %>% dplyr::select(gene_name, Gene_pSite, mod_sequence, phospho_site_STY, mean_delta_meltPoint, significant) delta_tm <- left_join( mean_deltaTm_df %>% dplyr::select(Gene_pSite, mean_delta_meltPoint, significant), huang_phospho_hq_potel_anno %>% mutate(mean_delta_meltPoint_huang = average_tm_phospho - average_tm_nbf) %>% dplyr::select(Gene_pSite, mean_delta_meltPoint_huang, p_value), by = "Gene_pSite" ) %>% na.omit() ggplot(delta_tm, aes(mean_delta_meltPoint_huang, mean_delta_meltPoint)) + geom_point(alpha = 0.5, shape = 16) + geom_point(data = filter(delta_tm, significant), color = "#b2182b", shape = 16) + geom_point(data = filter(delta_tm, p_value < 0.05), color = "#1f78b4", shape = 16) + geom_point(data = filter(delta_tm, p_value < 0.05, significant), color = "#ff7f00", shape = 16) + geom_smooth(method = "lm", color = "darkgray", size = 1) + ggrepel::geom_text_repel(aes(label = Gene_pSite), data = filter(delta_tm, p_value < 0.05, significant)) + ggpmisc::stat_poly_eq(formula = y ~ x, parse = TRUE) + labs(x = expression(Delta*' T'[m]^Potel*' '* '('*~degree*C*')'), y = expression(Delta* ' T'[m]^Huang*' '* '('*~degree*C*')')) + theme_bw()
Blue dots correspond to phosphosites found significantly shifted by Huang et al., red dots were found significantly shifted by Potel et al., salmon coloured dots were found by both.
sessionInfo()
References
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.