Figure4.R

######################################
## Figure 4                         ##
## author : Pascale Lemieux         ##
## Date : 2023-06-06                ##
######################################
library(tidyverse)
library(ggplot2)
library(ggpubr)
library(ggrepel)
library(viridis)
library(dplyr)
library(gtools)
library(magrittr)
library(readxl)

firstup <- function(x) {
  substring(x, 2) <- tolower(substring(x, 2))
  x
}

#Read Table S7

dockingData <- 
  read.csv('~/ancSH3_paper/Reviews/SupplementaryMaterial/TableS7.csv')

# Confirmed motif in vivo
conf <- 
data.frame(prey = c('Lsb3', 'Myo5', 'Osh2', 'Pkh2', 'Sla1', 'Ste20', 'Bzz1', 'Lsb3', 'Osh2', 'Sla1', 'Ste20'), 
           para = rep(c('Myo3', 'Myo5'), each = 6)[-11], 
           conf = T)

# Merge in vivo confirmed motif information with the docking dat

dockingData <- 
merge( 
      dockingData, 
      conf,
      by.x = c('SH3.Standard_name', 'Motif.Standard_name'), 
      by.y = c('para', 'prey'), 
      all = T)

# replace NAs generated by the merge with FALSE (no confirmation in vivo)
dockingData$conf <- na.replace(dockingData$conf, replace = FALSE)

# Read Table S1
PCA_complete <- 
  read_xlsx('~/ancSH3_paper/Reviews/SupplementaryMaterial/TableS1.xlsx')

# Select PCA data only for the extantSH3
PCA_subset <-  
 unique(PCA_complete[PCA_complete$sh3_sequence %in% c('extantMyo3', 'extantMyo5'), c(1:4, 16)])


# Select results for the first-ranked cluster
rank1 <- dockingData[dockingData$cluster_rank == 1, ]

# Select information for the docking which used the target SH3 
# for each motif (ex. remove the docking that used a target Myo3
# SH3 with a Myo5 predicted motif)
x <- vector('logical', length = nrow(rank1))

for(i in 1:nrow(rank1)){
  
  if(rank1[i, "target_SH3"] == 'both'){
    x[i] <- T
  }else if(rank1[i, "target_SH3"] != 'both'){
    
    if(rank1[i, "target_SH3"] == rank1[i, "SH3.Standard_name"])
    x[i] <- T
  } else if(rank1[i, "target_SH3"] != rank1[i, "SH3.Standard_name"]){
    x[i] <- F
  }
  
} 

rank1 <-rank1[x,]

# Merge docking data with PCA data
PCA_dock <- 
merge(rank1, 
      PCA_subset,
      by.x = c('SH3.Standard_name', 'Motif.Standard_name'), 
      by.y = c('Bait.Standard_name', 'Prey.Standard_name'))


# Figure 4B : Correlation between PPI score & Energy Interaction

colnames(PCA_dock)[6] <- 'Interact'

Fig4B <- 
ggplot(PCA_dock, aes(x =med.PPI_score , y = Interact))+
  facet_grid(cols = vars(SH3.Standard_name))+
  geom_point(aes(color = conf), size = 3)+
  scale_color_manual(values = c('grey20', '#e5ca28ff'), 
                    labels = c('predicted motif', 'validated motif'))+
  xlim(0,1)+
  xlab('med. PPI score')+
  ylab(expression(paste('med. ' , Delta, 'G (kcal/mol)')))+
  theme_bw()+
    geom_text_repel(
    aes(label = Motif.Standard_name),
    size = 3.5,
    box.padding = unit(0.4, "lines"),
    point.padding = unit(0.3, "lines"),
    min.segment.length = 0,
    segment.curvature = -0.1,
    segment.ncp = 3,
    segment.angle = 20
  )+
  ylim(-13, 0)+
  stat_cor( method = 'spearman', size = 4.5, cor.coef.name = c('r'), label.y = -0.4, color = 'grey20')+
  theme(strip.text = element_text(size =16),  
        axis.text = element_text(size =12),
          legend.text = element_text(size = 14), 
          axis.title = element_text(size =14), 
            legend.position = 'bottom', 
        panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
guides(color = guide_legend(title = "",
                            override.aes = aes(label = ""), nrow = 1), 
       shape = guide_legend(title = 'DHFR[1,2] strain', 
                            override.aes = aes(label=''), nrow = 2))


# Read raw docking result from Data S4
raw_dG <- 
  read.csv('~/ancSH3_paper/Reviews/SupplementaryMaterial/SupplementaryData4/Docking_RD.csv',
           header = T, sep = ',')


raw_dG$preys <- as.factor(raw_dG$preys)

# Compute minimal energy per cluster of each docking
raw_dG %>%
  dplyr::group_by(preys, SH3, cluster) %>%
  dplyr::summarise(min_ener =  min(Interaction_Energy)) -> min_dG

min_dG$cluster <- as.factor(min_dG$cluster)

# Comput minimal energy for each docking
min_dG%>%
  dplyr::group_by(preys, SH3)%>%
  dplyr::summarise(min_cluster = which.min(min_ener)) -> min_cluster

# Times that the minimal energy is in the first-rank cluster
tab_clus <- table(min_cluster[, c(2,3)])
apply(tab_clus, 1, sum)
tab_clus[,2]/28
# The minimal energy per docking is not necessarily in the first-ranked cluster

# Verify if the median of the first-ranked cluster is the lowest 
raw_dG %>%
  dplyr::group_by(preys, SH3, cluster) %>%
  dplyr::summarise(med_ener =  median(Interaction_Energy)) -> med_dG


med_dG%>%
  dplyr::group_by(preys, SH3)%>%
  dplyr::summarise(min_cluster = which.min(med_ener)) -> min_med_cluster

table(min_med_cluster[, c(2,3)])/28
# Median energy is a better indicator than the lowest energy  for each docking

library(olsrr)
library(ggpmisc)
library(reshape2)

pre1 <-setdiff(names(med_dG), c("med_ener", "SH3"))
mod_dres<-NULL
for (j in pre1) {
  model  <- lm( med_ener ~ cluster+ get(j), data = med_dG)
  bmodel <- broom::tidy(model)
  bmodel$term[3]<-j
  bmodel<-bmodel[3,]
  mod_dres<-rbind(mod_dres,bmodel)
}
mod_dres

pre1.plot = melt(med_dG, id.vars = c('med_ener', 'cluster')) %>% 
  dplyr::filter(variable == 'SH3')


FigSupp4B <-  
ggplot(pre1.plot, aes(x = cluster, y = med_ener))+
  geom_violin(aes(x = cluster, y = med_ener, group = cluster), fill = 'white', color = 'grey20')+
  #facet_grid(~SH3, scales = 'free')+
  #geom_boxplot(width=0.1, fill ='grey80', outlier.colour = 'transparent')+
  geom_jitter(aes(color = value), width = 0.2, alpha = 0.5, na.rm = T) +#, aes(color = preys))+
  scale_color_manual(values = c('#551A8B', '#3366CC', 'orangered'))+
  geom_smooth(aes(cluster, med_ener, colour=variable),
            formula = y ~ x, method="glm", se=FALSE,
            colour="grey20", na.rm = T)+
  scale_x_discrete(limits = c(1,2,3,4,5,6,7,8,9))+
  stat_poly_eq(formula = y ~ x, 
                        aes(x = cluster, y = med_ener, label = paste(after_stat(eq.label),
                                                                     after_stat(p.value.label),
                                                            sep = "~~~~")), 
                        parse = TRUE, label.y = 0.98, small.p = T)+
  xlab('Docking clusters')+
  ylab(expression(paste('med. ', Delta, 'G (kcal/mol)')))+
  theme_bw()+
  ylim(-15, 1)+
  theme(legend.position = 'bottom',
        strip.text = element_text(size =16),
        axis.text = element_text(size =12),
        legend.text = element_text(size = 14),
        axis.title = element_text(size =14),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  guides(color = guide_legend(title =''))


# Figure 4C : Compare the distribution of median Energy between the SH3s

library(rstatix)

# Prepare the data for plotting
c1_ener <- subset(raw_dG, subset = rank ==1)
table(c1_ener[, c(5,6)])
c1_ener$bait <- 
  matrix(unlist(strsplit(as.character(c1_ener$preys), split = '_', fixed = T)), ncol = 2, byrow = T)[,2]
c1_ener$bait <- as.factor(c1_ener$bait)
c1_ener$preys <- as.factor(c1_ener$preys)
c1_ener$SH3 <- as.factor(c1_ener$SH3)

# Is there a difference between the energy distribution between the SH3s
kruskal.test(Interaction_Energy ~ SH3, data = c1_ener)
#not significantly,  p-value = 0.069, close enough to check the next test
d_res <- 
  dunn_test(Interaction_Energy ~ SH3, data = c1_ener, p.adjust.method= 'bonferroni')

#.y.          group1 group2    n1    n2 stati…¹      p  p.adj
#* <chr>        <chr>  <chr>  <int> <int>   <dbl>  <dbl>  <dbl>
#  1 Interaction… AncC   Myo3     275   276   0.471 0.638  1     
#2 Interaction… AncC   Myo5     275   274  -1.72  0.0852 0.256 
#3 Interaction… Myo3   Myo5     276   274  -2.19  0.0283 0.0849

kruskal.test(Interaction_Energy ~ SH3, data = c1_ener[c1_ener$SH3 %in% c('Myo3', 'Myo5'), ])


# Statistical tests with median Energy
my_comparisons <- list( c("AncD", "Myo3"), c("AncD", "Myo5"), c("Myo3", "Myo5") )

med_c1_clust <- 
  aggregate(Interaction_Energy~SH3+preys, data = c1_ener, FUN = 'median')

pval <- 
  compare_means(Interaction_Energy ~ SH3, med_c1_clust, 
                method = 'wilcox', p.adjust.method = 'BH', paired = T)

pval$my.pval <- 
  paste0('p = ', pval$p.adj)

confi <- c('YHL007C', 'YMR109W', 'YBL007C', 'YDL019C', 'YFR024C-A', 'YOL100W', 'YHR114W')
med_c1_clust$conf <- 
matrix(unlist(strsplit(x = as.character(med_c1_clust$preys), split = '_', fixed = T)), ncol = 2, byrow = T)[,1 ] %in% confi

# Figure 4
Fig4C <- 
ggplot(med_c1_clust, aes(x = SH3, y = Interaction_Energy))+
  geom_violin( alpha = 0.8, fill = 'grey90', size = 1, width = 0.3, linewidth = 0.75)+
  xlim('Myo3', 'AncD', 'Myo5')+
  ylim(-12, 2)+
  #scale_fill_manual(values = c('#551A8B', '#3366CC', 'orangered'))+
  ylab(expression(paste('med. ',  Delta, 'G (kcal/mol)')))+
  #stat_cor(method = 'spearman', fontface = 'plain', size = 4, label.x = -7.5)+
  geom_jitter(aes(color = conf), size = 2, alpha = 0.7, width = 0.1)+
  #geom_line(aes(group = preys, color = conf), size =0.7, alpha = 0.7)+
  scale_color_manual(values = c('grey20', '#e5ca28ff'))+
  #geom_line(data = subset(med_c1_clust, preys == 'YHL007C_both'), aes(group = preys) ,color = 'orangered', size = 1)+
  #geom_point(data = subset(med_c1_clust, preys == 'YHL007C_both') ,color = 'orangered', size = 2)+
  stat_pvalue_manual(pval, label = "my.pval", y.position= c(-0.2, 0.5, 1.4), size = 4.5)+
  theme_bw()+
  theme(legend.position = 'none',
        strip.text = element_text(size =14),  
        axis.text = element_text(size =12),
        legend.text = element_text(size = 16), 
        axis.title = element_text(size =14), 
        panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ggsave('~/ancSH3_paper/Reviews/Figure4C.png', width = 4, height = 4)

# Arrange the Figure 4
library(cowplot)
library(svglite)

struc <- ggdraw()+
  draw_image('~/ancSH3_paper/Reviews/FigurePanels/Fig4A2.tif')

plot_grid(struc, Fig4B, Fig4C, nrow = 1,
          labels = 'AUTO', label_fontface = 'plain', label_size = 16, 
          align = 'h', axis = 'b', rel_widths = c(1,2,1))

ggsave('~/ancSH3_paper/Reviews/Figure4.png', width = 14, height = 4)



# Arrange Supplementary Figure 4

airs <- 
ggdraw()+
  draw_image('~/ancSH3_paper/Reviews/FigurePanels/FigSupp4A.tif')


plot_grid(airs, FigSupp4B, 
          nrow = 1, ncol = 2, labels = c('A','B'), 
          label_fontface = 'plain', label_size = 16, rel_widths = c(1.2, 2))

           
    
ggsave('~/ancSH3_paper/Reviews/SupplementaryMaterial/FigureS4.svg', 
       width = 12, height = 5)