---
title: "RMS RNA seq analysis"
output: html_document
date: "2022-11-18"
author: "Christina Putnam"
editor_options: 
  chunk_output_type: console
---

include=FALSE:

```{r, setup, include=FALSE}
knitr::opts_knit$set(root.dir = ("C:/Users/Christina/Documents/RMS RNA seq/Zoes data/zoes_raw_reads"))
library(knitr)
```

#load libraries
```{r}
#library(WGCNA)
library(readr)
library(ComplexHeatmap)
library(plyr)
library(factoextra)
library(scales)
library("factoextra")
library("FactoMineR")
  library(PCAtools)
library(magrittr)
library(GSVA)
library(circlize)
library(msigdbr)
library(ggpubr)
library(stringr)
library(clusterProfiler)
library(RColorBrewer)
library(edgeR)
library(org.Hs.eg.db)
library('biomaRt')
library(msigdb)
library(tidyverse)
library("gplots")
library(dplyr)
library("ggrepel")
library("readxl")
library("ggplot2")
library('knitr')
```


## Set up project
Load libraries, and read in data, check data dimensions
```{r}

setwd("~/RMS RNA seq/Zoes data/zoes_raw_reads")

#create count matrix by joining RawReadCountsOnly.out files
files <- list.files(path = "~/RMS RNA seq/Zoes data/zoes_raw_reads", pattern = "RawReadCounts")
counts <- map(files, 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)

#set data as numeric
rawData[1:24] = lapply(rawData[1:24], FUN = function(y){as.numeric(y)})

```

#Describe experimental groups
```{r}

########## DESCRIBE EXPERIMENTAL GROUPS #########################

#Make a variable to describe the experimental groups, timepoints and replicates for comparison

#for sample treatment

sampleID <- colnames(rawData)
treatment <- c()
for (i in 1:ncol(rawData)) {
  samplename <- sampleID[i]
  group <- if(grepl("RH30_P", samplename)) {print("RH30_P")
  } else if(grepl("RH_HD", samplename)) {print("RH_HD")}
  else if(grepl("RH_C7", samplename)) {print("RH_C7")
  } else if(grepl("RH_C13", samplename)) {print("RH_C13")
  }
  treatment[i] <- group
}

#for replicate
sampleID <- colnames(rawData)
replicate <- c()
for (i in 1:ncol(rawData)) {
  samplename <- sampleID[i]
  group <- if(grepl("1", samplename)) {print("1")
  } else if(grepl("2", samplename)) {print("2")}
  else if(grepl("3", samplename)) {print("3")}
   else if(grepl("4", samplename)) {print("4")}
   else if(grepl("5", samplename)) {print("5")}
   else if(grepl("6", samplename)) {print("6")}
  replicate[i] <- group
}

```

# Filtering to remove genes with low expression

Genes with very low counts across all libraries provide little evidence for differential expression and they interfere with some of the statistical approximations that are used later in the pipeline. They also add to the multiple testing burden when estimating false discovery rates, reducing power to detect differentially expressed genes. These genes should be filtered out prior to further analysis.

There are a few ways to filter out lowly expressed genes. When there are biological replicates in each group, in this case we have a sample size of 6 in each group, we favour filtering on a minimum counts per million threshold present in at least 6 samples. 6 represents the smallest sample size for each group in our experiment. In this dataset, we choose to retain genes if they are expressed at a counts-per-million (CPM) above 0.3 in at least 6 samples.

We’ll use the cpm function from the edgeR library (M D Robinson, McCarthy, and Smyth 2010) to generate the CPM values and then filter. Note that by converting to CPMs we are normalising for the different sequencing depths for each sample.

```{r}

# Obtain CPMs
myCPM <- cpm(rawData)
# Have a look at the output
head(myCPM)
# Let's have a look and see whether a threshold of 0.5 does indeed correspond to a count of about 10-15
# We will look at the first sample
plot(myCPM[,1],rawData[,1])
# Let us limit the x and y-axis so we can actually look to see what is happening at the smaller counts
plot(myCPM[,1],rawData[,1],ylim=c(0,50),xlim=c(0,3), ylab = "Raw counts", xlab= "CPM", main = "CPM vs raw counts")
# Add a vertical line at 0.5 CPM
abline(v=0.35)

#changing thresh to 0.35 to correspond to around 10-15 counts 


# Which values in myCPM are greater than 0.3?
thresh <- myCPM > 0.35
# This produces a logical matrix with TRUEs and FALSEs
head(thresh)

# Summary of how many TRUEs there are in each row
table(rowSums(thresh))

# we would like to keep genes that have at least 6 TRUES in each row of thresh
keep <- rowSums(thresh) >= 6
# Subset the rows of countdata to keep the more highly expressed genes
counts.keep <- rawData[keep,]
summary(keep)
dim(counts.keep)


#compare number of raw counts to number of filtered counts
dim(rawData)
dim(counts.keep)

```


# Create DGEList object
```{r}

edgeRData <- DGEList(counts=counts.keep, genes=rownames(counts.keep), group = treatment)
edgeRData
```

#Density plots
Desnity plot can be used to visually compare the filtered counts with the raw counts. In the raw counts there will be a high first peak which contains genes with no-low expression and a second peak containing expressed genes. The aim of filtering is to reduce the first peak and remove genes with low expression. 

```{r}
#first check density plot for rawData 
p <- cpm(rawData, log = TRUE) %>%
    as.data.frame() %>%
    pivot_longer(
        cols = everything(),
        names_to = "sample",
        values_to = "logCPM"
    ) %>%
    split(f = .$sample) %>%
    lapply(function(x){
        d <- density(x$logCPM)
        tibble(
            sample = unique(x$sample),
            x = d$x,
            y = d$y
        )
    }) %>%
    bind_rows()  %>%
    ggplot(aes(x, y, group = sample)) +
    geom_line()  +
    labs(
        x = "logCPM",
        y = "Density",
        colour = "Genotype"
    ) + ggtitle("Raw data")

print(p)

#can see the first peak containg genes with low expression

```

#Density plot filtered counts
```{r}
#check density plot for filtered counts
p <- cpm(counts.keep, log = TRUE) %>%
    as.data.frame() %>%
    pivot_longer(
        cols = everything(),
        names_to = "sample",
        values_to = "logCPM"
    ) %>%
    split(f = .$sample) %>%
    lapply(function(x){
        d <- density(x$logCPM)
        tibble(
            sample = unique(x$sample),
            x = d$x,
            y = d$y
        )
    }) %>%
    bind_rows()  %>%
    ggplot(aes(x, y, group = sample)) +
    geom_line()  +
    labs(
        x = "logCPM",
        y = "Density",
        colour = "Genotype"
    )+ ggtitle("Filtered data")

print(p)


#first peak has been greatly reduced
```


QC Plots before normalisation

Barplot of raw reads

```{r}
########## QC PLOTS ##################################
#make a bar plot of the number of reads, i.e. library size, in each sample


edgeRData$samples <- rownames_to_column(edgeRData$samples)
data <- as.tibble(edgeRData$samples )

#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(11,5,4,1)+.1)
barplot(data$lib.size,names=data$rowname, ylim=c(0,60000000),las=2,   cex.names = 0.8, ylab="Library size", xlab = "", cex.axis = 0.8);
title("Barplot of library sizes")
title(xlab="Sample", line=4, cex.lab=1)
dev.off()



```


#Boxplot of logCPM
```{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 

# Get log2 counts per million
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=5, cex.lab=1.1)
title("Boxplots of log-2 CPMs before normalisation")
```

#MDS plot

```{r}

treatment_names <- c("RH_C13", "RH_C7", "RH_HD", "RH30_P")

#MDS plot
col.group <- c( "red", "green",  "cyan",    "purple")[edgeRData$samples$group]

shps = c(15, 16, 17)
data.frame(edgeRData$samples$group, col.group)
plotMDS(edgeRData, col=col.group, cex=1.5, pch= 16)
legend(fill=c("red", "green","cyan","purple"),legend=treatment_names, cex=1, "topright")
  title("MDS plot")

  plotMDS(edgeRData)
  
```



#PCA
```{r}

# https://bioconductor.org/packages/devel/bioc/vignettes/PCAtools/inst/doc/PCAtools.html

#get normalised logged counts
logCPM <- cpm(edgeRData, log = TRUE)
CPM <- logCPM

#remove HD samples if needed
#CPM <- CPM[, !colnames(CPM) %in% c("RH_HD_1", "RH_HD_2", "RH_HD_3", "RH_HD_4", "RH_HD_5", "RH_HD_6")]

#sclae data to reduce influence of highly expressed genes
scaledata <- t(scale(t(CPM)))

#get metadata
   metadata <- data.frame(row.names = colnames(CPM))
  metadata$Treatment_condition <- rep(NA, ncol(CPM))

  #get sample group as extra col in df
sampleID <- row.names(metadata)
metadata$Treatment_condition <- c()
for (i in 1:nrow(metadata)) {
  samplename <- sampleID[i]
  Treatment_condition <- if(grepl("RH30_P", samplename)) {print("RH30_P")
  } else if(grepl("RH_HD", samplename)) {print("RH_HD")}
  else if(grepl("RH_C7", samplename)) {print("RH_C7")
  } else if(grepl("RH_C13", samplename)) {print("RH_C13")
  }
  metadata$Treatment_condition[i] <- Treatment_condition
}


#run PCA
 p <- pca(scaledata, metadata = metadata, removeVar = 0.1)

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

  screeplot(p,
    components = getComponents(p, 1:20),
    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))
 
 # biplot-----------------------------------------------------------------------------
   
  tiff(file="biplot_pc1vspc2_Zoes.tiff",
width=8, height=5, units="in", res=800)
  biplot(p,
          lab = NULL,
    colby = 'Treatment_condition', colkey = c("RH30_P"= "#000000", "RH_C13"= "#E69F00", "RH_C7"= "#56B4E9", "RH_HD"= "#CC79A7"),
    hline = 0, vline = 0,
    legendPosition = 'right', 
    ellipse = FALSE,
      ellipseType = 't',
      ellipseLevel = 0.95,
   
      ellipseAlpha = 1/4,
      ellipseLineSize = 1.0,
    colLegendTitle = "Treatment condition",
    pointSize = 5,
    title = 'PCA bi-plot',
    subtitle = 'PC1 versus PC2',)
  dev.off()
  

 
plotloadings(p,
    rangeRetain = 0.01,
    labSize = 4.0,
    title = 'Loadings plot',
    subtitle = 'PC1, PC2, PC3, PC4, PC5',
    caption = 'Top 1% variables',
    shape = 24,
    col = c('limegreen', 'black', 'red3', 'blue'),
    drawConnectors = TRUE)


  
  tiff(file="pairs_plot.tiff",
width=8, height=5, units="in", res=800)
  pairsplot(p,
      components = getComponents(p, c(1:4)),
      triangle = TRUE, trianglelabSize = 12,
      hline = 0, vline = 0,
      titleLabSize= 20,
      pointSize = 2,
      gridlines.major = FALSE, gridlines.minor = FALSE,
     colby = 'Treatment_condition', colkey = c("RH30_P"= "#000000", "RH_C13"= "#E69F00", "RH_C7"= "#56B4E9", "RH_HD"= "#CC79A7"),
      legendPosition = "none",
      title = 'Pairs plot', plotaxes = FALSE,
      margingaps = unit(c(-0.01, -0.01, -0.01, -0.01), 'cm'))
  dev.off()
```

#Normalisation 

> We will now adjust the raw read count data to account for any differences in sequencing depth between samples using the weighted trimmed mean of M-values method (TMM) (Robinson and Oshlack 2010). This calculates a global fold change between samples based on gene-wise log-fold-changes (M) and absolute expression levels (A) under the assumption that the majority of genes are not differentially expressed. These M and A values are trimmed by removing the upper and lower x% of the data before calculating the weighted average, hence a weighted trimmed mean. The weights are from the delta method for Binomial data. TMM normalisation (and most scaling normalisation methods) scale relative to one sample. 

```{r}

##################### NORMALISATION ###############################################

#use calcNormFactors function to calculate the normalization factors for each library. The product of these factors and the library sizes defines the effective library size. View the normalisation factors for each sample (norm.factors) by entering the object name.
edgeRData <- calcNormFactors(edgeRData, method="TMM")

```


#Boxplot of logCPM
```{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
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=5, cex.lab=1.1)
title("Boxplots of log-2 CPMs after normalisation")
```


> Stand-alone programs such as edgeR use log normalization and the negative binomial model to correct for differences due to both library size and composition. For the negative binomial model, we need to estimate a dispersion parameter but this is difficult because there are not enough samples to get an accurate estimate. Therefore, information is ‘shared’ between genes to compute a genewise dispersion (tagwise dispersion). In edgeR, an empirical Bayes method is then used to ’shrink’ the genewise dispersion estimates towards the common dispersion.

>When estimating the common dispersion parameter, edgeR also reports the biological coefficient of variation (BCV) which is given by the square root of the dispersion estimate. BCV is the relative variability of expression between biological replicates. For example, a BCV=0.24 means that the expression values vary up and down by 24% between replicates. Importantly, BCV measures the relative variability of true expression levels, not the variability of measured expression levels. It represents the variability that remains after the Poisson variability from sequencing has been removed.

>For human replicates, common dispersion values are usually between 0.2-0.4 which means expression varies between 20-40% (BCV) between human replicates. Dispersion is expected to be lower for cell lines (~0.1) because they are genetically identical.

```{r}
#Create design matrix 
design <- model.matrix(~0+group, data=edgeRData$samples)
colnames(design) <- levels(edgeRData$samples$group)
design


# create contrast matrix
my.contrasts <- makeContrasts(RH30_P_vs_RH_HD=RH_HD-RH30_P,
                              RH30_P_vs_RH_C7=RH_C7-RH30_P,
                              RH30_P_vs_RH_C13=RH_C13-RH30_P,
                              levels=design)

#Estimate dispersion
 edgeRData <- estimateDisp(edgeRData, design, robust=TRUE)


# note values of common dispersion and BCV
 edgeRData$common.dispersion
```

#BCV plot

> Plot the genewise biological coefficient of variation (BCV) against gene abundance (in log2 counts per million). The red and blue lines on this plot show the BCV and trended dispersion.

```{r}
plotBCV(edgeRData, xlab="Average log CPM", ylab="Biological coefficient of variation", pch=14, cex=0.2, col.common="red", col.trend="blue", col.tagwise="black")

 
  edgeRData$common.dispersion
plotBCV(edgeRData)

```


#DGE test

> All comparisons ran pairwise, number of significantly up and down DE with FDR <0.05 summarized, and top 10 DE presented.

```{r}

fit <- glmQLFit(edgeRData, design)
RH30_P_vs_RH_HD <- glmQLFTest(fit, contrast=my.contrasts[,"RH30_P_vs_RH_HD"])
RH30_P_vs_RH_C7 <- glmQLFTest(fit, contrast=my.contrasts[,"RH30_P_vs_RH_C7"])
RH30_P_vs_RH_C13 <- glmQLFTest(fit, contrast=my.contrasts[,"RH30_P_vs_RH_C13"])


### add FDR to each sample ########
RH30_P_vs_RH_HD$table$FDR <- p.adjust(RH30_P_vs_RH_HD$table$PValue, method = "BH")
RH30_P_vs_RH_C7$table$FDR <- p.adjust(RH30_P_vs_RH_C7$table$PValue, method = "BH")
RH30_P_vs_RH_C13$table$FDR <- p.adjust(RH30_P_vs_RH_C13$table$PValue, method = "BH")




```



#DGE visualisation


```{r}

dge_files <- list(RH30_P_vs_RH_HD, RH30_P_vs_RH_C7, RH30_P_vs_RH_C13)

filenames <- c("RH30_P_vs_RH_HD", "RH30_P_vs_RH_C7", "RH30_P_vs_RH_C13")

for (i in 1:length(dge_files)) {
  
  #Save number of up/down significant genes at FDR = 0.05
  j <- filenames[i]
  d <- get(j)
#  dge_sum <-summary(de <- decideTestsDGE(d, adjust.method = "BH", p.value = 0.05 ,lfc = 1))
#  print(dge_sum)
#  dge_sum <- data.frame(rbind(dge_sum))
  # Print top genes
#  print(topTags(d, 10, adjust.method = "BH", p.value = 0.05, sort.by = "PValue"))
#  detags <- rownames(edgeRData)[as.logical(de)]
#  plotSmear(d, de.tags=detags, cex=1, xlab="Average logCPM", ylab="logFC")
  # make subset of data
  data <- cbind(d$table)
  # add FDR column to dge object
  data$FDR <- p.adjust(data$PValue, method = "BH")
  head(data)
 
  
  
#Volcano Plots ----------------------------------------------------------
  

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



#save list of unfiltered DEGs
#write.csv(data, file=paste0("unfiltered_DEGs_", j,".csv"))


#get DEGs
#data <- data[data$Expression != "Unchanged", ] 
 
#save list of unfiltered DEGs
#write.csv(data, file=paste0("DEGs_LFC_1_", j,".csv"))

  

  # 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)
  

 #plot vol plot highlighting genelist and adding labels
 p3 <- ggplot(data = data, # Original data  
       aes(x = logFC, y = -log(FDR,10))) + 
  geom_hline(yintercept = -log10(0.05),
             linetype = "dashed") + 
  geom_vline(xintercept = c(-1, 1),
             linetype = "dashed") +  
 scale_colour_manual(values = cols)+
    ggtitle(filenames[[i]]) +
    geom_point(aes(colour = Expression), 
             alpha = 0.8, 
             shape = 16,
             size = 1)  
   top <- 10
  top_genes <- bind_rows(
    data %>%
      filter(Expression == 'Up-regulated') %>%
      arrange(FDR, desc(abs(logFC))) %>%
      head(top),
    data %>%
     filter(Expression == 'Down-regulated') %>%
      arrange(FDR, desc(abs(logFC))) %>%
      head(top)

  )
  top_genes %>%
    kable()

  p3 <- p3 +
    geom_label_repel(data = top_genes, label.size = 0.5, force = 0.5, max.overlaps = 15, 
                     mapping = aes(logFC, -log(FDR,10), label = row.names(top_genes)),
                     size = 4.5)
 
 #ggsave(p3, file=paste0("Volcano_plot_", j,".png"), width = 25, height = 25, units = "cm")
 
 
  #vol plot with labels of specific genes---------------------------------------------------------------

   #add expression
   vol_plot_data <- data %>% 
  mutate(
    Expression = case_when(logFC > 1.5 & FDR < 0.05 ~ "Up-regulated",
                           logFC < -1.5 & FDR < 0.05 ~ "Down-regulated",
                           TRUE ~ "Unchanged")
    )
head(data) %>% 
 kable()


  # 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)
  
# create gene list of genes to highlist

  
  gene_list <- c("ABCB1", "ABCB4", "ABCC2", "ABCB5")
#  gene_list <- diff_geneset
 #subset data selecting for genes in genelist
 ils <- vol_plot_data[row.names(vol_plot_data) %in% gene_list,]
 
 #tiff("Volplot_NB_diff_Isotretinoin_and_TAZvsDMSO_PT3.tiff", units = "in", width=8, height=6, res=300)

 #plot vol plot highlighting genelist and adding labels
 p3 <- ggplot(data = vol_plot_data, # Original data  
       aes(x = logFC, y = -log(FDR,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(filenames[[i]]) +
    geom_point(aes(colour = Expression), 
             alpha = 0.8, 
             shape = 16,
             size = 1)  + geom_point(data = ils, # New layer containing data subset il_genes
             size = 2,
             shape = 21,
             fill = "yellow",
             colour = "black") +  geom_label_repel(data = ils, label.size = 0.5, force = 0.5, max.overlaps = 15, 
                     mapping = aes(logFC, -log(FDR,10), label = rownames(ils)),
                     size = 3.8) # Add labels last to appear as the toplayer
                  
 print(p3)
  
 ggsave(p3, file=paste0("Volcano_plot_MDR_genes", j,".png"), width = 15, height = 10, units = "cm")

 # dev.off()
  
}

```


#Number of DEGs

```{r}

df <- read_xlsx("C:/Users/Christina/Documents/RMS RNA seq/Zoes data/No_DEGs.xlsx")


# Stacked
p <- ggplot(df, aes(fill=Direction, y=Value, x=Treatment)) + 
    geom_bar(position="stack", stat="identity") + ylab("Number of genes") + ggtitle("Number of DEGs intrinsic resistance") + theme(axis.text=element_text(size=14),
        axis.title=element_text(size=14))+
  scale_fill_manual(values=c('deepskyblue1', 'indianred1'))
  
  ggsave(p, file=paste0("No_DEGs_Zoes.png"), width = 25, height = 25, units = "cm")
```
  
#boxplot CPM
```{r}

myCPM <- cpm(edgeRData$counts)
 x <- myCPM['DDX1',] 
  m <- list(counts = as.numeric(x), group = as.factor(edgeRData$samples$group))
   m <- as.tibble(m)
 
#check reorder
   
 #  m$group <- sub("COMB1", "RA+TAZ", m$group)
  # m$group <- sub("COMB2", "RA+GSK", m$group)

 
 # add new column with varaible to seperate by timepoint 
  #names <- m$group
 # for (i in 1:54) {
#    string <- names[i]
#m$timepoint[i] <- if(grepl("PT3", string)) {print("Post-treatment day 3")} else if (grepl("PT7", string)) {print("Post-treatment day 7")} else {print("Post-treatment day 10")}}
  
   
   m$group <- factor(m$group, levels=c('DMSO', "RA", "TAZ", "GSK", "RA+TAZ", "RA+GSK") )
 
                        
  q <- ggplot(m, aes(group, counts)) + geom_boxplot() + geom_jitter(width = 0.1)
  q <- q + labs(x = "", y = "Normalized Counts ", title = "RARB") 
  q <- q + theme(axis.text.x = element_text(angle = 90, vjust = 1, hjust = 1, size=15) ,axis.title.y  = element_text(size = 15))

             
   print(q)
```


  
#MYCN
```{r}

myCPM <- cpm(edgeRData$counts)
 x <- myCPM['MYCN',] 
  m <- list(counts = as.numeric(x), group = as.factor(edgeRData$samples$group))
   m <- as.tibble(m)
 
   m$group <- factor(m$group, levels=c('RH30_P', "RH_C7", "RH_C13", "RH_HD") )
 
                        
  q <- ggplot(m, aes(group, counts)) + geom_boxplot() + geom_jitter(width = 0.1)
  q <- q + labs(x = "", y = "Normalized Counts ", title = "MYCN expression") 
  q <- q + theme(axis.text.x = element_text(angle = 90, vjust = 1, hjust = 1, size=15) ,axis.title.y  = element_text(size = 15))

   print(q)
   
```


#Scatterplot function
```{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("magenta", "hotpink1", 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))   + 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 ", 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")
print(scatterplot)

}

```




#Scatterplots of logfc 
```{r}


######### for C13 and C7 #######################
RH30_P_vs_RH_C13$table$FDR <- p.adjust(RH30_P_vs_RH_C13$table$PValue, method = "BH")
 RH30_P_vs_RH_C13$table <- RH30_P_vs_RH_C13$table %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "significant",
                           logFC < -1 & FDR < 0.05 ~ "significant",
                           TRUE ~ "Unchanged") )
 
RH30_P_vs_RH_C7$table$FDR <- p.adjust(RH30_P_vs_RH_C7$table$PValue, method = "BH")
 RH30_P_vs_RH_C7$table <- RH30_P_vs_RH_C7$table %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "significant",
                           logFC < -1 & FDR < 0.05 ~ "significant",
                           TRUE ~ "Unchanged") )
 
df <- merge(RH30_P_vs_RH_C7$table, RH30_P_vs_RH_C13$table, by = 'row.names', all = TRUE)

df <- df %>%
   mutate(sig = case_when(Expression.x == "significant" & Expression.y == "significant" ~ "Both",
                          Expression.x == "significant" ~ "C7",
                           Expression.y == "significant" ~ "C13",
                          
                           TRUE ~ "Unchanged") )
  


scatterplot <- ggplot(df, aes(x=logFC.x, y=logFC.y, color=sig)) + geom_point()   + geom_smooth(method=lm , color="blue", fill="#69b3a2", se=TRUE) + scale_color_manual(values=c("mediumorchid1", "steelblue1", "hotpink1", "grey")) + ggtitle("Differentially expressed genes for C7 vs P and C13 vs P") +
  xlab("Log-2 fold change C13 vs P") + ylab("Log-2 fold change C7 vs P") +theme(axis.text=element_text(size=12),
        axis.title=element_text(size=12)) + guides(color=guide_legend(title="Treatment differentially expressed in")) 

ggplot(df, aes(x=logFC.x, y=logFC.y)) +
    geom_point() +
  geom_smooth(method=lm , color="red", fill="#69b3a2", se=TRUE) +
  stat_regline_equation(label.y = 10) +
  stat_cor(label.y = 6)

#ggsave(scatterplot, file=paste0("Scatterplot_TAZ_iso_PT10.png"), width = 25, height = 25, units = "cm")
# + scale_color_viridis(discrete = TRUE, option = "D")  + stat_cor(method = "pearson",  label.x = -5, label.y = 10)





########## for HD and C13 ###########


RH30_P_vs_RH_HD$table$FDR <- p.adjust(RH30_P_vs_RH_HD$table$PValue, method = "BH")
 RH30_P_vs_RH_HD$table <- RH30_P_vs_RH_HD$table %>% 
  mutate(
    Expression = case_when(logFC > 1.5 & FDR < 0.05 ~ "significant",
                           logFC < -1.5 & FDR < 0.05 ~ "significant",
                           TRUE ~ "Unchanged") )
 
df <- merge(RH30_P_vs_RH_HD$table, RH30_P_vs_RH_C13$table, by = 'row.names', all = TRUE)

df <- df %>%
   mutate(sig = case_when(Expression.x == "significant" & Expression.y == "significant" ~ "Both",
                          Expression.x == "significant" ~ "HD",
                           Expression.y == "significant" ~ "C13",
                          
                           TRUE ~ "Unchanged") )
  
 #df <- df %>%
#   mutate(Significant = case_when(Expression.y == "significant" ~ "RA",TRUE ~ "Unchanged") )
  
   
ggplot(df, aes(x=logFC.x, y=logFC.y)) +
    geom_point() +
  geom_smooth(method=lm , color="red", fill="#69b3a2", se=TRUE) #+ scale_color_viridis(discrete = TRUE, option = "D")

ggplot(df, aes(x=logFC.x, y=logFC.y, color=sig)) + geom_point() + geom_smooth(method=lm , color="red", fill="#69b3a2", se=TRUE) + ggtitle("RH30 C13 and RH30 HD")# + scale_color_viridis(discrete = TRUE, option = "D")

############## for C7 and HD #####################
df <- merge(RH30_P_vs_RH_HD$table, RH30_P_vs_RH_C7$table, by = 'row.names', all = TRUE)

df <- df %>%
   mutate(sig = case_when(Expression.x == "significant" & Expression.y == "significant" ~ "Both",
                          Expression.x == "significant" ~ "HD",
                           Expression.y == "significant" ~ "C7",
                          
                           TRUE ~ "Unchanged") )

ggplot(df, aes(x=logFC.x, y=logFC.y, color=sig)) + geom_point() + geom_smooth(method=lm , color="red", fill="#69b3a2", se=TRUE) + ggtitle("RH30 C7 and RH30 HD")# + scale_color_viridis(discrete = TRUE, option = "D")
```


# ORA on overlapping downregulated DEGs in clones
```{r}
#find genes that are downregulated in C7 and C13 
data <- RH30_P_vs_RH_C7$table
   data$FDR <- p.adjust(data$PValue, method = "BH")
   data$geneID <- row.names(data)
  data <- data %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "Up-regulated",
                           logFC < -1 & FDR < 0.05 ~ "Down-regulated",
                           TRUE ~ "Unchanged"))
  

DOWN_DEG_list_RH30_P_vs_RH_C7 <- data[data$Expression == 'Down-regulated',]

data <- RH30_P_vs_RH_C13$table
   data$FDR <- p.adjust(data$PValue, method = "BH")
   data$geneID <- row.names(data)
  data <- data %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "Up-regulated",
                           logFC < -1 & FDR < 0.05 ~ "Down-regulated",
                           TRUE ~ "Unchanged"))
  

DOWN_DEG_list_RH30_P_vs_RH_C13 <- data[data$Expression == 'Down-regulated',]

#put all data frames into list
df_list <- list(DOWN_DEG_list_RH30_P_vs_RH_C13, DOWN_DEG_list_RH30_P_vs_RH_C7)
#merge all data frames in list
overlapping_DEGS <- df_list %>% reduce(merge, by="geneID")
overlapping_DEGS <- overlapping_DEGS$geneID

#perform ORA
geneList <- overlapping_DEGS

#set wd
setwd("C:/Users/Christina/Documents/RMS RNA seq/Zoes data/ORA on overlapping downregulated genes in clones")

#GO terms-----------------------------------------------------------------------
options(enrichplot.colours = c("red","blue"))
ego2 <- enrichGO(gene         = geneList,
                 OrgDb         = org.Hs.eg.db,
                 keyType       = 'SYMBOL',
                 ont           = "ALL",
                 pAdjustMethod = "BH",
                 pvalueCutoff  = 0.05,
                 qvalueCutoff  = 0.05)
head(ego2, 10)   
## Dotplot 
ego2 <- ego2 %>%
  dplyr::arrange(p.adjust)

GO <- try(barplot(ego2,  font.size = 16 ,showCategory=20, x="Count"))
try(ggsave(GO, file=paste0("GO_downregulated_overlapping_clones_png"), width = 25, height = 33, units = "cm"))


# Hallmark ORA --------------------------------------------------------------------------
m_t2g <- msigdbr(species = "Homo sapiens", category = "H") %>% dplyr::select(gs_name, gene_symbol)
head(m_t2g)
em <- enricher(geneList, pvalueCutoff = 0.05,  pAdjustMethod = "BH", TERM2GENE=m_t2g)
head(em)
write.csv(em, file=paste0("Hallmark_hub_genes_", i,".csv"))

#barplot
H <- print(ggplot(em,aes(x=Count,y=Description,fill=p.adjust))+
  scale_fill_continuous(low = "red", high = "blue", na.value = NA)+geom_bar(stat='identity')+theme_classic()+
                 scale_x_continuous(expand=c(0,0))+labs(fill = 'p.adjust')+xlab('Count')+ylab('')+theme(axis.text=element_text(size=16),axis.title=element_text(size=12)))



try(ggsave(H, file=paste0("Hallmark_.png"), width = 25, height = 10, units = "cm"))

#KEGG -----------------------------------------------------------------------

#convert to symbols
options(enrichplot.colours = c("red","blue"))
symbols <- mapIds(org.Hs.eg.db, geneList, 'ENTREZID', 'SYMBOL') %>% as.data.frame()
symbols <- symbols$.
kk <- enrichKEGG(gene         = symbols,
                 organism     = 'hsa',
                 keyType       = 'ncbi-geneid',
                 pAdjustMethod = "BH",
                 pvalueCutoff = 0.05, 
                 qvalueCutoff = 0.05)
head(kk)
write.csv(kk, file=paste0("Keg_hub_genes_", i,".csv"))
keg <- barplot(kk, font.size = 16, showCategory=20,x="Count")

try(ggsave(keg, file=paste0("KEGG_.png"), width = 25, height = 20, units = "cm"))

#Reactome ---------------------------------------------------------------------
options(enrichplot.colours = c("red","blue"))
reactome <- enrichPathway(gene = symbols, pvalueCutoff = 0.05, readable=TRUE, qvalueCutoff =0.05, pAdjustMethod = "BH", organism = "human")
head(reactome)
write.csv(reactome, file=paste0("Reactome_hub_genes_", i,".csv"))
reactome <- barplot(reactome, showCategory=20, font.size = 16,x="Count")

try(ggsave(reactome, file=paste0("Reactome_hub_genes_", i,".png"), width = 25, height = 35, units = "cm"))



```





#Compare lists of DEGs to find overlapping DEGs in ALL samples
```{r}



#get list of DEGS

#for RH30_P_vs_RH_HD

data <- RH30_P_vs_RH_HD$table
   data$FDR <- p.adjust(data$PValue, method = "BH")
   data$geneID <- row.names(data)
  data <- data %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "Up-regulated",
                           logFC < -1 & FDR < 0.05 ~ "Down-regulated",
                           TRUE ~ "Unchanged"))
DEG_list_RH30_P_vs_RH_HD<- data[data$Expression == 'Up-regulated' | data$Expression == 'Down-regulated',]
UP_DEG_list_RH30_P_vs_RH_HD <- data[data$Expression == 'Up-regulated',]
DOWN_DEG_list_RH30_P_vs_RH_HD <- data[data$Expression == 'Down-regulated',]

#write.csv(DEG_list_RH30_P_vs_RH_HD, "C:/Users/Christina/Documents/RMS RNA seq/Zoes data/DEG_list_RH30_P_vs_RH_HD.csv")


# for RH30_P_vs_RH_C7
data <- RH30_P_vs_RH_C7$table
   data$FDR <- p.adjust(data$PValue, method = "BH")
   data$geneID <- row.names(data)
  data <- data %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "Up-regulated",
                           logFC < -1 & FDR < 0.05 ~ "Down-regulated",
                           TRUE ~ "Unchanged"))
DEG_list_RH30_P_vs_RH_C7 <- data[data$Expression == 'Up-regulated' | data$Expression == 'Down-regulated',]
UP_DEG_list_RH30_P_vs_RH_C7 <- data[data$Expression == 'Up-regulated',]
DOWN_DEG_list_RH30_P_vs_RH_C7 <- data[data$Expression == 'Down-regulated',]

#write.csv(DEG_list_RH30_P_vs_RH_C7, "C:/Users/Christina/Documents/RMS RNA seq/Zoes data/DEG_list_RH30_P_vs_RH_C7.csv")




#for RH30_P_vs_RH_C13
data <- RH30_P_vs_RH_C13$table
   data$FDR <- p.adjust(data$PValue, method = "BH")
   data$geneID <- row.names(data)
  data <- data %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "Up-regulated",
                           logFC < -1 & FDR < 0.05 ~ "Down-regulated",
                           TRUE ~ "Unchanged"))
DEG_list_RH30_P_vs_RH_C13 <- data[data$Expression == 'Up-regulated' | data$Expression == 'Down-regulated',]
UP_DEG_list_RH30_P_vs_RH_C13 <- data[data$Expression == 'Up-regulated',]
DOWN_DEG_list_RH30_P_vs_RH_C13 <- data[data$Expression == 'Down-regulated',]

#write.csv(DEG_list_RH30_P_vs_RH_C13, "C:/Users/Christina/Documents/RMS RNA seq/Zoes data/DEG_list_RH30_P_vs_RH_C13.csv")


########## find overlapping upregulated DEGS ########################################
#put all data frames into list
df_list <- list(UP_DEG_list_RH30_P_vs_RH_C13, UP_DEG_list_RH30_P_vs_RH_C7, UP_DEG_list_RH30_P_vs_RH_HD)
#merge all data frames in list
overlapping_DEGS <- df_list %>% reduce(merge, by="geneID")
print(overlapping_DEGS)

########## find overlapping downregulated DEGS ########################################
#put all data frames into list
df_list <- list(DOWN_DEG_list_RH4ctrl_RH4vcr, DOWN_DEG_list_RH4ctrl_RH4ifo, DOWN_DEG_list_RMSYMctrl_RMSYMvcr, DOWN_DEG_list_RMSYMctrl_RMSYMifo)
#merge all data frames in list
overlapping_DEGS <- df_list %>% reduce(merge, by="geneID")



##### find all overlapping DEGS regardless of direction ################################
#put all data frames into list
df_list <- list(DEG_list_RH4ctrl_RH4vcr, DEG_list_RH4ctrl_RH4ifo, DEG_list_RMSYMctrl_RMSYMvcr, DEG_list_RMSYMctrl_RMSYMifo)
#merge all data frames in list
overlapping_DEGS <- df_list %>% reduce(merge, by="geneID")


```





Sense check genes


```{r}

myCPM <- cpm(edgeRData$counts)
 x <- myCPM['ST8SIA1',] 
  m <- list(counts = as.numeric(x), group = as.factor(edgeRData$samples$group))
   m <- as.tibble(m)
 
#check reorder
   
   m$group <- sub("COMB1", "RA+TAZ", m$group)
   m$group <- sub("COMB2", "RA+GSK", m$group)

 
 # add new column with varaible to seperate by timepoint 
  names <- m$group
  for (i in 1:54) {
    string <- names[i]
m$timepoint[i] <- if(grepl("PT3", string)) {print("Post-treatment day 3")} else if (grepl("PT7", string)) {print("Post-treatment day 7")} else {print("Post-treatment day 10")}}
  
  
     m$temp_timepoint = factor(m$timepoint, levels=c('Post-treatment day 3','Post-treatment day 7','Post-treatment day 10'))

   m$group <- sub("_PT3", "", m$group)
  m$group <- sub("_PT7", "", m$group)
   m$group <- sub("_PT10", "", m$group)
   
   m$group <- factor(m$group, levels=c('DMSO', "RA", "TAZ", "GSK", "RA+TAZ", "RA+GSK") )
 
q <- ggplot(data = m) +
    geom_dotplot(aes(y = counts, x = group, fill= group),dotsize = 2, stroke = NA,
                 binaxis = "y",         # which axis to bin along
                 binwidth = 0.05,       # Minimal difference considered diffeerent
                 stackdir = "center"    # Centered 
                 )   + stat_summary(aes(y = counts, x = group), fun.y = median, fun.ymin = median, fun.ymax = median,
                 geom = "crossbar", width = 0.5)             

 q <- q + facet_wrap(~temp_timepoint, scales = "free_x" )  +
       theme(strip.text.x = element_text(size = 15))
             
   print(q)
   
```
    
    

  
 


Number of DEGs

```{r}

#PT3 
df <- data.frame(
  treatment = c('RH4ctrl' ,'RH4ifo', 'RH4vcr','RMSYMctrl','RMSYMifo', 'RMSYMvcr') ,
  Direction = c('Down','Up'),
  Number = c(274, 390, 0, 654, 111, 263, 0, 496, 117, 569)
)

head(df)
p <- ggplot(df, aes(x = treatment, y = Number))+
  geom_col(aes(fill = Direction), width = 0.7)  + ylab("Number of genes") + ggtitle("Number of DEGs post-treatment day 3")

p


```

#Cluster profiler
GSEA using cluster profiler
```{r}
#function to do GSEA and dotplots

setwd("~/RMS RNA seq/Zoes data/Cluster profiler")
cluster_prolifer_plots <- function(msigdb_category){
  if(msigdb_category == "H"){
t2g <- msigdbr(species = "Homo sapiens", category = "H") %>% 
  dplyr::select(gs_name, gene_symbol)
head(t2g)
} else if(msigdb_category == "Reactome"){
t2g <- msigdbr(species = "Homo sapiens", category = "C2", subcategory = "CP:REACTOME") %>% 
  dplyr::select(gs_name, gene_symbol)
head(t2g)
}else if(msigdb_category == "GO"){
all_gene_sets = msigdbr(species = "Homo sapiens")
t2g <- all_gene_sets %>%
  dplyr::filter(gs_subcat %in% c("GO:BP","GO:CC","GO:MF")) %>% dplyr::select(gs_name, gene_symbol)
}

  set.seed(42)
em2 <- try(GSEA(geneList, TERM2GENE = t2g, nPermSimple = 100000))

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

#get top ten terms
top_ten <- em2@result$Description[1:10]

options(enrichplot.colours = c("red","blue"))
dotplot <- try(dotplot(em2, font.size=14, label_format = 60,x = "Count",
              split = ".sign",)+facet_grid(~.sign))

ggsave(dotplot, file=paste0(msigdb_category, "_", j ,".png"), width = 20, height = 20, units = "cm")
print(dotplot)

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

ggsave(core_plot, file=paste0(msigdb_category, "_", j ,"_coreplot.png"), width = 60, height = 60, units = "cm", bg = "white")

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

#do GSEA on samples--------------------------------------


dge_files <- list( RH30_P_vs_RH_HD)
filenames <- c("RH30_P_vs_RH_HD")

setwd("~/RMS RNA seq/Zoes data/Cluster profiler")

for (i in 1:length(dge_files)) {
  
  #get DGE object
  d <- filenames[i]
  d <- get(d)
  #get filename
  j <- filenames[i]

#get gene list
## feature 1: numeric vector (LFC)
geneList = d$table[,1]
## feature 2: named vector
names(geneList) = as.character(rownames(d$table))
## feature 3: decreasing orde
geneList = sort(geneList, decreasing = TRUE)

#perform GSEA function for GO
cluster_prolifer_plots(msigdb_category = "H")
#perform GSEA function for Hallmarks
#cluster_prolifer_plots(msigdb_category = "Reactome")
#perform GSEA function for reactome
#cluster_prolifer_plots(msigdb_category = "GO")

}

```



#Metaprogram scores
```{r}

#test for metaprogram scores -------------------------------------------------------

#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.

filenames <- c("RH30_P_vs_RH_HD", "RH30_P_vs_RH_C13", "RH30_P_vs_RH_C7")

for (i in filenames){
 diff_geneset <-  read_xlsx("~/RMS RNA seq/metaprogram_signatures_danielli.xlsx")

diff_geneset <- diff_geneset$Differentiated
   diff_geneset.ind <- ids2indices( diff_geneset, rownames(edgeRData$counts))
  # Run analysis:
  gst.camera <- camera.DGEList(edgeRData, design= design, index=diff_geneset,contrast=my.contrasts[, i], inter.gene.cor=0.05)

 #find signif gene sets at FDR 0.05
  gst.camera$FDR <-  p.adjust(gst.camera$PValue, method = "BH")
  sig_genesets <- gst.camera[gst.camera$PValue < 0.05, ]
  print(i)
  print(sig_genesets)
  try(rownames(sig_genesets) <- paste0(i))

  #output results to table
  write.table(sig_genesets, "RMS_metaprogram_sig_Differentiatedresults.csv", append = TRUE,sep = ",",col.names = FALSE)

#(read.csv(file = "test.csv"))


  
  
}

```


#Metaprogam GSVA
```{r}

filenames <- c("RH30_P_vs_RH_HD", "RH30_P_vs_RH_C13", "RH30_P_vs_RH_C7")

for (i in filenames){}
genes <-  read_xlsx("~/RMS RNA seq/metaprogram_signatures_danielli.xlsx")
genes <- genes$Progenitor
genes <- genes[1:171] %>% list()

gsvaPar <- gsvaParam(filtered_mapped_matrix, genes)
gsva_results <- gsva(gsvaPar)


gbmPar <- gsvaParam(CPM, genes, maxDiff=FALSE)
gbm_es <- gsva(gbmPar)

rownames(edgeRData) <- paste0("X", seq_len(nrow(edgeRData)))
edgeRData <- DGEList(counts=counts.keep, genes=rownames(counts.keep), group = treatment)
edgeRData <- calcNormFactors(edgeRData, method="TMM")
CPM <- cpm(edgeRData, log = TRUE)

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



gsvamethod <- gsvaParam(filtered_mapped_matrix,
  genes)
gsva_results <- gsva(gsvamethod)

#add phenotype data to df
pheno_data <- row.names(pd) %>% as.data.frame()
pheno_data$disease_state <- pd[,'Disease State'] #add to pheno data
pheno_data$Histology <- pd[,'Histology']
pheno_data$FFS <- pd[,'FFS']
pheno_data$Survival <- pd[,'Survival']
pheno_data$Grade <- pd[,'Disease Grade']
pheno_data$Stage <- pd[,'Disease Stage']
pheno_data$RG <- pd[,'RG']
pheno_data$Size <- pd[,'Size']
pheno_data$Node <- pd[,'Node']
pheno_data$Translocation <- pd[,'Translocation']
pheno_data$Invasion <- pd[,'Invasion']
pheno_data <- pheno_data %>% remove_rownames %>% column_to_rownames(var=".")


ARMS_samples <- subset(pd, Histology == "Alveolar" | Histology == "Mixed Alv/Emb") %>% rownames()
ERMS_samples <- subset(pd, Histology == "Embryonal") %>% rownames()


if(subtype=="ARMS"){
  #keep ARMS samples
pheno_data <- pheno_data[rownames(pheno_data) %in% ARMS_samples, ]
histology_colour <- c("violetred1","slateblue1")
}
if(subtype=="ERMS"){
  #keep ERMS samples only 
  pheno_data <- pheno_data[ rownames(pheno_data) %in% ERMS_samples, ]
  histology_colour <- c("navy")
} 

#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, pheno_data, by = "row.names")
colnames(gsva_results)[1] <- "sample"

  #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$Histology, 'Alveolar|Embryonal|Mixed Alv/Emb'))){histology_colour <- c("violetred1","slateblue1", "navy")}
  #  if(all(!str_detect(gsva_results$Histology, "Alveolar"))){histology_colour <- c("slateblue1")}
  #     if(all(!str_detect(gsva_results$Histology, "Embryonal"))){histology_colour <-  c("violetred1", "navy")}

#then selects options for clinical variable picked
  print(clinical_variable)
  
  #clinical_variable <- c("Stage")
    if (clinical_variable == "Stage") {
        variable_levels <- c("Stage 1", "Stage 2", "Stage 3", "Stage 4", "NA")
    readable_clinical_variable <- c("stage")
       my_comparisons <- list(c("Stage 1", "Stage 2"), c("Stage 1", "Stage 3"), c("Stage 1", "Stage 4"),c("Stage 2", "Stage 3"), c("Stage 2", "Stage 4"), c("Stage 3", "Stage 4"))
    
          #remove comparisons where no enoguh observations
       #in ARMS samples there is only one stage 1
    if(subtype=="ARMS"){
      my_comparisons <- list(c("Stage 2", "Stage 3"), c("Stage 2", "Stage 4"), c("Stage 3", "Stage 4"))
    }
       
   } else if(clinical_variable == "Grade"){
      variable_levels<-c("IA", "IB", "IIA", "IIB", "III", "IV", "NA")
  readable_clinical_variable <- c("grade")
   my_comparisons <- list(c("IA", "IB"), c("IA", "IIA"), c("IA", "IIB"),c("IA", "III"), c("IA", "IV"), c("IB", "IIA"), c("IB", "IIB"), c("IB", "III"), c("IB", "IV"), c("IIA", "IIB"), c("IIA", "III"), c("IIA", "IV"), c("IIB", "III"), c("IIB", "IV"), c("III", "IV"))
  
    #remove comparisons where no enoguh observations
  #in ARMS there is only 1 IB and IIB so can't compare that
    if(subtype=="ARMS"){
      my_comparisons <- list( c("IA", "IIA"),c("IA", "III"), c("IA", "IV"), c("IIA", "III"), c("IIA", "IV"), c("III", "IV"))
    }
   if(subtype == "ERMS"){
      my_comparisons <- list( c("IA", "IIA"), c("IA", "III"), c("IA", "IV"),   c("IIA", "III"), c("IIA", "IV"), c("III", "IV"))
  }
   
} else if(clinical_variable == "disease_state"){
  
    my_comparisons <- list( c("Alive", "Dead"))
  variable_levels <- c("Alive", "Dead", "NA")
    readable_clinical_variable <- c("disease state")
    
}  else if(clinical_variable == "RG"){
  
    my_comparisons <- list( c("LOW RISK", "INTERMEDIATE RISK"), c("LOW RISK", "HIGH RISK"), c("INTERMEDIATE RISK", "HIGH RISK"))
  variable_levels <- c("LOW RISK", "INTERMEDIATE RISK", "HIGH RISK", "NA")
    readable_clinical_variable <- c("risk group")
}else if(clinical_variable == "Size"){
  
    my_comparisons <- list( c("<= 5 cm", "> 5 cm"))
  variable_levels <-  c("<= 5 cm", "> 5 cm", "NA")
    readable_clinical_variable <- c("tumour size")
}else if(clinical_variable == "Node"){
  
    my_comparisons <- list( c("N-0", "N-1"))
  variable_levels <-   c("N-0", "N-1", "NA")
    readable_clinical_variable <- c("node")
}else if(clinical_variable == "Invasion"){
  
    my_comparisons <- list( c("T-1", "T-2"))
  variable_levels <-  c("T-1", "T-2", "NA")
    readable_clinical_variable <- c("invasion")
} 
 
#reorder plot
try(gsva_results[[clinical_variable]] <- factor(gsva_results[[clinical_variable]], variable_levels))
 try(fct_relevel(gsva_results[[clinical_variable]]))
 
    #boxplot
    p <- try(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")+ theme( axis.text.x = element_text(angle = 90,  size = 12))+ scale_fill_brewer(palette="BuPu", guide = "none")+geom_jitter(aes(col=Histology))+
  scale_color_manual(values = histology_colour)+ theme(legend.position = "right")+ stat_compare_means(comparisons = my_comparisons,label = "p.signif", method = paste0(correct_method)))
#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 = 5, height = 5)


if(clinical_variable == "FFS"){
  
    readable_clinical_variable <- c("FFS")
    gsva_results$FFS <- as.numeric(gsva_results$FFS)
  #make scatterplot
s <- ggplot(gsva_results,aes(x=FFS, y=gsva_results)) + 
    geom_point(aes( colour = Histology) )+scale_x_continuous(breaks = seq(0, 20))+ scale_color_manual(values=histology_colour) + ggtitle(paste0(resistance_type, " resistance score and ", readable_clinical_variable))+ylab("GSVA score")+
        geom_smooth(method = "lm")+
  stat_regline_equation(label.y = 1) +
  stat_cor(label.y = 1.1)+ xlab(paste0(clinical_variable, " (years)"))
 ggsave(s, file= paste0(resistance_type, "_scatterplot_", clinical_variable,"_", subtype, ".png")) 
  

}  else if(clinical_variable == "Survival"){
  
      readable_clinical_variable <- c("survival")
    gsva_results$Survival <- as.numeric(gsva_results$Survival)
  #make scatterplot
s <- ggplot(gsva_results,aes(x=Survival, y=gsva_results)) + 
    geom_point(aes( colour = Histology) )+scale_x_continuous(breaks = seq(0, 20))+ scale_color_manual(values=histology_colour) + ggtitle(paste0(resistance_type, " resistance score and ", readable_clinical_variable))+ylab("GSVA score")+
        geom_smooth(method = "lm")+
  stat_regline_equation(label.y = 1) +
  stat_cor(label.y = 1.1)+ xlab(paste0(clinical_variable, " (years)"))
 ggsave(s, file= paste0(resistance_type, "_scatterplot_", clinical_variable,"_", subtype, ".png")) 
  
} 

}
```








#GSEA REACTOME



```{r}

setwd("~/RMS RNA seq/Zoes data/GSEA results")
################################ REACTOME #################################


# 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.

filenames <- c( "RH30_P_vs_RH_HD", "RH30_P_vs_RH_C7" , "RH30_P_vs_RH_C13")
filenames <- c( "RH30_P_vs_RH_HD")

for (i in filenames){
# 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, rownames(edgeRData$counts))
# Run analysis:
gst.camera <- camera.DGEList(edgeRData, design= design, index=reactome.ind,contrast=my.contrasts[, i], 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$pathway_readable <- 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(pathway_readable, - 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(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(pathway_readable, - 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("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 
#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 <- rownames(edgeRData)[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") )
  
}
 
```

#GSEA HALLMARKS

```{r}

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

#use design matrix from GLM
# 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.

filenames <- c( "RH30_P_vs_RH_HD", "RH30_P_vs_RH_C13", "RH30_P_vs_RH_C7")
filenames <- c( "RH30_P_vs_RH_HD")
for (i in filenames){
# 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, rownames(edgeRData$counts))
# Run analysis:
gst.camera <- camera.DGEList(edgeRData, design= design, index=hallmarks.ind,contrast=my.contrasts[, i], 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$pathway_readable <- gsub("_", " ", gsea_df$pathway )

#create barplot of enriched pathways:
 #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(pathway_readable, - 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("Hallmark_barplot_", i,".png"), width = 35, height = 20, units = "cm")) 
    }else{
      p <- ggplot(gsea_df, aes(x = fct_reorder(pathway_readable, - 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))}

try(ggsave(p, file=paste0("Hallmarks_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 <- hallmarks.ind[[j]]
  #subset edgerdata for gene set
 genes_in_geneset <- rownames(edgeRData)[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") )
}


```

#GSEA GO

```{r}
##########for GO ##################################

#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.

filenames <- c( "RH30_P_vs_RH_HD","RH30_P_vs_RH_C7" , "RH30_P_vs_RH_C13")
filenames <- c( "RH30_P_vs_RH_HD")


for (i in filenames){
# 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, rownames(edgeRData$counts))
# Run analysis:
gst.camera <- camera.DGEList(edgeRData, design= design, index=GO.ind,contrast=my.contrasts[, i], 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$pathway_readable <- gsub("_", " ", gsea_df$pathway )

 #if more than 20 significant terms only plot the 20 most significant
    if(nrow(gsea_df)>=30){
      p <- ggplot(gsea_df[((nrow(gsea_df)-29):nrow(gsea_df)),], aes(x = fct_reorder(pathway_readable, - 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(pathway_readable, - 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 <- rownames(edgeRData)[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") )}

```

#fgsea
```{r}

dge_files <- list(RH30_P_vs_RH_HD, RH30_P_vs_RH_C7, RH30_P_vs_RH_C13)

filenames <- c("RH30_P_vs_RH_HD", "RH30_P_vs_RH_C7", "RH30_P_vs_RH_C13")

for (i in 1:length(dge_files)) {
  
  #Save number of up/down significant genes at FDR = 0.05
  j <- filenames[i]
  d <- get(j)
#  dge_sum <-summary(de <- decideTestsDGE(d, adjust.method = "BH", p.value = 0.05 ,lfc = 1))
#  print(dge_sum)
#  dge_sum <- data.frame(rbind(dge_sum))
  # Print top genes
#  print(topTags(d, 10, adjust.method = "BH", p.value = 0.05, sort.by = "PValue"))
#  detags <- rownames(edgeRData)[as.logical(de)]
#  plotSmear(d, de.tags=detags, cex=1, xlab="Average logCPM", ylab="logFC")
  # make subset of data
  data <- d$table

data$score <- -log10(data$PValue) * sign(data$logFC)

# Rank genes by score (higher absolute value means more important)
results_df <- data[order(abs(data$score), decreasing = TRUE), ]

# Assuming cox_results has your ranked genes
gene_ranks <- results_df$score
names(gene_ranks) <- rownames(results_df)

#H-----------------------------
# Download Hallmark gene sets for human
hallmark_genes <- msigdbr(species = "Homo sapiens", category = "H")

# Create a list of genes per gene set
hallmark_list <- split(hallmark_genes$gene_symbol, hallmark_genes$gs_name)

# Run GSEA using fgsea
gsea_results <- fgsea(pathways = hallmark_list,    # List of gene sets (Hallmark pathways)
                      stats = gene_ranks,        # Ranked gene list
                      minSize = 10,             # Minimum gene set size
                      maxSize = 5000)            # Number of permutations



gsea_results <- gsea_results[gsea_results$padj <0.05,]

# Sort the results by adjusted p-value
gsea_top <- gsea_results[order(gsea_results$padj), ][1:10, ]

# Create the bar plot and color bars by p.adj
p1 <- ggplot(gsea_top, aes(x = reorder(pathway, NES), y = NES, fill = padj)) +
  geom_bar(stat = "identity") +
  coord_flip() +  # Flip coordinates to make it horizontal
  scale_fill_gradient(low = "blue", high = "red") +  # Color scale from blue to red (adjust to your preference)
  labs(x = "Pathway", y = "Normalized Enrichment Score (NES)", fill = "Adjusted P-value") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
ggsave(p1, filename =paste0("~/RMS RNA seq/Zoes data/fgsea/H_", j, ".png"),bg = "white")

#GO-----------------------------------------------

all_gene_sets = msigdbr(species = "Homo sapiens")
all_gene_sets <- all_gene_sets[all_gene_sets$gs_subcollection == "GO:BP" | all_gene_sets$gs_subcollection == "GO:CC" |all_gene_sets$gs_subcollection == "GO:MF",  ]
# Create a list of genes per gene set
GO_list <- split(all_gene_sets$gene_symbol, all_gene_sets$gs_name)


# Run GSEA using fgsea
gsea_results <- fgsea(pathways = GO_list,    # List of gene sets (Hallmark pathways)
                      stats = gene_ranks,        # Ranked gene list
                      minSize = 10,             # Minimum gene set size
                      maxSize = 5000)            # Number of permutations



gsea_results <- gsea_results[gsea_results$padj <0.05,]

# Sort the results by adjusted p-value
gsea_top <- gsea_results[order(gsea_results$padj), ][1:10, ]

# Create the bar plot and color bars by p.adj
p2 <- ggplot(gsea_top, aes(x = reorder(pathway, NES), y = NES, fill = padj)) +
  geom_bar(stat = "identity") +
  coord_flip() +  # Flip coordinates to make it horizontal
  scale_fill_gradient(low = "blue", high = "red") +  # Color scale from blue to red (adjust to your preference)
  labs(x = "Pathway", y = "Normalized Enrichment Score (NES)", fill = "Adjusted P-value")  +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  theme_minimal()
ggsave(p2, filename =paste0("~/RMS RNA seq/Zoes data/fgsea/GO_", j, ".png"),bg = "white")

}

```





# Proliferation gene signature

```{r}


filenames <- c( "RH30_P_vs_RH_HD","RH30_P_vs_RH_C7" , "RH30_P_vs_RH_C13")

for (i in filenames){
  
  proliferation_gene_sig <- read_excel("~/RMS RNA seq/Zoes data/Proliferation_gene_signature.xlsx") 
  proliferation_gene_sig <- proliferation_gene_sig$Proliferation_gene_sig

   proliferation_gene_sig.ind <- ids2indices( proliferation_gene_sig, rownames(edgeRData$counts))
  # Run analysis:
  gst.camera <- camera.DGEList(edgeRData, design= design, index=proliferation_gene_sig.ind,contrast=my.contrasts[, i], inter.gene.cor=0.05)

 #find signif gene sets at FDR 0.05
  gst.camera$FDR <-  p.adjust(gst.camera$PValue, method = "BH")
  sig_genesets <- gst.camera[gst.camera$PValue < 0.05, ]
  print(i)
  print(sig_genesets)
  try(rownames(sig_genesets) <- paste0(i))

  #output results to table
  try(write.table(sig_genesets, "proliferation_signature_results.csv", append = TRUE,sep = ",",col.names = FALSE))

#try(read.csv(file = "proliferation_signature_results.csv"))

}

```


#Find overlapping GSEA GO terms
```{r}

# save GSEA results for each sample
############################## for GO ##################



############################### for RH30_P_vs_RH_HD ###################################

GO_pathways <- gmtPathways("C:/Users/Christina/Documents/Neuroblastoma RNAseq analysis/GSEA .gmt/c5.go.v2022.1.Hs.symbols.gmt")
GO.ind <- ids2indices( GO_pathways, rownames(edgeRData$counts))
# Run analysis:
gst.camera <- camera.DGEList(edgeRData, design= design, index=GO.ind,contrast=my.contrasts[, "RH30_P_vs_RH_HD"], inter.gene.cor=0.05)

#find signif gene sets at FDR 0.05
sig_genesets <- gst.camera[gst.camera$FDR < 0.05, ]

sig_genesets$pathway <- as.character(row.names(sig_genesets))
sig_GO_RH30_P_vs_RH_HD <- sig_genesets


###################################### for "RH30_P_vs_RH_C7"############################
# for GO
GO_pathways <- gmtPathways("C:/Users/Christina/Documents/Neuroblastoma RNAseq analysis/GSEA .gmt/c5.go.v2022.1.Hs.symbols.gmt")
GO.ind <- ids2indices( GO_pathways, rownames(edgeRData$counts))
# Run analysis:
gst.camera <- camera.DGEList(edgeRData, design= design, index=GO.ind,contrast=my.contrasts[, "RH30_P_vs_RH_C7"], inter.gene.cor=0.05)

#find signif gene sets at FDR 0.05
sig_genesets <- gst.camera[gst.camera$FDR < 0.05, ]

sig_genesets$pathway <- as.character(row.names(sig_genesets))
sig_GO_RH30_P_vs_RH_C7 <- sig_genesets

############################ for "RH30_P_vs_RH_C13" #############################################

# for GO
GO_pathways <- gmtPathways("C:/Users/Christina/Documents/Neuroblastoma RNAseq analysis/GSEA .gmt/c5.go.v2022.1.Hs.symbols.gmt")
GO.ind <- ids2indices( GO_pathways, rownames(edgeRData$counts))
# Run analysis:
gst.camera <- camera.DGEList(edgeRData, design= design, index=GO.ind,contrast=my.contrasts[, "RH30_P_vs_RH_C13"], inter.gene.cor=0.05)

#find signif gene sets at FDR 0.05
sig_genesets <- gst.camera[gst.camera$FDR < 0.05, ]

sig_genesets$pathway <- as.character(row.names(sig_genesets))
sig_GO_RH30_P_vs_RH_C13 <- sig_genesets


########## find overlapping GSEA GO terms ########################################
#put all data frames into list
df_list <- list(sig_GO_RH30_P_vs_RH_C13, sig_GO_RH30_P_vs_RH_C7, sig_GO_RH30_P_vs_RH_HD)
df_list <- list(sig_GO_RH30_P_vs_RH_C13, sig_GO_RH30_P_vs_RH_C7)
df_list <- list(sig_GO_RH30_P_vs_RH_C13, sig_GO_RH30_P_vs_RH_HD)
#merge all data frames in list
overlapping_GO_GSEA <- df_list %>% reduce(merge, by="pathway")



```


```{r}
########### FOR GO DO TOP 10  ###################
top_10_GO <- sig_genesets[1:10, ]


#dotplot
library(forcats) ## for reordering the factor

colourBy <- top_10_GO$FDR

#sig_genesets$Direction[sig_genesets$Direction == "Up"] = "upregulated"
#sig_genesets$Direction[sig_genesets$Direction == "Down"] = "downregulated"

#print dotplot
p <- ggplot(top_10_GO, aes(x = top_10_GO$NGenes, y = fct_reorder(row.names(top_10_GO), top_10_GO$NGenes))) + theme_bw(base_size = 5) +
  geom_point(aes(size = top_10_GO$NGenes, color = top_10_GO$FDR)) +
 theme(axis.text.y = element_text(size =6)) +
  
  scale_colour_gradient( low="red", high = "blue", name = "p.adjust")+
  ylab(NULL) + facet_grid(.~Direction) +  labs(size="Gene set size") +  xlab("Count") +ggtitle(i)

 print(p)
 
 #barplot
 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)+
  scale_x_discrete(expand = expansion(add = .5)) +
  coord_flip()+ facet_grid(~Direction)  + theme_bw() + labs(y="Count", x="")  +ggtitle(i)
  print(p)

}
```






Using fgsea for GSEA





```{r}

#load libraries
library(fgsea)
library(data.table)
library(ggplot2)


dge_files <- list(dge_COMB1_DMSO_PT3, dge_COMB1_DMSO_PT7, dge_COMB1_DMSO_PT10, dge_COMB2_DMSO_PT3, dge_COMB2_DMSO_PT7, dge_COMB2_DMSO_PT10, dge_GSK_DMSO_PT3, dge_GSK_DMSO_PT7, dge_GSK_DMSO_PT10, dge_RA_DMSO_PT3, dge_RA_DMSO_PT7, dge_RA_DMSO_PT10, dge_TAZ_DMSO_PT3, dge_TAZ_DMSO_PT7, dge_TAZ_DMSO_PT10)

filenames <- c("dge_COMB1_DMSO_PT3", "dge_COMB1_DMSO_PT7", "dge_COMB1_DMSO_PT10", "dge_COMB2_DMSO_PT3", "dge_COMB2_DMSO_PT7", "dge_COMB2_DMSO_PT10", "dge_GSK_DMSO_PT3", "dge_GSK_DMSO_PT7", "dge_GSK_DMSO_PT10", "dge_RA_DMSO_PT3", "dge_RA_DMSO_PT7", "dge_RA_DMSO_PT10", "dge_TAZ_DMSO_PT3", "dge_TAZ_DMSO_PT7", "dge_TAZ_DMSO_PT10")


for (i in 1:length(dge_files)) {
  
  #Save number of up/down significant genes at FDR = 0.05
  d <- filenames[i]
 
  d <- get(d)


a <- sub("....", "", filenames[i])
#set.seed(42)
#fgseaRes <- fgsea(examplePathways, exampleRanks,  maxSize=500, nproc=1)

#create vector for fgsea to rank data by lfc
  rankData <- d$table$logFC
  names(rankData) <- d$genes$genes
  head(rankData)

#load pathways
pathwaysR <- gmtPathways("~/GSEA .gmt/c2.cp.reactome.v2022.1.Hs.symbols.gmt")

#conduct analysis
set.seed(42)
fgseaRes <- fgsea(pathwaysR, 
                  rankData, 
                  minSize=15, 
                  maxSize = 500, nPermSimple = 10000,
                   nproc=1)

sig_fgsea <- fgseaRes[padj < 0.05] 
sig_fgsea <- sig_fgsea[order(sig_fgsea$padj),]
#print(sig_fgsea)

top_10_fgsea <- sig_fgsea[1:10,]
top_10_fgsea <- na.omit(top_10_fgsea)


library(ggplot2)
library(dplyr)
library(stringr)
library(forcats) ## for reordering the factor



gsea_df <- top_10_fgsea

colourBy <- gsea_df$padj

gsea_df$type = "upregulated"
gsea_df$type[gsea_df$NES < 0] = "downregulated"


p <- ggplot(gsea_df, aes(x = gsea_df$size, y = fct_reorder(gsea_df$pathway, gsea_df$size))) + 
  geom_point(aes(size = gsea_df$size, color = gsea_df$padj)) +
  theme_bw(base_size = 5) +  theme(axis.text.y = element_text(size =6))+
  scale_colour_gradient( low="red", high = "blue", name = "p.adjust")+
  ylab(NULL) + facet_grid(.~type) + labs(size="Gene set size") +  xlab("Count") + ggtitle(a)#+ ggtitle(filenames[[i]])

 print(p)
 }


```


#Venn
```{r}
#find overlapping DEGs between clones for the gene sig

  #find all overlappin DEGs----------------------------------------
#read in data
  data <- RH30_P_vs_RH_C7$table
  data <- data %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "Up-regulated",
                           logFC < -1 & FDR < 0.05 ~ "Down-regulated",
                           TRUE ~ "Unchanged")
    )
head(data) %>% 
  kable() 

data <- subset(data, Expression == 'Up-regulated' | Expression == 'Down-regulated')

#select for upregulated DEGS
all_C7 <- rownames(data)


#read in data
  data <- RH30_P_vs_RH_C13$table
  data <- data %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "Up-regulated",
                           logFC < -1 & FDR < 0.05 ~ "Down-regulated",
                           TRUE ~ "Unchanged")
    )
head(data) %>% 
  kable()

data <- subset(data, Expression == 'Up-regulated' | Expression == 'Down-regulated')

#select for upregulated DEGS
all_C13 <- rownames(data)


#read in data
  data <- RH30_P_vs_RH_HD$table
  data <- data %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "Up-regulated",
                           logFC < -1 & FDR < 0.05 ~ "Down-regulated",
                           TRUE ~ "Unchanged")
    )
head(data) %>% 
  kable()

data <- subset(data, Expression == 'Up-regulated' | Expression == 'Down-regulated')

#select for upregulated DEGS
all_HD <- rownames(data)

library(VennDiagram)
 
# Chart
venn.diagram(
  x = list(all_C7, all_C13, all_HD),
  imagetype = "png",
  category.names = c("C7 DEGs" , "C13 DEGS", "HD DEGs" ),
  filename = 'overlapping_all_DEGs_clones.png',
  output = TRUE ,
          resolution = 500 ,
  col=c("hotpink1", 'steelblue1', 'springgreen1'),
  fill = c(alpha("hotpink1",0.3), alpha('steelblue1',0.3), alpha('springgreen1',0.3)) ,
  cat.default.pos = "outer"  ,
 cat.pos = c(-1, 1, 150),
	inverted = TRUE,
 cat.dist = c(0.055, 0.055, 0.085),
	cat.cex = 1.4,	ext.line.lwd = 2,
 main.cex = 2,
 cex = 2,
 margin = 0.05,
  main = "All overlapping DEGs for clones")


  #find overlappin upregulated DEGs----------------------------------------
#read in data
  data <- RH30_P_vs_RH_C7$table
  data <- data %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "Up-regulated",
                           logFC < -1 & FDR < 0.05 ~ "Down-regulated",
                           TRUE ~ "Unchanged")
    )
head(data) %>% 
  kable() 

data <- subset(data, Expression == 'Up-regulated')

#select for upregulated DEGS
up_C7 <- rownames(data)


#read in data
  data <- RH30_P_vs_RH_C13$table
  data <- data %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "Up-regulated",
                           logFC < -1 & FDR < 0.05 ~ "Down-regulated",
                           TRUE ~ "Unchanged")
    )
head(data) %>% 
  kable()

#select for upregulated DEGS
data <- subset(data, Expression == 'Up-regulated')

#select for upregulated DEGS
up_C13 <- rownames(data)


#read in data
  data <- RH30_P_vs_RH_HD$table
  data <- data %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "Up-regulated",
                           logFC < -1 & FDR < 0.05 ~ "Down-regulated",
                           TRUE ~ "Unchanged")
    )
head(data) %>% 
  kable()

data <- subset(data, Expression == 'Up-regulated' )

#select for upregulated DEGS
up_HD <- rownames(data)


library(VennDiagram)
 
# Chart
venn.diagram(
  x = list(up_C7, up_C13, up_HD),
  imagetype = "png",
  category.names = c("C7
  up-regulated
  DEGs" , "C13
  up-regulated
  DEGs" , "HD up-regulated DEGs" ),
  filename = 'overlapping_up_DEGs_clones.png',
  output = TRUE ,
  height= 3000,
  width = 3000,
          resolution = 500 ,
  col=c("hotpink1", 'steelblue1', 'springgreen1'),
  fill = c(alpha("hotpink1",0.3), alpha('steelblue1',0.3), alpha('springgreen1',0.3)) ,
  cat.default.pos = "outer"  ,
 cat.pos = c(-2, 2, 150),
	inverted = TRUE,
 cat.dist = c(0.125, 0.125, 0.085),
	cat.cex = 1.4,
 main.cex = 2,
  margin = 0.05,
 cex = 2,
  main = "Overlapping up-regulated DEGs for clones")
         

#find overlappin downregulated DEGs----------------------------------------
#read in data
  data <- RH30_P_vs_RH_C7$table
  data <- data %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "Up-regulated",
                           logFC < -1 & FDR < 0.05 ~ "Down-regulated",
                           TRUE ~ "Unchanged")
    )
head(data) %>% 
  kable() 

data <- subset(data, Expression == 'Down-regulated')

#select for upregulated DEGS
down_C7 <- rownames(data)


#read in data
  data <- RH30_P_vs_RH_C13$table
  data <- data %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "Up-regulated",
                           logFC < -1 & FDR < 0.05 ~ "Down-regulated",
                           TRUE ~ "Unchanged")
    )
head(data) %>% 
  kable()

#select for upregulated DEGS
data <- subset(data, Expression == 'Down-regulated')

#select for upregulated DEGS
down_C13 <- rownames(data)


#read in data
  data <- RH30_P_vs_RH_HD$table
  data <- data %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "Up-regulated",
                           logFC < -1 & FDR < 0.05 ~ "Down-regulated",
                           TRUE ~ "Unchanged")
    )
head(data) %>% 
  kable()

data <- subset(data, Expression == 'Down-regulated' )

#select for upregulated DEGS
down_HD <- rownames(data)


library(VennDiagram)
 
# Chart
venn.diagram(
  x = list(down_C7, down_C13, down_HD),
  imagetype = "png",
  category.names = c("C7
  down-regulated
  DEGs" , "C13
  down-regulated
  DEGs" , "HD down-regulated DEGs"),
  filename = 'overlapping_down_DEGs_clones.png',
  output = TRUE ,
  height= 3100,
  width = 3100,
          resolution = 500 ,
  col=c("hotpink1", 'steelblue1', 'springgreen1'),
  fill = c(alpha("hotpink1",0.3), alpha('steelblue1',0.3), alpha('springgreen1',0.3)) ,
   cat.default.pos = "outer"  ,
 cat.pos = c(-2, 2, 150),
	inverted = TRUE,
 cat.dist = c(0.125, 0.125, 0.085),
	cat.cex = 1.4,
 main.cex = 2,
  margin = 0.05,
 cex = 2,
  main = "Overlapping down-regulated DEGs for clones")

```


# WGCNA
```{r}

#####Load in pacakges
library(WGCNA)
options(stringsAsFactors = FALSE)
library(dplyr)
library(ggplot2)
library(biomaRt)
library(colorspace)
library(igraph)
library(ggraph)
library(tidygraph)
library(visNetwork)
library(viridis)
library(viridisLite)

#need normalised counts in matrix with rows as samples and columns as genes
#get normalised counts
rownames(edgeRData) <- paste0("X", seq_len(nrow(edgeRData)))
edgeRData <- DGEList(counts=counts.keep, genes=rownames(counts.keep), group = treatment)
edgeRData <- calcNormFactors(edgeRData, method="TMM")
CPM <- cpm(edgeRData, log = TRUE)

#example removing samples
#logCPM <- logCPM[, !colnames(logCPM) %in% c("RH30_P_1", "RH30_P_2", "RH30_P_3", "RH30_P_4", "RH30_P_5", "RH30_P_6", "RH30_P_6")]

#transpose data so rows are samples and columns are genes
input_mat = t(CPM)
dim(input_mat) #check dimensions of data


datET <- input_mat #set the expression matrix with colums as genes and rows as samples
#setting important options for WGCNA
options(stringsAsFactors = FALSE)
##Checking for missing values
gsg <- goodSamplesGenes(datET, verbose = 3)
gsg$allOK

#Take a look to see how your samples compare: 
sampleTree = hclust(dist(datET), method ="average")
par(cex = 0.6); 
par(mar = c(0,4,2,0))
sizeGrWindow(16,9)
plot(sampleTree, main = "outlier detection", sub="",xlab="", cex.lab=1.5, cex.axis = 1.5, cex.main = 2)


# create trait data 
datTraits = read_excel("~/RMS RNA seq/Zoes data/WGCNA on all samples (incl HD)/zoes_trait_data.xlsx")
 datTraits <- as.matrix(datTraits)
row.names(datTraits) = row.names(input_mat)


# Re-cluster samples
sampleTree2 = hclust(dist(datET), method = "average")
# Convert traits to a color representation: white means low, red means high, grey means missing entry
traitColors = numbers2colors(datTraits, signed = FALSE);
# Plot the sample dendrogram and the colors underneath.
plotDendroAndColors(sampleTree2, traitColors,
                    groupLabels = names(datTraits),
                    main = "Sample dendrogram and trait heatmap")





#picking soft threshold
sft <- pickSoftThreshold(datET,
  dataIsExpr = TRUE,
  corFnc = cor,
  networkType = "signed"
)

sft_df <- data.frame(sft$fitIndices) %>%
  dplyr::mutate(model_fit = -sign(slope) * SFT.R.sq)

ggplot(sft_df, aes(x = Power, y = model_fit, label = Power)) +
  # Plot the points
  geom_point() +
  # We'll put the Power labels slightly above the data points
  geom_text(nudge_y = 0.1) +
  # We will plot what WGCNA recommends as an R^2 cutoff
  geom_hline(yintercept = 0.80, col = "red") +
  # Just in case our values are low, we want to make sure we can still see the 0.80 level
  ylim(c(min(sft_df$model_fit), 1.05)) +
  # We can add more sensible labels for our axis
  xlab("Soft Threshold (power)") +
  ylab("Scale Free Topology Model Fit, signed R^2") +
  ggtitle("Scale independence") +
  # This adds some nicer aesthetics to our plot
  theme_classic()

#Assess scale free topology of data
k <- softConnectivity(datE=datET, corFnc="bicor", type="signed")
##assessing distribution of k, and looking at scale-free topology
hist(k)
scaleFreePlot(k)

#Select or set soft threshold power 
#beta = 12

#
adjacency <- adjacency(datET, corFnc ="cor",corOptions = list(use = 'p', method ='kendall'), type="signed",power=12)
TOM <- TOMsimilarity(adjacency, TOMType = "signed", TOMDenom = "min") 
dissTOM <- 1 - TOM
geneTree <- hclust(as.dist(dissTOM), method = "average")
minModuleSize = 40

dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM, deepSplit = 4, pamRespectsDendro = FALSE, minClusterSize = minModuleSize) 
moduleColors <- labels2colors(dynamicMods)
length(unique(moduleColors)) 
plotDendroAndColors(geneTree, moduleColors, dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05, main = "Gene dendrogram and module colors")



#Calculating the module dissimilarity eigengenes
MEDiss = 1-cor(MEs)

#Clustering the eigengenes modules
METree = hclust(as.dist(MEDiss), method = "average")
#Plotting the result
#sizeGrWindow(7, 6)
plot(METree, main = "Clustering of module eigengenes",
     xlab = "", sub = "")


#Grouping the clusters from a cut-off
MEDissThres = 0.25
#Plotting a cut-off line
abline(h=MEDissThres, col = "red")


# Calculate eigengenes
datME <- moduleEigengenes(datET, colors = moduleColors)$eigengenes
MEs <- moduleEigengenes(datET, colors = moduleColors)$eigengenes

# adj <- (1 + cor(datME, use="p", method = "spearman"))/2
# diag(adj) <- 0


# code to merge modules
# Calculate dissimilarity of module eigengenes
MEDiss = 1-cor(MEs);
# Cluster module eigengenes
METree = hclust(as.dist(MEDiss), method = "average");
# Plot the result
sizeGrWindow(7, 6)
plot(METree, main = "Clustering of module eigengenes",
xlab = "", sub = "")
MEDissThres = 0.25
# Plot the cut line into the dendrogram
abline(h=MEDissThres, col = "red")
# Call an automatic merging function
merge = mergeCloseModules(datET, moduleColors, cutHeight = MEDissThres, verbose = 3)
# The merged module colors
mergedColors = merge$colors;
# Eigengenes of the new ed modules:
mergedMEs = merge$newMEs;

#plot merged vs non-merged modules
sizeGrWindow(12, 9)
#pdf(file = "Plots/geneDendro-3.pdf", wi = 9, he = 6)
plotDendroAndColors(geneTree, cbind(moduleColors, mergedColors),
c("Dynamic Tree Cut", "Merged dynamic"),
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05)
#dev.off()


# Rename to moduleColors
moduleColors = mergedColors
# Construct numerical labels corresponding to the colors
colorOrder = c("grey", standardColors(50));
moduleLabels = match(moduleColors, colorOrder)-1;
MEs = mergedMEs;


vis <- exportNetworkToVisANT(adj, weighted = TRUE, threshold = 0.75)


g <- graph_from_data_frame(vis, directed = FALSE)
g <- as_tbl_graph(g)
g <- g %>% mutate(centrality = centrality_degree())

nodes <- g %>% activate(nodes) %>% as_tibble()
edges <- g %>% activate(edges) %>% as_tibble()

nodes_vis <- data.frame(id=nodes$name, label=nodes$name, value = nodes$centrality, font.size = 10)
edges_vis <- data.frame(from=vis$from, to=vis$to, length=vis$weight)



net <- visNetwork(nodes_vis, edges_vis, width = "100%") %>%
  visPhysics(stabilization = F, timestep = .1,barnesHut = list(
    centralGravity = 0.5, 
    gravitationalConstant = -10000, 
    springLength = 100, 
    springConstant = 0.05, 
    avoidOverlap = 0)) %>%
  visNodes(shape = "dot",
           color = list(
             background = "#0085AF",
             border = "#013848",
             highlight = "#FF8000"
           ),
           shadow = list(enabled = FALSE)
  ) %>%
  visEdges(
    smooth = T,
    shadow = FALSE,
    color = list(highlight = "#C62F4B")
  ) %>%
  visOptions(highlightNearest = list(enabled = T, degree = 1, hover = T)) %>% 
  visLayout(randomSeed = 11)
net



#trait correlation
# In this analysis we would like to identify modules that are significantly associated with the measured clinical traits.Since we already have a summary profile (eigengene) for each module, we simply correlate eigengenes with external traits and look for the most significant associations:



# 
# # Form a data frame analogous to expression data that will hold the clinical traits.
# femaleSamples = rownames(datExpr0);
# traitRows = match(femaleSamples, allTraits$Mice);
# datTraits = allTraits[traitRows, -1];
# rownames(datTraits) = allTraits[traitRows, 1];
# collectGarbage();


# Define numbers of genes and samples
nGenes = ncol(datET);
nSamples = nrow(datET);
# Recalculate MEs with color labels
MEs0 = moduleEigengenes(datET, moduleColors)$eigengenes
MEs = orderMEs(MEs0)
moduleTraitCor = cor(MEs, datTraits, use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);



#combine moduletraitcor and moduletraitpval
moduleTraitCorDF <- as.data.frame(moduleTraitCor) 
moduleTraitPvalueDF <- as.data.frame(moduleTraitPvalue) 
moduleTraitCorAndPvalue <- merge(moduleTraitCorDF, moduleTraitPvalueDF,
                          by = 'row.names', all = TRUE)
#get only signficant module traits
SignficantModuleTraitCor <- subset(moduleTraitCorAndPvalue, Resistant.y < 0.05)



#Displaying correlations and its p-values
textMatrix =  paste(signif(moduleTraitCor, 2), "\n(",
                    signif(moduleTraitPvalue, 1), ")", sep = "")
dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(6, 8.5, 3, 3))

labeledHeatmap(Matrix = moduleTraitCor,
               xLabels = names(datTraits),
               yLabels = names(MEs),
               ySymbols = names(MEs),
               colorLabels = FALSE,
               colors = blueWhiteRed(50),
               textMatrix = textMatrix,
               setStdMargins = FALSE,
               cex.text = 0.5,
               zlim = c(-1,1),
               main = paste("Module-trait relationships"))

# Define variable weight containing the resistant column of datTrait
resistance = as.data.frame(datTraits);
names(resistance) = "resistance"
# names (colors) of the modules
modNames = substring(names(MEs), 3)
geneModuleMembership = as.data.frame(cor(datET, MEs, use = "p"));
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));

names(geneModuleMembership) = paste("MM", modNames, sep="");
names(MMPvalue) = paste("p.MM", modNames, sep="");





# binarytraitresults <- standardScreeningBinaryTrait(
#      datET, traitData,
#      corFnc = cor, corOptions = list(use = 'p'),
#      kruskalTest = FALSE, qValues = FALSE,
#      var.equal=FALSE, na.action="na.exclude",
#      getAreaUnderROC = TRUE)
# 



geneTraitSignificance = as.data.frame(cor(datET, datTraits, use = "p"));
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));
names(geneTraitSignificance) = paste("GS.", names(datTraits), sep="");
names(GSPvalue) = paste("p.GS.", names(datTraits), sep="");


# #Intramodular analysis: identifying genes with high GS and MM
# Using the GS and MM measures, we can identify genes that have a high significance for weight as well as high module
# membership in interesting modules. As an example, we look at the brown module that has the highest association
# with weight. We plot a scatterplot of Gene Significance vs. Module Membership in the brown module:



module = "orange"
column = match(module, modNames);
moduleGenes = moduleColors==module;
sizeGrWindow(7, 7);
par(mfrow = c(1,1));
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
abs(geneTraitSignificance[moduleGenes, 1]),
xlab = paste("Module Membership in", module, "module"),
ylab = "Gene significance for resistance",
main = paste("Module membership vs. gene significance\n"),
cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = "black")


module = "lightcyan"
column = match(module, modNames);
moduleGenes = moduleColors==module;
sizeGrWindow(7, 7);
par(mfrow = c(1,1));
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
abs(geneTraitSignificance[moduleGenes, 1]),
xlab = paste("Module Membership in", module, "module"),
ylab = "Gene significance for resistance",
main = paste("Module membership vs. gene significance\n"),
cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = "sienna3")


module_size <- colnames(datET)[moduleColors=="orange"]

#summary output 

module_list <- colnames(datME)
module_list <- sub("ME", "", module_list)

#get module sizes
list_module_sizes <- c()
for (i in module_list ) {
 
 module_size <- colnames(datET)[moduleColors==i] %>% length()
# print(i)
print(module_size)

list_module_sizes[i] <- module_size
}

#put module size info as df
list_module_sizes <- as.data.frame(list_module_sizes)
list_module_sizes



#get sum of module sizes
average_module_size <- sum(list_module_sizes$list_module_sizes) 



#get number of modules
number_of_modules <- unique(moduleColors) %>% length()

#find average module size
average_module_size/number_of_modules

# Create the starting data frame
geneInfo0 = data.frame(genes = colnames(datET),
moduleColor = moduleColors,
geneTraitSignificance,
GSPvalue)
# Order modules by their significance for weight
modOrder = order(-abs(cor(MEs, resistance, use = "p")));
# Add module membership information in the chosen order
for (mod in 1:ncol(geneModuleMembership))
{
oldNames = names(geneInfo0)
geneInfo0 = data.frame(geneInfo0, geneModuleMembership[, modOrder[mod]],
MMPvalue[, modOrder[mod]]);
names(geneInfo0) = c(oldNames, paste("MM.", modNames[modOrder[mod]], sep=""),
paste("p.MM.", modNames[modOrder[mod]], sep=""))
}
# Order the genes in the geneInfo variable first by module color, then by geneTraitSignificance
geneOrder = order(geneInfo0$moduleColor, -abs(geneInfo0$GS.));
geneInfo = geneInfo0[geneOrder, ]

#save geneinfo to file
write.csv(geneInfo, file = "geneInfo.csv")

#view geneinfo file
 fix(geneInfo)

 
 
 #getting genes from interesting modules

  module_genes <- colnames(datET)[moduleColors=="orange"] 
  write.csv(module_genes, "module_genes_orange.csv")

 #get gene info for module selected
 gene_info_orange <- geneInfo[geneInfo$genes %in% module_genes,] 

 genes_orange <- row.names(gene_info_orange)
 
#select genes with GS. >0.7 and MM > 0.7 as hub genes
hub_genes <-  subset(gene_info_orange, GS. >= 0.7 & MM.orange >= 0.7)
 

#find overlapping DEGS and hub genes

 test <- hub_genes[hub_genes$genes %in% DEG_list_RH30_P_vs_RH_HD$geneID,] 
  test <- hub_genes[hub_genes$genes %in% DEG_list_RH30_P_vs_RH_C13$geneID,] 
  test <- hub_genes[hub_genes$genes %in% DEG_list_RH30_P_vs_RH_C7$geneID,] 
 
 hub_gene_list <- hub_genes$genes




 #visualising networks
 
 
 
#  
#  # Calculate topological overlap anew: this could be done more efficiently by saving the TOM
# # calculated during module detection, but let us do it again here.
# dissTOM = 1-TOMsimilarityFromExpr(datET, power = 6);
# # Transform dissTOM with a power to make moderately strong connections more visible in the heatmap
# plotTOM = dissTOM^7;
# # Set diagonal to NA for a nicer plot
# diag(plotTOM) = NA;
# # Call the plot function
# sizeGrWindow(9,9)
# TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")
 
 # Recalculate module eigengenes
MEs = moduleEigengenes(datET, moduleColors)$eigengenes


# get clinical trait data

names(resistance) = "resistance"
# Add the weight to existing module eigengenes

 MET = orderMEs(cbind(MEs, resistance))
# Plot the relationships among the eigengenes and the trait
sizeGrWindow(5,7.5);
par(cex = 0.9)
plotEigengeneNetworks(MET, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8, xLabelsAngle
= 90)
 
 # Plot the dendrogram
sizeGrWindow(6,6);
par(cex = 1.0)
plotEigengeneNetworks(MET, "Eigengene dendrogram", marDendro = c(0,4,2,0),
plotHeatmaps = FALSE)
# Plot the heatmap matrix (note: this plot will overwrite the dendrogram plot)
par(cex = 1.0)
plotEigengeneNetworks(MET, "Eigengene adjacency heatmap", marHeatmap = c(3,4,2,2),
plotDendrograms = FALSE, xLabelsAngle = 90)

 
#  
# #with cortraits
# corTraits <- as.data.frame(cor(geneModuleMembership, geneTraitSignificance, use = "p", method="spearman")) #Finding correlation between the kME (module membership) and gene significance score
# 
# colour_df <- as.data.frame(numbers2colors(corTraits,signed = TRUE, centered = FALSE)) #using a WGCNA function to convert the numbers (cor coefficient) to a colour value (hex code) 
# colnames(colour_df) <- c("color") #labelling the name of the column on the dataframe
# colour_df$label <- colnames(datME) #adding in datME label
# 
# nodes_vis <- data.frame(id=nodes$name, label=nodes$name, value = nodes$centrality, font.size = 10) # setting up node vis - you have the script for this
# edges_vis <- data.frame(from=vis$from, to=vis$to, length=vis$weight) # ditto as above
# 
# nodes_vis <- left_join(nodes_vis, colour_df, by = "label") #using the left_join() function to merge the nodes_vis and the colour_df, visnetwork will then recognise the colour value
#  



#with module trait 
moduleTraitCor


test_colour_df <- as.data.frame(numbers2colors(moduleTraitCor,signed = TRUE, centered = FALSE))


colnames(test_colour_df) <- c("color") #labelling the name of the column on the dataframe
test_colour_df$label <- colnames(datME) #adding in datME label

nodes_vis <- data.frame(id=nodes$name, label=nodes$name, value = nodes$centrality, font.size = 10) # setting up node vis - you have the script for this
edges_vis <- data.frame(from=vis$from, to=vis$to, length=vis$weight) # ditto as above

nodes_vis <- left_join(nodes_vis, test_colour_df, by = "label") #using the left_join() function to merge the nodes_vis and the colour_df, visnetwork will then recognise the colour value
 
 
net <- visNetwork(nodes_vis, edges_vis, width = "100%") %>%
  visPhysics(stabilization = F, timestep = .1,barnesHut = list(
    centralGravity = 0.5, 
    gravitationalConstant = -10000, 
    springLength = 100, 
    springConstant = 0.05, 
    avoidOverlap = 0)) %>%
  visNodes(shape = "dot"
           , color = nodes_vis$color, 
           shadow = list(enabled = FALSE)
  ) %>%
  visEdges(
    smooth = T,
    shadow = FALSE,
    color = list(highlight = "#C62F4B")
  ) %>%
  visOptions(highlightNearest = list(enabled = T, degree = 1, hover = T)) %>% 
  visLayout(randomSeed = 11)
net

 
```


#GSEA on hub genes
```{r}

#use clusterprofiler to do GSEA on module associated with resistance

library(clusterProfiler)
library(org.Hs.eg.db)
library(ReactomePA)


hub_genes <- read.csv("hub_genes.csv") 
geneList <- hub_genes$x

#hub genes
ego2 <- enrichGO(gene         = geneList,
                OrgDb         = org.Hs.eg.db,
                keyType       = 'SYMBOL',
                ont           = "ALL",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.05,
                qvalueCutoff  = 0.05)
head(ego2, 3)   

## Dotplot 
dotplot(ego2, showCategory=50)


symbols <- hub_genes
symbols <- mapIds(org.Hs.eg.db, symbols, 'ENTREZID', 'SYMBOL')
symbols <- as.data.frame(symbols)
symbols <- symbols$symbols

kk <- enrichKEGG(gene         = symbols,
                 organism     = 'hsa',
                 keyType       = 'ncbi-geneid',
                 pvalueCutoff = 0.05)
head(kk)
dotplot(kk, showCategory=50)

x <- enrichPathway(gene=symbols, pvalueCutoff = 0.05, readable=TRUE)
head(x)
dotplot(x, showCategory=50)


```


#WGCNA CLONES ONLY

```{r}
#####Load in pacakges
library(WGCNA)
options(stringsAsFactors = FALSE)
library(dplyr)
library(ggplot2)
library(fgsea)
library(biomaRt)
library(colorspace)
library(igraph)
library(ggraph)
library(tidygraph)
library(visNetwork)
library(viridis)
library(viridisLite)

# #need normalised counts in matrix with rows as samples and columns as genes
# #get normalised counts
# rownames(edgeRData) <- paste0("X", seq_len(nrow(edgeRData)))
# edgeRData <- DGEList(counts=counts.keep, genes=rownames(counts.keep), group = treatment)
# edgeRData <- calcNormFactors(edgeRData, method="TMM")
# CPM <- cpm(edgeRData, log = TRUE)
# 
# #remove HD samples
# CPM <- CPM[, !colnames(CPM) %in% c("RH_HD_1", "RH_HD_2", "RH_HD_3", "RH_HD_4", "RH_HD_5", "RH_HD_6")]

#remove embryonal samples from raw data
rawData <- rawData[,!grepl("HD", colnames(rawData))]


#we have to refilter to remove genes not expressed in any Clone samples

#filter to remove lowly expressed genes
# Obtain CPMs
myCPM <- cpm(rawData)
# Have a look at the output
head(myCPM)
# Let's have a look and see whether a threshold of 0.5 does indeed correspond to a count of about 10-15
# We will look at the first sample
plot(myCPM[,1],rawData[,1])
# Let us limit the x and y-axis so we can actually look to see what is happening at the smaller counts
plot(myCPM[,1],rawData[,1],ylim=c(0,50),xlim=c(0,3))
# Add a vertical line at 0.5 CPM
abline(v=0.25)

# Which values in myCPM are greater than 0.3?
thresh <- myCPM > 0.25
# This produces a logical matrix with TRUEs and FALSEs
head(thresh)

# Summary of how many TRUEs there are in each row
table(rowSums(thresh))

# we would like to keep genes that have at least 6 TRUES in each row of thresh
keep <- rowSums(thresh) >= 6
# Subset the rows of countdata to keep the more highly expressed genes
counts.keep <- rawData[keep,]
summary(keep)
dim(counts.keep)


#compare number of raw counts to number of filtered counts
dim(rawData)
dim(counts.keep)

#get treatments just for rh4 remove rh30
#for sample treatment

sampleID <- colnames(rawData)
treatment <- c()
for (i in 1:ncol(rawData)) {
  samplename <- sampleID[i]
  group <- if(grepl("RH30_P", samplename)) {print("RH30_P")
  }
  else if(grepl("RH_C7", samplename)) {print("RH_C7")
  } else if(grepl("RH_C13", samplename)) {print("RH_C13")
  }
  treatment[i] <- group
}


edgeRData <- DGEList(counts=counts.keep, genes=rownames(counts.keep), group = treatment)
edgeRData



#need normalised counts in matrix with rows as samples and columns as genes
#get normalised counts
rownames(edgeRData) <- paste0("X", seq_len(nrow(edgeRData)))
edgeRData <- DGEList(counts=counts.keep, genes=rownames(counts.keep), group = treatment)
edgeRData <- calcNormFactors(edgeRData, method="TMM")
CPM <- cpm(edgeRData, log = TRUE)


#transpose data so rows are samples and columns are genes
input_mat = t(CPM)
dim(input_mat) #check dimensions of data


datET <- input_mat #set the expression matrix with colums as genes and rows as samples
#setting important options for WGCNA
options(stringsAsFactors = FALSE)
##Checking for missing values
gsg <- goodSamplesGenes(datET, verbose = 3)
gsg$allOK

#Take a look to see how your samples compare: 
sampleTree = hclust(dist(datET), method ="average")
par(cex = 0.6); 
par(mar = c(0,4,2,0))
sizeGrWindow(16,9)
plot(sampleTree, main = "outlier detection", sub="",xlab="", cex.lab=1.5, cex.axis = 1.5, cex.main = 2)


# create trait data 
datTraits = read_excel("~/RMS RNA seq/Zoes data/CLONES ONLY ANALYSIS/trait_data_clones_only.xlsx")
 datTraits <- as.matrix(datTraits)
row.names(datTraits) = row.names(input_mat)


#check trait data matches input matrix 
# Re-cluster samples
sampleTree2 = hclust(dist(datET), method = "average")
# Convert traits to a color representation: white means low, red means high, grey means missing entry
traitColors = numbers2colors(datTraits, signed = FALSE);
# Plot the sample dendrogram and the colors underneath.
plotDendroAndColors(sampleTree2, traitColors,
                    groupLabels = names(datTraits),
                    main = "Sample dendrogram and trait heatmap")





#picking soft threshold
sft <- pickSoftThreshold(datET,
  dataIsExpr = TRUE,
  corFnc = cor,
  networkType = "signed"
)

sft_df <- data.frame(sft$fitIndices) %>%
  dplyr::mutate(model_fit = -sign(slope) * SFT.R.sq)

ggplot(sft_df, aes(x = Power, y = model_fit, label = Power)) +
  # Plot the points
  geom_point() +
  # We'll put the Power labels slightly above the data points
  geom_text(nudge_y = 0.1) +
  # We will plot what WGCNA recommends as an R^2 cutoff
  geom_hline(yintercept = 0.80, col = "red") +
  # Just in case our values are low, we want to make sure we can still see the 0.80 level
  ylim(c(min(sft_df$model_fit), 1.05)) +
  # We can add more sensible labels for our axis
  xlab("Soft Threshold (power)") +
  ylab("Scale Free Topology Model Fit, signed R^2") +
  ggtitle("Scale independence") +
  # This adds some nicer aesthetics to our plot
  theme_classic()

#Assess scale free topology of data
k <- softConnectivity(datE=datET, corFnc="bicor", type="signed")
##assessing distribution of k, and looking at scale-free topology
hist(k)
#scaleFreePlot(k)

#Select or set soft threshold power 
#beta = 10

#
adjacency <- adjacency(datET, corFnc ="cor",corOptions = list(use = 'p', method ='kendall'), type="signed",power=10)
TOM <- TOMsimilarity(adjacency, TOMType = "signed", TOMDenom = "min") 
dissTOM <- 1 - TOM
geneTree <- hclust(as.dist(dissTOM), method = "average")
minModuleSize = 40

dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM, deepSplit = 4, pamRespectsDendro = FALSE, minClusterSize = minModuleSize) 
moduleColors <- labels2colors(dynamicMods)
length(unique(moduleColors)) 
plotDendroAndColors(geneTree, moduleColors, dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05, main = "Gene dendrogram and module colors")


# Calculate eigengenes
datME <- moduleEigengenes(datET, colors = moduleColors)$eigengenes
MEs = datME


# Save module colors and labels for use in subsequent parts
#save(MEs, moduleLabels, moduleColors, geneTree, file = "Zoes_data_WGCNA.RData")

# code to merge modules
# Calculate dissimilarity of module eigengenes
MEDiss = 1-cor(MEs);
# Cluster module eigengenes
METree = hclust(as.dist(MEDiss), method = "average");
# Plot the result
sizeGrWindow(7, 6)
plot(METree, main = "Clustering of module eigengenes",
xlab = "", sub = "")
MEDissThres = 0.25
# Plot the cut line into the dendrogram
abline(h=MEDissThres, col = "red")
# Call an automatic merging function
merge = mergeCloseModules(datET, moduleColors, cutHeight = MEDissThres, verbose = 3)
# The merged module colors
mergedColors = merge$colors;
# Eigengenes of the new ed modules:
mergedMEs = merge$newMEs;

#plot merged vs non-merged modules
sizeGrWindow(12, 9)
#pdf(file = "Plots/geneDendro-3.pdf", wi = 9, he = 6)
plotDendroAndColors(geneTree, cbind(moduleColors, mergedColors),
c("Dynamic Tree Cut", "Merged dynamic"),
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05)
#dev.off()


# Rename to moduleColors
moduleColors = mergedColors
# Construct numerical labels corresponding to the colors
colorOrder = c("grey", standardColors(50));
moduleLabels = match(moduleColors, colorOrder)-1;
MEs = mergedMEs;


vis <- exportNetworkToVisANT(adj, weighted = TRUE, threshold = 0.75)


g <- graph_from_data_frame(vis, directed = FALSE)
g <- as_tbl_graph(g)
g <- g %>% mutate(centrality = centrality_degree())

nodes <- g %>% activate(nodes) %>% as_tibble()
edges <- g %>% activate(edges) %>% as_tibble()

nodes_vis <- data.frame(id=nodes$name, label=nodes$name, value = nodes$centrality, font.size = 10)
edges_vis <- data.frame(from=vis$from, to=vis$to, length=vis$weight)



net <- visNetwork(nodes_vis, edges_vis, width = "100%") %>%
  visPhysics(stabilization = F, timestep = .1,barnesHut = list(
    centralGravity = 0.5, 
    gravitationalConstant = -10000, 
    springLength = 100, 
    springConstant = 0.05, 
    avoidOverlap = 0)) %>%
  visNodes(shape = "dot",
           color = list(
             background = "#0085AF",
             border = "#013848",
             highlight = "#FF8000"
           ),
           shadow = list(enabled = FALSE)
  ) %>%
  visEdges(
    smooth = T,
    shadow = FALSE,
    color = list(highlight = "#C62F4B")
  ) %>%
  visOptions(highlightNearest = list(enabled = T, degree = 1, hover = T)) %>% 
  visLayout(randomSeed = 11)
net



#trait correlation
# In this analysis we would like to identify modules that are significantly associated with the measured clinical traits.Since we already have a summary profile (eigengene) for each module, we simply correlate eigengenes with external traits and look for the most significant associations:


# # Form a data frame analogous to expression data that will hold the clinical traits.
# femaleSamples = rownames(datExpr0);
# traitRows = match(femaleSamples, allTraits$Mice);
# datTraits = allTraits[traitRows, -1];
# rownames(datTraits) = allTraits[traitRows, 1];
# collectGarbage();

# Define numbers of genes and samples
nGenes = ncol(datET);
nSamples = nrow(datET);
# Recalculate MEs with color labels
MEs0 = moduleEigengenes(datET, moduleColors)$eigengenes
MEs = orderMEs(MEs0)
moduleTraitCor = cor(MEs, datTraits, use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);

#combine moduletraitcor and moduletraitpval
moduleTraitCorDF <- as.data.frame(moduleTraitCor) 
moduleTraitPvalueDF <- as.data.frame(moduleTraitPvalue) 
moduleTraitCorAndPvalue <- merge(moduleTraitCorDF, moduleTraitPvalueDF,
                          by = 'row.names', all = TRUE)
#get only signficant module traits
SignficantModuleTraitCor <- subset(moduleTraitCorAndPvalue, Resistant.y < 0.05)

#Print correlation heatmap between modules and traits
textMatrix= paste(signif(moduleTraitCor, 2), "\n(",
                    signif(moduleTraitPvalue, 1), ")", sep = "")
dim(textMatrix)= dim(moduleTraitCor)
par(mar= c(5,10, 4, 3))
#display the corelation values with a heatmap plot
#INCLUE THE NEXT LINE TO SAVE TO FILE
#pdf(file="heatmap.pdf")
labeledHeatmap(Matrix= moduleTraitCor,
            xLabels= colnames(datTraits),
            yLabels= names(MEs),
            ySymbols= names(MEs),
            colorLabels= FALSE,
            colors= blueWhiteRed(50),
            textMatrix= textMatrix,
            setStdMargins= FALSE,
            cex.text= 0.5,
            zlim= c(-1,1),
            main= paste("Module-trait relationships"))

# Define variable weight containing the resistant column of datTrait
resistance = as.data.frame(datTraits);
names(resistance) = "resistance"
# names (colors) of the modules
modNames = substring(names(MEs), 3)
geneModuleMembership = as.data.frame(cor(datET, MEs, use = "p"));
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));

names(geneModuleMembership) = paste("MM", modNames, sep="");
names(MMPvalue) = paste("p.MM", modNames, sep="");

geneTraitSignificance = as.data.frame(cor(datET, datTraits, use = "p"));
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));
names(geneTraitSignificance) = paste("GS.", names(datTraits), sep="");
names(GSPvalue) = paste("p.GS.", names(datTraits), sep="");



#see what genes are contributing to the modules on a heatmap

neg_mods <- datME[, c("MEdarkred", "MEdarkorange")]
pos_mods <- datME[, c("MEdarkgreen", "MEdarkgrey")]

t(datME) %>% heatmap(col = brewer.pal(11,"RdBu"))
#plot heatmap of genes of interest
h1<-Heatmap(as.matrix(t(datME)),cluster_rows=TRUE,cluster_columns=FALSE,use_raster=FALSE,row_dend_reorder = TRUE,show_row_dend = TRUE,show_row_names=TRUE,show_column_names = TRUE, heatmap_legend_param = list(title = "Z-score"))
my_palette <- colorRampPalette(c("blue", "black", "red"))





# #Intramodular analysis: identifying genes with high GS and MM
# Using the GS and MM measures, we can identify genes that have a high significance for weight as well as high module
# membership in interesting modules. As an example, we look at the darkgreen/darkgrey module that has the highest association
# with weight. We plot a scatterplot of Gene Significance vs. Module Membership in the darkgreen/darkgrey module:



modules <- c("darkorange", "darkred")

for (i in modules){
  
module = i
column = match(module, modNames);
moduleGenes = moduleColors==module;
sizeGrWindow(7, 7);
par(mfrow = c(1,1));
png(filename=paste0("C:/Users/Christina/Documents/RMS RNA seq/Zoes data/CLONES ONLY ANALYSIS/Neg modules/MM_GS/MM_GS_", i,".png"), width=600, bg="white");
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
abs(geneTraitSignificance[moduleGenes, 1]),
xlab = paste("Module Membership in", module, "module"),
ylab = "Gene significance for resistance",
main = paste("Module membership vs. gene significance\n"),
cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = i)
}


module = "darkred"
column = match(module, modNames);
moduleGenes = moduleColors==module;
sizeGrWindow(7, 7);
par(mfrow = c(1,1));
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
abs(geneTraitSignificance[moduleGenes, 1]),
xlab = paste("Module Membership in", module, "module"),
ylab = "Gene significance for resistance",
main = paste("Module membership vs. gene significance\n"),
cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = "darkred")


module_size <- colnames(datET)[moduleColors=="darkgrey"]


#summary output of module size------------------------
module_list <- colnames(MEs)
module_list <- sub("ME", "", module_list)

#get module sizes
list_module_sizes <- c()
for (i in module_list ) {
 
 module_size <- colnames(datME)[moduleColors==i] %>% length()
# print(i)
print(module_size)

list_module_sizes[i] <- module_size
}

#put module size info as df
list_module_sizes <- as.data.frame(list_module_sizes)
list_module_sizes
#put module size info as df
list_module_sizes <- as.data.frame(list_module_sizes)
list_module_sizes
write.csv(list_module_sizes, "intrinsic_module_stats.csv")


#get sum of module sizes
average_module_size <- sum(list_module_sizes$list_module_sizes) 
#get number of modules
number_of_modules <- unique(moduleColors) %>% length()

#find average module size
average_module_size/number_of_modules


list_module_sizes <- read.csv("intrinsic_module_stats.csv")

# Create the data for the chart 
name   <- module_list
value=list_module_sizes$list_module_sizes 
GFG_Name <- c("lightgreen", "darkgreen" , "darkgrey"  , "darkred"   , "yellow"   ,  "skyblue",    "darkorange", "turquoise" )
X <- data.frame(name,GFG_Name,value) 
 
  
  
# Plot the bar chart  

par(mar=c(11,6,4,1)+.1)
p <- barplot(value,las=2,names.arg=GFG_Name,ylab="Size", ylim = c(0,5000)
        ,col=c("lightgreen", "darkgreen" , "darkgrey" ,  "darkred",    "yellow"   ,  "skyblue" ,   "darkorange", "turquoise" ), 
main="Module sizes of the intrinsic resistance data") 
title(xlab="Module name", line=6, cex.lab=1)
jpeg(file="module_size_intrinsic_zoes.jpeg")



#get sum of module sizes
average_module_size <- sum(list_module_sizes$list_module_sizes) 



#get number of modules
number_of_modules <- unique(moduleColors) %>% length()

#find average module size
average_module_size/number_of_modules

# Create the starting data frame
geneInfo0 = data.frame(genes = colnames(datET),
moduleColor = moduleColors,
geneTraitSignificance,
GSPvalue)
# Order modules by their significance for weight
modOrder = order(-abs(cor(MEs, resistance, use = "p")));
# Add module membership information in the chosen order
for (mod in 1:ncol(geneModuleMembership))
{
oldNames = names(geneInfo0)
geneInfo0 = data.frame(geneInfo0, geneModuleMembership[, modOrder[mod]],
MMPvalue[, modOrder[mod]]);
names(geneInfo0) = c(oldNames, paste("MM.", modNames[modOrder[mod]], sep=""),
paste("p.MM.", modNames[modOrder[mod]], sep=""))
}
# Order the genes in the geneInfo variable first by module color, then by geneTraitSignificance
geneOrder = order(geneInfo0$moduleColor, -abs(geneInfo0$GS.));
geneInfo = geneInfo0[geneOrder, ]

#save geneinfo to file
write.csv(geneInfo, file = "geneInfo_clones.csv")

#view geneinfo file
 fix(geneInfo)

 geneInfo <- read.csv("geneInfo_clones.csv")
 
 
 #getting genes from interesting modules--------------------------------------------------
 #for dark grey


 module_genes <- subset(geneInfo, moduleColor == 'darkgrey')
 
module_genes <-  module_genes$genes
 
 # write.csv(module_genes, "module_genes_darkgrey.csv")

 #get gene info for module selected
 gene_info <- geneInfo[geneInfo$genes %in% module_genes,] 

 genes <- gene_info$genes
 
#select genes with GS. >0.6 and MM > 0.6 as hub genes
hub_genes <-  subset(gene_info, GS. >0.6 & MM.darkgrey >0.6)
 
write.csv(hub_genes, file = "hub_genes_darkgrey_0.6.csv")



#for dark green


 module_genes <- subset(geneInfo, moduleColor == 'darkgreen')
 
module_genes <-  module_genes$genes
 
 # write.csv(module_genes, "module_genes_darkgreen.csv")

 #get gene info for module selected
 gene_info <- geneInfo[geneInfo$genes %in% module_genes,] 

 genes <- gene_info$genes
 
#select genes with GS. >0.6 and MM > 0.6 as hub genes
hub_genes <-  subset(gene_info, GS. >0.8 & MM.darkgreen >0.8)
 
write.csv(hub_genes, file = "hub_genes_darkgreen_0.8.csv")

#for light green


 module_genes <- subset(geneInfo, moduleColor == 'lightgreen')
 
module_genes <-  module_genes$genes
 
 # write.csv(module_genes, "module_genes_darkgreen.csv")

 #get gene info for module selected
 gene_info <- geneInfo[geneInfo$genes %in% module_genes,] 

 genes <- gene_info$genes
 
#select genes with GS. >0.6 and MM > 0.6 as hub genes
hub_genes <-  subset(gene_info, GS. >0.8 & MM.lightgreen >0.8)
 
write.csv(hub_genes, file = "hub_genes_lightgreen_0.8.csv")


# neg modules-----------------------------------------------------

#for dark orange

 module_genes <- subset(geneInfo, moduleColor == 'darkorange')
 
module_genes <-  module_genes$genes
 
  write.csv(module_genes, "module_genes_darkorange.csv")

 #get gene info for module selected
 gene_info <- geneInfo[geneInfo$genes %in% module_genes,] 

 genes <- gene_info$genes
 
#select genes with GS. >0.6 and MM > 0.6 as hub genes
hub_genes <-  subset(gene_info, GS. <-0.6 & MM.darkorange >0.6)
 
write.csv(hub_genes, file = "hub_genes_darkorange_0.6.csv")

#for darkred

 module_genes <- subset(geneInfo, moduleColor == 'darkred')
 
module_genes <-  module_genes$genes
 
  write.csv(module_genes, "module_genes_darkred.csv")

 #get gene info for module selected
 gene_info <- geneInfo[geneInfo$genes %in% module_genes,] 

 genes <- gene_info$genes
 
#select genes with GS. >0.6 and MM > 0.6 as hub genes
hub_genes <-  subset(gene_info, GS. <-0.6 & MM.darkred >0.6)
 
write.csv(hub_genes, file = "hub_genes_darkred_0.6.csv")
 
#find overlapping DEGS and hub genes-----------------------------------------------

 test <- hub_genes[hub_genes$genes %in% DEG_list_RH30_P_vs_RH_HD$geneID,] 
  test <- hub_genes[hub_genes$genes %in% DEG_list_RH30_P_vs_RH_C13$geneID,] 
  test <- hub_genes[hub_genes$genes %in% DEG_list_RH30_P_vs_RH_C7$geneID,] 

 
 
 
 #visualising networks

#  
#  # Calculate topological overlap anew: this could be done more efficiently by saving the TOM
# # calculated during module detection, but let us do it again here.
# dissTOM = 1-TOMsimilarityFromExpr(datET, power = 6);
# # Transform dissTOM with a power to make moderately strong connections more visible in the heatmap
# plotTOM = dissTOM^7;
# # Set diagonal to NA for a nicer plot
# diag(plotTOM) = NA;
# # Call the plot function
# sizeGrWindow(9,9)
# TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")
 
 # Recalculate module eigengenes
MEs = moduleEigengenes(datET, moduleColors)$eigengenes


# get clinical trait data

names(resistance) = "resistance"
# Add the weight to existing module eigengenes

 MET = orderMEs(cbind(MEs, resistance))
# Plot the relationships among the eigengenes and the trait
sizeGrWindow(5,7.5);
par(cex = 0.9)
plotEigengeneNetworks(MET, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8, xLabelsAngle
= 90)
 
 # Plot the dendrogram
sizeGrWindow(6,6);
par(cex = 1.0)
plotEigengeneNetworks(MET, "Eigengene dendrogram", marDendro = c(0,4,2,0),
plotHeatmaps = FALSE)
# Plot the heatmap matrix (note: this plot will overwrite the dendrogram plot)
par(cex = 1.0)
plotEigengeneNetworks(MET, "Eigengene adjacency heatmap", marHeatmap = c(3,4,2,2),
plotDendrograms = FALSE, xLabelsAngle = 90)

 
#  
# #with cortraits
# corTraits <- as.data.frame(cor(geneModuleMembership, geneTraitSignificance, use = "p", method="spearman")) #Finding correlation between the kME (module membership) and gene significance score
# 
# colour_df <- as.data.frame(numbers2colors(corTraits,signed = TRUE, centered = FALSE)) #using a WGCNA function to convert the numbers (cor coefficient) to a colour value (hex code) 
# colnames(colour_df) <- c("color") #labelling the name of the column on the dataframe
# colour_df$label <- colnames(datME) #adding in datME label
# 
# nodes_vis <- data.frame(id=nodes$name, label=nodes$name, value = nodes$centrality, font.size = 10) # setting up node vis - you have the script for this
# edges_vis <- data.frame(from=vis$from, to=vis$to, length=vis$weight) # ditto as above
# 
# nodes_vis <- left_join(nodes_vis, colour_df, by = "label") #using the left_join() function to merge the nodes_vis and the colour_df, visnetwork will then recognise the colour value
#  



#with module trait 
moduleTraitCor


test_colour_df <- as.data.frame(numbers2colors(moduleTraitCor,signed = TRUE, centered = FALSE))


colnames(test_colour_df) <- c("color") #labelling the name of the column on the dataframe
test_colour_df$label <- colnames(datME) #adding in datME label

nodes_vis <- data.frame(id=nodes$name, label=nodes$name, value = nodes$centrality, font.size = 10) # setting up node vis - you have the script for this
edges_vis <- data.frame(from=vis$from, to=vis$to, length=vis$weight) # ditto as above

nodes_vis <- left_join(nodes_vis, test_colour_df, by = "label") #using the left_join() function to merge the nodes_vis and the colour_df, visnetwork will then recognise the colour value
 
 
net <- visNetwork(nodes_vis, edges_vis, width = "100%") %>%
  visPhysics(stabilization = F, timestep = .1,barnesHut = list(
    centralGravity = 0.5, 
    gravitationalConstant = -10000, 
    springLength = 100, 
    springConstant = 0.05, 
    avoidOverlap = 0)) %>%
  visNodes(shape = "dot"
           , color = nodes_vis$color, 
           shadow = list(enabled = FALSE)
  ) %>%
  visEdges(
    smooth = T,
    shadow = FALSE,
    color = list(highlight = "#C62F4B")
  ) %>%
  visOptions(highlightNearest = list(enabled = T, degree = 1, hover = T)) %>% 
  visLayout(randomSeed = 11)
net

 
```

#function to get hub genes
```{r}
#function to get hub genes changing the MM and GS values

setwd("C:/Users/Christina/Documents/RMS RNA seq/Zoes data/CLONES ONLY ANALYSIS/Hub genes")

#read in gene info
gene_info <- read.csv("C:/Users/Christina/Documents/RMS RNA seq/Zoes data/CLONES ONLY ANALYSIS/geneInfo_clones.csv")

get_hub_genes <- function(module_colour, threshold){

module_genes <- subset(gene_info, moduleColor == module_colour)
module_genes <-  module_genes$genes

#get gene info for module selected
gene_info <- gene_info[gene_info$genes %in% module_genes,] %>% as_tibble()

#select columns that match module colour
subset_module_colour <- gene_info %>% dplyr::select(X, GS., contains(module_colour))

#select genes with GS. >0.6 and MM > 0.6 as hub genes
hub_genes <-  subset_module_colour %>% subset(.[[2]] >threshold & .[[3]] >threshold)

write.csv(hub_genes, file = paste0("hub_genes_",module_colour, "_", threshold,".csv"))
}
module_colours <- c("darkgreen", "darkgrey")

for (i in module_colours){
get_hub_genes(i, 0.9)}

```



#GSEA on hub genes
```{r}

#use clusterprofiler to do GSEA on module associated with resistance

library(clusterProfiler)
library(org.Hs.eg.db)
library(ReactomePA)
#read in geneInfo

gene_info <- read.csv("C:/Users/Christina/Documents/RMS RNA seq/Zoes data/CLONES ONLY ANALYSIS/geneInfo_clones.csv")

options(enrichplot.colours = c("red","blue"))
#modules <- c("darkgrey", "darkgreen")#pos modules
modules <- c( "darkgreen", "lightgreen")#neg modules


for (i in modules){
  #get hub genes from module
hub_genes <- read.csv(file=paste0("C:/Users/Christina/Documents/RMS RNA seq/Zoes data/CLONES ONLY ANALYSIS/New analysis/hub_genes_", i,"_0.8.csv") )
#get just genes 
geneList <- hub_genes$genes


#Enrichment analysis on hub genes 

#GO terms-----------------------------------------------------------------------
options(enrichplot.colours = c("red","blue"))
ego2 <- enrichGO(gene         = geneList,
                 OrgDb         = org.Hs.eg.db,
                 keyType       = 'SYMBOL',
                 ont           = "ALL",
                 pAdjustMethod = "BH",
                 pvalueCutoff  = 0.05,
                 qvalueCutoff  = 0.05)
head(ego2, 10)   
## Dotplot 
ego2 <- ego2 %>%
  dplyr::arrange(p.adjust)

GO <- try(barplot(ego2,  font.size = 16 ,showCategory=20, x="Count"))
try(ggsave(GO, file=paste0("GO_hub_genes_", i,".png"), width = 25, height = 33, units = "cm"))


# Hallmark ORA --------------------------------------------------------------------------
m_t2g <- msigdbr(species = "Homo sapiens", category = "H") %>% dplyr::select(gs_name, gene_symbol)
head(m_t2g)
em <- enricher(geneList, pvalueCutoff = 0.05,  pAdjustMethod = "BH", TERM2GENE=m_t2g)
head(em)
write.csv(em, file=paste0("Hallmark_hub_genes_", i,".csv"))

#barplot
H <- print(ggplot(em,aes(x=Count,y=Description,fill=p.adjust))+
  scale_fill_continuous(low = "red", high = "blue", na.value = NA)+geom_bar(stat='identity')+theme_classic()+
                 scale_x_continuous(expand=c(0,0))+labs(fill = 'p.adjust')+xlab('Count')+ylab('')+theme(axis.text=element_text(size=16),axis.title=element_text(size=12)))



try(ggsave(H, file=paste0("Hallmark_hub_genes_", i,".png"), width = 25, height = 10, units = "cm"))

#KEGG -----------------------------------------------------------------------

#convert to symbols
options(enrichplot.colours = c("red","blue"))
symbols <- mapIds(org.Hs.eg.db, geneList, 'ENTREZID', 'SYMBOL') %>% as.data.frame()
symbols <- symbols$.
kk <- enrichKEGG(gene         = symbols,
                 organism     = 'hsa',
                 keyType       = 'ncbi-geneid',
                 pAdjustMethod = "BH",
                 pvalueCutoff = 0.05, 
                 qvalueCutoff = 0.05)
head(kk)
write.csv(kk, file=paste0("Keg_hub_genes_", i,".csv"))
keg <- barplot(kk, font.size = 16, showCategory=20,x="Count")

try(ggsave(keg, file=paste0("KEGG_hub_genes_", i,".png"), width = 25, height = 20, units = "cm"))

#Reactome ---------------------------------------------------------------------
options(enrichplot.colours = c("red","blue"))
reactome <- enrichPathway(gene = symbols, pvalueCutoff = 0.05, readable=TRUE, qvalueCutoff =0.05, pAdjustMethod = "BH", organism = "human")
head(reactome)
write.csv(reactome, file=paste0("Reactome_hub_genes_", i,".csv"))
reactome <- barplot(reactome, showCategory=20, font.size = 16,x="Count")

try(ggsave(reactome, file=paste0("Reactome_hub_genes_", i,".png"), width = 25, height = 35, units = "cm"))

}


```


#Get gene sig 

```{r}

#for orange ------------------------------------------------------
 #read in hub proteins 

 hub_proteins <- read.csv("C:/Users/Christina/Documents/RMS RNA seq/Luciles data/RH4 WGCNA/selected_nodes_orange.csv")
hub_proteins <- hub_proteins$X

#read in WGCNA hub genes from the module
module_genes <- read.csv("C:/Users/Christina/Documents/RMS RNA seq/Luciles data/RH4 WGCNA/hub_genes_orange_0.6.csv")


 gene_sig_RH4 <- module_genes[module_genes$genes %in% hub_proteins,]
  gene_sig_RH4_orange <- gene_sig_RH4$X
  
  
write.csv(gene_sig_RH4_orange, "gene_sig_RH4_orange.csv")
  
#for lightgreen-----------------------------------------------------------------
 #read in hub proteins 

 hub_proteins <- read.csv("C:/Users/Christina/Documents/RMS RNA seq/Luciles data/RH4 WGCNA/selected_nodes_lightgreen.csv")
hub_proteins <- hub_proteins$X

#read in WGCNA hub genes from the module
module_genes <- read.csv("C:/Users/Christina/Documents/RMS RNA seq/Luciles data/RH4 WGCNA/hub_genes_lightgreen_0.6.csv")


 gene_sig_RH4 <- module_genes[module_genes$genes %in% hub_proteins,]
  gene_sig_RH4_lightgreen <- gene_sig_RH4$X
  
write.csv(gene_sig_RH4_lightgreen, "gene_sig_RH4_lightgreen.csv")



```




# Score
```{r}

#using GSVA to calculate a chemo-resistant score


#need normalised counts in matrix with rows as samples and columns as genes
#get normalised counts
rownames(edgeRData) <- paste0("X", seq_len(nrow(edgeRData)))
edgeRData <- DGEList(counts=counts.keep, genes=rownames(counts.keep), group = treatment)
edgeRData <- calcNormFactors(edgeRData, method="TMM")
CPM <- cpm(edgeRData, log = TRUE)

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

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

gsvamethod <- gsvaParam(filtered_mapped_matrix,
  hub_genes)
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)

#get sample group as extra col in df
sampleID <- row.names(gsva_results)
gsva_results$sample <- c()
for (i in 1:nrow(gsva_results)) {
  samplename <- sampleID[i]
  group <- if(grepl("RH30_P", samplename)) {print("RH30_P")
  } else if(grepl("RH_HD", samplename)) {print("RH_HD")}
  else if(grepl("RH_C7", samplename)) {print("RH_C7")
  } else if(grepl("RH_C13", samplename)) {print("RH_C13")
  }
  gsva_results$sample[i] <- group
}




#make barplot
intrinsic_resistance_score_plot=ggplot(data=gsva_results, mapping=aes(x=sample, y=gsva_results, fill=sample))+geom_boxplot() +
  stat_summary(fun = "mean", geom = "point", shape = 8,
               size = 2, color = "white") +labs(title="Intrinsic resistance score: Zoe's data",
        x ="Sample", y = "Intrinsic resistance score GSVA score")

intrinsic_resistance_score_plot


#statistical tests
#set factor
gsva_factors <- as.factor(gsva_results$sample)

#first test for normality- From the output, the p-value > 0.05 implying that the distribution of the data are not significantly different from normal distribution. In other words, we can assume the normality.

shapiro.test(gsva_results$gsva_results)


#data not normally distributed 

#we do non-parametric anova or kruskal test

kruskal.test(gsva_results ~ sample, data = gsva_results)

#As the p-value is less than the significance level 0.05, we can conclude that there are significant differences between the treatment groups.


#From the output of the Kruskal-Wallis test, we know that there is a significant difference between groups, but we don’t know which pairs of groups are different. It’s possible to use the function pairwise.wilcox.test() to calculate pairwise comparisons between group levels with corrections for multiple testing.


pairwise.wilcox.test(gsva_results$gsva_results, gsva_factors,
                 p.adjust.method = "BH")


#or ---------------------------------------------------------------

# First do nromality test then a global test
qqPlot(gsva_results$gsva_results)
qqPlot(gsva_results$gsva_results, main = "Q-Q plot of the intrinsic resistance GSVA score data", ylab = "GSVA score", distribution = "norm", envelope = FALSE)

compare_means(gsva_results ~ sample,  data = gsva_results, method = "kruskal.test")



# Default method = "kruskal.test" for multiple groups
ggboxplot(gsva_results, x = "sample", y = "gsva_results",
          color = "sample")+
  stat_compare_means()


# Perorm pairwise comparisons
compare_means(gsva_results ~ sample,  data = gsva_results, method = "wilcox.test", p.adjust.method = "BH" )


# Visualize: Specify the comparisons you want
my_comparisons <- list( c("RH4 vincristine-resistant", "RH4 ifosfamide-resistant"), c("RH4 vincristine-resistant", "RH4 control"), c("RH4 ifosfamide-resistant", "RH4 control"), c("RMSYM vincristine-resistant", "RMSYM ifosfamide-resistant"), c("RMSYM vincristine-resistant", "RMSYM control"), c("RMSYM ifosfamide-resistant", "RMSYM control") )

#reorder plot
gsva_results <- gsva_results %>%
  arrange(sample) %>%
  mutate(sample = factor(sample, levels=c("RH4 control", "RH4 vincristine-resistant", "RH4 ifosfamide-resistant", "RMSYM control", "RMSYM vincristine-resistant", "RMSYM ifosfamide-resistant")))


#boxplot of results for Luciles data only 
p <- ggboxplot(gsva_results  , x = "sample", y = "gsva_results",
          color = "treatment")+labs(title="ARMS acquired resistance score",
        x ="Treatment condition", y = "GSVA score for ARMS acquired resistance signature")+ 
  stat_compare_means(comparisons = my_comparisons,label = "p.signif")+ # Add pairwise comparisons p-value
 stat_compare_means(label.y = 1.5) + theme(legend.position = "none", axis.text.x = element_text(angle = 90,  size = 12))   # Add global p-value
ggsave("ARMS_RH4_acquired_resistance_score_GSVA.png")






#get plot of Zoes and Luciles data---------------------------------------------------

#read in GSVA scores for Zoes data
#get hub genes 
gsva_results_Luciles_data <- read.csv("C:/Users/Christina/Documents/RMS RNA seq/Zoes data/CLONES ONLY ANALYSIS/Luciles_intrinsic_score.csv") 

#drop treatment column to make number of cols match
gsva_results_Luciles_data = subset(gsva_results_Luciles_data, select = -c(treatment) )
#convert rownames to columns in luciles data
#gsva_results <- tibble::rownames_to_column(gsva_results, "X")


#merge dataframes
gsva_scores <- rbind(gsva_results_Luciles_data, gsva_results)


#reorder for plot
gsva_scores <- gsva_scores %>%
  arrange(sample) %>%
  mutate(sample = factor(sample, levels=c("RH4 control", "RH4 vincristine-resistant", "RH4 ifosfamide-resistant", "RMSYM control", "RMSYM vincristine-resistant", "RMSYM ifosfamide-resistant", "RH30_P", "RH_C7", "RH_C13", "RH_HD")))


#do stats test -------------------------

#frist normality test 
shapiro.test(gsva_scores$gsva_results) #not normally disturbuted 

qqPlot(gsva_scores$gsva_results, main = "Q-Q plot of the intrinsic resistance GSVA score data", ylab = "GSVA score")

#then global test
compare_means(gsva_results ~ sample,  data = gsva_scores, method = "kruskal.test", p.adjust.method = "BH")

# Default method = "kruskal.test" for multiple groups
ggboxplot(gsva_scores, x = "sample", y = "gsva_results",
          color = "sample")+
  stat_compare_means()


# Perorm pairwise comparisons
compare_means(gsva_results ~ sample,  data = gsva_scores, method = "wilcox.test" , p.adjust.method = "BH")


# Visualize: Specify the comparisons you want
my_comparisons <- list( c("RH4 vincristine-resistant", "RH4 ifosfamide-resistant"), c("RH4 vincristine-resistant", "RH4 control"), c("RH4 ifosfamide-resistant", "RH4 control"), c("RMSYM vincristine-resistant", "RMSYM ifosfamide-resistant"), c("RMSYM vincristine-resistant", "RMSYM control"), c("RMSYM ifosfamide-resistant", "RMSYM control"), c("RH30_P", "RH_C7"), c("RH30_P", "RH_C13"), c("RH30_P", "RH_HD"), c("RH_C7", "RH_C13"), c("RH_C7", "RH_HD"), c("RH_C13", "RH_HD") )



#boxplot of gsva scores for ARMS acquired resistance signarure for all data
p <- ggboxplot(gsva_scores  , x = "sample", y = "gsva_results",
          color = "sample")+labs(title="Intrinsic resistance score",
        x ="Treatment condition", y = "GSVA score for intrinsic resistance signature")+ 
  stat_compare_means(comparisons = my_comparisons,label = "p.signif")+ # Add pairwise comparisons p-value
 stat_compare_means(label.y = 1.5) + theme(legend.position = "none", axis.text.x = element_text(angle = 90,  size = 12))   # Add global p-value


#save gsva score results as csv
write_excel_csv(gsva_scores, "gsva_scores_intrinsic_res_all_samples.csv")


```

#PPI
```{r}



library(dplyr)
library(tidygraph)
library(visNetwork)
library(igraph)
library(RandomWalkRestartMH)
library(WGCNA)


gene_info <- read.csv('geneInfo_clones.csv')

"darkgrey"
modules <- c("darkgrey" )
# geneList <- all_clone_hub_genes
# i <- "all_clone_hub_genes"

for (i in modules)
{
 #get all genes from module
all_genes <- read.csv(file=paste0("module_genes_", i,".csv") )
#get just genes 
geneList <- all_genes$x

  #get hub genes from module
hub_genes <- read.csv(file=paste0("hub_genes_", i,"_0.7.csv") )
#get just genes 
hub_genes <- hub_genes$genes


geneList <- hub_genes

# #for C7 DEGs
# x <- 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")
# geneList <- x$X
# 
# #for C13 DEGs
# x <- 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")
# geneList <- x$X

# 
# library(dplyr)
# library(tidygraph)
# library(visNetwork)
# library(igraph)
# library(RandomWalkRestartMH)
# library(WGCNA)


pnet <- read.table("C:/Users/Christina/Documents/STRING/9606.protein.physical.links.full.v11.5.txt.gz", sep="", header = T, colClasses = "character")
pnet$experiments <- as.numeric(pnet$experiments)
filt <- pnet %>% filter(., experiments >= 400) #start with 700 change to 400 or higher. Just keep in mind that 400 is low confidence but very very few interactions are high confidence. 
filt <- filt[, c(1,2, 4, 10)] 
filt <- filt[complete.cases(filt), ]
ids <- read.table("C:/Users/Christina/Documents/STRING/9606.protein.info.v11.5.txt.gz", sep="\t", colClasses = "character", header = F)

idx <- match(filt$protein1, ids$V1)
symbol1 <- as.character(ids$V2[idx])

idx <- match(filt$protein2, ids$V1)
symbol2 <- as.character(ids$V2[idx])

filt$symbol1 <- symbol1
filt$symbol2 <- symbol2

df <- data.frame("from"=filt$symbol1, "to"=filt$symbol2)
df <- df %>% tidyr::drop_na()
df1 <- df[df$from %in% geneList,] # filter for genes in genelsit 
df2 <- df1[df1$to %in% geneList,] # filter for genes in genelsit 
g0 <- graph_from_data_frame(df2, directed = FALSE)
PNet <- g0 

hub_scores <- hub_score(PNet, scale = TRUE, weights = NULL, options = arpack_defaults)
hub_scores <- sort(hub_scores$vector, decreasing = TRUE) 
hub_scores <- as.data.frame(hub_scores)

write.csv(hub_scores, file=paste0("hub_gene_scores_PPI", i,".csv"))

###### export to visANT but don't filter on the threshold settting they proivde in WGCNA as it eliminates nodes (not just the edge)

#take protein interactio network and convert to tbl graph (a df to make a network)
Pnet <- as_tbl_graph(df2, directed = FALSE)
#activate the notes and can filter based on centraility ect
Pnet <- Pnet %>% activate(nodes) %>% mutate(centrality = centrality_degree()) %>% filter(centrality > 0)
#pulls out nodes and edges as tibbles
nodes <- Pnet %>% activate(nodes) %>% as_tibble()
edges <- Pnet %>% activate(edges) %>% as.tibble()
#gets df for vis of network
nodes_vis <- data.frame(id=1:nrow(nodes), label=nodes$name, value = nodes$centrality)
edges_vis <- data.frame(from=edges$from, to=edges$to, width = 5)


#colour df by connectivity 
colour_df <- as.data.frame(numbers2colors(nodes$centrality,signed = FALSE, centered = FALSE)) #using a WGCNA function to convert the numbers (cor coefficient) to a colour value (hex code)
colnames(colour_df) <- c("color") #labelling the name of the column on the dataframe
colour_df$label <- nodes$name #adding in datME label

nodes_vis <- left_join(nodes_vis, colour_df, by = "label") #using the left_join() function to merge the nodes_vis and the colour_df, visnetwork will then recognise the colour value


#vis of network
 projection <- visNetwork(nodes_vis, edges_vis, width = "100%") %>%
    visPhysics(stabilization = FALSE, timestep = .1,barnesHut = list(
      centralGravity = 0.5, 
      gravitationalConstant = -30000, 
      springLength = 1000, 
      springConstant = 0.04, 
      avoidOverlap = 0)) %>% 
    visNodes(shape = "dot",
             color = list(
               background = "#0085AF",
               border = "#013848",
               highlight = "#FF8000"
             ),
             scaling = list(min = 20, max = 100),
             shadow = list(enabled = FALSE)
    ) %>%
    visEdges(
      smooth = T,
      shadow = FALSE,
      color = list(highlight = "#C62F4B")
    ) %>%
    visOptions(highlightNearest = list(enabled = TRUE) ) %>% 
    visLayout(randomSeed = 123) %>% visIgraphLayout() 

  projection



#top nodes
top_nodes <- sort(nodes$centrality, decreasing= TRUE)
  

#select genes with score greater than or equal to 6
top_nodes <- as.data.frame(top_nodes)
   table(top_nodes$top_nodes >= 8) #get number for which this is true
  selected_nodes <- subset(top_nodes, top_nodes >= 8) #subset data to select for score >=6
  write.csv(selected_nodes, file=paste0("selected_nodes_", i,".csv")) #save selected nodes
  
  #find overlap between selected nodes and hub genes
  selected_nodes <- tibble::rownames_to_column(selected_nodes, "genes")#chnage rownames to col so all data is retained when subset
  gene_sig_using_PPI <- selected_nodes[selected_nodes$genes %in% hub_genes,]
  #save results
 write.csv(gene_sig_using_PPI, file=paste0("gene_sig_using_PPI_", i,".csv"))
  
}
  
  
  
  
#Run the below function first before you do so: 
function(edges, nodes)
  
{
  suppressMessages(convert_tbl_graph_to_iGraph)
  convert <- data.frame("from" = edges$from, "to" = edges$to, "direction" = 0)
  
  
  idtonode <- data.frame()
  for (j in 1:length(nodes$label)) {
    nm <- nodes$label
    nodename <- nm[j]
    ni <- nodes$id
    nodeid <- ni[j]
    for (i in 1:length(convert$to)) {
      to <- convert$to
      toid <- to[i]
      if (toid == nodeid) {idtonode[i, 1] <- nodename}
    }}
  
  idfromnode <- data.frame()
  for (j in 1:length(nodes$label)) {
    nm <- nodes$label
    nodename <- nm[j]
    ni <- nodes$id
    nodeid <- ni[j]
    for (i in 1:length(convert$from)) {
      from <- convert$from
      fromid <- from[i]
      if (fromid == nodeid) {idfromnode[i, 1] <- nodename}
    }}
  graph_df2 <- data.frame("from" = idfromnode$V1 , "to" = idtonode$V1, "direction" = 0, "method" = "M0039", "weight" = edges$weight )
  
  return(graph_df2)
  
}


Ggraph <- graph_from_data_frame(convert_tbl_graph_to_iGraph(edges = edges_vis, nodes = nodes_vis), directed = FALSE) 

#If the random walk with repeat on the multiplex network does not like the objects going in due to matching columns. Make sure both objects going in have the same columns (e.g., $to and $from)

#Below is the randomwalk with repeat
multiplexNet <- create.multiplex(list( GCN = Ggraph, PPIN = PNet ))
matrix <- compute.adjacency.matrix(multiplexNet)
Normalised <- normalize.multiplex.adjacency(matrix)
seed <- #Select a gene that is important to your network
Res <- Random.Walk.Restart.Multiplex(Normalised, multiplexNet, seed, r = 0.70)
RWRstatprob <- data.frame(Res$RWRM_Results)
TopResults <- create.multiplexNetwork.topResults(Res, multiplexNet, k=199)


##### Projection ########
nodes <- gENet %>% activate(nodes) %>% as_tibble()

edges <- gENet %>% activate(edges) %>% as_tibble()
hubs <- slice_max(nodes, order_by = nodes$centrality, n = 3) # change the n to highlight different numbers of nodes
nodes_vis <- data.frame(id=1:nrow(nodes), label=nodes$name, value = value, font.size = ifelse(value > 10, 120, 90 ), color = ifelse(nodes$name %in% hubs$name, "red", "steelblue"))
edges_vis <- data.frame(from=edges$from, to=edges$to, type =edges$type, width = 5)
edges_vis <- edges_vis %>% mutate(color = case_when(type == "GCN" ~ "red", type == "PPIN" ~ "blue"))


  projection <- visNetwork(nodes_vis, edges_vis, width = "100%") %>%
    visPhysics(stabilization = FALSE, timestep = .1,barnesHut = list(
      centralGravity = 0.5, 
      gravitationalConstant = -30000, 
      springLength = 1000, 
      springConstant = 0.04, 
      avoidOverlap = 0)) %>% 
    visNodes(shape = "dot",
             color = list(
               background = "#0085AF",
               border = "#013848",
               highlight = "#FF8000"
             ),
             scaling = list(min = 20, max = 100),
             shadow = list(enabled = FALSE)
    ) %>%
    visEdges(
      smooth = T,
      shadow = FALSE,
      color = list(highlight = "#C62F4B")
    ) %>%
    visOptions(highlightNearest = list(enabled = TRUE) ) %>% 
    visLayout(randomSeed = 123) 
  projection


```





#Get gene sig 

```{r}

#for darkgreen------------------------------------------------------
 #read in hub proteins 

 hub_proteins <- read.csv("C:/Users/Christina/Documents/RMS RNA seq/Zoes data/CLONES ONLY ANALYSIS/PPI/selected_nodes_darkgreen.csv")
hub_proteins <- hub_proteins$X

#read in WGCNA hub genes from the module
module_genes <- read.csv("C:/Users/Christina/Documents/RMS RNA seq/Zoes data/CLONES ONLY ANALYSIS/hub_genes_darkgreen_0.6.csv")


 gene_sig_intrinsic <- module_genes[module_genes$genes %in% hub_proteins,]
  gene_sig_intrinsic_darkgreen <- gene_sig_intrinsic$X
  
  
write.csv(gene_sig_intrinsic_darkgreen, "gene_sig_intrinsic_darkgreen.csv")
  
#for darkgrey------------------------------------------------------
 #read in hub proteins 

 hub_proteins <- read.csv("C:/Users/Christina/Documents/RMS RNA seq/Zoes data/CLONES ONLY ANALYSIS/PPI/selected_nodes_darkgrey.csv")
hub_proteins <- hub_proteins$X

#read in WGCNA hub genes from the module
module_genes <- read.csv("C:/Users/Christina/Documents/RMS RNA seq/Zoes data/CLONES ONLY ANALYSIS/hub_genes_darkgrey_0.6.csv")


 gene_sig_intrinsic <- module_genes[module_genes$genes %in% hub_proteins,]
  gene_sig_intrinsic_darkgrey <- gene_sig_intrinsic$X
  
  
write.csv(gene_sig_intrinsic_darkgrey, "gene_sig_intrinsic_darkgrey.csv")

```





#heatmap of genes of interest
```{r}


genes_of_interest <- c("MYH9", "CLTCL1", "PKP4", "ARHGEF7", "FLCN", "NIN", "MICAL3", "CEP104", "NCOR1", "PARP4", "MAP7D3", "PKD2", "KIF16B", "LATS2", "CHAMP1", "CDC16", "TUBGCP3", "NUMA1", "HAUS7", "CENPF", "BOD1L1", "KNTC1", "SPECC1L", "CLTC")

genes_of_interest <- read.csv("combined_DEGs_clones.csv")
genes_of_interest <- genes_of_interest$x

#get normalised count matrix
rownames(edgeRData) <- paste0("X", seq_len(nrow(edgeRData)))
edgeRData <- DGEList(counts=counts.keep, genes=rownames(counts.keep), group = treatment)
edgeRData <- calcNormFactors(edgeRData, method="TMM")
CPM <- cpm(edgeRData, log = FALSE)

CPM <- as.matrix(CPM)
scaledata <- t(scale(t(CPM)))

#subset for genes of interest
expression_of_genes_of_interest <- scaledata[row.names(scaledata) %in% genes_of_interest ,]

h1<-Heatmap(as.matrix(expression_of_genes_of_interest[,1:(ncol(expression_of_genes_of_interest))]),cluster_rows=TRUE,cluster_columns=FALSE,use_raster=FALSE,row_dend_reorder = TRUE,show_row_dend = FALSE,column_split=factor(c(rep('RH_C13',6),rep('RH_C7',6),rep('RH_HD',6), rep("RH30_P", 6)),levels=c('RH_C13','RH_C7','RH_HD', 'RH30_P')),show_row_names=TRUE,show_column_names = FALSE, heatmap_legend_param = list(title = "Z-score"))

hlist <- h1
tiff(paste0('hallamrks_immune_only_Hmap_kmeans_k','_clusters.tiff'),width=2000,height=3000,res=300)
draw(hlist,cluster_rows=FALSE, row_dend_reorder = TRUE,show_row_dend = FALSE)
dev.off()

```




#RH4 acquired resistance score ZOes data
```{r}
#using GSVA to calculate a chemo-resistant score

#get normalised filtered count matrix
#get normalised counts
rownames(edgeRData) <- paste0("X", seq_len(nrow(edgeRData)))
edgeRData <- DGEList(counts=counts.keep, genes=rownames(counts.keep), group = treatment)
edgeRData <- calcNormFactors(edgeRData, method="TMM")
CPM <- cpm(edgeRData, log = TRUE)
library(GSVA)

#get genes 
hub_genes <- read_xlsx("C:/Users/Christina/Documents/RMS RNA seq/Luciles data/RH4 WGCNA/RH4_acquired_res_signature.xlsx") 
hub_genes <- hub_genes$Gene
hub_genes <- list(hub_genes)

filtered_mapped_matrix <- CPM %>%
  as.matrix()

gsva_results <- gsva(
  filtered_mapped_matrix,
  hub_genes,
  method = "gsva",
  # Minimum gene set size
  min.sz = 1,
  # Maximum gene set size
  max.sz = 500
)


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


#get sample group as extra col in df
sampleID <- row.names(gsva_results)
gsva_results$sample <- c()
for (i in 1:nrow(gsva_results)) {
  samplename <- sampleID[i]
  group <- if(grepl("RH30_P", samplename)) {print("RH30_P")
  } else if(grepl("RH_HD", samplename)) {print("RH_HD")}
  else if(grepl("RH_C7", samplename)) {print("RH_C7")
  } else if(grepl("RH_C13", samplename)) {print("RH_C13")
  }
  gsva_results$sample[i] <- group
}


#save resutls to csv to join to plot in luciles data
write_excel_csv(gsva_results, "Zoes_RH4_acquired_score")



#make barplot
intrinsic_resistance_score_plot=ggplot(data=gsva_results, mapping=aes(x=sample, y=gsva_results, fill=sample))+geom_boxplot() +
  stat_summary(fun = "mean", geom = "point", shape = 8,
               size = 2, color = "white") +labs(title="Intrinsic chemo-resistance score",
        x ="Condition", y = "Intrinsic resistance GSVA score") +geom_signif(stat="signif",position="identity",
                 comparisons=list(c("RH_C13","RH30_P")),map_signif_level = TRUE,annotations="*")+geom_signif(stat="signif",position="identity",
                 comparisons=list(c("RH_C7","RH30_P")),map_signif_level = TRUE,annotations="*")+geom_signif(stat="signif",position="identity",
                 comparisons=list(c("RH_C13","RH_HD")),map_signif_level = TRUE,annotations="*")+geom_signif(stat="signif",position="identity",
                 comparisons=list(c("RH_C7","RH_HD")),map_signif_level = TRUE,annotations="*")


intrinsic_resistance_score_plot

stat.test <- gsva_results$gsva_results %>% pairwise.wilcox.test( gsva_factors,
                 p.adjust.method = "BH")
stat.test


# Box plot
stat.test <- stat.test %>% add_xy_position(x = "sample")
bxp + stat_pvalue_manual(stat.test, label = "p.adj.signif", tip.length = 0.01)



#statistical tests
#set factor
gsva_factors <- as.factor(gsva_results$sample)

#first test for normality- From the output, the p-value > 0.05 implying that the distribution of the data are not significantly different from normal distribution. In other words, we can assume the normality.

shapiro.test(gsva_results$gsva_results)

qqPlot(gsva_results$gsva_results)

#data not normally distributed 

#we do non-parametric anova or kruskal test

kruskal.test(gsva_results ~ gsva_factors, data = gsva_results)

#As the p-value is less than the significance level 0.05, we can conclude that there are significant differences between the treatment groups.


#From the output of the Kruskal-Wallis test, we know that there is a significant difference between groups, but we don’t know which pairs of groups are different. It’s possible to use the function pairwise.wilcox.test() to calculate pairwise comparisons between group levels with corrections for multiple testing.


pairwise.wilcox.test(gsva_results$gsva_results, gsva_factors,
                 p.adjust.method = "BH")

```



#RMSYM acquired resistance score ZOes data
```{r}
#using GSVA to calculate a chemo-resistant score

#get normalised filtered count matrix
#get normalised counts
rownames(edgeRData) <- paste0("X", seq_len(nrow(edgeRData)))
edgeRData <- DGEList(counts=counts.keep, genes=rownames(counts.keep), group = treatment)
edgeRData <- calcNormFactors(edgeRData, method="TMM")
CPM <- cpm(edgeRData, log = TRUE)
library(GSVA)

#get genes 
hub_genes <- read_xlsx("C:/Users/Christina/Documents/RMS RNA seq/Luciles data/RMSYM WGCNA/RMSYM_acquired_resistance_gene_signature.xlsx") 
hub_genes <- hub_genes$Gene
hub_genes <- list(hub_genes)

filtered_mapped_matrix <- CPM %>%
  as.matrix()

gsva_results <- gsva(
  filtered_mapped_matrix,
  hub_genes,
  method = "gsva",
  # Minimum gene set size
  min.sz = 1,
  # Maximum gene set size
  max.sz = 500
)


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


#get sample group as extra col in df
sampleID <- row.names(gsva_results)
gsva_results$sample <- c()
for (i in 1:nrow(gsva_results)) {
  samplename <- sampleID[i]
  group <- if(grepl("RH30_P", samplename)) {print("RH30_P")
  } else if(grepl("RH_HD", samplename)) {print("RH_HD")}
  else if(grepl("RH_C7", samplename)) {print("RH_C7")
  } else if(grepl("RH_C13", samplename)) {print("RH_C13")
  }
  gsva_results$sample[i] <- group
}


#save resutls to csv to join to plot in luciles data
write_excel_csv(gsva_results, "Zoes_RMSYM_acquired_score.csv")



#make barplot
intrinsic_resistance_score_plot=ggplot(data=gsva_results, mapping=aes(x=sample, y=gsva_results, fill=sample))+geom_boxplot() +
  stat_summary(fun = "mean", geom = "point", shape = 8,
               size = 2, color = "white") +labs(title="Intrinsic chemo-resistance score",
        x ="Condition", y = "Intrinsic resistance GSVA score") +geom_signif(stat="signif",position="identity",
                 comparisons=list(c("RH_C13","RH30_P")),map_signif_level = TRUE,annotations="*")+geom_signif(stat="signif",position="identity",
                 comparisons=list(c("RH_C7","RH30_P")),map_signif_level = TRUE,annotations="*")+geom_signif(stat="signif",position="identity",
                 comparisons=list(c("RH_C13","RH_HD")),map_signif_level = TRUE,annotations="*")+geom_signif(stat="signif",position="identity",
                 comparisons=list(c("RH_C7","RH_HD")),map_signif_level = TRUE,annotations="*")


intrinsic_resistance_score_plot

stat.test <- gsva_results$gsva_results %>% pairwise.wilcox.test( gsva_factors,
                 p.adjust.method = "BH")
stat.test


# Box plot
stat.test <- stat.test %>% add_xy_position(x = "sample")
bxp + stat_pvalue_manual(stat.test, label = "p.adj.signif", tip.length = 0.01)



#statistical tests
#set factor
gsva_factors <- as.factor(gsva_results$sample)

#first test for normality- From the output, the p-value > 0.05 implying that the distribution of the data are not significantly different from normal distribution. In other words, we can assume the normality.

shapiro.test(gsva_results$gsva_results)

qqPlot(gsva_results$gsva_results)

#data not normally distributed 

#we do non-parametric anova or kruskal test

kruskal.test(gsva_results ~ gsva_factors, data = gsva_results)

#As the p-value is less than the significance level 0.05, we can conclude that there are significant differences between the treatment groups.


#From the output of the Kruskal-Wallis test, we know that there is a significant difference between groups, but we don’t know which pairs of groups are different. It’s possible to use the function pairwise.wilcox.test() to calculate pairwise comparisons between group levels with corrections for multiple testing.


pairwise.wilcox.test(gsva_results$gsva_results, gsva_factors,
                 p.adjust.method = "BH")

```



#HEATMAP

```{r}

genes_of_interest = c("TUBGCP3","NUMA1","CENPF","KNTC1","CDC42BPA","FARP1","SPTAN1","ITSN1","RAB3GAP1","PCNT","CKAP5","DOCK2","MID1","ESPL1","TAOK2","NF1","CSNK1D","LMNB1","NCK2","MARK4","MAP3K11","CEP250","ARHGEF12","KIF4A")

genes_of_interest = c("TUBGCP3","NUMA1","CENPF","KNTC1","CDC42BPA","FARP1","SPTAN1","ITSN1","RAB3GAP1","PCNT","CKAP5","DOCK2","MID1","ESPL1","TAOK2","NF1","CSNK1D","LMNB1","NCK2","MARK4","MAP3K11","CEP250","ARHGEF12","KIF4A", "MYH10","FLNB","MYH9","ARHGEF7","RASAL2","BCAR1","RFC1","NIN","FLNA","PKD2","RABGAP1","TRIO","NOTCH2","ACTN4","ANLN","TIAM1","CLIP1","SMC4","PAFAH1B1","RAPGEF6")

genes_of_interest <- unique(genes_of_interest)


#get normalised count matrix
rownames(edgeRData) <- paste0("X", seq_len(nrow(edgeRData)))
edgeRData <- DGEList(counts=counts.keep, genes=rownames(counts.keep), group = treatment)
edgeRData <- calcNormFactors(edgeRData, method="TMM")
CPM <- cpm(edgeRData, log = TRUE)

 #subset raw data for the list
 expression_of_genes_of_interest <- CPM[row.names(CPM) %in% genes_of_interest,]
 
#Since we are comparing gene expression patterns we need to scale the data otherwise all of the highly expressed genes will cluster together even if they have different patterns among the samples. The function scale scales collumns. So to scale per gene (rows) data needs to be transposed. Z-scores are calculated for each row (each gene) and these are plotted instead of the normalized expression values; this ensures that the expression patterns/trends that we want to visualize are not overwhelmed by the expression values.Z-scores are computed on a gene-by-gene basis by subtracting the mean and then dividing by the standard deviation. The Z-scores are computed after the clustering, so that it only affects the graphical aesthetics and the color visualization is improved.
#scale data. https://scienceparkstudygroup.github.io/rna-seq-lesson/08-cluster-analysis/index.html
DEgenes <- as.matrix(expression_of_genes_of_interest)
scaledata <- t(scale(t(DEgenes)))


#reoorder cols
col_order <- c("DMSO_PT3_1", "DMSO_PT3_2", "DMSO_PT3_3",
                "DMSO_PT7_1", "DMSO_PT7_2", "DMSO_PT7_3", "DMSO_PT10_1", "DMSO_PT10_2", "DMSO_PT10_3",
                "Isotretinoin_PT3_1", "Isotretinoin_PT3_2", "Isotretinoin_PT3_3", "Isotretinoin_PT7_1", "Isotretinoin_PT7_2", "Isotretinoin_PT7_3", "Isotretinoin_PT10_1", "Isotretinoin_PT10_2", "Isotretinoin_PT10_3",
                "TAZ_PT3_1", "TAZ_PT3_2", "TAZ_PT3_3",    "TAZ_PT7_1", "TAZ_PT7_2", "TAZ_PT7_3",  "TAZ_PT10_1", "TAZ_PT10_2", "TAZ_PT10_3",
                "GSK_PT3_1", "GSK_PT3_2", "GSK_PT3_3",    "GSK_PT7_1", "GSK_PT7_2", "GSK_PT7_3",  "GSK_PT10_1", "GSK_PT10_2", "GSK_PT10_3",
                "Isotretinoin_and_TAZ_PT3_1", "Isotretinoin_and_TAZ_PT3_2", "Isotretinoin_and_TAZ_PT3_3",    "Isotretinoin_and_TAZ_PT7_1", "Isotretinoin_and_TAZ_PT7_2", "Isotretinoin_and_TAZ_PT7_3",  "Isotretinoin_and_TAZ_PT10_1", "Isotretinoin_and_TAZ_PT10_2", "Isotretinoin_and_TAZ_PT10_3",
                "Isotretinoin_and_GSK_PT3_1", "Isotretinoin_and_GSK_PT3_2", "Isotretinoin_and_GSK_PT3_3",    "Isotretinoin_and_GSK_PT7_1", "Isotretinoin_and_GSK_PT7_2", "Isotretinoin_and_GSK_PT7_3",  "Isotretinoin_and_GSK_PT10_1", "Isotretinoin_and_GSK_PT10_2", "Isotretinoin_and_GSK_PT10_3")
scaledata <- scaledata[, col_order]

#plot heatmap of genes of interest
h1<-Heatmap(as.matrix(scaledata[,1:(ncol(scaledata))]),cluster_rows=TRUE,cluster_columns=FALSE,use_raster=FALSE,row_dend_reorder = TRUE,show_row_dend = FALSE,column_split=factor(c(rep('DMSO',9),rep('Isotretinoin',9),rep('TAZ',9), rep("GSK", 9), rep("Isotretinoin_and_TAZ", 9), rep("Isotretinoin_and_GSK",9)),levels=c('DMSO','Isotretinoin','TAZ', 'GSK','Isotretinoin_and_TAZ','Isotretinoin_and_GSK')),show_row_names=FALSE,show_column_names = FALSE, heatmap_legend_param = list(title = "Z-score"))

hlist <- h1
tiff(paste0('Immune_response_TAZ.tiff'),width=3600,height=3000,res=300)
draw(hlist,  cluster_row_slices = TRUE,cluster_rows=FALSE, row_dend_reorder = TRUE,show_row_dend = FALSE)
dev.off()


colors = c(seq(-3,-2,length=100),seq(-2,0.5,length=100),seq(0.5,6,length=100))

my_palette <- colorRampPalette(c("blue", "black", "red"))



#heatmap ----------------------------------------------------------------------
heatmap.2(expression_of_genes_of_interest,
          col=my_palette(100),
          scale="row",
          margins = c(10, 7),
          cexCol = 0.7,
          labRow = F,
          main = "Mitotic spindle response genes",
          trace = "none")
```


