setwd("~/UiO Dropbox/Ankush Sharma/LymphomaBiology/2024/MyklebustLab/FL_Multiomic_Revision/MichaelGreenPaper")

library(Seurat)

library(dplyr)

library(ggplot2)

library(EnhancedVolcano)

currentjob="Green_FL"

#Single Cell RNASeq Data from Han et al, Cancer Discovery article, GEO ID GSE203610_FL-gex-meta  

GreenFL1 <- readRDS("../Data/FL_MultiomicRevision/Data_MGreenPaper/Green_FL_CZIbioscience.rds")

#Extract gene names

gene_names <- GreenFL1@assays[["RNA"]]@meta.features[["feature_name"]]

head(feature_names)

# Ensure that the length of gene_names matches the number of rows in your counts matrix

if(length(gene_names) == nrow(GreenFL1@assays[["RNA"]]@counts)) {

  # Step 2: Replace Ensembl IDs with gene names in the counts matrix

  rownames(GreenFL1@assays[["RNA"]]@counts) <- gene_names

} else {

  stop("The number of gene names does not match the number of Ensembl IDs in counts.")

}

rownames(GreenFL1@assays[["RNA"]]@data) <- gene_names

#Define the group based on sample_id

#Corresponding Sample name "FL-013", "FL-044", "FL-030", "FL-039") mapped for samples (FL03,FL11,FL06,FL10) that are EZH2 Mutated in Metadata for HAN et al  manuscript GSE203610_FL-gex-meta 

# Subset to include only 'Malignant' cells based on the 'cell_type' metadata

GreenFL1_Malignant <- subset(GreenFL1, subset = cell_type == "malignant cell")

fl_group_ids <- c("FL-013", "FL-044", "FL-030", "FL-039")

GreenFL1_Malignant$Group <- ifelse(GreenFL1_Malignant$sample_id %in% fl_group_ids, "EZH2Mutated", "EZH2WT")

table(GreenFL1_Malignant$Group)

GreenFL1_Malignant <- SetIdent(GreenFL1_Malignant, value = "Group")

# Perform FindMarkers between activated Tregs and naive Tregs

de_results <- FindMarkers(

  GreenFL1_Malignant,

  ident.1 = "EZH2Mutated", ident.2 = "EZH2WT",

  min.pct = 0.1

)

write.table(de_results,"GreenPaper_MalignantCells_FL_EZH2Mutated_WT_updated20240508.csv")

Histone <- c("H3C8","H2AC17", "H2AC7", "H2AC8", "H2BC10", "H2BC5", "H3C4", "H2AC19", "H2BC4", "H3C1", "H2BC14", "H4C4", "H2BC9", "H4C13", "H4C16", "H1-2", "H1-3", "H2AC4", "H2BC13", "H2AC16", "H2BC11", "H2AC6", "H4C12", "H4C8", "H4C2", "H2AC18", "H3C13", "H2AC20", "H2AC15", "H1-5", "H2AC12", "H2BC8", "H4C1", "H2BC3", "H2BC17", "H2AC13", "H2AC11", "H4C3", "H3C12", "H2BC7", "H2BC12", "H4C6", "H2BC18", "H4C14", "H2AC21", "H2BC15", "HLA-DQB1", "HLA-DQB2", "HLA-DRA", "TCL1A", "HLA-DRB1", "HLA-DOB", "HLA-B", "HLA-F", "LGALS1", "FOS", "JUN", "HLA_DQA", "SOX5", "CLDN5", "H2AC25", "SYT11", "CLEC4A", "IL2RA", "SLC25A23", "HLA-DQA2", "ALKAL2", "LMO7", "HLA-DMB", "HLA-E","HLA-DMB")

  # Assigning basic colors based on significance and fold change

  keyvals <- ifelse(de_results$p_val_adj > 0.000005, "#434c54", 

                    ifelse(de_results$avg_log2FC < -0.2 , '#5870B3',

                           ifelse(de_results$avg_log2FC > 0.2, '#C86065', 'lightgray')))

  keyvals[is.na(keyvals)] <- '#292524' # Assigning 'darkgrey' to NA values

  # Updated gene names based on color coding

  names(keyvals)[keyvals == '#5870B3' & de_results$avg_log2FC < -0.2] <- 'downregulated'

  names(keyvals)[keyvals == '#C86065' & de_results$avg_log2FC > 0.2] <- 'upregulated'

  names(keyvals)[keyvals == '#434c54'] <- 'non significant'

  names(keyvals)[keyvals == 'lightgray'] <- 'significant'

  Green_Volc <- EnhancedVolcano(de_results,

                                lab = rownames(de_results),

                                x = 'avg_log2FC',

                                y = 'p_val_adj',

                                selectLab = Histone,

                                xlim = c(-4,2),

                                xlab = bquote(~Log[2]~ 'fold change'),

                                title = 'EZH2 Mutated vs WT',

                                subtitle = 'Han et al.(2022) Cancer Discovery',

                                pCutoff = 0.0005,

                                FCcutoff = 0.2,

                                pointSize = 0.5,

                                labSize = 1.5,

                                axisLabSize=6,

                                titleLabSize = 8,  # Smaller font size for the title

                                subtitleLabSize = 6,

                                colCustom = keyvals,

                                colAlpha = 1,

                                legendLabSize =6,

                                legendIconSize = 3.0,

                                legendPosition = 'right',

                                legendLabels=c('Not significant','Avg log2FC','Adj.Pval','Adj.Pval & Avg. log2FC'),

                                drawConnectors = TRUE,

                                widthConnectors = 0.1,

                                colConnectors = 'grey10',

                                gridlines.major = FALSE,

                                gridlines.minor = FALSE,

                                border = 'full',

                                borderWidth = 0.3,

                                borderColour = 'black')

Green_Volc

# Close the PDF device to finalize and save the file

dev.copy(pdf,paste0(currentjob, "_Volcano_EZH2Mutated_vsWT_MalignantCells_default_20240508_v1.pdf"), width=6, height=7,paper='special')

dev.off()

dev.set(dev.next())

saveRDS(GreenFL1,"GREENFL1_group_AS.rds") 

#########################