PCA_analysis_complete_paralog.R

######################################
## Full-lenght paralog PCA analysis ##
## author : Pascale Lemieux         ##
## Date : 2023-01-10                ##
######################################

library(tidyverse)
library(magrittr)
library(gtools)
library(mixtools)

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

##############################################
## Import raw data generated by pyphe-quant ##
##############################################
## Recheck with rawdata from pyphe

#  For prey array plates, final diploid selection plates and final MTX plates #
# Read array plan in 1536 format
plate1536 <-
  read.csv('~/ancSH3_paper/SupplementaryMaterial/Data/plate1536_complete_paralog.csv')

# Set working directory with raw data from pyphe (for diploid S2 and MTX2)
setwd("~/ancSH3_paper/SupplementaryMaterial/Data/complete_paralog_RD/S2_diploid/")
file = dir(pattern = ".csv")

# Assign names to each data frame object
j = 1
for (i in file) {
  assign(x = paste0('plate', j), read.csv(i))
  j = j + 1
  
}

# merge each data frame with array plan
n = 48 #number of data frame
  for (i in 1:n) {
    p <- get(paste0('plate', i))
    assign('p',
           merge(
             p,
             plate1536,
             by.x = c('row', 'column'),
             by.y = c('row', 'col'),
             all = T
           ))
    
    p$plate <- i
    assign(paste0('plate', i), p)
  }

# assemble all data frame
all_plate <- 
  do.call(rbind.data.frame, mget(ls(pattern = "^plate.{1,2}$")))

# write file with data of each selection step
# MTX2_data.csv , diploid2_data.csv or DHFR3_data.csv
filename = '~/ancSH3_paper/SupplementaryMaterial/Data/complete_paralog_data/diploid2_data.csv'
  
write.csv(all_plate,
             file = filename,
             row.names = T,
             quote = T)

######################################
## Data verification  and filtering ##
######################################

#  Removal of non-growing position in the prey array and
#  diploid selection for the MTX data  #

# Set working directory where all the output files will be generated
setwd('~/ancSH3_paper/SupplementaryMaterial/Data/complete_paralog_data/')

# Read assembled data 
MTX2_data <- 
  read.csv(file = './MTX2_data.csv')

diploid2_data <- 
  read.csv(file = './diploid2_data.csv')

DHFR3_data <- 
  read.csv(file = './DHFR3_data.csv')


# Verification that all absent positions in plate DHFR3 96/384  
# are also absent 1536 DHFR3 plate

ngrowth <- DHFR3_data[is.na(DHFR3_data$area),c('row', 'column', 'orf')]

# Histogram of colony area distribution of prey array
ggplot(data = DHFR3_data)+
  geom_histogram(mapping = aes(x = area),
                 binwidth = 20)+
  theme_minimal()

table(ngrowth$orf)
# YDR480W YLR078C YPL049C 
# 5       5       4 


# Verification of no growth positions in diploid2 selection plates
# Histogram of colony area distribution of diploid selection
ggplot(data = diploid2_data)+
  geom_histogram(mapping = aes(x = area),
                 binwidth = 20)+
  xlim(0, 8000)+
  geom_vline(xintercept = 500)+
  theme_minimal()


# Determination of positions to withdraw for each plate. Each orf has to be
# present in two replicates on each plates to be conserved in the analysis
# Important to verify that the absence of growth is consistent with the prey 
# tested and not technical errors. Remove all non-growing positions and remove
# all prey wich has less than 2 replicates per plate

diploid2.data.na <- diploid2_data
diploid2.data.na[, 'area'] <- na.replace(diploid2.data.na[ , 'area'], 0)
diploid2.data.na$growth <- diploid2.data.na$area > 500


# Frequence of total non-growing diploid for each orf per plate
pos <- table(diploid2.data.na[!diploid2.data.na$growth, c('orf', 'plate')])
pos <- as.data.frame(pos)

# Frequence of total expected diploid for each orf per plate
orf<- diploid2.data.na[!diploid2.data.na$growth, 'orf']
  
orf.to.with <-  
  table(diploid2.data.na[diploid2.data.na$orf %in% orf , c('orf', 'plate')])
orf.to.with <- as.data.frame(orf.to.with)

# Determine which orf has less than 2 replicate by plate
orf.to.with$ngrowth <- pos$Freq
orf.to.with$with <- orf.to.with$Freq - orf.to.with$ngrowth < 2

pos.to.with <- 
  diploid2.data.na[diploid2.data.na$area < 500, c('row', 'column', 'plate')]
pos.to.with$with <- TRUE

# Withdrawal of the nongrowing orf in each diploid2 plate

diploid2_data.tresh <- 
  merge(x = diploid2_data,
        y = orf.to.with, 
        by = c('orf', 'plate'), 
        all = TRUE)

diploid2_data.tresh <- 
  merge(x = diploid2_data.tresh, 
        y = pos.to.with, 
        by = c('row', 'column', 'plate'), 
        all = TRUE)

diploid2.data.clean <- 
  subset(diploid2_data.tresh, 
         subset = (diploid2_data.tresh$with.x | diploid2_data.tresh$with.y) %in% c(FALSE, NA), 
         select = colnames(diploid2_data))

nrow(diploid2_data.tresh) - nrow(diploid2.data.clean)
#2426

# Comparison 
ggplot(data = diploid2.data.clean)+
  geom_histogram(mapping = aes(x = area),
                 binwidth = 20, alpha = 0.7, fill = 'red')+
  geom_histogram(data = diploid2_data, mapping = aes(x = area),
                 binwidth = 20, alpha = 0.7)+
  xlim(0,8000)+
  theme_minimal() 


# Withdrawal of nongrowing orf on diploid2 plates in each MTX2 plate 
MTX2.data.tresh <- merge(x = MTX2_data, 
                         y = orf.to.with, 
                         by = c('orf', 'plate'), 
                         all = TRUE)

MTX2.data.tresh <- merge(x = MTX2.data.tresh,
                         y = pos.to.with, 
                         by = c('row', 'column', 'plate'), 
                         all = TRUE)

MTX2.data.clean <- 
  subset(MTX2.data.tresh, 
         subset = (MTX2.data.tresh$with.x | MTX2.data.tresh$with.y) %in% c(FALSE, NA), 
         select = colnames(diploid2_data))

MTX2.data.clean$area <- 
  na.replace(MTX2.data.clean$area, 0)

# Histrogram of the distribution of MTX colony area
ggplot(data = MTX2.data.clean)+
  geom_histogram(mapping = aes(x = area),
                 binwidth = 20)+
  theme_minimal()

### Export clean MTX2 data
write.csv(x = MTX2.data.clean, 
          file = '~/ancSH3_paper/SupplementaryMaterial/Data/complete_paralog_data/MTX2_data_clean.csv', 
          quote = TRUE, 
          fileEncoding = 'UTF-8')

#############################
## MTX data transformation ##
#############################
library(viridis)

MTX2.data <- 
  read.csv(file = '~/ancSH3_paper/SupplementaryMaterial/Data/complete_paralog_data/MTX2_data_clean.csv', 
           fileEncoding = 'UTF-8', )

# Columns removal and transformation
MTX2.data$X <- NULL
MTX2.data$X.1 <- NULL
MTX2.data$X.2 <- NULL
MTX2.data$plate <- as.factor(MTX2.data$plate)

# Removal of unwanted plates (growth problem or other strains)
MTX2.data <- subset(MTX2.data, 
                    subset = !(MTX2.data$plate %in% c(2, 10, 26, 34)))

# Adjust area = 0 to area = 1 for logarithm transformation
MTX2.data %<>% mutate(adjust = area + 1)

# Transform colony adjusted area in log2 scale
MTX2.data$log2.area <- log2(MTX2.data$adjust)

#Distribution of log2(colony area) per plate (Supplementary Figure 2A)
sf2a <- 
ggplot(data = MTX2.data) +
  geom_density(mapping = aes(x = log2.area, color = plate))+ 
  theme_minimal()+
  theme(legend.position = 'none')+
  theme_minimal()+
  geom_vline(xintercept = 3, color = 'grey20', linewidth = 1)+
  geom_vline(xintercept = 13.13, color = 'grey20', linewidth = 1)+
  geom_vline(xintercept = 8.60733, color = 'orangered', linewidth = 1)+
  theme(legend.position = 'none')+
  xlim(0, 14)+
  theme_bw()+
  xlab(expression(paste(log[2], '(colony area)')))+
  ylab('Density')+
  scale_color_viridis(option = 'G', discrete = T)+
  theme(legend.position = 'none',
        legend.text = element_text(size = 14),
        legend.title = element_text(size = 16),
        axis.title = element_text(size =14), 
        axis.text = element_text(size =12), 
        strip.text.x = element_text(size = 16), 
        panel.grid.major = element_blank(), panel.grid.minor = element_blank())

saveRDS(sf2a, file = '~/ancSH3_paper/SupplementaryMaterial/FigurePanels/SuppFig2A.rds')

# Remove the empty positions
MTX2.data <-  subset(MTX2.data, 
                     MTX2.data$log2.area >= 3)

# Remove aberrant data
MTX2.data <-  subset(MTX2.data, 
                     MTX2.data$area < 9000)


# Import bait description
bait <- read.csv2('~/ancSH3_paper/SupplementaryMaterial/Data/complete_paralog_data/bait_description.csv')
colnames(bait) <- c('plate', 'sh3_sequence', 'Bait.Standard_name', 'tech_rep')

MTX2.data <- 
merge(bait, 
      MTX2.data, 
      by ='plate')

MTX2.data$strain <- str_c(MTX2.data$sh3_sequence, MTX2.data$Bait.Standard_name, sep = '.')

saveRDS(MTX2.data, file = '~/ancSH3_paper/SupplementaryMaterial/Data/complete_paralog_data/MTX2_trans.rds')

###############################
## Normalization : PCA score ##
###############################
library(ggpubr)

# Important to compare the plates with each other
MTX2.list <- readRDS('~/ancSH3_paper/SupplementaryMaterial/Data/complete_paralog_data/MTX2_trans.rds')

# Withdrawal of the border
#MTX2.list %<>% filter(orf != 'border')  


# Verification if there is bias between the baits  
exbiais <- 
  MTX2.list[, c(2,3,7,13,15)]

# Withdrawal of the border for biais check
exbiais %<>% filter(orf != 'border')  

exbiais <- 
  exbiais[exbiais$sh3_sequence %in% c('MYO3_WT', 'MYO5_WT'), ]

exbiais$Bait.Standard_name <- 
  firstup(gsub(exbiais$Bait.Standard_name, pattern = 'p', replacement= ''))

exbiais$sh3_sequence <- firstup(gsub(exbiais$sh3_sequence, pattern = '_WT', replacement = ''))

compare_means(area~Bait.Standard_name, exbiais)
# bias confirmed

sf2b <- 
  ggplot(exbiais)+
  geom_violin(aes(y = log2.area, x = Bait.Standard_name, 
                 fill = Bait.Standard_name), 
              fill = 'grey60',alpha = 0.4, color = 'grey20')+
  stat_compare_means(aes(x = Bait.Standard_name, y = log2.area), 
                     method = 'wilcox.test', size = 4.5, label.x = 1)+
  ylim(0,15)+
  ylab(expression(paste(log[2], '(colony area)')))+
  xlab('WT bait')+
  theme_bw()+
  theme(
    #legend.position = 'none', 
    legend.title = element_blank(), 
        legend.text = element_text(size = 14),
        #legend.title = element_text(size = 16),
        axis.title = element_text(size =14), 
        axis.text = element_text(size =12), 
        strip.text.x = element_text(size = 16), 
        panel.grid.major = element_blank(), panel.grid.minor = element_blank())

saveRDS(sf2b, file = '~/ancSH3_paper/SupplementaryMaterial/FigurePanels/SuppFig2B.rds')   


# Calculation of min/max for each plate
  MTX2.list %<>% 
    group_by(plate) %>% 
    mutate(max = max(log2.area, na.rm = TRUE))%>%
    ungroup()
  
  MTX2.list %<>% 
    group_by(plate) %>% 
    mutate(min = min(log2.area, na.rm = TRUE))%>%
    ungroup() 
  
x <- MTX2.list[MTX2.list$max == MTX2.list$log2.area, ]
  
# Calculation of standardized data using the min/max calculated 
 MTX2.list %<>% 
  mutate(norm = ((log2.area - min)/(max-min)))

 # Withdrawal of the border
 MTX2.list %<>% filter(orf != 'border')  
 
# Verification of the correlation between technical replicates

   dfrep1 = MTX2.list %>% filter(tech_rep == 1)
   dfrep2 = MTX2.list %>% filter(tech_rep == 2)
   
   dfrep1 %<>% mutate(baitXprey_1 = interaction(strain, orf)) %>% select(baitXprey_1, norm, log2.area)
   
   dfrep2 %<>% mutate(baitXprey_2 = interaction(strain, orf)) %>% select(baitXprey_2, norm, log2.area)
   
   crep <-  merge(dfrep1,
                  dfrep2,
                  by.x="baitXprey_1",
                  by.y="baitXprey_2")

# Supplementary Figure 2C 
sf2c <-    
   ggplot(data = crep, aes(x = norm.x, y = norm.y)) +
     geom_point(aes(x = norm.x, y = norm.y), alpha = 0.2)+
     xlab('PCA score replicate 1')+
     ylab('PCA score replicate 2')+
     lims(y = c(0,1), x = c(0,1))+
     stat_cor(size = 4.5, cor.coef.name = c('r'), label.y = 1)+
     theme_bw()+
     theme(legend.position = 'none', legend.title = element_blank(), 
           legend.text = element_text(size = 14),
           axis.title = element_text(size =14), 
           axis.text = element_text(size =12), 
           strip.text.x = element_text(size = 16), 
           panel.grid.major = element_blank(), panel.grid.minor = element_blank())
   
saveRDS(sf2c, file = '~/ancSH3_paper/SupplementaryMaterial/FigurePanels/SuppFig2C.rds')   

# Export normalized data
saveRDS(MTX2.list, file = '~/ancSH3_paper/SupplementaryMaterial/Data/complete_paralog_data/MTX2_norm.rds')   

###########################
## PPI score computation ##
###########################

MTX2.list <-  readRDS(file = '~/ancSH3_paper/SupplementaryMaterial/Data/complete_paralog_data/MTX2_norm.rds')
MTX2.list.rep <- vector(mode = 'list', length = 1)

# Compute the median of biological replicates (unique prey-bait combination) (PPI_score)  
  MTX2.list.rep <-  
    MTX2.list%>%
    group_by(sh3_sequence, Bait.Standard_name, orf) %>%
    dplyr::summarise(n(), medsize = median(adjust ,na.rm = TRUE ), 
                     medrep = median(norm, na.rm = TRUE ))

# Transform columns
MTX2.list.rep %<>% mutate(strain= str_c(sh3_sequence, Bait.Standard_name, sep = '.'))
MTX2.list.rep %<>% mutate(log2medsize = log2(medsize))

# Create object with PCA and PPI scores
MTX2.list.rep <-   
  merge(MTX2.list ,
        MTX2.list.rep[, c('orf', 'strain', 'medsize','log2medsize','medrep','n()')], 
        by = c('orf', 'strain'))

# Establishing the threshold for median colony area size 
# Bellow the threshold is noise 

x = 390 #colony size threshold

ggplot(data = MTX2.list.rep ) +
  geom_density(aes(x = medsize, colour = strain)) +
  theme_minimal()+
  geom_vline(xintercept = x)+
  theme(legend.position = 'none')


# Threshold the data & remove preys with less than 20 data points

 for(j in unique(MTX2.list.rep$orf)){
    test <- 
      subset(MTX2.list.rep, 
             subset = orf == j)
    if(all(test$medsize < x)){
      MTX2.list.rep %<>% subset(subset = orf != j)
    } else if(nrow(test) < 20){
      MTX2.list.rep %<>% subset(subset = orf != j)
    }
  }  

# Export the data (used for PCA_SH3_dependency_test.R)
saveRDS(MTX2.list.rep , 
        file = '~/ancSH3_paper/SupplementaryMaterial/Data/complete_paralog_data/MTX2data_complete_paralog.rds')