library(Seurat)
library(dplyr)

BM.combo <- readRDS("/DATA/31/molteni/vexas_bm/results/integration/pc15/BM.combo_anno.rds")

BM.combo@meta.data$mutations <- with(BM.combo@meta.data,
                                     ifelse(BM.combo@meta.data$orig.ident == "BM-01","THR",
                                            ifelse(BM.combo@meta.data$orig.ident == "BM-02","VAL",
                                                   ifelse(BM.combo@meta.data$orig.ident == "BM-03","VAL", 
                                                          ifelse(BM.combo@meta.data$orig.ident == "BM-04","LEU",
                                                                 ifelse(BM.combo@meta.data$orig.ident == "BM-08","VAL",
                                                                        ifelse(BM.combo@meta.data$orig.ident == "BM-09","LEU",
                                                                               ifelse(BM.combo@meta.data$orig.ident == "PT1","VAL",
                                                                                      ifelse(BM.combo@meta.data$orig.ident == "PT2","THR",
                                                                                             ifelse(BM.combo@meta.data$orig.ident == "PT3","THR",
                                                                                                    ifelse(BM.combo@meta.data$orig.ident == "PT4","THR",
                                                                                                           ifelse(BM.combo@meta.data$orig.ident == "PT5","THR",
                                                                                                                  ifelse(BM.combo@meta.data$orig.ident == "PT6","LEU",
                                                                                                                         ifelse(BM.combo@meta.data$orig.ident == "PT7","LEU",
                                                                                                                                ifelse(BM.combo@meta.data$orig.ident == "PT8","THR",
                                                                                                                                       ifelse(BM.combo@meta.data$orig.ident == "PT9","THR",
                                                                                                                                                     ifelse(BM.combo@meta.data$orig.ident == "HD1","CTRL",
                                                                                                                                                            ifelse(BM.combo@meta.data$orig.ident == "HD2","CTRL",
                                                                                                                                                                   ifelse(BM.combo@meta.data$orig.ident == "HD3","CTRL",
                                                                                                                                                                          ifelse(BM.combo@meta.data$orig.ident == "HD4","CTRL","NA"))))))))))))))))))))
#"BM-01","BM-02","BM-03","BM-04", "BM-08","BM-09",

#COMPARISON THR VEXAS PATIENTS vs HEALTHY CONTROLS WU

#select only the subjects involved in the comparison by mutations
Idents(BM.combo) <- "orig.ident"
BM.combo_subset <- subset(BM.combo, idents = c("BM-01","HD1", "HD2", "HD3", "HD4"))
clusters <- unique(BM.combo_subset$celltypes_annot)
Idents(BM.combo_subset) <- "celltypes_annot"

setwd("/DATA/31/molteni/vexas_bm/results/integration/DE_markers_OSR_Wu/Thr_CTRL/")

for(i in clusters) {
  cells_group1 = sum(BM.combo_subset$mutations == "THR" & BM.combo_subset$celltypes_annot == i)
  cells_group2 = sum(BM.combo_subset$mutations == "CTRL" & BM.combo_subset$celltypes_annot == i)
  if(cells_group1 >= 3 & cells_group2 >= 3 ) {
    markers <- FindMarkers(BM.combo_subset, ident.1 = "THR", group.by = "mutations", subset.ident = i)
    markers['cluster'] <- i
    filename=paste("markers_cluster_",i,".csv", sep="")
    write.csv(markers, file=filename)
  }
  else {
    print(paste("Skipping comparison for", i, "due to insufficient cell count in one or both groups."))
  }
}

setwd("/DATA/31/molteni/vexas_bm/results/integration/DE_markers/Thr_CTRL")

df <- list.files(path="/DATA/31/molteni/vexas_bm/results/integration/DE_markers/Thr_CTRL") %>% 
  lapply(read.csv) %>% 
  bind_rows

df <- df[df$p_val_adj<0.05,]
setwd("/DATA/31/molteni/vexas_bm/results/integration/DE_markers_OSR_Wu/summary")
write.csv(df, file="THRvsCTRL_sigMarkers.csv")

#COMPARISON VAL VEXAS PATIENTS vs HEALTHY CONTROLS WU

#select only the subjects involved in the comparison
Idents(BM.combo) <- "orig.ident"
BM.combo_subset <- subset(BM.combo, idents = c("BM-02","BM-03","BM-08","HD1", "HD2", "HD3", "HD4"))
clusters <- unique(BM.combo_subset$celltypes_annot)
Idents(BM.combo_subset) <- "celltypes_annot"
setwd("/DATA/31/molteni/vexas_bm/results/integration/DE_markers_OSR_Wu/Val_CTRL/")

for(i in clusters) {
  cells_group1 = sum(BM.combo_subset$mutations == "VAL" & BM.combo_subset$celltypes_annot == i)
  cells_group2 = sum(BM.combo_subset$mutations == "CTRL" & BM.combo_subset$celltypes_annot == i)
  if(cells_group1 >= 3 & cells_group2 >= 3 ) {
    markers <- FindMarkers(BM.combo_subset, ident.1 = "VAL", group.by = "mutations", subset.ident = i)
    markers['cluster'] <- i
    filename=paste("markers_cluster_",i,".csv", sep="")
    write.csv(markers, file=filename)
  }
  else {
    print(paste("Skipping comparison for", i, "due to insufficient cell count in one or both groups."))
  }
}

setwd("/DATA/31/molteni/vexas_bm/results/integration/DE_markers_OSR_Wu/Val_CTRL")

df <- list.files(path="/DATA/31/molteni/vexas_bm/results/integration/DE_markers_OSR_Wu/Val_CTRL") %>% 
  lapply(read.csv) %>% 
  bind_rows

df <- df[df$p_val_adj<0.05,]
setwd("/DATA/31/molteni/vexas_bm/results/integration/DE_markers_OSR_Wu/summary")
write.csv(df, file="VALvsCTRL_sigMarkers.csv")

#COMPARISON LEU VEXAS PATIENTS vs HEALTHY CONTROLS WU

#select only the subjects involved in the comparison
Idents(BM.combo) <- "orig.ident"
BM.combo_subset <- subset(BM.combo, idents = c("BM-04", "BM-09","HD1", "HD2", "HD3", "HD4"))
clusters <- unique(BM.combo_subset$celltypes_annot)
Idents(BM.combo_subset) <- "celltypes_annot"

setwd("/DATA/31/molteni/vexas_bm/results/integration/DE_markers_OSR_Wu/Leu_CTRL/")

for(i in clusters) {
  cells_group1 = sum(BM.combo_subset$mutations == "LEU" & BM.combo_subset$celltypes_annot == i)
  cells_group2 = sum(BM.combo_subset$mutations == "CTRL" & BM.combo_subset$celltypes_annot == i)
  if(cells_group1 >= 3 & cells_group2 >= 3 ) {
    markers <- FindMarkers(BM.combo_subset, ident.1 = "LEU", group.by = "mutations", subset.ident = i)
    markers['cluster'] <- i
    filename=paste("markers_cluster_",i,".csv", sep="")
    write.csv(markers, file=filename)
  }
  else {
    print(paste("Skipping comparison for", i, "due to insufficient cell count in one or both groups."))
  }
}

setwd("/DATA/31/molteni/vexas_bm/results/integration/DE_markers_OSR_Wu/Leu_CTRL")


df <- list.files(path="/DATA/31/molteni/vexas_bm/results/integration/DE_markers_OSR_Wu/Leu_CTRL") %>% 
  lapply(read.csv) %>% 
  bind_rows

df <- df[df$p_val_adj<0.05,]
setwd("/DATA/31/molteni/vexas_bm/results/integration/DE_markers_OSR_Wu/summary")
write.csv(df, file="LEUvsCTRL_sigMarkers.csv")


#COMPARISON VEXAS PATIENTS vs HEALTHY CONTROLS WU

#select only the subjects involved in the comparison
Idents(BM.combo) <- "orig.ident"
BM.combo_subset <- subset(BM.combo, idents = c("BM-01", "BM-02", "BM-03", "BM-04", "BM-08","BM-09","HD1", "HD2", "HD3", "HD4"))
                          
clusters <- unique(BM.combo_subset$celltypes_annot)
Idents(BM.combo_subset) <- "celltypes_annot"
setwd("/DATA/31/molteni/vexas_bm/results/integration/DE_markers_OSR_Wu/Pz_CTRL/")

for(i in clusters) {
  cells_group1 = sum(BM.combo_subset$condition == "PZ" & BM.combo_subset$celltypes_annot == i)
  cells_group2 = sum(BM.combo_subset$condition == "HD" & BM.combo_subset$celltypes_annot == i)
  if(cells_group1 >= 3 & cells_group2 >= 3 ) {
    markers <- FindMarkers(BM.combo_subset, ident.1 = "PZ", group.by = "condition", subset.ident = i)
    markers['cluster'] <- i
    filename=paste("markers_cluster_",i,".csv", sep="")
    write.csv(markers, file=filename)
  }
  else {
    print(paste("Skipping comparison for", i, "due to insufficient cell count in one or both groups."))
  }
}

setwd("/DATA/31/molteni/vexas_bm/results/integration/DE_markers_OSR_Wu/Pz_CTRL")
df <- list.files(path="/DATA/31/molteni/vexas_bm/results/integration/DE_markers_OSR_Wu/Pz_CTRL") %>% 
  lapply(read.csv) %>% 
  bind_rows

df <- df[df$p_val_adj<0.05,]

setwd("/DATA/31/molteni/vexas_bm/results/integration/DE_markers_OSR_Wu/summary")

write.csv(df, file="PZvsCTRL_sigMarkers.csv")