---
title: "Relapse_patient_RNA_Seq"
author: "Christina Putnam"
date: "2024-10-02"
output: html_document
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

#load librarys
```{r}
library(GEOquery)
library(oligo)
library(stringr)
library(tidyverse)
library(clusterProfiler)
library(pheatmap)
#library(ReactomePA)
library(knitr)
library(fgsea)
library(HYPEtools)
library(ggpubr)
library(hgu133plus2.db)
library(affycoretools)
library(viridis)
library("ArrayExpress")
library(edgeR)
library(tidyverse)
library(PCAtools)
library(GSVA)
#library("arrayQualityMetrics")
library("readxl")
library(jetset)
library("affyio")
library(dplyr)
library(oligoClasses)
```

Project description
RNA sequencing data from patients pre and mid treatment (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE240287). Will be analysed to score gene signatures in pre and mid treatment samples to determine if resistance scores are higher in mid treatment samples. 

#get data
```{r}

#untar series file
fileList <- untar("C:/Users/Christina/Documents/Relapse RNA seq/GSE240287_RAW.tar",exdir = "C:/Users/Christina/Documents/Relapse RNA seq/New dir", list = T)

#read in each file and extract info needed to create count matrix
for (i in fileList){
  
 #read in file
  raw_count_file <- read.delim(file = paste0("C:/Users/Christina/Documents/Relapse RNA seq/New dir/", i)) 
  
  #extract geneSymbol and count columns
  raw_count_file <- raw_count_file[,c("geneID","count")]
  mode(raw_count_file$count) <- "integer"
  colnames(raw_count_file) <- c("GeneID", str_sub(i, 1, 10)) 

  
  write_delim(raw_count_file, file = paste0("C:/Users/Christina/Documents/Relapse RNA seq/Count files/", i), delim = "\t")
}

path = "C:/Users/Christina/Documents/Relapse RNA seq/Count files/"
setwd(path)
#create count matrix by joining RawReadCountsOnly.out files
counts <- list.files(pattern = "genes.counts.txt.gz", path = path) %>% map(read.table, header=TRUE) %>%
 purrr::reduce(inner_join, by = "GeneID")

#set the first column as row names
rawData <- counts[,-1]
rownames(rawData) <- counts[,1]
dim(rawData)

write.csv(rawData, file = "C:/Users/Christina/Documents/Relapse RNA seq/Count files/EM.csv")
#write.csv(counts.keep, "~/Neuroblastoma RNAseq analysis/filtered_EM.csv")


#clinical data
data <- getGEO(GEO = "GSE240287")
clindata <- data[["GSE240287_series_matrix.txt.gz"]]@phenoData@data
clindata <- clindata[,c("pax-foxo1 fusion_status:ch1","title")]
clindata <- clindata %>%
      tidyr::extract(title, c('Patient', 'Treatment'), '(\\w+\\s\\w+)\\s(.*)')

```

#design matrix
```{r}

Subject <- gsub(" ", "_", clindata$Patient) %>% factor()
Treat <- gsub("-", "_", clindata$Treatment)%>% factor(levels=c("pre_treatment","mid_treatment"))

 design <- model.matrix(~0+Treat+Subject)
 #design <- model.matrix(~Subject+Treat)
                      
```


#DGEList
```{r}

#get ensembl ids as df to add to DGElist$gene 
gene_annotation <- as.data.frame(rownames(rawData))
colnames(gene_annotation) <- c("ensembl_id")

#get DGElist
edgeRData <- DGEList(counts=rawData, genes=gene_annotation)
dim(edgeRData)

#remove transcript number from the DGEList$gene (as it can't map with the transcript number)
edgeRData$genes$ensembl_id = sapply(strsplit(edgeRData$genes$ensembl_id, ".", fixed=T), function(x) x[1]) #remove anythign after . in rowname

#get gene symbols for ensembl ids
genes <-  mapIds(
    # Replace with annotation package for the organism relevant to your data
    org.Hs.eg.db,
    keys = edgeRData$genes$ensembl_id,
    # Replace with the type of gene identifiers in your data
    keytype = "ENSEMBL",
    # Replace with the type of gene identifiers you would like to map to
    column = "SYMBOL",
    # This will keep only the first mapped value for each Ensembl ID
    multiVals = "first"
  ) %>% as.data.frame()

#add gene symbol to DGEList$genes
edgeRData$genes$symbol <- genes$.

#filter using edger with design matrix
keep <- filterByExpr(edgeRData, design)
edgeRData <- edgeRData[keep,keep.lib.sizes=FALSE]
dim(edgeRData)


```

#Boxplot of logCPM before normalisation
```{r}

#Sample outliers can be detected based on the distribution of read counts between samples which are visualised using a box and whisker plot of log2 counts per milliom (CPM) for each sample. CPM is used to correct for differences in library size (ie depth of sequencing per sample) which are logged because the count data are not normally distributed. We will use the cpm function in edgeR to get log2(CPM) which also adds a small offset to avoid taking a log of zero.

# Get log2 counts per million
par(mar=c(7,7,4,1)+.1)
logcounts <- cpm(edgeRData,log=TRUE)
# Check distributions of samples using boxplots
boxplot(logcounts, xlab="", ylab="Log2 counts per million",las=2)
# Let's add a blue horizontal line that corresponds to the median logCPM
abline(h=median(logcounts),col="blue")
title(xlab="Sample", line=6, cex.lab=1)
title("Boxplots of log-2 CPMs before normalisation")
```


#Normalisation
```{r}

edgeRData <- calcNormFactors(edgeRData, method="TMM")

#get CPM
logCPM <- cpm(edgeRData, log = TRUE) %>% as.data.frame()

```


#boxplot of counts
```{r}

#make a bar plot of the number of reads, i.e. library size, in each sample
png(filename="lib_size.png", width=900, bg="white")
par(mar=c(7,7,4,1)+.1)
barplot(edgeRData$samples$lib.size,names=rownames(edgeRData$samples),las=3, ylim=range(pretty(c(edgeRData$samples$lib.size))), cex.names = 0.8, ylab="Library size", cex.axis = 0.8);
title(xlab="Sample", line=5.5, cex.lab=1.1)
title("Barplot of library sizes");
dev.off()

#after normalisation
png(filename="logCPM_after_norm.png", width=900, bg="white")
par(mar=c(7,7,4,1)+.1)
for_median <- cpm(edgeRData, log=TRUE) #creating seperate df to caluclate median as after I reorder cannot calc median
boxplot(logCPM, xlab="", ylab="Log2 counts per million", las=2, cex.names = 0.8, cex.axis = 0.8)
abline(h=median(for_median), col="blue")
title("Boxplots of log-2 CPMs after normalisation")
title(xlab="Sample", line=6, cex.lab=1.1, family="Calibri Light")
dev.off()

par(mar=c(7,7,4,1)+.1)
logcounts <- cpm(edgeRData,log=TRUE)
# Check distributions of samples using boxplots
boxplot(logcounts, xlab="", ylab="Log2 counts per million",las=2)
# Let's add a blue horizontal line that corresponds to the median logCPM
abline(h=median(logcounts),col="blue")
title(xlab="Sample", line=6, cex.lab=1)
title("Boxplots of log-2 CPMs before normalisation")

```


#PCA
```{r}
#get normalised counts and remove dups/NAs
expression_data <- logCPM
p <- pca(expression_data, metadata = clindata, removeVar = 0.1)


###determoine how many pc to retain----------------------------------------------------
    horn <- parallelPCA(expression_data)
  horn$n 
  elbow <- findElbowPoint(p$variance)
  elbow

  screeplot(p,
    components = getComponents(p, 1:23),
    vline = c(horn$n, elbow)) +

    geom_label(aes(x = horn$n + 0, y = 50,
      label = 'Horn\'s', vjust = -2, size = 8)) +
    geom_label(aes(x = elbow + 1, y = 50,
      label = 'Elbow method', vjust = -1, size = 8))
  

 #biplots ----------------------------------------------- 
    
 
   variables <- c("Patient", "pax-foxo1 fusion_status:ch1", "Treatment")
  for (i in variables) {   
    if(i == "Patient"){legend_title <- "Patient" }
    if(i == "pax-foxo1 fusion_status:ch1"){legend_title <- "Fusion gene status" }
     if(i == "Treatment"){legend_title <- "Treatment stage" }
    
    #find length of variables
  n = length(unique(p$metadata[[i]]))
 p$metadata[[i]] <-  as.factor(p$metadata[[i]])
#set colours
  myCols <- setNames(palette.colors(palette = "Okabe-Ito"), unique(p$metadata[[i]]))
  #myCols <- viridis(n, option = "G", direction = -1, end = 0.9)
 # legend_title <- sapply(strsplit(i,":"), `[`, 1)
  #make filename title
   j <- sapply(strsplit(i,":"), `[`, 1)
   filename <- paste0("PCA_biplot",j,".png")
   #label <- subset(p$metadata, `histology:ch1` == "Mixed Alv/Emb") %>% rownames() #label only mixed samples
 g <- biplot(p,
         lab = NULL,
         #max.overlaps = 200,
    colby = paste0(i),
    colkey = myCols,
    hline = 0, vline = 0,
    legendPosition = 'right', 
    ellipse = FALSE,
      ellipseType = 't',
      ellipseLevel = 0.95,
      ellipseFill = TRUE,
      ellipseAlpha = 1/4,
    colLegendTitle = legend_title,
      ellipseLineSize = 1.0,
   shape = "Treatment",
 #  shapeLegendTitle = "Histology",
    pointSize = 5,
    title = 'PCA bi-plot',
    subtitle = 'PC1 versus PC2',)
 print(g)
ggsave(g, file = filename, width = 6, height = 4, dpi = 300)
  }
  

  
```

#GSVA plot function
```{r}

# This function plots GSVA results with statistics
plot_resistance_score_with_statistics <- function(genes, subtype, clinical_variable, resistance_type, test_variable, histology) {

#hub_genes <- hub_genes$Gene
genes <- data.frame(
  "entrez_id" = mapIds(
    # Replace with annotation package for the organism relevant to your data
    org.Hs.eg.db,
    keys = genes,
    # Replace with the type of gene identifiers in your data
    keytype = "SYMBOL",
    # Replace with the type of gene identifiers you would like to map to
    column = "ENSEMBL",
    # This will keep only the first mapped value for each Ensembl ID
    multiVals = "first"
  )
) 

genes <-list(genes = genes$entrez_id)

#remove transcript number from the rowname
logCPM$geneID = sapply(strsplit(rownames(logCPM), ".", fixed=T), function(x) x[1]) #remove anythign after . in rowname
logCPM <- logCPM[!duplicated(logCPM$geneID), ] #remove duplicate transcripts
#set the 19th column as row names
rownames(logCPM) <- logCPM[,ncol(logCPM)]
logCPM  = subset(logCPM, select = -c(geneID) ) 
dim(logCPM)

#change to matrix
filtered_mapped_matrix <- logCPM %>%
  as.matrix()

gsvamethod <- gsvaParam(filtered_mapped_matrix,
  genes, minSize = 1)
gsva_results <- gsva(gsvamethod)

#order samples by score
gsva_results <- gsva_results[,order(gsva_results[1,], decreasing = T)]
gsva_results <- as.data.frame(gsva_results)
#View(gsva_results)

#merge gsva results with pheno data
gsva_results <- merge(gsva_results, clindata, by = "row.names")
colnames(gsva_results)[1] <- "sample"


ARMS_samples <- subset(gsva_results, `pax-foxo1 fusion_status:ch1` == "positive") %>% rownames()
ERMS_samples <- subset(gsva_results, `pax-foxo1 fusion_status:ch1` == "negative") %>% rownames()

if(subtype=="FP"){
  #keep FP samples
gsva_results <- gsva_results[rownames(gsva_results) %in% ARMS_samples, ]
}
if(subtype=="FN"){
  #keep FN samples only
  gsva_results <- gsva_results[ rownames(gsva_results) %in% ERMS_samples, ]
}


  #first performs normality test to decide if parametric or non parametric pairwise tests should be used
  
  shapiro_results <- shapiro.test(gsva_results$gsva_results)
  
  if (shapiro_results$p.value > 0.05) {
    print("Parametric")
    correct_method <- c("t.test")
    #correct_label <- "T-test, p = {p}" #add this if you want to label test used
  } else {
    print("Non parametric")
   correct_method <- c("wilcox.test")
 #  correct_label <- "Wilcoxon test, p = {p}"
  }

  #select correct colour for histology
     if(all(str_detect(gsva_results$`pax-foxo1 fusion_status:ch1`, 'negative|positive'))){histology_colour <- c("#0072B2", "#CC79A7")}
   if(all(!str_detect(gsva_results$`pax-foxo1 fusion_status:ch1`, "negative"))){histology_colour <- c("#CC79A7")}
      if(all(!str_detect(gsva_results$`pax-foxo1 fusion_status:ch1`, "positive"))){histology_colour <-  c("#0072B2")}

#then selects options for clinical variable picked
  print(clinical_variable)
  
    if (clinical_variable == "Treatment") {
        variable_levels <- c("pre-treatment", "mid-treatment")
    readable_clinical_variable <- c("treatment stage")
       my_comparisons <- list(c("pre-treatment", "mid-treatment"))
    }
  
#reorder plot
gsva_results[[clinical_variable]] <- factor(gsva_results[[clinical_variable]], variable_levels)
 fct_relevel(gsva_results[[clinical_variable]])

okabe_palette <- palette.colors(palette = "Okabe-Ito")
patient_colors <- setNames(okabe_palette[1:length(unique(gsva_results$Patient))],
                           unique(gsva_results$Patient))
    #boxplot
    p <- ggboxplot(gsva_results  , x =  paste0(clinical_variable), y = "gsva_results", fill = paste0(clinical_variable), outlier.shape = NA)+labs(title= paste0(subtype, "-", resistance_type," score and\n", readable_clinical_variable),
        x =paste0(toupper(substr(readable_clinical_variable, 1, 1)), substr(readable_clinical_variable, 2, nchar(readable_clinical_variable)), sep=""), y = "GSVA score", colour = "Patient")+ theme( axis.text.x = element_text(angle = 90,  size = 12))+ scale_fill_brewer(palette="BuPu", guide = "none")+geom_jitter(aes(col=Patient))+ theme(legend.position = "right")+ stat_compare_means(comparisons = my_comparisons,label = "p.signif", method = paste0(correct_method),paired = TRUE)+
  scale_color_manual(values = patient_colors)
#stat_pwc(method = paste0(correct_method), label = "p.adj.signif", p.adjust.method = "BH",na.rm = TRUE,) 
    print(p)
ggsave(file= paste0(resistance_type, "_", clinical_variable, "_", subtype,".png"), width = 3.5, height = 4)
}
    

# p <- ggboxplot(gsva_results  , x =  paste0(clinical_variable), y = "gsva_results", fill = paste0(clinical_variable), outlier.shape = NA)+labs(title= paste0(resistance_type, " resistance score and ", readable_clinical_variable),
#         x =paste0(toupper(substr(readable_clinical_variable, 1, 1)), substr(readable_clinical_variable, 2, nchar(readable_clinical_variable)), sep=""), y = "GSVA score", colour = "Fusion gene status")+ theme( axis.text.x = element_text(angle = 90,  size = 12))+ scale_fill_brewer(palette="BuPu", guide = "none")+geom_jitter(aes(col=`pax-foxo1 fusion_status:ch1`))+
#   scale_color_manual(values = histology_colour)+ theme(legend.position = "right")+ stat_compare_means(comparisons = my_comparisons,label = "p.adjust", method = paste0(correct_method),paired = TRUE)
# #stat_pwc(method = paste0(correct_method), label = "p.adj.signif", p.adjust.method = "BH",na.rm = TRUE,)
    
```

#GSVA score plots of signature
```{r}

#for intrinsic all samples -----------------------------------------------------------
setwd("C:/Users/Christina/Documents/Relapse RNA seq/GSVA scores with jitter as patient")

#get genes 
genes <- read_xlsx("C:/Users/Christina/Documents/RMS RNA seq/Zoes data/CLONES ONLY ANALYSIS/intrinsic_resistance_signature.xlsx") #intrinsic sig
genes <- genes$Gene

  # plot_resistance_score_with_statistics(genes = genes, subtype = "both",  clinical_variable = "Treatment", resistance_type = "Intrinsic")
  
  #for intrinsic ARMS only -----------------------------------------------------------------
    plot_resistance_score_with_statistics(genes = genes, subtype = "FP",  clinical_variable = "Treatment", resistance_type = "intrinsic", histology = "ARMS")
    
    #acquired RH4-------------------------------------------------------------------------------------
 genes <- read_xlsx("C:/Users/Christina/Documents/RMS RNA seq/Luciles data/RH4 WGCNA/RH4_acquired_res_signature.xlsx") #acquired rh4 sig
 genes <- genes$Gene
 
  plot_resistance_score_with_statistics(genes = genes, subtype = "FP",  clinical_variable = "Treatment", resistance_type = "acquired", histology = "ARMS")
 
 #acquired RMSYM---------------------------------------------------------
 genes <- read_xlsx("C:/Users/Christina/Documents/RMS RNA seq/Luciles data/RMSYM WGCNA/RMSYM_acquired_resistance_gene_signature.xlsx") #RMSYM acquired sig
  genes <- genes$Gene
  
  plot_resistance_score_with_statistics(genes = genes, subtype = "FN",  clinical_variable = "Treatment", resistance_type = "acquired", histology = "ERMS")
```

#GSVA score plots of DEGs
Score samples using DEGs as resistance signature
```{r}

#using GSVA to calculate a chemo-resistant score
#intrinsic all samples -----------------------------------------------------------
setwd("C:/Users/Christina/Documents/Relapse RNA seq/DEGs as signature")
#get normalised counts

#using overlapping DEGs upregulated in C7 and C13 
DEGS_C7 <- read.csv("C:/Users/Christina/Documents/RMS RNA seq/Zoes data/DEG lists/Significant DEGs (LFC 1)/DEGs_LFC_1_RH30_P_vs_RH_C7.csv") 
DEGS_C7_UP <- DEGS_C7[DEGS_C7$Expression == 'Up-regulated',]

DEGS_C13 <- read.csv("C:/Users/Christina/Documents/RMS RNA seq/Zoes data/DEG lists/Significant DEGs (LFC 1)/DEGs_LFC_1_RH30_P_vs_RH_C13.csv")
DEGS_C13_UP <- DEGS_C13[DEGS_C13$Expression == 'Up-regulated',]

genes <- intersect(DEGS_C13_UP$X,DEGS_C7_UP$X)

write.csv(genes, "intrinsic_overlapping_DEGs.csv")
  # plot_resistance_score_with_statistics(genes = genes, subtype = "both",  clinical_variable = "Treatment", resistance_type = "Intrinsic" )

#for intrinsic ARMS only --------------------------------------------------------------------
  genes <- intersect(DEGS_C13_UP$X,DEGS_C7_UP$X)
  plot_resistance_score_with_statistics(genes = genes, subtype = "FP",  clinical_variable = "Treatment", resistance_type = "intrinsic", histology = "ARMS")

#acquired RH4-----------------------------------------------------------------
#using overlapping DEGs upregulated in RH4 IFO and VCR resistant
DEGS_RH4_IFO <- read.csv("C:/Users/Christina/Documents/RMS RNA seq/Luciles data/DEG lists/DEGs_LFC_1_RH4ctrl_RH4ifo.csv") 
DEGS_RH4_IFO_UP <- DEGS_RH4_IFO[DEGS_RH4_IFO$Expression == 'Up-regulated',]

DEGS_RH4_VCR <- read.csv("C:/Users/Christina/Documents/RMS RNA seq/Luciles data/DEG lists/DEGs_LFC_1_RH4ctrl_RH4vcr.csv") 
DEGS_RH4_VCR_UP <- DEGS_RH4_VCR[DEGS_RH4_VCR$Expression == 'Up-regulated',]

genes <- intersect(DEGS_RH4_IFO_UP$X,DEGS_RH4_VCR_UP$X)
write.csv(genes, "ARMS_acquired_overlapping_DEGs.csv")

  plot_resistance_score_with_statistics(genes = genes, subtype = "FP",  clinical_variable = "Treatment", resistance_type = "acquired", histology = "ARMS")

#acquired RMSYM -----------------------------------------------------------------
#using overlapping DEGs upregulated in RH4 IFO and VCR resistant
DEGS_RMSYM_IFO <- read.csv("C:/Users/Christina/Documents/RMS RNA seq/Luciles data/DEG lists/DEGs_LFC_1_RMSYMctrl_RMSYMifo.csv") 
DEGS_RMSYM_IFO_UP <- DEGS_RMSYM_IFO[DEGS_RMSYM_IFO$Expression == 'Up-regulated',]

DEGS_RMSYM_VCR <- read.csv("C:/Users/Christina/Documents/RMS RNA seq/Luciles data/DEG lists/DEGs_LFC_1_RMSYMctrl_RMSYMvcr.csv") 
DEGS_RMSYM_VCR_UP <- DEGS_RMSYM_VCR[DEGS_RMSYM_VCR$Expression == 'Up-regulated',]

genes <- intersect(DEGS_RMSYM_IFO_UP$X,DEGS_RMSYM_VCR_UP$X)
write.csv(genes, "ERMS_acquired_overlapping_DEGs.csv")
  plot_resistance_score_with_statistics(genes = genes, subtype = "FN",  clinical_variable = "Treatment", resistance_type = "acquired", histology = "ERMS")
```

#BCV
```{r}

  edgeRData <- estimateDisp(edgeRData, design, robust=TRUE)
  edgeRData$common.dispersion
plotBCV(edgeRData)

```

#DGE
```{r}

con <- makeContrasts(Treatmid_treatment-Treatpre_treatment, levels=design)
fit <- glmQLFit(edgeRData, design, robust=TRUE)
qlf <- glmQLFTest(fit, contrast=con)
 topTags(qlf)
 
 #get DEGs
qlf$table$p.adjust <-p.adjust(qlf$table$PValue,method="BH")
#DEGs_object <- qlf[qlf$table$p.adjust < 0.05,]

 
```


#Vol plot
```{r}
setwd("~/Relapse RNA seq/DEGs")
  data <- qlf$table
data$gene <- qlf$genes$symbol
  head(data)

  ############## Volcano Plots ################
 
  # Make a basic volcano plot for up and down genes
  
  # p1 <- ggplot(data, aes(logFC, -log(FDR,10))) +
  #   geom_point(size = 2/5) +
  #   xlab(expression("log"[2]*"FC")) +
  #   ylab(expression("-log"[10]*"FDR"))

  data <- data %>% 
  mutate(
    Expression = case_when(logFC > 1.5 & p.adjust < 0.05 ~ "Up-regulated",
                           logFC < -1.5 & p.adjust < 0.05 ~ "Down-regulated",
                           TRUE ~ "Unchanged")
    )
head(data)  %>% 
   kable()

write.csv(data, "unfiltered_DEGs_mid_vs_pre.csv")

 # plot volcano plot with up and downregulated genes coloured
  p2 <- ggplot(data, aes(logFC, -log(p.adjust,10))) +
    geom_point(aes(color = Expression), size = 2/5) +
    xlab(expression("log"[2]*"FC")) +
    ylab(expression("-log"[10]*"p.adjust")) +
    scale_color_manual(values = c("dodgerblue3", "gray50", "firebrick3")) +
    guides(colour = guide_legend(override.aes = list(size=1.5)))
   print(p2)
  
   
   

   ############# vol plot for highlighting specific genes  in gene_list ############
  #add expression
   vol_plot_data <- data %>% 
  mutate(
    Expression = case_when(logFC > 1.5 & p.adjust < 0.05 ~ "Up-regulated",
                           logFC < -1.5 & p.adjust < 0.05 ~ "Down-regulated",
                           TRUE ~ "Unchanged")
    )
head(data) %>% 
 kable()

dim(data[data$Expression == "Up-regulated",])
dim(data[data$Expression == "Down-regulated",])

  # Add colour, size and alpha (transparency) to volcano plot --------------------
cols <- c("Up-regulated" = "#ffad73", "Down-regulated" = "#26b3ff", "Unchanged" = "grey") 
sizes <- c("Up-regulated" = 2, "Down-regulated" = 2, "Unchanged" = 1) 
alphas <- c("Up-regulated" = 1, "Down-regulated" = 1, "Unchanged" = 0.5)

 top <- 10
  top_genes <- bind_rows(
    data %>%
      filter(Expression == 'Up-regulated') %>%
      dplyr::arrange(p.adjust, desc(abs(logFC))) %>%
      head(top),
    data %>%
     filter(Expression == 'Down-regulated') %>%
      dplyr::arrange(p.adjust, desc(abs(logFC))) %>%
      head(top)

  )
  top_genes %>%
    kable()
  
 #plot vol plot highlighting genelist and adding labels
 p3 <- ggplot(data = vol_plot_data, # Original data  
       aes(x = logFC, y = -log(p.adjust,10))) + 
  geom_hline(yintercept = -log10(0.05),
             linetype = "dashed") + 
  geom_vline(xintercept = c(-1.5, 1.5),
             linetype = "dashed") +  
 scale_colour_manual(values = cols)+
    ggtitle("Mid- vs pre-treatment") +
    geom_point(aes(colour = Expression), 
             alpha = 0.8, 
             shape = 16,
             size = 1)  
                  
 print(p3)
  p3 <- p3 +
    geom_label_repel(data = top_genes, label.size = 0.5, force = 1.5, max.overlaps = 25, 
                     mapping = aes(logFC, -log(p.adjust,10), label = gene),
                     size = 4.5)
  print(p3)
  
  
  
```

#GSEA using cluster profiler
```{r}
## feature 1: numeric vector
geneList = DEGs_object$table[,1]
## feature 2: named vector
names(geneList) = as.character(DEGs_object$genes$symbol)
## feature 3: decreasing orde
geneList = sort(geneList, decreasing = TRUE)


#for GO------------------------------------------------------
ego3 <- gseGO(geneList     = geneList,
              OrgDb        = org.Hs.eg.db,
              ont          = "ALL",
              minGSSize    = 10,
              maxGSSize    = 600,
              nPermSimple = 10000,
              pvalueCutoff = 0.05,
              verbose      = FALSE, 
               keyType = "SYMBOL", 
              eps = 0)
## Dotplot 
ego3 <- ego3 %>%
  dplyr::arrange(p.adjust)


GO <- dotplot(ego3, showCategory = 10, font.size=14, x = "Count",
              split = ".sign",)+facet_grid(~.sign)
print

ggsave(GO, file=paste0("GO_mid_vs_pre.png"), width = 20, height = 30, units = "cm")


#for reactome

y <- gsePathway(geneList, 
                pvalueCutoff = 0.05,
                pAdjustMethod = "BH", 
                organism = "human",
                verbose = FALSE)
head(y)


```

#MSigdb using cluster profiler
```{r}
setwd("~/Relapse RNA seq/GSEA cluster profiler")

cluster_prolifer_plots <- function(msigdb_category, msigdb_subcategory){
  
  if(length(msigdb_subcategory) <= 0){
t2g <- msigdbr(species = "Homo sapiens", category = paste0(msigdb_category)) %>% 
  dplyr::select(gs_name, gene_symbol)
head(t2g)}
  else if(length(msigdb_subcategory) > 0)
t2g <- msigdbr(species = "Homo sapiens", category = paste0(msigdb_category), subcategory = paste0(msigdb_subcategory)) %>% 
  dplyr::select(gs_name, gene_symbol)
head(t2g)

set.seed(42)
em2 <- GSEA(geneList, TERM2GENE = t2g)
head(em2)

## Dotplot 
em2 <- em2 %>%
  dplyr::arrange(p.adjust)

H <- dotplot(em2, font.size=14,label_format = 60, x = "Count",
              split = ".sign",)+facet_grid(~.sign)
print(H)

ggsave(H, file=paste0(msigdb_category, "_mid_vs_pre.png"), width = 20, height = 20, units = "cm")

#core enriched genes 
core_plot <- cnetplot(em2, foldChange=geneList, categorySize="count",
                                           showCategory = 5,
                                           vertex.label.font=5,
                                           colorEdge = TRUE)


ggsave(core_plot, file=paste0(msigdb_category, "_coreplot_mid_vs_pre.png"), width = 30, height = 30, units = "cm", bg = "white")



# #for just upregulated
# up <- em2$Description[order(em2$NES, decreasing=TRUE)][1:6]
# cnetplot(em2, foldChange=geneList, showCategory=up)
}

cluster_prolifer_plots(msigdb_category = "H", msigdb_subcategory = NULL)
cluster_prolifer_plots(msigdb_category = "C5", "GO")

```

#Scatterplot function
LFC of RMSYM VCR-res and mid treatment
```{r}


make_scatterplot_LFC <- function(filepath_1, filepath_2, treatment_condition_1, treatment_condition_2, LFC_TC1, LFC_TC2, colour_TC1, title_name){

#get TC1 DEGs
TC1_DEGs <- read.csv(paste0(filepath_1))
TC2_DEGs <- read.csv(paste0(filepath_2))

TC1_DEGs <- TC1_DEGs %>% mutate(
    Expression = case_when(logFC > as.numeric(paste0(LFC_TC1)) & FDR < 0.05 ~ "significant",
                           logFC < as.numeric(paste0(-LFC_TC1)) & FDR < 0.05 ~ "significant",
                           TRUE ~ "Unchanged") )

#get TC2 DEGs
TC2_DEGs <- TC2_DEGs %>% mutate(
    Expression = case_when(logFC > as.numeric(paste0(LFC_TC2)) & p.adjust < 0.05 ~ "significant",
                           logFC < as.numeric(paste0(-LFC_TC2)) & p.adjust < 0.05 ~ "significant",
                           TRUE ~ "Unchanged") )
 
df <- merge(TC1_DEGs, TC2_DEGs, by.x = 'X', by.y = "gene")

df <- df %>%
   mutate(sig = case_when(Expression.x == "significant" & Expression.y == "significant" ~ "Both",
                          Expression.x == "significant" ~ paste0(treatment_condition_1),
                           Expression.y == "significant" ~ paste0(treatment_condition_2), Expression.x == "Unchanged" & Expression.y == "Unchanged" ~ "Unchanged" ))

levels(df$sig) <- c("Both", treatment_condition_2, treatment_condition_1, "Unchanged")
myColors <- c("#009E73", "#000000", colour_TC1, "grey")
names(myColors) <- levels(df$sig)
colScale <- scale_colour_manual(name = "grp",values = myColors)

  
scatterplot <- ggplot(df, aes(x=logFC.x, y=logFC.y)) + geom_point(aes(color=sig))   +theme_bw()+ geom_smooth(method=lm , color="blue", fill="#69b3a2", se=TRUE) + scale_color_manual(name = "Sig", values=myColors) + ggtitle(title_name) +
  xlab(paste0("Log-2 fold change\n", treatment_condition_1)) + ylab(paste0("Log-2 fold change ", treatment_condition_2)) +theme(axis.text=element_text(size=12),
        axis.title=element_text(size=12)) + guides(color=guide_legend(title="Differentially expressed in:")) +
  stat_regline_equation(label.y = 7) +
  stat_cor(label.y = 6)
ggsave(scatterplot, filename = paste0("Scatterplot_", treatment_condition_1, "_", treatment_condition_2, ".png"),bg= "white", width = 5.5, height=4)
print(scatterplot)

}

```

#Make scatterplots
Using function
```{r}
setwd("~/Relapse RNA seq/scatterplots")

#RMSYM -------------------------------------
  make_scatterplot_LFC(filepath_1 = "~/RMS RNA seq/Luciles data/Unfiltered_DEG_lists/unfiltered_DEGs_RMSYMctrl_RMSYMvcr.csv", treatment_condition_1 = "RMSYM VCR-resistant vs control", LFC_TC1 = 1, colour_TC1 = "#56B4E9", filepath_2 = "~/Relapse RNA seq/DEGs/unfiltered_DEGs_mid_vs_pre.csv", treatment_condition_2 = "Mid- vs pre-treatment", LFC_TC2 = 1.5, title_name = "Correlation of LFC between RMSYM VCR-resistant and\nmid-treatment samples")
   
  make_scatterplot_LFC(filepath_1 = "~/RMS RNA seq/Luciles data/Unfiltered_DEG_lists/unfiltered_DEGs_RMSYMctrl_RMSYMifo.csv", treatment_condition_1 = "RMSYM IFO-resistant vs control", LFC_TC1 = 1, colour_TC1 = "#D55E00", filepath_2 = "~/Relapse RNA seq/DEGs/unfiltered_DEGs_mid_vs_pre.csv", treatment_condition_2 = "Mid- vs pre-treatment", LFC_TC2 = 1.5, title_name = "Correlation of LFC between RMSYM IFO-resistant and\nmid-treatment samples")
  
  #RH4 -------------------------------------
  make_scatterplot_LFC(filepath_1 = "~/RMS RNA seq/Luciles data/Unfiltered_DEG_lists/unfiltered_DEGs_RH4ctrl_RH4vcr.csv", treatment_condition_1 = "RH4 VCR-resistant vs control", LFC_TC1 = 1, colour_TC1 = "#CC79A7", filepath_2 = "~/Relapse RNA seq/DEGs/unfiltered_DEGs_mid_vs_pre.csv", treatment_condition_2 = "Mid- vs pre-treatment", LFC_TC2 = 1.5, title_name = "Correlation of LFC between RH4 VCR-resistant and\nmid-treatment samples")
   
  make_scatterplot_LFC(filepath_1 = "~/RMS RNA seq/Luciles data/Unfiltered_DEG_lists/unfiltered_DEGs_RH4ctrl_RH4ifo.csv", treatment_condition_1 = "RH4 IFO-resistant vs control", LFC_TC1 = 1, colour_TC1 = "#F0E442", filepath_2 = "~/Relapse RNA seq/DEGs/unfiltered_DEGs_mid_vs_pre.csv", treatment_condition_2 = "Mid- vs pre-treatment", LFC_TC2 = 1.5, title_name = "Correlation of LFC between RH4 IFO-resistant and\nmid-treatment samples")
  
    #Clones -------------------------------------
  make_scatterplot_LFC(filepath_1 = "~/RMS RNA seq/Zoes data/UNFILTERED DEGs/unfiltered_DEGs_RH30_P_vs_RH_C7.csv", treatment_condition_1 = "Clone 7 vs parent", LFC_TC1 = 1, colour_TC1 = "skyblue1", filepath_2 = "~/Relapse RNA seq/DEGs/unfiltered_DEGs_mid_vs_pre.csv", treatment_condition_2 = "Mid- vs pre-treatment", LFC_TC2 = 1.5, title_name = "Correlation of LFC between clone 7 and\nmid-treatment samples")
   
  make_scatterplot_LFC(filepath_1 =  "~/RMS RNA seq/Zoes data/UNFILTERED DEGs/unfiltered_DEGs_RH30_P_vs_RH_C13.csv", treatment_condition_1 = "Clone 13 vs parent", LFC_TC1 = 1, colour_TC1 = "#0072B2", filepath_2 = "~/Relapse RNA seq/DEGs/unfiltered_DEGs_mid_vs_pre.csv", treatment_condition_2 = "Mid- vs pre-treatment", LFC_TC2 = 1.5, title_name = "Correlation of LFC between clone 13 and\nmid-treatment samples")
  
     
  make_scatterplot_LFC(filepath_1 =  "~/RMS RNA seq/Zoes data/UNFILTERED DEGs/unfiltered_DEGs_RH30_P_vs_RH_HD.csv", treatment_condition_1 = "HD vs parent", LFC_TC1 = 1, colour_TC1 = "#E69F00", filepath_2 = "~/Relapse RNA seq/DEGs/unfiltered_DEGs_mid_vs_pre.csv", treatment_condition_2 = "Mid- vs pre-treatment", LFC_TC2 = 1.5, title_name = "Correlation of LFC between high dose and\nmid-treatment samples")

```




#REACTOME CAMERA
```{r}

setwd("~/Relapse RNA seq/GSEA")
# Check contrasts:
print(colnames(design))

#CAMERA takes as input the DGEList object dgeObj, the indexed list of gene sets c2.ind, the design matrix, the contrast being tested, as well as some other arguments. By default, CAMERA can estimate the correlation for each gene set separately. However, in practise, it works well to set a small inter-gene correlation of about 0.05 using the inter.gene.cor argument.
i <- c("Pre vs mid-treatment")

# Get pathways from gmt file
Reactome_pathways <- gmtPathways("C:/Users/Christina/Documents/GSEA .gmt/c2.cp.reactome.v2024.1.Hs.symbols.gmt")
#create index 
reactome.ind <- ids2indices( Reactome_pathways, edgeRData$genes$symbol)
# Run analysis:
gst.camera <- camera.DGEList(edgeRData, design= design, index=reactome.ind,contrast=con, inter.gene.cor=0.05)

#get number of signif gene sets at FDR 0.05
table(gst.camera$FDR < 0.05)

#find signif gene sets at FDR 0.05
sig_genesets <- gst.camera[gst.camera$FDR < 0.05, ]
print(sig_genesets)
 
#ensure cols are correct class
gsea_df <- sig_genesets
gsea_df$padj <- as.numeric(gsea_df$FDR)
gsea_df$pathway <- as.character(row.names(gsea_df))
gsea_df<-gsea_df[order(gsea_df$NGenes),]
gsea_df$readable_pathway <- gsub("_", " ", gsea_df$pathway )


 #if more than 20 significant terms only plot the 20 most significant
    if(nrow(gsea_df)>=20){
      p <- ggplot(gsea_df[((nrow(gsea_df)-19):nrow(gsea_df)),], aes(x = fct_reorder(readable_pathway, - FDR), y = NGenes))+
  geom_bar(aes(fill = padj), stat = "identity", width = 0.8) +
  scale_fill_continuous(low = "red", high = "blue", na.value = NA, name = "p.adjust")+
  scale_x_discrete(expand = expansion(add = .5)) +
  coord_flip()+ try(facet_grid(~Direction))  + theme_bw() + labs(y="Count", x="")  +ggtitle(i) + 
theme(axis.text=element_text(size=17))
      try(print(p))
 try(ggsave(p, file=paste0("Reactome_barplot_", i,".png"), width = 35, height = 20, units = "cm")) 
    }else{
      p <- ggplot(gsea_df, aes(x = fct_reorder(readable_pathway, - FDR), y = NGenes))+
  geom_bar(aes(fill = padj), stat = "identity", width = 0.8) +
  scale_fill_continuous(low = "red", high = "blue", na.value = NA, name = "p.adjust")+
  scale_x_discrete(expand = expansion(add = .5)) +
  coord_flip()+ try(facet_grid(~Direction))  + theme_bw() + labs(y="Count", x="")  +ggtitle(i) + 
theme(axis.text=element_text(size=17
                             ))
try(ggsave(p, file=paste0("Reactome_barplot_", i,".png"), width = 35, height = 20, units = "cm"))
try(print(p))
    }

# #create barplot of enriched pathways:
# p <- ggplot(gsea_df, aes(x = fct_reorder(pathway, - FDR), y = NGenes))+
#   geom_bar(aes(fill = gsea_df$padj), stat = "identity", width = 0.8) +
#   scale_fill_continuous(low = "red", high = "blue", na.value = NA, name = "p.adjust")+
#   scale_x_discrete(expand = expansion(add = .5)) +
#   coord_flip()+ try(facet_grid(~Direction))  + theme_bw() + labs(y="Count", x="")  +ggtitle(i) + 
# theme(axis.text=element_text(size=12))

#try(ggsave(p, file=paste0("Reactome_barplot_", i,".png"), width = 30, height = 20, units = "cm"))

  ### add gene IDs to the dataframe and save to csv
#create empty list
my_list <- list()
for (j in gsea_df$pathway){
#get the reactome indices for i
  pathway <- reactome.ind[[j]]
  #subset edgerdata for gene set
 genes_in_geneset <- edgeRData$genes$symbol[pathway]
 #append into a list
  my_list[[length(my_list) + 1]] <- genes_in_geneset    # Append new list element
}

#name list elements as gsea pathways
names(my_list) <- c(gsea_df$pathway)
#combine list with gsea dataframe
gsea_df$gene_ids <- my_list[match(gsea_df$pathway, names(my_list))]
#chnage to character so can save as excel file
gsea_df <- apply(gsea_df,2,as.character)

#save results:
write.csv(gsea_df, file=paste0("Reactome_", i,".csv") )


 
```

#HALLMARKS CAMERA

```{r}

#################### HALLMARKS ###############################

#use design matrix from GLM

#CAMERA takes as input the DGEList object dgeObj, the indexed list of gene sets c2.ind, the design matrix, the contrast being tested, as well as some other arguments. By default, CAMERA can estimate the correlation for each gene set separately. However, in practise, it works well to set a small inter-gene correlation of about 0.05 using the inter.gene.cor argument.

# Get pathways from gmt file
Hallmarks_pathways <- gmtPathways("C:/Users/Christina/Documents/GSEA .gmt/h.all.v2024.1.Hs.symbols.gmt")
#create index 
hallmarks.ind <- ids2indices( Hallmarks_pathways, edgeRData$genes$symbol)
# Run analysis:
gst.camera <- camera.DGEList(edgeRData, design= design, index=hallmarks.ind,contrast=con, inter.gene.cor=0.05)

#get number of signif gene sets at FDR 0.05
table(gst.camera$FDR < 0.05)

#find signif gene sets at FDR 0.05
sig_genesets <- gst.camera[gst.camera$FDR < 0.05, ]
print(sig_genesets)
 
#ensure cols are correct class
gsea_df <- sig_genesets
gsea_df$padj <- as.numeric(gsea_df$FDR)
gsea_df$pathway <- as.character(row.names(gsea_df))
gsea_df$readable_pathway <- gsub("_", " ", gsea_df$pathway )

 #if more than 20 significant terms only plot the 20 most significant
    if(nrow(gsea_df)>=20){
      p <- ggplot(gsea_df[((nrow(gsea_df)-19):nrow(gsea_df)),], aes(x = fct_reorder(readable_pathway, - FDR), y = NGenes))+
  geom_bar(aes(fill = padj), stat = "identity", width = 0.8) +
  scale_fill_continuous(low = "red", high = "blue", na.value = NA, name = "p.adjust")+
  scale_x_discrete(expand = expansion(add = .5)) +
  coord_flip()+ try(facet_grid(~Direction))  + theme_bw() + labs(y="Count", x="")  +ggtitle(i) + 
theme(axis.text=element_text(size=17))
      try(print(p))
 try(ggsave(p, file=paste0("Reactome_barplot_", i,".png"), width = 35, height = 20, units = "cm")) 
    }else{
      p <- ggplot(gsea_df, aes(x = fct_reorder(readable_pathway, - FDR), y = NGenes))+
  geom_bar(aes(fill = padj), stat = "identity", width = 0.8) +
  scale_fill_continuous(low = "red", high = "blue", na.value = NA, name = "p.adjust")+
  scale_x_discrete(expand = expansion(add = .5)) +
  coord_flip()+ try(facet_grid(~Direction))  + theme_bw() + labs(y="Count", x="")  +ggtitle(i) + 
theme(axis.text=element_text(size=17
                             ))
try(ggsave(p, file=paste0("Reactome_barplot_", i,".png"), width = 35, height = 20, units = "cm"))
try(print(p))}

### add gene IDs to the dataframe 
#create empty list
my_list <- list()
for (j in gsea_df$pathway){
#get the indices for i
  pathway <- hallmarks.ind[[j]]
  #subset edgerdata for gene set
 genes_in_geneset <- edgeRData$genes$symbol[pathway]
 #append into a list
  my_list[[length(my_list) + 1]] <- genes_in_geneset    # Append new list element
}

#name list elements as gsea pathways
names(my_list) <- c(gsea_df$pathway)
#combine list with gsea dataframe
gsea_df$gene_ids <- my_list[match(gsea_df$pathway, names(my_list))]
#chnage to character so can save as excel file
gsea_df <- apply(gsea_df,2,as.character)

#save results:
write.csv(gsea_df, file=paste0("Hallmarks_", i,".csv") )



```

#GO CAMERA
```{r}
# Get pathways from gmt file
GO_pathways <- gmtPathways("C:/Users/Christina/Documents/GSEA .gmt/c5.all.v2024.1.Hs.symbols.gmt")
#create index 
GO.ind <- ids2indices( GO_pathways, edgeRData$genes$symbol)
# Run analysis:
gst.camera <- camera.DGEList(edgeRData, design= design, index=GO.ind,contrast=con, inter.gene.cor=0.05)

#get number of signif gene sets at FDR 0.05
table(gst.camera$FDR < 0.05)

#find signif gene sets at FDR 0.05
sig_genesets <- gst.camera[gst.camera$FDR < 0.05, ]
print(sig_genesets)
 
#ensure cols are correct class
gsea_df <- sig_genesets
gsea_df$padj <- as.numeric(gsea_df$FDR)
gsea_df$pathway <- as.character(row.names(gsea_df))
gsea_df<-gsea_df[order(gsea_df$NGenes),]
gsea_df$readable_pathway <- gsub("_", " ", gsea_df$pathway )

 #if more than 20 significant terms only plot the 20 most significant
    if(nrow(gsea_df)>=20){
      p <- ggplot(gsea_df[((nrow(gsea_df)-19):nrow(gsea_df)),], aes(x = fct_reorder(readable_pathway, - FDR), y = NGenes))+
  geom_bar(aes(fill = padj), stat = "identity", width = 0.8) +
  scale_fill_continuous(low = "red", high = "blue", na.value = NA, name = "p.adjust")+
  scale_x_discrete(expand = expansion(add = .5)) +
  coord_flip()+ try(facet_grid(~Direction))  + theme_bw() + labs(y="Count", x="")  +ggtitle(i) + 
theme(axis.text=element_text(size=16))
 try(ggsave(p, file=paste0("GO_barplot_", i,".png"), width = 35, height = 20, units = "cm")) 
      try(print(p))
    }else{
      p <- ggplot(gsea_df, aes(x = fct_reorder(readable_pathway, - FDR), y = NGenes))+
  geom_bar(aes(fill = padj), stat = "identity", width = 0.8) +
  scale_fill_continuous(low = "red", high = "blue", na.value = NA, name = "p.adjust")+
  scale_x_discrete(expand = expansion(add = .5)) +
  coord_flip()+ try(facet_grid(~Direction))  + theme_bw() + labs(y="Count", x="")  +ggtitle(i) + 
theme(axis.text=element_text(size=16))
try(ggsave(p, file=paste0("GO_barplot_", i,".png"), width = 35, height = 20, units = "cm"))
try(print(p))
    }
# #create barplot of enriched pathways:
# p <- ggplot(gsea_df, aes(x = fct_reorder(pathway, - FDR), y = NGenes))+
#   geom_bar(aes(fill = gsea_df$padj), stat = "identity", width = 0.8) +
#   scale_fill_continuous(low = "red", high = "blue", na.value = NA, name = "p.adjust")+
#   scale_x_discrete(expand = expansion(add = .5)) +
#   coord_flip()+ try(facet_grid(~Direction))  + theme_bw() + labs(y="Count", x="")  +ggtitle(i) + 
# theme(axis.text=element_text(size=10))
# 
# try(ggsave(p, file=paste0("GO_barplot_", i,".png"), width = 30, height = 20, units = "cm"))

### add gene IDs to the dataframe 
#create empty list
my_list <- list()
for (j in gsea_df$pathway){
#get the indices for i
  pathway <- GO.ind[[j]]
  #subset edgerdata for gene set
 genes_in_geneset <- edgeRData$genes$symbol[pathway]
 #append into a list
  my_list[[length(my_list) + 1]] <- genes_in_geneset    # Append new list element
}

#name list elements as gsea pathways
names(my_list) <- c(gsea_df$pathway)
#combine list with gsea dataframe
gsea_df$gene_ids <- my_list[match(gsea_df$pathway, names(my_list))]
#chnage to character so can save as excel file
gsea_df <- apply(gsea_df,2,as.character)

#save results:
write.csv(gsea_df, file=paste0("GO", i, ".csv") )

```

