LMNB1 stainings.Rmd

title: "R Notebook"
output: html_notebook
library(ggplot2)
library(pheatmap)
library(RColorBrewer)
#display.brewer.all()
library(tidyverse)
library(ggpubr)
library(ggsignif)
library(openxlsx)
library(viridis)
library(data.table)
library(readxl)
library(janitor)
library(dplyr)
library(jcolors)
library(stringr)


setwd("//atlas.uni.lux/users/isabel.rosety/GBA/Stainings/Plots/")

data1<-read.csv("//atlas.uni.lux/users/isabel.rosety/GBA/Stainings/Plots/488TUJ1_568LMNB1_Sox2647/IR_20211208_488TUJ1_568LMNB1_647Sox2_RS_e41e42e45d30_20211216_094400.csv",sep = ",")
data2<-read.csv("//atlas.uni.lux/users/isabel.rosety/GBA/Stainings/Plots/488TUJ1_568LMNB1_Sox2647/IR_20211208_488TUJ1_568LMNB1_647Sox2_RS_e43e35to36d30_20211208_110531.csv",sep = ";")

data3<-read.csv("//atlas.uni.lux/users/isabel.rosety/GBA/Stainings/Plots/488TUJ1_568LMNB1_Sox2647/IR_20211208_488TUJ1_568LMNB1_647Sox2_RS_e46e48d30_20220730_160215.csv",sep = ";")


data = bind_rows(data1,data2,data3)


tmp=do.call(rbind,strsplit(data$AreaName,"_")) #stringsplit will take the column Areaname and will split this depending on the underscore, then we put them otgeter one ofter the other


data$Condition = tmp[,1]

data %>%
   mutate(Condition = str_replace_all(Condition,
            pattern = "MUT", replacement = "GBA-PD")) %>%
     mutate(Condition = str_replace_all(Condition,
            pattern = "WT", replacement = "CTRL")) ->data
data %>%
  mutate(CellLine = AreaName) ->data

data %>%
   mutate(CellLine = str_replace_all(CellLine,
            pattern = "WT_56", replacement = "CTRL1")) %>%
    mutate(CellLine = str_replace_all(CellLine,
            pattern = "WT_39", replacement = "CTRL2"))  %>%
    mutate(CellLine = str_replace_all(CellLine,
            pattern = "WT_68", replacement = "CTRL3")) %>%
    mutate(CellLine = str_replace_all(CellLine,
            pattern = "MUT_309", replacement = "PD1")) %>%
    mutate(CellLine = str_replace_all(CellLine,
            pattern = "MUT_KTI6", replacement = "PD2")) %>%
    mutate(CellLine = str_replace_all(CellLine,
            pattern = "MUT_SGO1", replacement = "PD3")) ->data

toselect = c("Barcode","AreaName","Condition","Batch","OrganoidIdx","Day","Batch", "LMNB1hMFIinSox2")


data %>% #same
  dplyr::select(toselect)->data
#feature_names<-colnames(data[,7:ncol(data)])
feature_names<-"LMNB1hMFIinSox2"
#Mean of replicates
data%>%
    pivot_longer(feature_names,names_to = "Features", values_to = "Measure") -> data_all_long
data2 = drop_na(data_all_long)

data2 %>%
  group_by(Condition,AreaName,Batch,Features, Barcode,Day) %>%
  summarize(MeanFeatures = mean(Measure))-> data_Mean

#Normalizing to mean of controls
data_Mean %>%
    dplyr::filter(Condition=="CTRL") %>%
    #group_by(Condition, Day, Features,Batch) %>%
    group_by(Condition, Features,Batch,Day) %>%
    summarise(baseline = mean(MeanFeatures)) %>%
    ungroup() %>%
    dplyr::select(-Condition) %>%
    #full_join(data_Mean, by=c("Day","Features","Batch")) %>%
    full_join(data_Mean, by=c("Features","Batch","Day")) %>%
    mutate(Foldchange = MeanFeatures /baseline) ->data_all_based
    included_vars =c("AreaName", "Condition","Batch","Barcode","Day") # here we already have the mean of each section
    data_all_based %>%
    pivot_wider(all_of(included_vars),names_from = Features,values_from = Foldchange)->data_all_based_wide2

#write.csv(data_all_based_wide2, file = 'All_data_normalized.csv')

#If I want to get the wide table without normalizing to mean of controls:
data_Mean %>%
   full_join(data_Mean, by=c("Features","MeanFeatures","Batch","AreaName", "Barcode", "Condition","Day")) ->data_Mean
included_vars =c("AreaName", "Condition","Batch","Barcode","Day")
data_Mean %>%
    pivot_wider(all_of(included_vars),names_from = Features,values_from = MeanFeatures)->data_Mean_wide

#write.csv(data_Mean_wide, file = 'All_data_not_normalized.csv')


for (i in 1:length(feature_names)) {
data_all_based_wide2 %>%
  pivot_longer(cols=feature_names, names_to = "feature", values_to = "value") %>%
  filter(feature %in% feature_names[i]) %>%
  ggplot( aes(x=factor(Condition, level = c("CTRL", "GBA-PD")), y=value)) +
  #geom_violin( aes(fill=Condition),show.legend = T, trim=T),
  geom_violin( aes(fill=Condition),show.legend = T,scale = "width", trim=F)+
geom_dotplot(binaxis = "y",stackdir = "center",dotsize=0.8)+
  #scale_fill_manual(values= c("#bdd7e7","#2171b5"),name = "Condition", guide = FALSE)+
    scale_fill_manual(values= alpha(c("#1565C0","#CC0000"),0.75),name = "Condition",guide = "none")+
    #geom_point(aes(color=Batch),shape=18, size=3,show.legend = T)+
    #scale_color_manual(values = rev(brewer.pal(n=6, name="OrRd")))+
    scale_color_jcolors("pal7")+
    #scale_color_viridis(option = "D", discrete=TRUE)+
    #geom_point(shape = 1,size = 3,colour = "black")+
    theme(legend.key=element_blank()) +

   geom_signif(comparisons = list(c("CTRL", "GBA-PD")), test='wilcox.test',
              vjust=0.6, size=0.5, margin_top=0.5, textsize=9, map_signif_level=c("***"=0.001, "**"=0.01, "*"=0.05,  " "=2) ) +
  #facet_grid(~fct_relevel(Day, "d30","d60"), scales="free") +
      labs(x     = "",
       y     = paste(feature_names[i]),
       #y     = paste(names[i]),
       fill  = "Condition",
       #title = paste(feature_names[i])) +
       title = "" )+
  theme_bw() +
  theme(
    axis.line = element_line(colour = 'black', size = 0.5) ,
    axis.title.x = element_blank(),
    axis.text.x = element_text(size=12, color="black"),
    axis.title.y = element_text(size = 12),
    axis.text.y = element_text(size=10, color="black"),
    axis.ticks.y = element_line(),
    axis.ticks.length=unit(.25, "cm"),
    #change legend text font size)
    #legend.key.size = unit(0.7, "cm"),
    #legend.key.width = unit(0.6,"cm"),
    legend.key=element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    plot.title = element_text(size = 20, hjust=0.5, vjust= 1, face = "bold"),
    plot.subtitle = element_blank(),#element_text(size = 2, hjust=0.5)
    strip.text = element_text(size=12, vjust=0.5),
    strip.background = element_rect(fill="lightgray"),
   # panel.border = element_rect(fill = NA, color = "black"),
    panel.spacing.y = unit(0.8, "lines"),
    strip.switch.pad.wrap=unit(20, "lines"),
    legend.position="right",
    legend.text = element_text(size=17),
    legend.title = element_text(size=19)

  )  -> p
  #t<- cowplot::ggdraw(cowplot::add_sub(p, "Wilcox-test, ***p=0.001, **p=0.01, *p=0.05",hjust=-0.2, size=13))
  print(p)
  ##ggsave(paste0(Sys.Date(),"_", names[i], ".pdf"), plot=t)
#ggsave(paste0(Sys.Date(),(sprintf("Plot_%s.pdf",feature_names[i]))),height=3.5,width=3)
}