---
title: "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/All_Raw_Reads')
library(knitr)
```

#load libraries
```{r}
library(dplyr)
library(readr)
library(clusterProfiler)
library(ComplexHeatmap)
#library("Rsubread")
library(GSVA)
library(plyr)
#library(rrvgo)
library(factoextra)
library(scales)
library(circlize)
library(stringr)
library(msigdbr)
library(forcats)
library(tidyverse)
library(ggpubr)
library(readxl)
library(RColorBrewer)
library(org.Hs.eg.db)
#library(ReactomePA)
library(fgsea)
library(msigdbr)
library(stringr)
library(edgeR)
library("gplots")
library(affycoretools)
library("ggrepel")
library("ggplot2")
library('knitr')
library(ggplot2)
library(factoextra)
```

# Set up project
Load libraries, and read in data, check data dimensions
```{r}
#set wd
setwd("~/NB/All_Raw_Reads")

######### COUNT MATRIX #####################################

# #create count matrix by joining RawReadCountsOnly.out files
# files <- list.files(pattern = "RawReadCountsOnly.out" )
# counts <- map(files, read.table, header=TRUE) %>% 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:54] = lapply(rawData[1:54], FUN = function(y){as.numeric(y)})
# 
# #write.csv(rawData, "~/Neuroblastoma RNAseq analysis/EM.csv")
# #write.csv(counts.keep, "~/Neuroblastoma RNAseq analysis/filtered_EM.csv")


#read in count matrix
counts <- read.csv("EM.csv")
counts <- counts %>% column_to_rownames(var="X")
rawData <- counts


colnames(rawData) = gsub("COMB1", "Isotretinoin_and_TAZ", colnames(rawData))

colnames(rawData) = gsub("COMB2", "Isotretinoin_and_GSK", colnames(rawData))


colnames(rawData) = gsub("RA", "Isotretinoin", colnames(rawData))

#rawData = subset(rawData, select = c(TAZ_PT3_1 , TAZ_PT3_2, TAZ_PT3_3, DMSO_PT3_1, DMSO_PT3_2, DMSO_PT3_3))

```



* genes = `r dim(rawData)[1]`
* samples = `r dim(rawData)[2]`

```{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("_and_TAZ", samplename)) {print("Isotretinoin_and_TAZ")
  } else if(grepl("_and_GSK", samplename)) {print("Isotretinoin_and_GSK")
    } else if (grepl("GSK", samplename)) {print("GSK")
  } else if(grepl("Isotretinoin", samplename)) {print("Isotretinoin")}
  else if(grepl("DMSO", samplename)) {print("DMSO")
  } 
  else {print("TAZ")}
  treatment[i] <- group
}

#for timepoint

sampleID <- colnames(rawData)
timepoint <- c()
for (i in 1:ncol(rawData)) {
  samplename <- sampleID[i]
  group <- if(grepl("PT3", samplename)) {print("PT3")
  } else if(grepl("PT7", samplename)) {print("PT7")}
  else if(grepl("PT10", samplename)) {print("PT10")}
  timepoint[i] <- group
}

#set timepoint as numeric
#timepoint[1:54] = lapply(timepoint, FUN = function(y){as.numeric(y)})


#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")}
  replicate[i] <- group
}


#set replicate as numeric
#replicate[1:54] = lapply(replicate[1:54], FUN = function(y){as.numeric(y)})


#for both treatment and timepoint
sampleID <- colnames(rawData)
treatment_timepoint <- c()
for (i in 1:ncol(rawData)) {
  samplename <- sampleID[i]
  group <- if(grepl("Isotretinoin\\_and_TAZ_PT3", samplename)) {print("Isotretinoin_and_TAZ_PT3")} 
  else if(grepl("Isotretinoin\\_and_TAZ_PT7", samplename)) {print("Isotretinoin_and_TAZ_PT7")}
  else if(grepl("Isotretinoin\\_and_TAZ_PT10", samplename)) {print("Isotretinoin_and_TAZ_PT10")}
   else if(grepl("Isotretinoin\\_and_GSK_PT3", samplename)) {print("Isotretinoin_and_GSK_PT3")}
   else if(grepl("Isotretinoin\\_and_GSK_PT7", samplename)) {print("Isotretinoin_and_GSK_PT7")}
    else if(grepl("Isotretinoin\\_and_GSK_PT10", samplename)) {print("Isotretinoin_and_GSK_PT10")}
      else if(grepl("DMSO_PT3", samplename)) {print("DMSO_PT3")}
  else if(grepl("DMSO_PT7", samplename)) {print("DMSO_PT7")}
  else if(grepl("DMSO_PT10", samplename)) {print("DMSO_PT10")}
  else if(grepl("GSK_PT3", samplename)) {print("GSK_PT3")}
  else if(grepl("GSK_PT7", samplename)) {print("GSK_PT7")}
  else if(grepl("GSK_PT10", samplename)) {print("GSK_PT10")}
  else if(grepl("Isotretinoin_PT3", samplename)) {print("Isotretinoin_PT3")}
  else if(grepl("Isotretinoin_PT7", samplename)) {print("Isotretinoin_PT7")}
  else if(grepl("Isotretinoin_PT10", samplename)) {print("Isotretinoin_PT10")}
  else if(grepl("TAZ_PT3", samplename)) {print("TAZ_PT3")}
  else if(grepl("TAZ_PT7", samplename)) {print("TAZ_PT7")}
  else if(grepl("TAZ_PT10", samplename)) {print("TAZ_PT10")
  }
  treatment_timepoint[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 3 in each group, we favour filtering on a minimum counts per million threshold present in at least 3 samples. 3 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.5 in at least 3 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[,54],rawData[,54])
# Let us limit the x and y-axis so we can actually look to see what is happening at the smaller counts
plot(myCPM[,54],rawData[,54],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.4)


#thresh to 0.5 to correspond to around 15 counts 


# Which values in myCPM are greater than 0.5?
thresh <- myCPM > 0.4
# 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 3 TRUES in each row of thresh
keep <- rowSums(thresh) >= 3
# 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)

```



#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

```


```{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
```



#Create DGEList object
```{r}

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

```


#add sample data to edgeR object

```{r}


#add replicate to edger
edgeRData$samples$replicate
#for replicate
sampleID <- colnames(rawData)
edgeRData$samples$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")}
  edgeRData$samples$replicate[i] <- group
}



#timepoint 
names <- colnames(edgeRData$counts)
for (i in 1:54) {
  string <- names[i]
  edgeRData$samples$timepoint[i] <- if(grepl("PT3", string)) {print("PT3")} else if (grepl("PT7", string)) {print("PT7")} else {print("PT10")}}

#timepint treatment



 Group <- factor(paste(edgeRData$samples$group,edgeRData$samples$timepoint,sep="."))
edgeRData$samples <- cbind(edgeRData$samples,Group=Group)
levels(edgeRData$samples$Group)


levels(edgeRData$samples$replicate) <- c("1", "3")

```


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 )

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

png(filename="logCPM_after_norm.png", width=900, bg="white")
par(mar=c(11,6,4,1)+.1)
logCountsPerM <- cpm(edgeRData, log=TRUE)
for_median <- cpm(edgeRData, log=TRUE) #creating seperate df to caluclate median as after I reorder cannot calc median

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")
logCountsPerM <- as.data.frame(logCountsPerM)
logCountsPerM <- logCountsPerM[, col_order]
#logCountsPerM <- as.numeric(unlist(logCountsPerM))

boxplot(logCountsPerM, xlab="", ylab="Log2 counts per million", las=2, cex.names = 0.8, cex.axis = 0.8)
abline(h=median(for_median), col="blue")
title("Boxplots of log-2 CPMs after normalisation")
title(xlab="Sample", line=6, cex.lab=1.1, family="Calibri Light")
dev.off()

```

#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("Boxplots of logCPMs ")
```

#MDS plot

```{r}


  #plotting by treatment and timepoint


names <- colnames(edgeRData$counts)
for (i in 1:54) {
  string <- names[i]
  edgeRData$samples$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")}}

levels = c("DMSO", "Isotretinoin", "TAZ", "GSK", "Isotretinoin_and_TAZ", "Isotretinoin_and_GSK")


png(filename="MDS.png", width=900, bg="white")
par(mar=c(5, 4, 4, 12), xpd=TRUE)
  #plotting by treatment 
treatment_names <- c("DMSO", "Isotretinoin", "TAZ", "GSK", "Isotretinoin_and_TAZ", "Isotretinoin_and_GSK")
pchs<- c(0:3)
gr.moi <- factor(edgeRData$samples$timepoint, c("Post-treatment day 3", "Post-treatment day 7", "Post-treatment day 10") )
edgeRData$samples$group<- factor(edgeRData$samples$group, c("DMSO", "Isotretinoin", "TAZ", "GSK", "Isotretinoin_and_TAZ", "Isotretinoin_and_GSK"))

#MDS plot
col.group <- c("red",  "darkgoldenrod1", "springgreen1",  "turquoise3",  "dodgerblue1",  "magenta1")[edgeRData$samples$group]
shps = c(15, 16, 17)
data.frame(edgeRData$samples$group, levels = c("DMSO", "Isotretinoin", "TAZ", "GSK", "Isotretinoin_and_TAZ", "Isotretinoin_and_GSK") ,col.group)
plotMDS(edgeRData, levels = c("DMSO", "Isotretinoin", "TAZ", "GSK", "Isotretinoin_and_TAZ", "Isotretinoin_and_GSK"), col=col.group, cex=2, pch= pchs[gr.moi])
legend(fill=c("red",  "darkgoldenrod1", "springgreen1",  "turquoise3",  "dodgerblue1",  "magenta1"),legend=treatment_names, title = "Sample", "topright",   inset=c(-0.25, 0))
  title("MDS plot")
  
  legend( "right", legend=levels(gr.moi), bty = "n", col= c("black"), pch = pchs, cex = 1,  inset=c(-0.25, 0) ,title = "Timepoint")
  dev.off()
  
  

  #by replicate

names <- colnames(edgeRData$counts)
for (i in 1:54) {
  string <- names[i]
  edgeRData$samples$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")}}

levels = c("DMSO", "Isotretinoin", "TAZ", "GSK", "Isotretinoin_and_TAZ", "Isotretinoin_and_GSK")


png(filename="MDS_rep.png", width=900, bg="white")
par(mar=c(5, 4, 4, 12), xpd=TRUE)
  #plotting by treatment 
treatment_names <- c("DMSO", "Isotretinoin", "TAZ", "GSK", "Isotretinoin_and_TAZ", "Isotretinoin_and_GSK")
pchs<- c(0:3)
gr.moi <- factor(edgeRData$samples$replicate, c("1", "2", "3") )
edgeRData$samples$group<- factor(edgeRData$samples$group, c("DMSO", "Isotretinoin", "TAZ", "GSK", "Isotretinoin_and_TAZ", "Isotretinoin_and_GSK"))

#MDS plot
col.group <- c("red",  "darkgoldenrod1", "springgreen1",  "turquoise3",  "dodgerblue1",  "magenta1")[edgeRData$samples$group]
shps = c(15, 16, 17)
data.frame(edgeRData$samples$group, levels = c("DMSO", "Isotretinoin", "TAZ", "GSK", "Isotretinoin_and_TAZ", "Isotretinoin_and_GSK") ,col.group)
plotMDS(edgeRData, levels = c("DMSO", "Isotretinoin", "TAZ", "GSK", "Isotretinoin_and_TAZ", "Isotretinoin_and_GSK"), col=col.group, cex=2, pch= pchs[gr.moi])
legend(fill=c("red",  "darkgoldenrod1", "springgreen1",  "turquoise3",  "dodgerblue1",  "magenta1"),legend=treatment_names, title = "Sample", "topright",   inset=c(-0.25, 0))
  title("MDS plot")
  
  legend( "right", legend=levels(gr.moi), bty = "n", col= c("black"), pch = pchs, cex = 1,  inset=c(-0.25, 0) ,title = "Replicate")
  dev.off()
  
  levels(edgeRData$samples$timepoint) <- c("Post-treatment day 3", "Post-treatment day 7", "Post-treatment day 10")
  

#just one treatment and timepoint
  
png(filename="MDS_rep_TAZvsDMSO_PT3.png", width=900, bg="white")
par(mar=c(5, 4, 4, 12), xpd=TRUE)
  #plotting by treatment 
pchs<- c(0:3)
gr.moi <- factor(edgeRData$samples$replicate, c("1", "2", "3") )
  
#MDS plot
shps = c(15, 16, 17)
col.group <- c("red", "springgreen1")[edgeRData$sample$group]
plotMDS(edgeRData, levels = c("DMSO", "Isotretinoin", "TAZ", "GSK", "Isotretinoin_and_TAZ", "Isotretinoin_and_GSK"), col=col.group, cex=2, pch= pchs[gr.moi])
  title("MDS plot TAZ and DMSO Post-treatment day 3")
  legend( "right", legend=levels(gr.moi), bty = "n", col= c("black"), pch = pchs, cex = 1,  inset=c(-0.25, 0) ,title = "Replicate")
  legend(fill=c("orange", "green"),legend=c("DMSO", "TAZ"), title = "Sample", "topright",   inset=c(-0.35, 0))
  dev.off()





```



#Batch effects
```{r}

rep <- edgeRData$samples$replicate
Group <- edgeRData$samples$Group

#no batch effect included:
# design <- model.matrix(~0+Group, data=edgeRData$samples)
# colnames(design) <- levels(edgeRData$samples$Group)
# design

#batch effect included:
 design <- model.matrix(~ 0 + Group + rep)
 rownames(design) <- colnames(edgeRData)
 design

 
 edgeRData <- estimateDisp(edgeRData, design, robust=TRUE)
 
 fit <- glmQLFit(edgeRData, design, robust=TRUE)
 qlf <- glmQLFTest(fit, coef=1:2)
top <- topTags(qlf)

#write.csv(top, file = "rep_LFC.csv")

# Before we fit GLMs, we need to define our design matrix based on the experimental design.
# We want to test for differential expressions between ∆hrcC challenged and mock-inoculated
# samples within batches, i.e. adjusting for differences between batches. In statistical terms,
# this is an additive linear model. So the design matrix is created as:


```



#Hierarchical clustering with heatmaps
An alternative to plotMDS for examining relationships between samples is using hierarchical clustering. Heatmaps are a nice visualisation to examine hierarchical clustering of your samples. We can do this using the heatmap.2 function from the gplots package. In this example heatmap.2 calculates a matrix of euclidean distances from the logCPM (logcounts object) for the 500 most variable genes. (Note this has more complicated code than plotting principle components using plotMDS.)

The RColorBrewer package has nicer colour schemes, accessed using the brewer.pal function. “RdYlBu” is a common choice, and “Spectral” is also nice.

Note:The png function will create a png file to save the plots created straight after, and will close this file when dev.off() is called. To see your plots interactively, simply omit those two lines.

Let’s select data for the 500 most variable genes and plot the heatmap

```{r}
# We estimate the variance for each row in the logcounts matrix
var_genes <- apply(logcounts, 1, var)
head(var_genes)
# Get the gene names for the top 500 most variable genes
select_var <- names(sort(var_genes, decreasing=TRUE))[1:500]
head(select_var)
# Subset logcounts matrix
highly_variable_lcpm <- logcounts[select_var,]
dim(highly_variable_lcpm)
head(highly_variable_lcpm)

## Get some nicer colours
mypalette <- brewer.pal(11, "RdYlBu")
morecols <- colorRampPalette(mypalette)
# Set up colour vector for celltype variable
col.cell <- c("purple","orange")[edgeRData$samples$group]

# Plot the heatmap
heatmap.2(highly_variable_lcpm, 
          col=rev(morecols(50)),
          trace="column", 
          main="Top 500 most variable genes across samples",
          ColSideColors=col.cell,scale="row")

# Save the heatmap
png(file="High_var_genes.heatmap.png")
heatmap.2(highly_variable_lcpm,col=rev(morecols(50)),trace="none", main="Top 500 most variable genes\nacross samples",ColSideColors=col.cell,scale="row")
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")
edgeRData



```

> 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}
#Estimate common, trended and tagwise dispersion.
# 
 # edgeRData <- estimateCommonDisp(edgeRData, verbose = TRUE)
 # edgeRData <- estimateTrendedDisp(edgeRData)
 # edgeRData <- estimateTagwiseDisp(edgeRData)
  #plotBCV(edgeRData)


# note values of common dispersion and BCV
```

#Design matrix
```{r}



```



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

  edgeRData <- estimateDisp(edgeRData, design, robust=TRUE)
  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}

###################### DGE glmfit ###############################################

#create design matrix


#for both treatment and timepoint
sampleID <- colnames(rawData)
treatment_timepoint <- c()
for (i in 1:ncol(rawData)) {
  samplename <- sampleID[i]
  group <- if(grepl("Isotretinoin\\_and_TAZ_PT3", samplename)) {print("Isotretinoin_and_TAZ_PT3")}
  else if(grepl("Isotretinoin\\_and_TAZ_PT7", samplename)) {print("Isotretinoin_and_TAZ_PT7")}
  else if(grepl("Isotretinoin\\_and_TAZ_PT10", samplename)) {print("Isotretinoin_and_TAZ_PT10")}
   else if(grepl("Isotretinoin\\_and_GSK_PT3", samplename)) {print("Isotretinoin_and_GSK_PT3")}
   else if(grepl("Isotretinoin\\_and_GSK_PT7", samplename)) {print("Isotretinoin_and_GSK_PT7")}
    else if(grepl("Isotretinoin\\_and_GSK_PT10", samplename)) {print("Isotretinoin_and_GSK_PT10")}
      else if(grepl("DMSO_PT3", samplename)) {print("DMSO_PT3")}
  else if(grepl("DMSO_PT7", samplename)) {print("DMSO_PT7")}
  else if(grepl("DMSO_PT10", samplename)) {print("DMSO_PT10")}
  else if(grepl("GSK_PT3", samplename)) {print("GSK_PT3")}
  else if(grepl("GSK_PT7", samplename)) {print("GSK_PT7")}
  else if(grepl("GSK_PT10", samplename)) {print("GSK_PT10")}
  else if(grepl("Isotretinoin_PT3", samplename)) {print("Isotretinoin_PT3")}
  else if(grepl("Isotretinoin_PT7", samplename)) {print("Isotretinoin_PT7")}
  else if(grepl("Isotretinoin_PT10", samplename)) {print("Isotretinoin_PT10")}
  else if(grepl("TAZ_PT3", samplename)) {print("TAZ_PT3")}
  else if(grepl("TAZ_PT7", samplename)) {print("TAZ_PT7")}
  else if(grepl("TAZ_PT10", samplename)) {print("TAZ_PT10")
  }
  treatment_timepoint[i] <- group
}



# create contrast matrix
my.contrasts <- makeContrasts(
  
  ##########treatment vs DMSO ------------------------------------------------------
  
  #iso_TAZ vs DMSO
  Isotretinoin_and_TAZvsDMSO_PT3 = GroupIsotretinoin_and_TAZ.PT3-GroupDMSO.PT3, 
  Isotretinoin_and_TAZvsDMSO_PT7=GroupIsotretinoin_and_TAZ.PT7-GroupDMSO.PT7, 
  Isotretinoin_and_TAZvsDMSO_PT10=GroupIsotretinoin_and_TAZ.PT10-GroupDMSO.PT10,

  #iso+GSK vs DMSO
   Isotretinoin_and_GSKvsDMSO_PT3 = GroupIsotretinoin_and_GSK.PT3-GroupDMSO.PT3, 
  Isotretinoin_and_GSKvsDMSO_PT7=GroupIsotretinoin_and_GSK.PT7-GroupDMSO.PT7, 
  Isotretinoin_and_GSKvsDMSO_PT10=GroupIsotretinoin_and_GSK.PT10-GroupDMSO.PT10,

#GSK vs DMSO
    GSKvsDMSO_PT3 = GroupGSK.PT3-GroupDMSO.PT3, 
  GSKvsDMSO_PT7=GroupGSK.PT7-GroupDMSO.PT7, 
  GSKvsDMSO_PT10=GroupGSK.PT10-GroupDMSO.PT10,

#TAZ vs DMSO
TAZvsDMSO_PT3 = GroupTAZ.PT3-GroupDMSO.PT3, 
  TAZvsDMSO_PT7=GroupTAZ.PT7-GroupDMSO.PT7, 
  TAZvsDMSO_PT10=GroupTAZ.PT10-GroupDMSO.PT10,

#isotretinoin vs DMSO

IsotretinoinvsDMSO_PT3 = GroupIsotretinoin.PT3-GroupDMSO.PT3, 
  IsotretinoinvsDMSO_PT7=GroupIsotretinoin.PT7-GroupDMSO.PT7, 
  IsotretinoinvsDMSO_PT10=GroupIsotretinoin.PT10-GroupDMSO.PT10,



# treatment vs Isotretinoin-----------------------------------------------------


  #iso_TAZ vs DMSO
  Isotretinoin_and_TAZvsIsotretinoin_PT3 = GroupIsotretinoin_and_TAZ.PT3-GroupIsotretinoin.PT3, 
  Isotretinoin_and_TAZvsIsotretinoin_PT7=GroupIsotretinoin_and_TAZ.PT7-GroupIsotretinoin.PT7, 
  Isotretinoin_and_TAZvsIsotretinoin_PT10=GroupIsotretinoin_and_TAZ.PT10-GroupIsotretinoin.PT10,

  #iso+GSK vs DMSO
   Isotretinoin_and_GSKvsIsotretinoin_PT3 = GroupIsotretinoin_and_GSK.PT3-GroupIsotretinoin.PT3, 
  Isotretinoin_and_GSKvsIsotretinoin_PT7=GroupIsotretinoin_and_GSK.PT7-GroupIsotretinoin.PT7, 
  Isotretinoin_and_GSKvsIsotretinoin_PT10=GroupIsotretinoin_and_GSK.PT10-GroupIsotretinoin.PT10,

#GSK vs DMSO
    GSKvsIsotretinoin_PT3 = GroupGSK.PT3-GroupIsotretinoin.PT3, 
  GSKvsIsotretinoin_PT7=GroupGSK.PT7-GroupIsotretinoin.PT7, 
  GSKvsIsotretinoin_PT10=GroupGSK.PT10-GroupIsotretinoin.PT10,

#TAZ vs DMSO
TAZvsIsotretinoin_PT3 = GroupTAZ.PT3-GroupIsotretinoin.PT3, 
  TAZvsIsotretinoin_PT7=GroupTAZ.PT7-GroupIsotretinoin.PT7, 
  TAZvsIsotretinoin_PT10=GroupTAZ.PT10-GroupIsotretinoin.PT10,


# treatment vs TAZ ------------------------------------------------------

 #iso_TAZ vs TAZ
  Isotretinoin_and_TAZvsTAZ_PT3 = GroupIsotretinoin_and_TAZ.PT3-GroupTAZ.PT3, 
  Isotretinoin_and_TAZvsTAZ_PT7=GroupIsotretinoin_and_TAZ.PT7-GroupTAZ.PT7, 
  Isotretinoin_and_TAZvsTAZ_PT10=GroupIsotretinoin_and_TAZ.PT10-GroupTAZ.PT10,

  #iso+GSK vs TAZ
   Isotretinoin_and_GSKvsTAZ_PT3 = GroupIsotretinoin_and_GSK.PT3-GroupTAZ.PT3, 
  Isotretinoin_and_GSKvsTAZ_PT7=GroupIsotretinoin_and_GSK.PT7-GroupTAZ.PT7, 
  Isotretinoin_and_GSKvsTAZ_PT10=GroupIsotretinoin_and_GSK.PT10-GroupTAZ.PT10,

#GSK vs TAZ
    GSKvsTAZ_PT3 = GroupGSK.PT3-GroupTAZ.PT3, 
  GSKvsTAZ_PT7=GroupGSK.PT7-GroupTAZ.PT7, 
  GSKvsTAZ_PT10=GroupGSK.PT10-GroupTAZ.PT10,

# treatment vs GSK ------------------------------------------------------

 #iso_TAZ vs TAZ
  Isotretinoin_and_TAZvsGSK_PT3 = GroupIsotretinoin_and_TAZ.PT3-GroupGSK.PT3, 
  Isotretinoin_and_TAZvsGSK_PT7=GroupIsotretinoin_and_TAZ.PT7-GroupGSK.PT7, 
  Isotretinoin_and_TAZvsGSK_PT10=GroupIsotretinoin_and_TAZ.PT10-GroupGSK.PT10,

  #iso+GSK vs TAZ
   Isotretinoin_and_GSKvsGSK_PT3 = GroupIsotretinoin_and_GSK.PT3-GroupGSK.PT3, 
  Isotretinoin_and_GSKvsGSK_PT7=GroupIsotretinoin_and_GSK.PT7-GroupGSK.PT7, 
  Isotretinoin_and_GSKvsGSK_PT10=GroupIsotretinoin_and_GSK.PT10-GroupGSK.PT10,



levels=colnames(design))
  rownames(my.contrasts) <- gsub("treatment_timepoint", "", rownames(my.contrasts))
my.contrasts


#run GLFit

fit <- glmQLFit(edgeRData, design ,robust=TRUE)

#Now we test for significant differential expression in each gene using the QL F-test.
# First we check whether there was a genuine need to adjust for the experimental times. We
# do this by testing for differential expression between the three times. There is considerable
# differential expression, justifying our decision to adjust for the batch effect:

# qlf <- glmQLFTest(fit, coef=2:3)
# topTags(qlf)
# 
#  FDR <- p.adjust(qlf$table$PValue, method="BH")
#  sum(FDR < 0.05)
#  
#  top <- rownames(topTags(qlf))
# cpm(edgeRData)[top,]

 
#treatment vs control
Isotretinoin_and_TAZvsDMSO_PT3 <- glmQLFTest(fit, contrast=my.contrasts[,"Isotretinoin_and_TAZvsDMSO_PT3"])
Isotretinoin_and_TAZvsDMSO_PT7 <- glmQLFTest(fit, contrast=my.contrasts[,"Isotretinoin_and_TAZvsDMSO_PT7"])
Isotretinoin_and_TAZvsDMSO_PT10 <- glmQLFTest(fit, contrast=my.contrasts[,"Isotretinoin_and_TAZvsDMSO_PT10"])

Isotretinoin_and_GSKvsDMSO_PT3 <- glmQLFTest(fit, contrast=my.contrasts[,"Isotretinoin_and_GSKvsDMSO_PT3"])
Isotretinoin_and_GSKvsDMSO_PT7 <- glmQLFTest(fit, contrast=my.contrasts[,"Isotretinoin_and_GSKvsDMSO_PT7"])
Isotretinoin_and_GSKvsDMSO_PT10 <- glmQLFTest(fit, contrast=my.contrasts[,"Isotretinoin_and_GSKvsDMSO_PT10"])

GSKvsDMSO_PT3 <- glmQLFTest(fit, contrast=my.contrasts[,"GSKvsDMSO_PT3"])
GSKvsDMSO_PT7 <- glmQLFTest(fit, contrast=my.contrasts[,"GSKvsDMSO_PT7"])
GSKvsDMSO_PT10 <- glmQLFTest(fit, contrast=my.contrasts[,"GSKvsDMSO_PT10"])

TAZvsDMSO_PT3 <- glmQLFTest(fit, contrast=my.contrasts[,"TAZvsDMSO_PT3"])
TAZvsDMSO_PT7 <- glmQLFTest(fit, contrast=my.contrasts[,"TAZvsDMSO_PT7"])
TAZvsDMSO_PT10 <- glmQLFTest(fit, contrast=my.contrasts[,"TAZvsDMSO_PT10"])

IsotretinoinvsDMSO_PT3 <- glmQLFTest(fit, contrast=my.contrasts[,"IsotretinoinvsDMSO_PT3"])
IsotretinoinvsDMSO_PT7 <- glmQLFTest(fit, contrast=my.contrasts[,"IsotretinoinvsDMSO_PT7"])
IsotretinoinvsDMSO_PT10 <- glmQLFTest(fit, contrast=my.contrasts[,"IsotretinoinvsDMSO_PT10"])


#treatment vs RA

Isotretinoin_and_TAZvsIsotretinoin_PT3 <- glmQLFTest(fit, contrast=my.contrasts[,"Isotretinoin_and_TAZvsIsotretinoin_PT3"])
Isotretinoin_and_TAZvsIsotretinoin_PT7 <- glmQLFTest(fit, contrast=my.contrasts[,"Isotretinoin_and_TAZvsIsotretinoin_PT7"])
Isotretinoin_and_TAZvsIsotretinoin_PT10 <- glmQLFTest(fit, contrast=my.contrasts[,"Isotretinoin_and_TAZvsIsotretinoin_PT10"])

Isotretinoin_and_GSKvsIsotretinoin_PT3 <- glmQLFTest(fit, contrast=my.contrasts[,"Isotretinoin_and_GSKvsIsotretinoin_PT3"])
Isotretinoin_and_GSKvsIsotretinoin_PT7 <- glmQLFTest(fit, contrast=my.contrasts[,"Isotretinoin_and_GSKvsIsotretinoin_PT7"])
Isotretinoin_and_GSKvsIsotretinoin_PT10 <- glmQLFTest(fit, contrast=my.contrasts[,"Isotretinoin_and_GSKvsIsotretinoin_PT10"])

TAZvsIsotretinoin_PT3 <- glmQLFTest(fit, contrast=my.contrasts[,"TAZvsIsotretinoin_PT3"])
TAZvsIsotretinoin_PT7 <- glmQLFTest(fit, contrast=my.contrasts[,"TAZvsIsotretinoin_PT7"])
TAZvsIsotretinoin_PT10 <- glmQLFTest(fit, contrast=my.contrasts[,"TAZvsIsotretinoin_PT10"])

GSKvsIsotretinoin_PT3 <- glmQLFTest(fit, contrast=my.contrasts[,"GSKvsIsotretinoin_PT3"])
GSKvsIsotretinoin_PT7 <- glmQLFTest(fit, contrast=my.contrasts[,"GSKvsIsotretinoin_PT7"])
GSKvsIsotretinoin_PT10 <- glmQLFTest(fit, contrast=my.contrasts[,"GSKvsIsotretinoin_PT10"])


## treatment vs TAZ

Isotretinoin_and_TAZvsTAZ_PT3 <- glmQLFTest(fit, contrast=my.contrasts[,"Isotretinoin_and_TAZvsTAZ_PT3"])
Isotretinoin_and_TAZvsTAZ_PT7 <- glmQLFTest(fit, contrast=my.contrasts[,"Isotretinoin_and_TAZvsTAZ_PT7"])
Isotretinoin_and_TAZvsTAZ_PT10 <- glmQLFTest(fit, contrast=my.contrasts[,"Isotretinoin_and_TAZvsTAZ_PT10"])

Isotretinoin_and_GSKvsTAZ_PT3 <- glmQLFTest(fit, contrast=my.contrasts[,"Isotretinoin_and_GSKvsTAZ_PT3"])
Isotretinoin_and_GSKvsTAZ_PT7 <- glmQLFTest(fit, contrast=my.contrasts[,"Isotretinoin_and_GSKvsTAZ_PT7"])
Isotretinoin_and_GSKvsTAZ_PT10 <- glmQLFTest(fit, contrast=my.contrasts[,"Isotretinoin_and_GSKvsTAZ_PT10"])

GSKvsTAZ_PT3 <- glmQLFTest(fit, contrast=my.contrasts[,"GSKvsTAZ_PT3"])
GSKvsTAZ_PT7 <- glmQLFTest(fit, contrast=my.contrasts[,"GSKvsTAZ_PT7"])
GSKvsTAZ_PT10 <- glmQLFTest(fit, contrast=my.contrasts[,"GSKvsTAZ_PT10"])

# treatment vs GSK


Isotretinoin_and_TAZvsGSK_PT3 <- glmQLFTest(fit, contrast=my.contrasts[,"Isotretinoin_and_TAZvsGSK_PT3"])
Isotretinoin_and_TAZvsGSK_PT7 <- glmQLFTest(fit, contrast=my.contrasts[,"Isotretinoin_and_TAZvsGSK_PT7"])
Isotretinoin_and_TAZvsGSK_PT10 <- glmQLFTest(fit, contrast=my.contrasts[,"Isotretinoin_and_TAZvsGSK_PT10"])

Isotretinoin_and_GSKvsGSK_PT3 <- glmQLFTest(fit, contrast=my.contrasts[,"Isotretinoin_and_GSKvsGSK_PT3"])
Isotretinoin_and_GSKvsGSK_PT7 <- glmQLFTest(fit, contrast=my.contrasts[,"Isotretinoin_and_GSKvsGSK_PT7"])
Isotretinoin_and_GSKvsGSK_PT10 <- glmQLFTest(fit, contrast=my.contrasts[,"Isotretinoin_and_GSKvsGSK_PT10"])
```






###################### DGE- exact tests ###############################################



######## Isotretinoin_and_TAZ VS DMSO #####################################

#Isotretinoin_and_TAZ VS DMSO PT3
Isotretinoin_and_TAZ_DMSO_PT3 <- exactTest(edgeRData, pair=c("DMSO_PT3", "Isotretinoin_and_TAZ_PT3"))

#Isotretinoin_and_TAZ VS DMSO PT7
Isotretinoin_and_TAZ_DMSO_PT7 <- exactTest(edgeRData, pair=c("DMSO_PT7", "Isotretinoin_and_TAZ_PT7"))

#Isotretinoin_and_TAZ VS DMSO PT10
Isotretinoin_and_TAZ_DMSO_PT10 <- exactTest(edgeRData, pair=c("DMSO_PT10", "Isotretinoin_and_TAZ_PT10"))



######### Isotretinoin_and_GSK VS DMSO ####################

#Isotretinoin_and_GSK VS DMSO PT3
dge_Isotretinoin_and_GSK_DMSO_PT3 <- exactTest(edgeRData, pair = c("DMSO_PT3","Isotretinoin_and_GSK_PT3"))

#Isotretinoin_and_GSK VS DMSO PT7
dge_Isotretinoin_and_GSK_DMSO_PT7 <- exactTest(edgeRData, pair = c("DMSO_PT7","Isotretinoin_and_GSK_PT7"))

#Isotretinoin_and_GSK VS DMSO PT10
dge_Isotretinoin_and_GSK_DMSO_PT10 <- exactTest(edgeRData, pair=c("DMSO_PT10", "Isotretinoin_and_GSK_PT10"))

########### GSK VS DMSO ###########################################


#GSK VS DMSO PT3
dge_GSK_DMSO_PT3 <- exactTest(edgeRData, pair = c("DMSO_PT3","GSK_PT3"))

#GSK VS DMSO PT7
dge_GSK_DMSO_PT7 <- exactTest(edgeRData, pair = c("DMSO_PT7","GSK_PT7"))

#GSK VS DMSO PT10
dge_GSK_DMSO_PT10 <- exactTest(edgeRData, pair = c("DMSO_PT10","GSK_PT10"))


############ RA VS DMSO ##########################

#RA VS DMSO PT3
dge_RA_DMSO_PT3 <- exactTest(edgeRData, pair = c("DMSO_PT3","RA_PT3"))

#RA VS DMSO PT7
dge_RA_DMSO_PT7 <- exactTest(edgeRData, pair = c("DMSO_PT7","RA_PT7"))

#RA VS DMSO PT10
dge_RA_DMSO_PT10 <- exactTest(edgeRData, pair = c("DMSO_PT10","RA_PT10"))

############## TAZ VS DMSO ################################


#TAZ VS DMSO PT3
dge_TAZ_DMSO_PT3 <- exactTest(edgeRData, pair = c("DMSO_PT3","TAZ_PT3"))

#TAZ VS DMSO PT7
dge_TAZ_DMSO_PT7 <- exactTest(edgeRData, pair = c("DMSO_PT7","TAZ_PT7"))

#TAZ VS DMSO PT10
dge_TAZ_DMSO_PT10 <- exactTest(edgeRData, pair = c("DMSO_PT10","TAZ_PT10"))



#DGE visualisation


```{r}

dge_files <- list(Isotretinoin_and_TAZvsDMSO_PT3, Isotretinoin_and_TAZvsDMSO_PT7, Isotretinoin_and_TAZvsDMSO_PT10, Isotretinoin_and_GSKvsDMSO_PT3, Isotretinoin_and_GSKvsDMSO_PT7, Isotretinoin_and_GSKvsDMSO_PT10, GSKvsDMSO_PT3, GSKvsDMSO_PT7, GSKvsDMSO_PT10, IsotretinoinvsDMSO_PT3, IsotretinoinvsDMSO_PT7, IsotretinoinvsDMSO_PT10, TAZvsDMSO_PT3, TAZvsDMSO_PT7, TAZvsDMSO_PT10, Isotretinoin_and_TAZvsIsotretinoin_PT3, Isotretinoin_and_TAZvsIsotretinoin_PT7, Isotretinoin_and_TAZvsIsotretinoin_PT10, Isotretinoin_and_GSKvsIsotretinoin_PT3, Isotretinoin_and_GSKvsIsotretinoin_PT7, Isotretinoin_and_GSKvsIsotretinoin_PT10, TAZvsIsotretinoin_PT3, TAZvsIsotretinoin_PT7, TAZvsIsotretinoin_PT10, GSKvsIsotretinoin_PT3, GSKvsIsotretinoin_PT7, GSKvsIsotretinoin_PT10, Isotretinoin_and_TAZvsTAZ_PT3, Isotretinoin_and_TAZvsTAZ_PT7, Isotretinoin_and_TAZvsTAZ_PT10, Isotretinoin_and_GSKvsTAZ_PT3, Isotretinoin_and_GSKvsTAZ_PT7, Isotretinoin_and_GSKvsTAZ_PT10, GSKvsTAZ_PT3, GSKvsTAZ_PT7, GSKvsTAZ_PT10, Isotretinoin_and_TAZvsGSK_PT3, Isotretinoin_and_TAZvsGSK_PT7, Isotretinoin_and_TAZvsGSK_PT10, Isotretinoin_and_TAZvsGSK_PT3, Isotretinoin_and_TAZvsGSK_PT7, Isotretinoin_and_TAZvsGSK_PT10
               )


filenames <- c("Isotretinoin_and_TAZvsDMSO_PT3", "Isotretinoin_and_TAZvsDMSO_PT7", "Isotretinoin_and_TAZvsDMSO_PT10", "Isotretinoin_and_GSKvsDMSO_PT3", "Isotretinoin_and_GSKvsDMSO_PT7", "Isotretinoin_and_GSKvsDMSO_PT10", "GSKvsDMSO_PT3", "GSKvsDMSO_PT7", "GSKvsDMSO_PT10", "IsotretinoinvsDMSO_PT3", "IsotretinoinvsDMSO_PT7", "IsotretinoinvsDMSO_PT10", "TAZvsDMSO_PT3", "TAZvsDMSO_PT7", "TAZvsDMSO_PT10", "Isotretinoin_and_TAZvsIsotretinoin_PT3", "Isotretinoin_and_TAZvsIsotretinoin_PT7", "Isotretinoin_and_TAZvsIsotretinoin_PT10", "Isotretinoin_and_GSKvsIsotretinoin_PT3", "Isotretinoin_and_GSKvsIsotretinoin_PT7", "Isotretinoin_and_GSKvsIsotretinoin_PT10", "TAZvsIsotretinoin_PT3", "TAZvsIsotretinoin_PT7", "TAZvsIsotretinoin_PT10", "GSKvsIsotretinoin_PT3", "GSKvsIsotretinoin_PT7", "GSKvsIsotretinoin_PT10", "Isotretinoin_and_TAZvsTAZ_PT3", "Isotretinoin_and_TAZvsTAZ_PT7", "Isotretinoin_and_TAZvsTAZ_PT10", "Isotretinoin_and_GSKvsTAZ_PT3", "Isotretinoin_and_GSKvsTAZ_PT7", "Isotretinoin_and_GSKvsTAZ_PT10", "GSKvsTAZ_PT3", "GSKvsTAZ_PT7", "GSKvsTAZ_PT10", "Isotretinoin_and_TAZvsGSK_PT3", "Isotretinoin_and_TAZvsGSK_PT7", "Isotretinoin_and_TAZvsGSK_PT10", "Isotretinoin_and_TAZvsGSK_PT3", "Isotretinoin_and_TAZvsGSK_PT7", "Isotretinoin_and_TAZvsGSK_PT10")


#  dge_files <- list(TAZvsDMSO_PT3, TAZvsDMSO_PT7, TAZvsDMSO_PT10)
# #
# #
#  filenames <- c( "TAZvsDMSO_PT3", "TAZvsDMSO_PT7", "TAZvsDMSO_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)
  j <- filenames[i]
  title <- filenames[[i]]
  title <- gsub("_", " ",title)
  title <- gsub("vs", " vs ",title)
  #dge_sum <-summary(de <- decideTestsDGE(d, adjust.method = "BH", p.value = 0.05 ,lfc = 1.5))
  #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 ################
 
  # Make a basic volcano plot for up and down genes
  
  # p1 <- ggplot(data, aes(logFC, -log(FDR,10))) +
  #   geom_point(size = 2/5) +
  #   xlab(expression("log"[2]*"FC")) +
  #   ylab(expression("-log"[10]*"FDR"))

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

  # check if number of up/down genes match previous
 # print(d) 
  #up_down <- data %>%
   #count(Expression) %>%
    #kable()
  #print(up_down)
#  print(dge_sum)
  
  
 #save list of all DEGs
    # DEG_list <- data
    # j <- filenames[i]
    # write.csv(DEG_list, file=paste0("unfiltered_DEGs_",j ,".csv") )
   #  


  #save list of only signfi DEGs
   # DEG_list <- data[data$Expression == 'Up-regulated' | data$Expression == 'Down-regulated',]
   # j <- filenames[i]
   # write.csv(DEG_list, file=paste0("significant_DEGs_",j ,".csv") )

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

  
   

   ############# vol plot for highlighting specific genes  in gene_list ############
  #add expression
   vol_plot_data <- data %>% 
  mutate(
    Expression = case_when(logFC > 1.5 & 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("APOBEC3B")
# diff_geneset <- read_excel("~/Neuroblastoma RNAseq analysis/diff_gene_set.xlsx")
 # diff_geneset <- as.list(diff_geneset)
 # gene_list <- diff_geneset$NEURONAL_DIFF
  
 #  NK_ligand_geneset <- read_excel("~/NK ligands.xlsx")
 #  NK_ligand_geneset <- as.list(NK_ligand_geneset)
 # gene_list <- NK_ligand_geneset$`HLA-C`
 #  

# adr_geneset <- read_excel("~/adr_sig.xlsx")
# adr_geneset <- as.list(adr_geneset$`Supplementary Table 2. List of ADRN and MES signature genes.`)
# 
# gene_list <- adr_geneset
#   
  
#   gene_list <- c("SLC11A2", "SCN1B", "IRF7", "TLR2", "HAS2", "IL1R1", "BEST1", "P2RX7", "CALCRL", "RELA", "GABBR1", "OSMR", "CSF1", "EIF2AK2", "LPAR1", "CD40", "ACVR2A", "NLRP3", "ATP2B1", "ITGB3", "AXL", "SELL", "ATP2A2", "SPHK1", "P2RX4", "CSF3R", "SRI", "ABI1", "LIF", "TPBG", "RAF1", "IRF1", "HRH1", "ITGA5", "PTPRE", "PSEN1", "VIP", "BDKRB1", "ADRM1", "HIF1A", "PDE4B", "TNFSF9", "KLF6", "SELENOS", "ICAM4", "SLC28A2", "SEMA4D", "FZD5", "LYN", "PIK3R5", "MYC", "PDPN", "KCNMB2", "PVR", "SLC31A1", "EDN1", "ATP2C1", 
# "PLAUR", "LCK", "ITGB8", "SERPINE1", "MSR1", "RTP4", "GPC3", "TACR1", "RASGRP1", "CMKLR1", "CD55", "KCNA3", "ACVR1B", "IL6", "GCH1", "GP1BA", "RNF144B", "MET", "PROK2", "TIMP1", "IL4R", "KIF1B", "PCDH7", "NFKB1", "OPRK1", "STAB1", "ADGRE1", "MXD1", "SLC4A4", "IL15", "CD82", "NOD2", "LDLR", "CCL2", "TNFRSF9", "NAMPT", "TNFRSF1B", "ABCA1", "CX3CL1", "IL10RA", "HBEGF", "CD70", "SLC7A2", "IFNAR1", "LAMP3", "IFNGR2", "TACR3", "PTGIR", "BTG2", "IL2RB", "EBI3", "GNA15", "EMP3", "GNAI3", "PTGER2", "CDKN1A", 
# "MEFV", "SLC1A2", "MMP14", "P2RY2", "AHR", "NFKBIA", "NMUR1", "RIPK2", "NPFFR2", "RHOG", "CXCL8", "F3", "CXCL6", "ADM", "NMI", "KCNJ2", "DCBLD2", "ADORA2B", "ICAM1", "PTGER4", "CHST2", "BST2", "SLC7A1", "C5AR1", "IRAK2", "INHBA", "RGS16", "SGMS2")
  
 #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(title) +
    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, # Add labels last to appear as the top layer
                  
 print(p3)
  
 #ggsave(p3, file=paste0("Volcano_plot_", j,".png"), width = 10, height = 10, units = "cm")
 # dev.off()
  
 #  write.csv(ils, file=paste0("NK_ligand_Andys_", j,".csv") )
  # add colour by significance

   top <- 10
  top_genes <- bind_rows(
    data %>%
      filter(Expression == 'Up-regulated') %>%
      dplyr::arrange(FDR, desc(abs(logFC))) %>%
      head(top),
    data %>%
     filter(Expression == 'Down-regulated') %>%
      dplyr::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 = 3)
 
ggsave(p3, file=paste0("Volcano_plot_", j,".png"), width = 15, height = 10, units = "cm")
 
 
  data <- data %>%
    mutate(
      Significance = case_when(
        abs(logFC) >= 1.5 & FDR <= 0.05 & FDR > 0.01 ~ "FDR 0.05",
        abs(logFC) >= 1.5 & FDR <= 0.01 & FDR > 0.001 ~ "FDR 0.01",
        abs(logFC) >= 1.5 & FDR <= 0.001 ~ "FDR 0.001",
        TRUE ~ "Unchanged")
    )
  head(data) %>%
    kable()
  
 
  p4 <- ggplot(data, aes(logFC, -log(FDR,10))) +
    geom_point(aes(color = Significance), size = 2/5) +
    xlab(expression("log"[2]*"FC")) +
    ylab(expression("-log"[10]*"FDR")) +
    scale_color_viridis_d() +
    guides(colour = guide_legend(override.aes = list(size=1.5))) +
    ggtitle(filenames[[i]])
  # print(p3)
  
  
 # data %>%
 #   count(Expression, Significance) %>%
  #  kable()
  
  top <- 10
  top_genes <- bind_rows(
    data %>%
      filter(Expression == 'Up-regulated') %>%
      dplyr::arrange(FDR, desc(abs(logFC))) %>%
      head(top),
    data %>%
     filter(Expression == 'Down-regulated') %>%
      dplyr::arrange(FDR, desc(abs(logFC))) %>%
      head(top)
   
  )
  top_genes %>%
    kable()
  
  
  p4 <-  p4 +
    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)
#print(p4)
  #ggsave(p4, file=paste0("Volcano_plot_", j,".png"), width = 25, height = 25, units = "cm")
}


```


# Gene sig score


```{r}
#averaged log2(TPM + 1) expression of all genes within each gene set. Mesenchymal cell lines were called based on MES score ≥4.5 and ADRN score ≤5.

# To convert the MES and ADRN gene signatures into a single value per sample, we created a signature score as follows: within each sample, we ranked all genes according to their expression level (using HugoOnce to represent each gene by only one reporter). For each gene, we established the rank order in this list and calculated the percentile relative to the full list. The average of the percentiles for all MES or for all ADRN signature genes was used as a signature score. These values (signature scores) can readily be compared across samples. The method is implemented and available on the R2 platform.


#get gene lengths

#get GTF file to use to get gene lengths
hg38.gtf <- Rsubread::flattenGTF("hg38.refGene.gtf")
#get gene lengths from GTF file
GeneLength <- rowsum(hg38.gtf$End-hg38.gtf$Start+1, hg38.gtf$GeneID)
#remove BAGE5
row_names_df_to_remove <- c('BAGE5')
GeneLength <- GeneLength[!(row.names(GeneLength) %in% row_names_df_to_remove),]
#set to data frame
GeneLength <- as.data.frame(GeneLength)
#get order of gene ids in my data
gene_order <- row.names(rawData)
#order the gene lengths to the same order as my data
GeneLength <- GeneLength[ order(match(rownames(GeneLength), gene_order)), ]


#get mean fragment length for each sample
#change dir to insert size metrics
#read in insert size metrics.txt files
files <- list.files(pattern = ".txt")
#create empty list
mean_frag_length_list <-  list()

for (i in files){

#get the 8th line (containing mean frag length)
line8_insert_metric_sizes <- read.table(i, header = FALSE, sep = '\t' , skip = 7, nrows = 1)
#get the 6hth column (mean frag length)
mean_frag_length <- line8_insert_metric_sizes$V6
mean_frag_length_list[[length(mean_frag_length_list) + 1]] <- mean_frag_length    # Append new list element

}


#Convert counts to transcripts per million (TPM).
#' 
#' Convert a numeric matrix of features (rows) and conditions (columns) with
#' raw feature counts to transcripts per million.
#' 
#'    Lior Pachter. Models for transcript quantification from RNA-Seq.
#'    arXiv:1104.3889v2 
#'    
#'    Wagner, et al. Measurement of mRNA abundance using RNA-seq data:
#'    RPKM measure is inconsistent among samples. Theory Biosci. 24 July 2012.
#'    doi:10.1007/s12064-012-0162-3
#'    
#' @param counts A numeric matrix of raw feature counts i.e.
#'  fragments assigned to each gene.
#' @param featureLength A numeric vector with feature lengths.
#' @param meanFragmentLength A numeric vector with mean fragment lengths.
#' @return tpm A numeric matrix normalized by library size and feature length.

#get counts
counts <- rawData
#get feature length
featureLength <- GeneLength
#get mean fragment length
 meanFragmentLength <- mean_frag_length_list %>% as.numeric()
 

 
 
   counts_to_tpm <- function(counts, featureLength, meanFragmentLength) {
  
  # Ensure valid arguments.
  stopifnot(length(featureLength) == nrow(counts))
  stopifnot(length(meanFragmentLength) == ncol(counts))
  
  # Compute effective lengths of features in each library.
  effLen <- do.call(cbind, lapply(1:ncol(counts), function(i) {
    featureLength - meanFragmentLength[i] + 1
  }))
  
  # Exclude genes with length less than the mean fragment length.
  idx <- apply(effLen, 1, function(x) min(x) > 1)
  counts <- counts[idx,]
  effLen <- effLen[idx,]
  featureLength <- featureLength[idx]
  
  # Process one column at a time.
  tpm <- do.call(cbind, lapply(1:ncol(counts), function(i) {
    rate = log(counts[,i]) - log(effLen[,i])
    denom = log(sum(exp(rate)))
    exp(rate - denom + log(1e6))
  }))

  # Copy the row and column names from the original matrix.
  colnames(tpm) <- colnames(counts)
  rownames(tpm) <- rownames(counts)
  return(tpm)
}

tpm <- counts_to_tpm(counts, featureLength, meanFragmentLength)
   
#get list of genes in the score
ADR_sig <- read_xlsx("C:/Users/Christina/Documents/Neuroblastoma RNAseq analysis/GSEA custom gene sets/ADR_signature.xlsx")
ADR_genes <- ADR_sig$`Supplementary Table 2. List of ADRN and MES signature genes.`
MES_sig <- read_xlsx("C:/Users/Christina/Documents/Neuroblastoma RNAseq analysis/GSEA custom gene sets/MES_signature.xlsx")
MES_genes <- MES_sig$`Mes signature genes`

TPM_add_1 <- tpm+1
log2_TPM_add_1 <- log2(TPM_add_1)

#for ADR signtaure
#pull out columns with ADR genes
ADR_genes_expression <- log2_TPM_add_1 [rownames(log2_TPM_add_1) %in% ADR_genes, ]
average_ADR_expression <- colMeans(ADR_genes_expression) %>% as.data.frame()



# for MES signature
MES_genes_expression <- log2_TPM_add_1 [rownames(log2_TPM_add_1) %in% MES_genes, ]
average_MES_expression <- colMeans(MES_genes_expression) %>% as.data.frame()

merged_ADR_MES_score <- merge(average_MES_expression, average_ADR_expression,
                          by = 'row.names', all = TRUE)


#calculated based on an averaged log2(TPM + 1) expression of all genes within each gene set. Mesenchymal cell lines #were called based on MES score 
#≥4.5 and ADRN score ≤5.



sampleID <- merged_ADR_MES_score$Row.names
treatment <- c()
for (i in 1:nrow(merged_ADR_MES_score)) {
  samplename <- sampleID[i]
  group <- if(grepl("Isotretinoin_and_TAZ", samplename)) {print("Isotretinoin_and_TAZ")
  } else if(grepl("Isotretinoin_and_GSK", samplename)) {print("Isotretinoin_and_GSK")}
  else if(grepl("DMSO", samplename)) {print("DMSO")
  } else if(grepl("GSK", samplename)) {print("GSK")
  } else if(grepl("Isotretinoin", samplename)) {print("Isotretinoin")
  } else {print("TAZ")}
  treatment[i] <- group
}

#add treatment group to split by colour in ggplot
merged_ADR_MES_score$group <- treatment
merged_ADR_MES_score$Timepoint <- timepoint


merged_ADR_MES_score$Timepoint <- order(merged_ADR_MES_score$Timepoint, "PT3", "PT7", "PT10" )
reorder.factor(merged_ADR_MES_score$Timepoint)

merged_ADR_MES_score <- merged_ADR_MES_score %>% arrange(Timepoint)

#plot ADRN/MES scores
ggplot(merged_ADR_MES_score, aes(x=..x, y=..y, color = group)) +
    geom_point(size = 4) + ylab("ADRN score") + xlab("MES score")  #+ scale_color_viridis(discrete = TRUE, option = "D",)

#with shape by timpeoint
ggplot(merged_ADR_MES_score, aes(x=..x, y=..y, color = group))+ labs(fill='NEW LEGEND TITLE')  +
    geom_point(aes(shape = Timepoint), size = 4)  +
  xlab("Mesenchymal score") + ylab("Adrenergic score") + 
  guides(color = guide_legend(title = "Treatment")) + ggtitle("Adrenergic and mesenchymal scores")  +
        geom_smooth(method = "lm",se = FALSE,)
#+ scale_color_viridis(discrete = TRUE, option = "D",)

 colnames(merged_ADR_MES_score)[2] <- c("MES")
  colnames(merged_ADR_MES_score)[3] <- c("ADRN")


```

#ADRN/MES
```{r}
#get the adrenergic and mesenchymal gene lists
ADR_sig <- read_xlsx("C:/Users/Christina/Documents/Neuroblastoma RNAseq analysis/GSEA custom gene sets/ADR_signature.xlsx") 
ADR_genes <- ADR_sig$`Supplementary Table 2. List of ADRN and MES signature genes.` %>% list() 
MES_sig <- read_xlsx("C:/Users/Christina/Documents/Neuroblastoma RNAseq analysis/GSEA custom gene sets/MES_signature.xlsx")
MES_genes <- MES_sig$`Mes signature genes`  %>% list() 

#do GSVA on the genes to account for differences in the sizes of the lists
#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)

#change to matrix
filtered_mapped_matrix <- CPM %>%
  as.matrix()
#GSVA sets kcdf to gaussian as default, suitable for normalised gene expression data


#for ADRN-----------------------------------------------------------------------------------

#conduct GSVA
ADRN_scores <- gsva(
  filtered_mapped_matrix,
  ADR_genes,
  method = "gsva",
  # Minimum gene set size
  min.sz = 1,
  # Maximum gene set size
  max.sz = 500
)

#change rowname
rownames(ADRN_scores) <- "ADRN_score"
#order samples by score
ADRN_scores <- ADRN_scores[,order(ADRN_scores[1,], decreasing = T)] %>% as.data.frame()
View(ADRN_scores)

#get treatment  as extra col in df
sampleID <- row.names(ADRN_scores)
ADRN_scores$treatment <- c()
for (i in 1:nrow(ADRN_scores)) {
  samplename <- sampleID[i]
  group <- if(grepl("_and_TAZ", samplename)) {print("Isotretinoin_and_TAZ")
  } else if(grepl("_and_GSK", samplename)) {print("Isotretinoin_and_GSK")
    } else if (grepl("GSK", samplename)) {print("GSK")
  } else if(grepl("Isotretinoin", samplename)) {print("Isotretinoin")}
  else if(grepl("DMSO", samplename)) {print("DMSO")
  } 
  else {print("TAZ")}
  ADRN_scores$treatment[i] <- group
}

p=ggplot(data=ADRN_scores, mapping=aes(x=treatment, y=., fill=treatment))+geom_boxplot()

#stats------------------------------------------------------------------------------
#first check if the data is normally distributed
shapiro.test(ADRN_scores$.)

ADRN_scores <- as.data.frame(ADRN_scores)

compare_means(. ~ treatment, data = ADRN_scores, method = "t.test", p.adjust.method = "BH")

my_comparisons <- list( c("DMSO", "Isotretinoin"), c("DMSO", "TAZ"), c("DMSO", "GSK"), c("DMSO", "Isotretinoin_and_TAZ"), c("DMSO", "Isotretinoin_and_GSK"),c("Isotretinoin", "TAZ"),c("Isotretinoin", "GSK"),c("Isotretinoin", "Isotretinoin_and_TAZ"),c("Isotretinoin", "Isotretinoin_and_GSK"),c("TAZ", "GSK"),c("TAZ", "Isotretinoin_and_TAZ"),c("TAZ", "Isotretinoin_and_GSK"),c("GSK", "Isotretinoin_and_TAZ"),c("GSK", "Isotretinoin_and_GSK"),c("Isotretinoin_and_TAZ", "Isotretinoin_and_GSK"))


      ADRN_scores$treatment <- factor(ADRN_scores$treatment, levels=c('DMSO', "Isotretinoin", "TAZ", "GSK", "Isotretinoin_and_TAZ", "Isotretinoin_and_GSK") )
      

#plot results with statistical data
p=ggplot(data=ADRN_scores, mapping=aes(x=treatment, y=., fill=treatment))+geom_boxplot()+ theme( axis.text.x = element_text(angle = 90,  size = 12)) + stat_compare_means(comparisons = my_comparisons, label = "p.signif", method = "t.test") + scale_fill_discrete(name = "Treatment") + ggtitle("Adrenergic scores")+ ylab("GSVA score")
ggsave("adr_score_stats.png", p, dpi = 300, width = 5, height = 6.5)


#for MES-----------------------------------------------------------------------------------

#conduct GSVA
MES_scores <- gsva(
  filtered_mapped_matrix,
  MES_genes,
  method = "gsva",
  # Minimum gene set size
  min.sz = 1,
  # Maximum gene set size
  max.sz = 500
)
#change rowname
rownames(MES_scores) <- "MES_score"
#order samples by score
MES_scores <- MES_scores[,order(MES_scores[1,], decreasing = T)] %>% as.data.frame()
View(MES_scores)

#get treatment  as extra col in df
sampleID <- row.names(MES_scores)
MES_scores$treatment <- c()
for (i in 1:nrow(MES_scores)) {
  samplename <- sampleID[i]
  group <- if(grepl("_and_TAZ", samplename)) {print("Isotretinoin_and_TAZ")
  } else if(grepl("_and_GSK", samplename)) {print("Isotretinoin_and_GSK")
    } else if (grepl("GSK", samplename)) {print("GSK")
  } else if(grepl("Isotretinoin", samplename)) {print("Isotretinoin")}
  else if(grepl("DMSO", samplename)) {print("DMSO")
  } 
  else {print("TAZ")}
  MES_scores$treatment[i] <- group
}


      MES_scores$treatment <- factor(MES_scores$treatment, levels=c('DMSO', "Isotretinoin", "TAZ", "GSK", "Isotretinoin_and_TAZ", "Isotretinoin_and_GSK") )
      

#plot results with statistical data
p=ggplot(data=MES_scores, mapping=aes(x=treatment, y=., fill=treatment))+geom_boxplot()+ theme( axis.text.x = element_text(angle = 90,  size = 12)) + stat_compare_means(comparisons = my_comparisons, label = "p.signif", method = "t.test") + scale_fill_discrete(name = "Treatment") + ggtitle("Mesenchymal scores")+ ylab("GSVA score")
ggsave("mes_score_stats.png", p, dpi = 300, width = 5, height = 6.5)

#stats------------------------------------------------------------------------------
#first check if the data is normally distributed
shapiro.test(MES_scores$.)

compare_means(. ~ treatment, data = MES_scores, method = "t.test", p.adjust.method = "BH")

my_comparisons <- list( c("DMSO", "Isotretinoin"), c("DMSO", "TAZ"), c("DMSO", "GSK"), c("DMSO", "Isotretinoin_and_TAZ"), c("DMSO", "Isotretinoin_and_GSK"),c("Isotretinoin", "TAZ"),c("Isotretinoin", "GSK"),c("Isotretinoin", "Isotretinoin_and_TAZ"),c("Isotretinoin", "Isotretinoin_and_GSK"),c("TAZ", "GSK"),c("TAZ", "Isotretinoin_and_TAZ"),c("TAZ", "Isotretinoin_and_GSK"),c("GSK", "Isotretinoin_and_TAZ"),c("GSK", "Isotretinoin_and_GSK"),c("Isotretinoin_and_TAZ", "Isotretinoin_and_GSK"))


#plot results with statistical data
p=ggplot(data=MES_scores, mapping=aes(x=treatment, y=., fill=treatment))+geom_boxplot()+ theme( axis.text.x = element_text(angle = 90,  size = 12)) + stat_compare_means(comparisons = my_comparisons, label = "p.signif", method = "t.test") + ggtitle("Mesenchymal scores")


#get ratio of scores
ADRN_MES_ratio <- ADRN_scores/MES_scores
#rename score column
colnames(ADRN_MES_ratio)[1] <- "ADRN_MES_ratio"

#add extra column containing treatment group info for plotting

#get treatment  as extra col in df
sampleID <- row.names(ADRN_MES_ratio)
ADRN_MES_ratio$treatment <- c()
for (i in 1:nrow(ADRN_MES_ratio)) {
  samplename <- sampleID[i]
  group <- if(grepl("_and_TAZ", samplename)) {print("Isotretinoin_and_TAZ")
  } else if(grepl("_and_GSK", samplename)) {print("Isotretinoin_and_GSK")
    } else if (grepl("GSK", samplename)) {print("GSK")
  } else if(grepl("Isotretinoin", samplename)) {print("Isotretinoin")}
  else if(grepl("DMSO", samplename)) {print("DMSO")
  } 
  else {print("TAZ")}
  ADRN_MES_ratio$treatment[i] <- group
}


#get treatment  as extra col in df
sampleID <- row.names(ADRN_MES_ratio)
ADRN_MES_ratio$time <- c()
for (i in 1:nrow(ADRN_MES_ratio)) {
  samplename <- sampleID[i]
  group <- if(grepl("PT3", samplename)) {print("PT3")
  } else if(grepl("PT7", samplename)) {print("PT7")
    } else if (grepl("PT10", samplename)) {print("PT10")
  }
  ADRN_MES_ratio$time[i] <- group
}


#quick plot to look at scores
p=ggplot(data=ADRN_MES_ratio, mapping=aes(x=treatment, y=ADRN_MES_ratio, fill=treatment))+geom_boxplot()

#only run the code ## once otherwise it will continue to concatnate into the csvs. Once it has been run once you can just read in the csv

#normalise scores to DMSO and timepoint
#for PT3---------------------------------------------------------------------------------
# PT3 <- ADRN_MES_ratio %>% filter(time == "PT3") #filter for only pT3
# DMSO <- PT3 %>% filter(treatment == "DMSO") #get DMSO scores
# 
# t <- c("DMSO", "Isotretinoin", "TAZ", "GSK", "Isotretinoin_and_TAZ"  , "Isotretinoin_and_GSK"         )
# for (i in t) {
#   #filter for treatment
#   treatment_scores <- PT3 %>% filter(treatment == i)
#   #normalise scores to DMSO by doing treatment scores/DMSO scores
#   treatment_scores$normalised_scores <- treatment_scores$ADRN_MES_ratio/DMSO$ADRN_MES_ratio
#    #output results to table
#   write.table(treatment_scores, "Normalised_PT3.csv", append = TRUE,sep = ",",col.names = FALSE)
# }
ADRN_MES_ratio_PT3 <- (read.csv(file = "C:/Users/Christina/Documents/Neuroblastoma RNAseq analysis/ADRN_MES scores/Normalised_PT3.csv", header = FALSE))



#for PT7-------------------------------------------------------------------------------
# PT7 <- ADRN_MES_ratio %>% filter(time == "PT7") #filter for only pT7
# DMSO <- PT7 %>% filter(treatment == "DMSO")#get DMSO scores
# 
# t <- c("DMSO", "Isotretinoin", "TAZ", "GSK", "Isotretinoin_and_TAZ"  , "Isotretinoin_and_GSK"         )
# for (i in t) {
#   #filter for treatment
#   treatment_scores <- PT7 %>% filter(treatment == i)
#   #normalise scores to DMSO by doing treatment scores/DMSO scores
#   treatment_scores$normalised_scores <- treatment_scores$ADRN_MES_ratio/DMSO$ADRN_MES_ratio
#    #output results to table
#   write.table(treatment_scores, "Normalised_PT7.csv", append = TRUE,sep = ",",col.names = FALSE)
# }
ADRN_MES_ratio_PT7 <- (read.csv(file = "C:/Users/Christina/Documents/Neuroblastoma RNAseq analysis/ADRN_MES scores/Normalised_PT7.csv", header = FALSE))


# #for PT10-----------------------------------------------------------------------------
# PT10 <- ADRN_MES_ratio %>% filter(time == "PT10") #filter for only PT10
# DMSO <- PT10 %>% filter(treatment == "DMSO")#get DMSO scores
# 
# t <- c("DMSO", "Isotretinoin", "TAZ", "GSK", "Isotretinoin_and_TAZ"  , "Isotretinoin_and_GSK"         )
# for (i in t) {
#   #filter for treatment
#   treatment_scores <- PT10 %>% filter(treatment == i)
#   #normalise scores to DMSO by doing treatment scores/DMSO scores
#   treatment_scores$normalised_scores <- treatment_scores$ADRN_MES_ratio/DMSO$ADRN_MES_ratio
#    #output results to table
#   write.table(treatment_scores, "Normalised_PT10.csv", append = TRUE,sep = ",",col.names = FALSE)
# }
ADRN_MES_ratio_PT10 <- (read.csv(file = "C:/Users/Christina/Documents/Neuroblastoma RNAseq analysis/ADRN_MES scores/Normalised_PT10.csv", header = FALSE))

#concatenate results
ADRN_MES_ratio <- rbind(ADRN_MES_ratio_PT3, ADRN_MES_ratio_PT7, ADRN_MES_ratio_PT10)

#change colnames
colnames(ADRN_MES_ratio)[5] <- "Normalised_score"
colnames(ADRN_MES_ratio)[3] <- "Treatment"
colnames(ADRN_MES_ratio)[2] <- "ADRN_MES_ratio"
colnames(ADRN_MES_ratio)[1] <- "Sample"
colnames(ADRN_MES_ratio)[4] <- "Time"


 #add levels and reorder
  
      ADRN_MES_ratio$Treatment <- factor(ADRN_MES_ratio$Treatment, levels=c('DMSO', "Isotretinoin", "TAZ", "GSK", "Isotretinoin_and_TAZ", "Isotretinoin_and_GSK") )
      
     ADRN_MES_ratio$Time = factor(ADRN_MES_ratio$Time, levels=c('PT3','PT7','PT10'))


#plot with normalised scores
p=ggplot(data=ADRN_MES_ratio, mapping=aes(x=Treatment, y=Normalised_score, fill=Treatment))+geom_boxplot()+ theme( axis.text.x = element_text(angle = 90,  size = 12)) + ggtitle("Adrenergic and mesenchymal scores")

p <- p + facet_wrap(~Time, scales = "free_x" )  +
       theme(strip.text.x = element_text(size = 15))+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

# #normalise scores to DMSO
# #get DMSO scores
# DMSO <- ADRN_MES_ratio %>% filter(treatment == "DMSO")
# 
# #code below will get scores normalised to DMSO and append them to a csv. I've only ran this once as if it is run again then the data will be appended to the file again
# 
# # treatment <- c("DMSO", "Isotretinoin", "TAZ", "GSK", "Isotretinoin_and_TAZ"  , "Isotretinoin_and_GSK"         )
# # for (i in treatment) {
# #   #filter for treatment
# #   treatment_scores <- ADRN_MES_ratio %>% filter(treatment == i)
# #   #normalise scores to DMSO by doing treatment scores/DMSO scores
# #   treatment_scores$normalised_scores <- treatment_scores$ADRN_MES_ratio/DMSO$ADRN_MES_ratio
# #    #output results to table
# #   write.table(treatment_scores, "test.csv", append = TRUE,sep = ",",col.names = FALSE)
# # } 
# 
# #read in csv containing normalised scores
# ADRN_MES_ratio <- (read.csv(file = "C:/Users/Christina/Documents/Neuroblastoma RNAseq analysis/ADRN_MES scores/test.csv", header = FALSE))
# 
# #change colnames
# colnames(ADRN_MES_ratio)[4] <- "Normalised_score"
# colnames(ADRN_MES_ratio)[3] <- "Treatment"
# colnames(ADRN_MES_ratio)[2] <- "ADRN_MES_score"
# 
# #plot with normalised scores
# p=ggplot(data=ADRN_MES_ratio, mapping=aes(x=Treatment, y=Normalised_score, fill=Treatment))+geom_boxplot()+ theme( axis.text.x = element_text(angle = 90,  size = 12)) + ggtitle("Adrenergic and mesenchymal scores")





#do statistical tests-----------------------------------------------------------------------------
#add treatment timepoint column
ADRN_MES_ratio <- ADRN_MES_ratio %>% remove_rownames %>% column_to_rownames(var="Sample")


#get treatment  as extra col in df
sampleID <- row.names(ADRN_MES_ratio)
ADRN_MES_ratio$treatment_time <- c()
for (i in 1:nrow(ADRN_MES_ratio)) {
  samplename <- sampleID[i]
  group <- if(grepl("Isotretinoin\\_and_TAZ_PT3", samplename)) {print("Isotretinoin_and_TAZ_PT3")}
  else if(grepl("Isotretinoin\\_and_TAZ_PT7", samplename)) {print("Isotretinoin_and_TAZ_PT7")}
  else if(grepl("Isotretinoin\\_and_TAZ_PT10", samplename)) {print("Isotretinoin_and_TAZ_PT10")}
   else if(grepl("Isotretinoin\\_and_GSK_PT3", samplename)) {print("Isotretinoin_and_GSK_PT3")}
   else if(grepl("Isotretinoin\\_and_GSK_PT7", samplename)) {print("Isotretinoin_and_GSK_PT7")}
    else if(grepl("Isotretinoin\\_and_GSK_PT10", samplename)) {print("Isotretinoin_and_GSK_PT10")}
      else if(grepl("DMSO_PT3", samplename)) {print("DMSO_PT3")}
  else if(grepl("DMSO_PT7", samplename)) {print("DMSO_PT7")}
  else if(grepl("DMSO_PT10", samplename)) {print("DMSO_PT10")}
  else if(grepl("GSK_PT3", samplename)) {print("GSK_PT3")}
  else if(grepl("GSK_PT7", samplename)) {print("GSK_PT7")}
  else if(grepl("GSK_PT10", samplename)) {print("GSK_PT10")}
  else if(grepl("Isotretinoin_PT3", samplename)) {print("Isotretinoin_PT3")}
  else if(grepl("Isotretinoin_PT7", samplename)) {print("Isotretinoin_PT7")}
  else if(grepl("Isotretinoin_PT10", samplename)) {print("Isotretinoin_PT10")}
  else if(grepl("TAZ_PT3", samplename)) {print("TAZ_PT3")}
  else if(grepl("TAZ_PT7", samplename)) {print("TAZ_PT7")}
  else if(grepl("TAZ_PT10", samplename)) {print("TAZ_PT10")}
  ADRN_MES_ratio$treatment_time[i] <- group
}



#first check if the data is normally distributed
shapiro.test(ADRN_MES_ratio$ADRN_MES_ratio)
#p is less than 0.05 suggesting the data is significantly different from a normal distribution therefor a non-parametric test needs to be used

#we do non-parametric anova or kruskal test
kruskal.test(ADRN_MES_ratio ~ treatment_time, data = ADRN_MES_ratio)
#p is greater than 0.05 indicating that there is no significant difference between any groups


compare_means(ADRN_MES_ratio ~ treatment_time, data = ADRN_MES_ratio, method = "kruskal.test")


#plot results with statistical data
p=ggplot(data=ADRN_MES_ratio, mapping=aes(x=Treatment, y=Normalised_score, fill=Treatment))+geom_boxplot()+ theme( axis.text.x = element_text(angle = 90,  size = 12)) + ggtitle("Adrenergic and mesenchymal scores")

p <- p + facet_wrap(~Time, scales = "free_x" )  +
       theme(strip.text.x = element_text(size = 15))+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
 stat_compare_means(label.y = 3.5, label.x = 2) 





















#get the adrenergic and mesenchymal gene lists
ADR_sig <- read_xlsx("C:/Users/Christina/Documents/Neuroblastoma RNAseq analysis/GSEA custom gene sets/ADR_signature.xlsx") 
ADR_genes <- ADR_sig$`Supplementary Table 2. List of ADRN and MES signature genes.` %>% list() 
MES_sig <- read_xlsx("C:/Users/Christina/Documents/Neuroblastoma RNAseq analysis/GSEA custom gene sets/MES_signature.xlsx")
MES_genes <- MES_sig$`Mes signature genes`  %>% list() 

#do GSVA on the genes to account for differences in the sizes of the lists
#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)

#change to matrix
filtered_mapped_matrix <- CPM %>%
  as.matrix()
#GSVA sets kcdf to gaussian as default, suitable for normalised gene expression data


#for ADRN-----------------------------------------------------------------------------------

#conduct GSVA
ADRN_scores <- gsva(
  filtered_mapped_matrix,
  ADR_genes,
  method = "gsva",
  # Minimum gene set size
  min.sz = 1,
  # Maximum gene set size
  max.sz = 500
)

#change rowname
rownames(ADRN_scores) <- "ADRN_score"
#order samples by score
ADRN_scores <- ADRN_scores[,order(ADRN_scores[1,], decreasing = T)] %>% as.data.frame()
View(ADRN_scores)

#for MES-----------------------------------------------------------------------------------

#conduct GSVA
MES_scores <- gsva(
  filtered_mapped_matrix,
  MES_genes,
  method = "gsva",
  # Minimum gene set size
  min.sz = 1,
  # Maximum gene set size
  max.sz = 500
)
#change rowname
rownames(MES_scores) <- "MES_score"
#order samples by score
MES_scores <- MES_scores[,order(MES_scores[1,], decreasing = T)] %>% as.data.frame()
View(MES_scores)

#get ratio of scores
ADRN_MES_ratio <- ADRN_scores/MES_scores
#rename score column
colnames(ADRN_MES_ratio)[1] <- "ADRN_MES_ratio"

#add extra column containing treatment group info for plotting

#get treatment  as extra col in df
sampleID <- row.names(ADRN_MES_ratio)
ADRN_MES_ratio$treatment <- c()
for (i in 1:nrow(ADRN_MES_ratio)) {
  samplename <- sampleID[i]
  group <- if(grepl("_and_TAZ", samplename)) {print("Isotretinoin_and_TAZ")
  } else if(grepl("_and_GSK", samplename)) {print("Isotretinoin_and_GSK")
    } else if (grepl("GSK", samplename)) {print("GSK")
  } else if(grepl("Isotretinoin", samplename)) {print("Isotretinoin")}
  else if(grepl("DMSO", samplename)) {print("DMSO")
  } 
  else {print("TAZ")}
  ADRN_MES_ratio$treatment[i] <- group
}

#plot in a boxplot

#quick plot to look at scores
p=ggplot(data=ADRN_MES_ratio, mapping=aes(x=treatment, y=ADRN_MES_ratio, fill=treatment))+geom_boxplot()

#normalise scores to DMSO
#get DMSO scores
DMSO <- ADRN_MES_ratio %>% filter(treatment == "DMSO")

#code below will get scores normalised to DMSO and append them to a csv. I've only ran this once as if it is run again then the data will be appended to the file again

# treatment <- c("DMSO", "Isotretinoin", "TAZ", "GSK", "Isotretinoin_and_TAZ"  , "Isotretinoin_and_GSK"         )
# for (i in treatment) {
#   #filter for treatment
#   treatment_scores <- ADRN_MES_ratio %>% filter(treatment == i)
#   #normalise scores to DMSO by doing treatment scores/DMSO scores
#   treatment_scores$normalised_scores <- treatment_scores$ADRN_MES_ratio/DMSO$ADRN_MES_ratio
#    #output results to table
#   write.table(treatment_scores, "test.csv", append = TRUE,sep = ",",col.names = FALSE)
# } 

#read in csv containing normalised scores
ADRN_MES_ratio <- (read.csv(file = "C:/Users/Christina/Documents/Neuroblastoma RNAseq analysis/ADRN_MES scores/test.csv", header = FALSE))

#change colnames
colnames(ADRN_MES_ratio)[4] <- "Normalised_score"
colnames(ADRN_MES_ratio)[3] <- "Treatment"
colnames(ADRN_MES_ratio)[2] <- "ADRN_MES_score"

#plot with normalised scores
p=ggplot(data=ADRN_MES_ratio, mapping=aes(x=Treatment, y=Normalised_score, fill=Treatment))+geom_boxplot()+ theme( axis.text.x = element_text(angle = 90,  size = 12)) + ggtitle("Adrenergic and mesenchymal scores")

#do statistical tests-----------------------------------------------------------------------------

#first check if the data is normally distributed
shapiro.test(ADRN_MES_ratio$ADRN_MES_score)
#p is less than 0.05 suggesting the data is significantly different from a normal distribution therefor a non-parametric test needs to be used

#we do non-parametric anova or kruskal test
kruskal.test(ADRN_MES_score ~ Treatment, data = ADRN_MES_ratio)
#p is greater than 0.05 indicating that there is no significant difference between any groups

```





#sense check genes

```{r}

myCPM <- cpm(edgeRData$counts, log = FALSE)
 x <- myCPM['MYCN',] 
  m <- list(counts = as.numeric(x), group = as.factor(edgeRData$samples$group), timepoint =  as.factor(timepoint))
   m <- as.tibble(m)
 
   #rename PT
   m$timepoint <- sub("PT3", "Post-treatment day 3", m$timepoint)
    m$timepoint <- sub("PT7", "Post-treatment day 7", m$timepoint)
      m$timepoint <- sub("PT10", "Post-treatment day 10", m$timepoint)
      
      #add levels and reorder
     m$timepoint = factor(m$timepoint, levels=c('Post-treatment day 3','Post-treatment day 7','Post-treatment day 10'))
     
      m$group <- factor(m$group, levels=c('DMSO', "Isotretinoin", "TAZ", "GSK", "Isotretinoin_and_TAZ", "Isotretinoin_and_GSK") )

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

q <- ggplot(data = m) +
    geom_boxplot(aes(y = counts, x = group, fill= group),dotsize =3, 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)  + 
  xlab("Sample") + ylab("Log-2 count per million")     + guides(fill=guide_legend(title="Treatment"))

 
 q <- q + facet_wrap(~timepoint, scales = "free_x" )  +
       theme(strip.text.x = element_text(size = 14))+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +ggtitle(my_title)
 
 my_title <- expression(paste("Expression of ", italic("MYCN")))
 + labs(y=my_y_title)
 
   #ggsave(q, file=paste0("Sense_check_PHOX2B.png"), width = 25, height = 25, units = "cm")
   
```
    
    
    
    

#sense check genes boxplot and dotplot merged

```{r}

myCPM <- cpm(edgeRData$counts, log = FALSE)
 x <- myCPM['MYCN',] 
  m <- list(counts = as.numeric(x), group = as.factor(edgeRData$samples$group), timepoint =  as.factor(timepoint))
   m <- as.tibble(m)
 
   #rename PT
   m$timepoint <- sub("PT3", "Post-treatment day 3", m$timepoint)
    m$timepoint <- sub("PT7", "Post-treatment day 7", m$timepoint)
      m$timepoint <- sub("PT10", "Post-treatment day 10", m$timepoint)
      
      #add levels and reorder
     m$timepoint = factor(m$timepoint, levels=c('Post-treatment day 3','Post-treatment day 7','Post-treatment day 10'))
     
      m$group <- factor(m$group, levels=c('DMSO', "Isotretinoin", "TAZ", "GSK", "Isotretinoin_and_TAZ", "Isotretinoin_and_GSK") )

      
       my_title <- expression(paste("Expression of ", italic("MYCN")))

dot1 <- ggplot(m, aes( x = group, y = counts))
q <- dot1 + geom_boxplot(outlier.size = 0) +
  geom_dotplot(aes(y = counts, x = group, fill= group),dotsize =1500, 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)  + 
  xlab("Sample") + ylab("Log-2 count per million")     + guides(fill=guide_legend(title="Treatment"))+ facet_wrap(~timepoint, scales = "free_x" )  +
       theme(strip.text.x = element_text(size = 13))+ggtitle(my_title) +
  theme(legend.key.size = unit(0.9, 'cm'), #change legend key size
        legend.title = element_text(size=13), #change legend title font size
        legend.text = element_text(size=13)) +theme(axis.text.y =element_text(size=13),
        axis.title=element_text(size=13), axis.text.x = element_text(size =13, angle = 90, vjust = 0.5, hjust = 1))

      
     # ggsave(q, file=paste0("Sense_check_MYCN.png"), width = 25, height = 25, units = "cm")
      
   
```
    
    
    
    
#Venn PT3

```{r}
#find overlapping DEGs PT3

#ISO VS TAZ -------------------------------------------------------
#read in data
data <- read.csv("~/Neuroblastoma RNAseq analysis/Significant DEGs LFC 1.5/significant_DEGs_IsotretinoinvsDMSO_PT3.csv")

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

#select for upregulated DEGS
all_isotretinoin <- data$X

#read in data
  data <- read.csv("~/Neuroblastoma RNAseq analysis/Significant DEGs LFC 1.5/significant_DEGs_TAZvsDMSO_PT3.csv")
data <- subset(data, Expression == 'Up-regulated' | Expression == 'Down-regulated')

#select for upregulated DEGS
all_TAZ <- data$X

#find overlap and save as csv
overlapping_DEGs_iso_TAZ <- Reduce(intersect, list(all_TAZ,all_isotretinoin))
write.csv(overlapping_DEGs_iso_TAZ, "overlapping_DEGs_isovsDMSO_and_TAZvsDMSO_PT3.csv")


library(VennDiagram)
 
# Chart
venn.diagram(
  x = list(all_isotretinoin, all_TAZ),
  imagetype = "png",
  category.names = c("Isotretinoin" , "TAZ"),
  filename = 'isovsTAZ_venn_PT3.png',
  
  height= 3300,
  width = 3300,
  output = TRUE ,
          resolution = 500 ,
  col=c("darkgoldenrod1", 'springgreen1'),
  fill = c(alpha("darkgoldenrod1",0.3), alpha('springgreen1',0.3)) ,
  cat.default.pos = "outer"  ,
 cat.pos = c(-1, 1),
	inverted = TRUE,
cat.dist = c(0.125, 0.125),
	cat.cex = 1.4,	ext.line.lwd = 2,
 cex = 2,
 margin = 0.05)


#ISO VS GSK---------------------------------------------------------------
#read in data

#read in data
  data <- read.csv("~/Neuroblastoma RNAseq analysis/Significant DEGs LFC 1.5/significant_DEGs_GSKvsDMSO_PT3.csv")
data <- subset(data, Expression == 'Up-regulated' | Expression == 'Down-regulated')

#select for upregulated DEGS
all_GSK <- data$X

#find overlap and save as csv
overlapping_DEGs_iso_GSK <- Reduce(intersect, list(all_GSK,all_isotretinoin))
write.csv(overlapping_DEGs_iso_GSK, "overlapping_DEGs_isovsDMSO_and_GSKvsDMSO_PT3.csv")

# Chart
venn.diagram(
  x = list(all_isotretinoin, all_GSK),
  imagetype = "png",
  category.names = c("Isotretinoin" , "GSK"),
  filename = 'isovsGSK_venn_PT3.png',
  
  height= 3300,
  width = 3300,
  output = TRUE ,
          resolution = 500 ,
  col=c("darkgoldenrod1", 'turquoise3'),
  fill = c(alpha("darkgoldenrod1",0.3), alpha('turquoise3',0.3)) ,
  cat.default.pos = "outer"  ,
 cat.pos = c(-1, 1),
	inverted = TRUE,
cat.dist = c(0.125, 0.125),
	cat.cex = 1.4,	ext.line.lwd = 2,
 cex = 2,
 margin = 0.05)

#TAZ VS GSK---------------------------------------------------------------
#read in data

#read in data


#find overlap and save as csv
overlapping_DEGs_TAZ_GSK <- Reduce(intersect, list(all_GSK,all_TAZ))
write.csv(overlapping_DEGs_TAZ_GSK, "overlapping_DEGs_TAZvsDMSO_and_GSKvsDMSO_PT3.csv")

 
# Chart
venn.diagram(
  x = list(all_TAZ, all_GSK),
  imagetype = "png",
  category.names = c("TAZ" , "GSK"),
  filename = 'TAZvsGSK_venn_PT3.png',
  
  height= 3300,
  width = 3300,
  output = TRUE ,
          resolution = 500 ,
  col=c("springgreen1", 'turquoise3'),
  fill = c(alpha("springgreen1",0.3), alpha('turquoise3',0.3)) ,
  cat.default.pos = "outer"  ,
cat.dist = c(0.125, 0.125),
	cat.cex = 1.4,	ext.line.lwd = 2,
 cex = 2,
 margin = 0.05)



#ISO TAZ and ISO+TAZ---------------------------------------------------------------
#read in data

#read in data
#read in data
data <- read.csv("~/Neuroblastoma RNAseq analysis/Significant DEGs LFC 1.5/significant_DEGs_Isotretinoin_and_TAZvsDMSO_PT3.csv")

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

#select for upregulated DEGS
all_isotretinoin_and_TAZ <- data$X

#find overlap and save as csv
overlapping_DEGs_TAZ_iso_isotretinoin_and_TAZ <- Reduce(intersect, list(all_isotretinoin,all_TAZ, all_isotretinoin_and_TAZ))
write.csv(overlapping_DEGs_TAZ_iso_isotretinoin_and_TAZ, "overlapping_DEGs_TAZvsDMSO_and_ISOvsDMSO_and_isotretinoin_and_TAZvsDMSO_PT3.csv")

 
# Chart
venn.diagram(
  x = list(all_TAZ, all_isotretinoin, all_isotretinoin_and_TAZ),
  imagetype = "png",
  category.names = c("TAZ" , "ISO", "Isotrinoin + TAZ" ),
  filename = 'TAZvsISOvsiso_andTAZ_venn_PT3.png',
  
  height= 3300,
  width = 3300,
  output = TRUE ,
          resolution = 500 ,
  col=c("springgreen1", 'darkgoldenrod1', "dodgerblue1"),
  fill = c(alpha("springgreen1",0.3), alpha('darkgoldenrod1',0.3), alpha('dodgerblue1',0.3)) ,
  cat.default.pos = "outer"  ,
	cat.cex = 1.4,	ext.line.lwd = 2,
 cex = 2,
 margin = 0.05)



#ISO GSK and ISO+GSK--------------------------------------------------------------
#read in data

#read in data
#read in data
data <- read.csv("~/Neuroblastoma RNAseq analysis/Significant DEGs LFC 1.5/significant_DEGs_Isotretinoin_and_GSKvsDMSO_PT3.csv")

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

#select for upregulated DEGS
all_isotretinoin_and_GSK <- data$X

#find overlap and save as csv
overlapping_DEGs_GSK_iso_isotretinoin_and_GSK<- Reduce(intersect, list(all_isotretinoin,all_GSK, all_isotretinoin_and_GSK))
write.csv(overlapping_DEGs_GSK_iso_isotretinoin_and_GSK, "overlapping_DEGs_GSKvsDMSO_and_ISOvsDMSO_and_isotretinoin_and_GSKvsDMSO_PT3.csv")

 
# Chart
venn.diagram(
  x = list(all_GSK, all_isotretinoin, all_isotretinoin_and_GSK),
  imagetype = "png",
  category.names = c("GSK" , "ISO", "Isotrinoin + GSK" ),
  filename = 'GSKvsISOvsiso_andGSK_venn_PT3.png',
  
  height= 3300,
  width = 3300,
  output = TRUE ,
          resolution = 500 ,
  col=c("turquoise3", "darkgoldenrod1", "magenta1"),
  fill = c( alpha("turquoise3",0.3), alpha('darkgoldenrod1',0.3), alpha('magenta1',0.3)) ,
  cat.default.pos = "outer"  ,

	cat.cex = 1.4,	ext.line.lwd = 2,
 cex = 2,
 margin = 0.05)


```

      
#Venn PT7

```{r}
#find overlapping DEGs PT7

#ISO VS TAZ -------------------------------------------------------
#read in data
data <- read.csv("~/Neuroblastoma RNAseq analysis/Significant DEGs LFC 1.5/significant_DEGs_IsotretinoinvsDMSO_PT7.csv")

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

#select for upregulated DEGS
all_isotretinoin <- data$X

#read in data
  data <- read.csv("~/Neuroblastoma RNAseq analysis/Significant DEGs LFC 1.5/significant_DEGs_TAZvsDMSO_PT7.csv")
data <- subset(data, Expression == 'Up-regulated' | Expression == 'Down-regulated')

#select for upregulated DEGS
all_TAZ <- data$X

#find overlap and save as csv
overlapping_DEGs_iso_TAZ <- Reduce(intersect, list(all_TAZ,all_isotretinoin))
write.csv(overlapping_DEGs_iso_TAZ, "overlapping_DEGs_isovsDMSO_and_TAZvsDMSO_PT7.csv")


library(VennDiagram)
 
# Chart
venn.diagram(
  x = list(all_isotretinoin, all_TAZ),
  imagetype = "png",
  category.names = c("Isotretinoin" , "TAZ"),
  filename = 'isovsTAZ_venn_PT7.png',
  
  height= 3300,
  width = 3300,
  output = TRUE ,
          resolution = 500 ,
  col=c("darkgoldenrod1", 'springgreen1'),
  fill = c(alpha("darkgoldenrod1",0.3), alpha('springgreen1',0.3)) ,
  cat.default.pos = "outer"  ,
 cat.pos = c(-1, 1),
	inverted = TRUE,
cat.dist = c(0.125, 0.125),
	cat.cex = 1.4,	ext.line.lwd = 2,
 cex = 2,
 margin = 0.05)


#ISO VS GSK---------------------------------------------------------------
#read in data

#read in data
  data <- read.csv("~/Neuroblastoma RNAseq analysis/Significant DEGs LFC 1.5/significant_DEGs_GSKvsDMSO_PT7.csv")
data <- subset(data, Expression == 'Up-regulated' | Expression == 'Down-regulated')

#select for upregulated DEGS
all_GSK <- data$X

#find overlap and save as csv
overlapping_DEGs_iso_GSK <- Reduce(intersect, list(all_GSK,all_isotretinoin))
write.csv(overlapping_DEGs_iso_GSK, "overlapping_DEGs_isovsDMSO_and_GSKvsDMSO_PT7.csv")

# Chart
venn.diagram(
  x = list(all_isotretinoin, all_GSK),
  imagetype = "png",
  category.names = c("Isotretinoin" , "GSK"),
  filename = 'isovsGSK_venn_PT7.png',
  
  height= 3300,
  width = 3300,
  output = TRUE ,
          resolution = 500 ,
  col=c("darkgoldenrod1", 'turquoise3'),
  fill = c(alpha("darkgoldenrod1",0.3), alpha('turquoise3',0.3)) ,
  cat.default.pos = "outer"  ,
 cat.pos = c(-1, 1),
	inverted = TRUE,
cat.dist = c(0.125, 0.125),
	cat.cex = 1.4,	ext.line.lwd = 2,
 cex = 2,
 margin = 0.05)

#TAZ VS GSK---------------------------------------------------------------
#read in data

#read in data


#find overlap and save as csv
overlapping_DEGs_TAZ_GSK <- Reduce(intersect, list(all_GSK,all_TAZ))
write.csv(overlapping_DEGs_TAZ_GSK, "overlapping_DEGs_TAZvsDMSO_and_GSKvsDMSO_PT7.csv")

 
# Chart
venn.diagram(
  x = list(all_TAZ, all_GSK),
  imagetype = "png",
  category.names = c("TAZ" , "GSK"),
  filename = 'TAZvsGSK_venn_PT7.png',
  
  height= 3300,
  width = 3300,
  output = TRUE ,
          resolution = 500 ,
  col=c("springgreen1", 'turquoise3'),
  fill = c(alpha("springgreen1",0.3), alpha('turquoise3',0.3)) ,
  cat.default.pos = "outer"  ,
cat.dist = c(0.125, 0.125),
	cat.cex = 1.4,	ext.line.lwd = 2,
 cex = 2,
 margin = 0.05)



#ISO TAZ and ISO+TAZ---------------------------------------------------------------
#read in data

#read in data
#read in data
data <- read.csv("~/Neuroblastoma RNAseq analysis/Significant DEGs LFC 1.5/significant_DEGs_Isotretinoin_and_TAZvsDMSO_PT7.csv")

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

#select for upregulated DEGS
all_isotretinoin_and_TAZ <- data$X

#find overlap and save as csv
overlapping_DEGs_TAZ_iso_isotretinoin_and_TAZ <- Reduce(intersect, list(all_isotretinoin,all_TAZ, all_isotretinoin_and_TAZ))
write.csv(overlapping_DEGs_TAZ_iso_isotretinoin_and_TAZ, "overlapping_DEGs_TAZvsDMSO_and_ISOvsDMSO_and_isotretinoin_and_TAZvsDMSO_PT7.csv")

 
# Chart
venn.diagram(
  x = list(all_TAZ, all_isotretinoin, all_isotretinoin_and_TAZ),
  imagetype = "png",
  category.names = c("TAZ" , "ISO", "Isotrinoin + TAZ" ),
  filename = 'TAZvsISOvsiso_andTAZ_venn_PT7.png',
  
  height= 3300,
  width = 3300,
  output = TRUE ,
          resolution = 500 ,
  col=c("springgreen1", 'darkgoldenrod1', "dodgerblue1"),
  fill = c(alpha("springgreen1",0.3), alpha('darkgoldenrod1',0.3), alpha('dodgerblue1',0.3)) ,
  cat.default.pos = "outer"  ,
	cat.cex = 1.4,	ext.line.lwd = 2,
 cex = 2,
 margin = 0.05)



#ISO GSK and ISO+GSK--------------------------------------------------------------
#read in data

#read in data
#read in data
data <- read.csv("~/Neuroblastoma RNAseq analysis/Significant DEGs LFC 1.5/significant_DEGs_Isotretinoin_and_GSKvsDMSO_PT7.csv")

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

#select for upregulated DEGS
all_isotretinoin_and_GSK <- data$X

#find overlap and save as csv
overlapping_DEGs_GSK_iso_isotretinoin_and_GSK<- Reduce(intersect, list(all_isotretinoin,all_GSK, all_isotretinoin_and_GSK))
write.csv(overlapping_DEGs_GSK_iso_isotretinoin_and_GSK, "overlapping_DEGs_GSKvsDMSO_and_ISOvsDMSO_and_isotretinoin_and_GSKvsDMSO_PT7.csv")

 
# Chart
venn.diagram(
  x = list(all_GSK, all_isotretinoin, all_isotretinoin_and_GSK),
  imagetype = "png",
  category.names = c("GSK" , "ISO", "Isotrinoin + GSK" ),
  filename = 'GSKvsISOvsiso_andGSK_venn_PT7.png',
  
  height= 3300,
  width = 3300,
  output = TRUE ,
          resolution = 500 ,
  col=c("turquoise3", "darkgoldenrod1", "magenta1"),
  fill = c( alpha("turquoise3",0.3), alpha('darkgoldenrod1',0.3), alpha('magenta1',0.3)) ,
  cat.default.pos = "outer"  ,

	cat.cex = 1.4,	ext.line.lwd = 2,
 cex = 2,
 margin = 0.05)


```

        
#Venn PT10

```{r}
#find overlapping DEGs PT10

#ISO VS TAZ -------------------------------------------------------
#read in data
data <- read.csv("~/Neuroblastoma RNAseq analysis/Significant DEGs LFC 1.5/significant_DEGs_IsotretinoinvsDMSO_PT10.csv")

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

#select for upregulated DEGS
all_isotretinoin <- data$X

#read in data
  data <- read.csv("~/Neuroblastoma RNAseq analysis/Significant DEGs LFC 1.5/significant_DEGs_TAZvsDMSO_PT10.csv")
data <- subset(data, Expression == 'Up-regulated' | Expression == 'Down-regulated')

#select for upregulated DEGS
all_TAZ <- data$X

#find overlap and save as csv
overlapping_DEGs_iso_TAZ <- Reduce(intersect, list(all_TAZ,all_isotretinoin))
write.csv(overlapping_DEGs_iso_TAZ, "overlapping_DEGs_isovsDMSO_and_TAZvsDMSO_PT10.csv")


library(VennDiagram)
 
# Chart
venn.diagram(
  x = list(all_isotretinoin, all_TAZ),
  imagetype = "png",
  category.names = c("Isotretinoin" , "TAZ"),
  filename = 'isovsTAZ_venn_PT10.png',
  
  height= 3300,
  width = 3300,
  output = TRUE ,
          resolution = 500 ,
  col=c("darkgoldenrod1", 'springgreen1'),
  fill = c(alpha("darkgoldenrod1",0.3), alpha('springgreen1',0.3)) ,
  cat.default.pos = "outer"  ,
 cat.pos = c(-1, 1),
	inverted = TRUE,
cat.dist = c(0.125, 0.125),
	cat.cex = 1.4,	ext.line.lwd = 2,
 cex = 2,
 margin = 0.05)


#ISO VS GSK---------------------------------------------------------------
#read in data

#read in data
  data <- read.csv("~/Neuroblastoma RNAseq analysis/Significant DEGs LFC 1.5/significant_DEGs_GSKvsDMSO_PT10.csv")
data <- subset(data, Expression == 'Up-regulated' | Expression == 'Down-regulated')

#select for upregulated DEGS
all_GSK <- data$X

#find overlap and save as csv
overlapping_DEGs_iso_GSK <- Reduce(intersect, list(all_GSK,all_isotretinoin))
write.csv(overlapping_DEGs_iso_GSK, "overlapping_DEGs_isovsDMSO_and_GSKvsDMSO_PT10.csv")

# Chart
venn.diagram(
  x = list(all_isotretinoin, all_GSK),
  imagetype = "png",
  category.names = c("Isotretinoin" , "GSK"),
  filename = 'isovsGSK_venn_PT10.png',
  
  height= 3300,
  width = 3300,
  output = TRUE ,
          resolution = 500 ,
  col=c("darkgoldenrod1", 'turquoise3'),
  fill = c(alpha("darkgoldenrod1",0.3), alpha('turquoise3',0.3)) ,
  cat.default.pos = "outer"  ,
 cat.pos = c(-1, 1),
	inverted = TRUE,
cat.dist = c(0.125, 0.125),
	cat.cex = 1.4,	ext.line.lwd = 2,
 cex = 2,
 margin = 0.05)

#TAZ VS GSK---------------------------------------------------------------
#read in data

#read in data


#find overlap and save as csv
overlapping_DEGs_TAZ_GSK <- Reduce(intersect, list(all_GSK,all_TAZ))
write.csv(overlapping_DEGs_TAZ_GSK, "overlapping_DEGs_TAZvsDMSO_and_GSKvsDMSO_PT10.csv")

 
# Chart
venn.diagram(
  x = list(all_TAZ, all_GSK),
  imagetype = "png",
  category.names = c("TAZ" , "GSK"),
  filename = 'TAZvsGSK_venn_PT10.png',
  
  height= 3300,
  width = 3300,
  output = TRUE ,
          resolution = 500 ,
  col=c("springgreen1", 'turquoise3'),
  fill = c(alpha("springgreen1",0.3), alpha('turquoise3',0.3)) ,
  cat.default.pos = "outer"  ,
cat.dist = c(0.125, 0.125),
	cat.cex = 1.4,	ext.line.lwd = 2,
 cex = 2,
 margin = 0.05)



#ISO TAZ and ISO+TAZ---------------------------------------------------------------
#read in data

#read in data
#read in data
data <- read.csv("~/Neuroblastoma RNAseq analysis/Significant DEGs LFC 1.5/significant_DEGs_Isotretinoin_and_TAZvsDMSO_PT10.csv")

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

#select for upregulated DEGS
all_isotretinoin_and_TAZ <- data$X

#find overlap and save as csv
overlapping_DEGs_TAZ_iso_isotretinoin_and_TAZ <- Reduce(intersect, list(all_isotretinoin,all_TAZ, all_isotretinoin_and_TAZ))
write.csv(overlapping_DEGs_TAZ_iso_isotretinoin_and_TAZ, "overlapping_DEGs_TAZvsDMSO_and_ISOvsDMSO_and_isotretinoin_and_TAZvsDMSO_PT10.csv")

 
# Chart
venn.diagram(
  x = list(all_TAZ, all_isotretinoin, all_isotretinoin_and_TAZ),
  imagetype = "png",
  category.names = c("TAZ" , "ISO", "Isotrinoin + TAZ" ),
  filename = 'TAZvsISOvsiso_andTAZ_venn_PT10.png',
  
  height= 3300,
  width = 3300,
  output = TRUE ,
          resolution = 500 ,
  col=c("springgreen1", 'darkgoldenrod1', "dodgerblue1"),
  fill = c(alpha("springgreen1",0.3), alpha('darkgoldenrod1',0.3), alpha('dodgerblue1',0.3)) ,
  cat.default.pos = "outer"  ,
	cat.cex = 1.4,	ext.line.lwd = 2,
 cex = 2,
 margin = 0.05)



#ISO GSK and ISO+GSK--------------------------------------------------------------
#read in data

#read in data
#read in data
data <- read.csv("~/Neuroblastoma RNAseq analysis/Significant DEGs LFC 1.5/significant_DEGs_Isotretinoin_and_GSKvsDMSO_PT10.csv")

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

#select for upregulated DEGS
all_isotretinoin_and_GSK <- data$X

#find overlap and save as csv
overlapping_DEGs_GSK_iso_isotretinoin_and_GSK<- Reduce(intersect, list(all_isotretinoin,all_GSK, all_isotretinoin_and_GSK))
write.csv(overlapping_DEGs_GSK_iso_isotretinoin_and_GSK, "overlapping_DEGs_GSKvsDMSO_and_ISOvsDMSO_and_isotretinoin_and_GSKvsDMSO_PT10.csv")

 
# Chart
venn.diagram(
  x = list(all_GSK, all_isotretinoin, all_isotretinoin_and_GSK),
  imagetype = "png",
  category.names = c("GSK" , "ISO", "Isotrinoin + GSK" ),
  filename = 'GSKvsISOvsiso_andGSK_venn_PT10.png',
  
  height= 3300,
  width = 3300,
  output = TRUE ,
          resolution = 500 ,
  col=c("turquoise3", "darkgoldenrod1", "magenta1"),
  fill = c( alpha("turquoise3",0.3), alpha('darkgoldenrod1',0.3), alpha('magenta1',0.3)) ,
  cat.default.pos = "outer"  ,

	cat.cex = 1.4,	ext.line.lwd = 2,
 cex = 2,
 margin = 0.05)


```



# NK ligands
```{r}

genes <- c( "MICB", "ULBP1", "ULBP2", "ULBP3","ICAM1", "CD58", 
             "HLA-C", "CD274", "PVR", "NECTIN2", "CD276", 
             "FAS", "RAE1" )

for (i in genes) {

myCPM <- cpm(edgeRData$counts, log = FALSE)
 x <- myCPM[i,] 
 
  m <- list(counts = as.numeric(x), group = as.factor(edgeRData$samples$group), timepoint =  as.factor(timepoint))
   m <- as.tibble(m)
 
   #rename PT
   m$timepoint <- sub("PT3", "Post-treatment day 3", m$timepoint)
    m$timepoint <- sub("PT7", "Post-treatment day 7", m$timepoint)
      m$timepoint <- sub("PT10", "Post-treatment day 10", m$timepoint)
      
      #add levels and reorder
     m$timepoint = factor(m$timepoint, levels=c('Post-treatment day 3','Post-treatment day 7','Post-treatment day 10'))
     
      m$group <- factor(m$group, levels=c('DMSO', "Isotretinoin", "TAZ", "Isotretinoin_and_TAZ", "GSK",  "Isotretinoin_and_GSK") )

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

q <- ggplot(data = m) +
    geom_boxplot(aes(y = counts, x = group, fill= group),dotsize =0.7, stroke = NA,
                 binaxis = "y",         # which axis to bin along
                      # 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)  + 
  xlab("Sample") + ylab("Log-2 count per million")+ guides(fill=guide_legend(title="Treatment")) +
    ggtitle(i)

 q <- q + facet_wrap(~timepoint, scales = "free_x" )  +
       theme(strip.text.x = element_text(size = 15))+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
             
   print(q)
   ggsave(q, file=paste0("NK_ligand_", i, ".png" ), width = 25, height = 25, units = "cm")
}
   
```
    
    

#Number of DEGs

```{r}

#PT3 --------------------------------------------------------------------------

df <- read_xlsx("C:/Users/Christina/Documents/Neuroblastoma RNAseq analysis/No_DEGs_PT3.xlsx")

#reorder
df$Treatment <- factor(df$Treatment, levels=c("Isotretinoin", "TAZ", "GSK", 'Isotretinoin + TAZ', 'Isotretinoin + GSK'))

# 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 post-treatment day 3") + theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14))+
  scale_fill_manual(values=c('deepskyblue1', 'indianred1'))+ theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  ggsave(p, file=paste0("No_DEGs_PT3.png"), width = 25, height = 25, units = "cm")

  
# PT7-----------------------------------------------------------------

  df <- read_xlsx("C:/Users/Christina/Documents/Neuroblastoma RNAseq analysis/No_DEGs_PT7.xlsx")

#reorder
df$Treatment <- factor(df$Treatment, levels=c("Isotretinoin", "TAZ", "GSK", 'Isotretinoin + TAZ', 'Isotretinoin + GSK'))

# 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 post-treatment day 7")  + theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14))+
  scale_fill_manual(values=c('deepskyblue1', 'indianred1'))+ theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  ggsave(p, file=paste0("No_DEGs_PT7.png"), width = 25, height = 25, units = "cm")

  # PT10-----------------------------------------------------------------

  df <- read_xlsx("C:/Users/Christina/Documents/Neuroblastoma RNAseq analysis/No_DEGs_PT10.xlsx")

#reorder
df$Treatment <- factor(df$Treatment, levels=c("Isotretinoin", "TAZ", "GSK", 'Isotretinoin + TAZ', 'Isotretinoin + GSK'))

# 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 post-treatment day 10")  + theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14))+
  scale_fill_manual(values=c('deepskyblue1', 'indianred1'))+ theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  ggsave(p, file=paste0("No_DEGs_PT10.png"), width = 25, height = 25, units = "cm")

  
```

   

#Number of DEGs treatment vs treatment 

```{r}

#PT3 --------------------------------------------------------------------------

df <- read_xlsx("C:/Users/Christina/Documents/Neuroblastoma RNAseq analysis/No_DEGs_PT3_treatmentvstreatment.xlsx")

#reorder
df$Treatment <- factor(df$Treatment, levels=c("Isotretinoin + TAZ vs isotretinoin", "Isotretinoin + GSK vs isotretinoin", "Isotretinoin + TAZ vs TAZ", 'Isotretinoin + GSK vs GSK'))

# 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 for combination therapy vs single agent post-treatment day 3")+ theme(plot.title = element_text(size=20)) + theme( axis.text.x = element_text(size=20), axis.text.y = element_text(size=20),
        axis.title.y = element_text(size=20), axis.title.x = element_text(size=20),  legend.text = element_text(size=20),legend.title=element_text(size=20))+
  scale_fill_manual(values=c('deepskyblue1', 'indianred1'))+ theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  ggsave(p, file=paste0("No_DEGs_PT3.png"), width = 17, height = 26, units = "cm")

  
# PT7-----------------------------------------------------------------

df <- read_xlsx("C:/Users/Christina/Documents/Neuroblastoma RNAseq analysis/No_DEGs_PT7_treatmentvstreatment.xlsx")

#reorder
df$Treatment <- factor(df$Treatment, levels=c("Isotretinoin + TAZ vs isotretinoin", "Isotretinoin + GSK vs isotretinoin", "Isotretinoin + TAZ vs TAZ", 'Isotretinoin + GSK vs GSK'))

# 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 for combination therapy vs single agent post-treatment day 7")+ theme(plot.title = element_text(size=20)) + theme( axis.text.x = element_text(size=20), axis.text.y = element_text(size=20),
        axis.title.y = element_text(size=20), axis.title.x = element_text(size=20),  legend.text = element_text(size=20),legend.title=element_text(size=20))+
  scale_fill_manual(values=c('deepskyblue1', 'indianred1'))+ theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  ggsave(p, file=paste0("No_DEGs_PT7.png"), width = 17, height = 26, units = "cm")


  # PT10-----------------------------------------------------------------

 
df <- read_xlsx("C:/Users/Christina/Documents/Neuroblastoma RNAseq analysis/No_DEGs_PT10_treatmentvstreatment.xlsx")

#reorder
df$Treatment <- factor(df$Treatment, levels=c("Isotretinoin + TAZ vs isotretinoin", "Isotretinoin + GSK vs isotretinoin", "Isotretinoin + TAZ vs TAZ", 'Isotretinoin + GSK vs GSK'))

# 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 for combination therapy vs single agent post-treatment day 10")+ theme(plot.title = element_text(size=20)) + theme( axis.text.x = element_text(size=20), axis.text.y = element_text(size=20),
        axis.title.y = element_text(size=20), axis.title.x = element_text(size=20),  legend.text = element_text(size=20),legend.title=element_text(size=20))+
  scale_fill_manual(values=c('deepskyblue1', 'indianred1'))+ theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  ggsave(p, file=paste0("No_DEGs_PT10.png"), width = 17, height = 26, units = "cm")

  
```


#Venn diagram
```{r}
#edit contrast matrix to treatment vs dmso only





venn <- makeVenn(fit, my.contrasts, design, lfc = "1.5", p.value = 0.05, method = "sameup", adj.meth = "BH", )


```





#Scatterplot comparing LFC for TAZ and ra


```{r}
#PT10=====================================================


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

df <- df %>%
   mutate(sig = case_when(Expression_TAZ == "significant" & Expression_RA == "significant" ~ "Both",
                          Expression_TAZ == "significant" ~ "TAZ",
                            Expression_RA == "significant" ~ "Isotretinoin",
                          
                           TRUE ~ "Unchanged") )
   
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)

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("yellow2", "darkgoldenrod1", "seagreen2", "grey")) + ggtitle("Comparision of DEGs for TAZ vs DMSO and isotretinoin vs DMSO PT10")+ 
  theme(plot.title = element_text(size = 15)) +
  xlab("Log-2 fold change TAZ vs DMSO") + ylab("Log-2 fold change isotretinoin vs DMSO PT10") +theme(axis.text=element_text(size=15),
        axis.title=element_text(size=15), legend.text = element_text(size=14)) + guides(color=guide_legend(title="Treatment differentially expressed in"))+ guides(colour = guide_legend(override.aes = list(size=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)


#PT3 --------------------------------------------------------------------

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

df <- df %>%
   mutate(sig = case_when(Expression_TAZ == "significant" & Expression_RA == "significant" ~ "Both",
                          Expression_TAZ == "significant" ~ "TAZ",
                            Expression_RA == "significant" ~ "Isotretinoin",
                          
                           TRUE ~ "Unchanged") )
   
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)

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("yellow2", "darkgoldenrod1", "seagreen2", "grey")) + ggtitle("Comparision of DEGs for TAZ vs DMSO and isotretinoin vs DMSO PT3") + 
  theme(plot.title = element_text(size = 15))+
  xlab("Log-2 fold change TAZ vs DMSO") + ylab("Log-2 fold change isotretinoin vs DMSO PT3") +theme(axis.text=element_text(size=15),
        axis.title=element_text(size=15)) + guides(color=guide_legend(title="Treatment differentially expressed in"))

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


#PT7 --------------------------------------------------------------------


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

df <- df %>%
   mutate(sig = case_when(Expression_TAZ == "significant" & Expression_RA == "significant" ~ "Both",
                          Expression_TAZ == "significant" ~ "TAZ",
                            Expression_RA == "significant" ~ "Isotretinoin",
                          
                           TRUE ~ "Unchanged") )
   
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)

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("yellow2", "darkgoldenrod1", "seagreen2", "grey")) + ggtitle("Comparision of DEGs for TAZ vs DMSO and isotretinoin vs DMSO PT7") + 
  theme(plot.title = element_text(size = 15)) +
  xlab("Log-2 fold change TAZ vs DMSO") + ylab("Log-2 fold change isotretinoin vs DMSO PT7") +theme(axis.text=element_text(size=15),
        axis.title=element_text(size=15)) + guides(color=guide_legend(title="Treatment differentially expressed in"))

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

```


#Scatterplot comparing LFC for TAZ and GSK


```{r}
#PT10=====================================================


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

df <- df %>%
   mutate(sig = case_when(Expression_TAZ == "significant" & Expression_RA == "significant" ~ "Both",
                          Expression_TAZ == "significant" ~ "TAZ",
                            Expression_RA == "significant" ~ "GSK",
                          
                           TRUE ~ "Unchanged") )
   
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)

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("yellow2", "turquoise3", "seagreen2", "grey")) + ggtitle("Comparision of DEGs for TAZ vs DMSO and GSK vs DMSO PT10") + 
  theme(plot.title = element_text(size = 15)) +
  xlab("Log-2 fold change TAZ vs DMSO") + ylab("Log-2 fold change GSK vs DMSO PT10") +theme(axis.text=element_text(size=15),
        axis.title=element_text(size=15)) + guides(color=guide_legend(title="Treatment differentially expressed in"))

#ggsave(scatterplot, file=paste0("Scatterplot_TAZ_GSK_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)


#PT3 --------------------------------------------------------------------

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

df <- df %>%
   mutate(sig = case_when(Expression_TAZ == "significant" & Expression_RA == "significant" ~ "Both",
                          Expression_TAZ == "significant" ~ "TAZ",
                            Expression_RA == "significant" ~ "GSK",
                          
                           TRUE ~ "Unchanged") )
   
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)

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("yellow2", "turquoise3", "seagreen2", "grey")) + ggtitle("Comparision of DEGs for TAZ vs DMSO and GSK vs DMSO PT3") + 
  theme(plot.title = element_text(size = 15))+
  xlab("Log-2 fold change TAZ vs DMSO") + ylab("Log-2 fold change GSK vs DMSO PT3") +theme(axis.text=element_text(size=15),
        axis.title=element_text(size=15)) + guides(color=guide_legend(title="Treatment differentially expressed in"))

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


#PT7 --------------------------------------------------------------------


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

df <- df %>%
   mutate(sig = case_when(Expression_TAZ == "significant" & Expression_RA == "significant" ~ "Both",
                          Expression_TAZ == "significant" ~ "TAZ",
                            Expression_RA == "significant" ~ "GSK",
                          
                           TRUE ~ "Unchanged") )
   
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)

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("yellow2", "turquoise3", "seagreen2", "grey")) + ggtitle("Comparision of DEGs for TAZ vs DMSO and GSK vs DMSO PT7")  + 
  theme(plot.title = element_text(size = 15))+
  xlab("Log-2 fold change TAZ vs DMSO") + ylab("Log-2 fold change GSK vs DMSO PT7") +theme(axis.text=element_text(size=15),
        axis.title=element_text(size=14), legend.text = element_text(size=15)) + guides(color=guide_legend(title="Treatment differentially expressed in")) + guides(colour = guide_legend(override.aes = list(size=6)))

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


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

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(1234)
em2 <- try(GSEA(geneList, TERM2GENE = t2g, nPermSimple = 100000, seed = TRUE))
head(em2)

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

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

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

print(dotplot)

ggsave(dotplot, file=paste0("~/Neuroblastoma RNAseq analysis/Cluster profiler/dotplots/", msigdb_category, "_", j ,".png"), width = 25, height = 23, units = "cm")

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


#ggsave(core_plot, file=paste0("~/Neuroblastoma RNAseq analysis/Cluster profiler/core plots/",msigdb_category, "_", j ,"_coreplot.png"), width = 30, height = 30, units = "cm", bg = "white")

#write results to csv
#try(write.csv(em2, file=paste0("~/Neuroblastoma RNAseq analysis/Cluster profiler/tables/",msigdb_category, "_", j ,".csv")))

# #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(Isotretinoin_and_TAZvsDMSO_PT3, Isotretinoin_and_TAZvsDMSO_PT7, Isotretinoin_and_TAZvsDMSO_PT10, Isotretinoin_and_GSKvsDMSO_PT3, Isotretinoin_and_GSKvsDMSO_PT7, Isotretinoin_and_GSKvsDMSO_PT10, GSKvsDMSO_PT3, GSKvsDMSO_PT7, GSKvsDMSO_PT10, IsotretinoinvsDMSO_PT3, IsotretinoinvsDMSO_PT7, IsotretinoinvsDMSO_PT10, TAZvsDMSO_PT3, TAZvsDMSO_PT7, TAZvsDMSO_PT10, Isotretinoin_and_TAZvsIsotretinoin_PT3, Isotretinoin_and_TAZvsIsotretinoin_PT7, Isotretinoin_and_TAZvsIsotretinoin_PT10, Isotretinoin_and_GSKvsIsotretinoin_PT3, Isotretinoin_and_GSKvsIsotretinoin_PT7, Isotretinoin_and_GSKvsIsotretinoin_PT10, TAZvsIsotretinoin_PT3, TAZvsIsotretinoin_PT7, TAZvsIsotretinoin_PT10, GSKvsIsotretinoin_PT3, GSKvsIsotretinoin_PT7, GSKvsIsotretinoin_PT10, Isotretinoin_and_TAZvsTAZ_PT3, Isotretinoin_and_TAZvsTAZ_PT7, Isotretinoin_and_TAZvsTAZ_PT10, Isotretinoin_and_GSKvsTAZ_PT3, Isotretinoin_and_GSKvsTAZ_PT7, Isotretinoin_and_GSKvsTAZ_PT10, GSKvsTAZ_PT3, GSKvsTAZ_PT7, GSKvsTAZ_PT10, Isotretinoin_and_TAZvsGSK_PT3, Isotretinoin_and_TAZvsGSK_PT7, Isotretinoin_and_TAZvsGSK_PT10, Isotretinoin_and_TAZvsGSK_PT3, Isotretinoin_and_TAZvsGSK_PT7, Isotretinoin_and_TAZvsGSK_PT10
#                )

# 
# filenames <- c("Isotretinoin_and_TAZvsDMSO_PT3", "Isotretinoin_and_TAZvsDMSO_PT7", "Isotretinoin_and_TAZvsDMSO_PT10", "Isotretinoin_and_GSKvsDMSO_PT3", "Isotretinoin_and_GSKvsDMSO_PT7", "Isotretinoin_and_GSKvsDMSO_PT10", "GSKvsDMSO_PT3", "GSKvsDMSO_PT7", "GSKvsDMSO_PT10", "IsotretinoinvsDMSO_PT3", "IsotretinoinvsDMSO_PT7", "IsotretinoinvsDMSO_PT10", "TAZvsDMSO_PT3", "TAZvsDMSO_PT7", "TAZvsDMSO_PT10", "Isotretinoin_and_TAZvsIsotretinoin_PT3", "Isotretinoin_and_TAZvsIsotretinoin_PT7", "Isotretinoin_and_TAZvsIsotretinoin_PT10", "Isotretinoin_and_GSKvsIsotretinoin_PT3", "Isotretinoin_and_GSKvsIsotretinoin_PT7", "Isotretinoin_and_GSKvsIsotretinoin_PT10", "TAZvsIsotretinoin_PT3", "TAZvsIsotretinoin_PT7", "TAZvsIsotretinoin_PT10", "GSKvsIsotretinoin_PT3", "GSKvsIsotretinoin_PT7", "GSKvsIsotretinoin_PT10", "Isotretinoin_and_TAZvsTAZ_PT3", "Isotretinoin_and_TAZvsTAZ_PT7", "Isotretinoin_and_TAZvsTAZ_PT10", "Isotretinoin_and_GSKvsTAZ_PT3", "Isotretinoin_and_GSKvsTAZ_PT7", "Isotretinoin_and_GSKvsTAZ_PT10", "GSKvsTAZ_PT3", "GSKvsTAZ_PT7", "GSKvsTAZ_PT10", "Isotretinoin_and_TAZvsGSK_PT3", "Isotretinoin_and_TAZvsGSK_PT7", "Isotretinoin_and_TAZvsGSK_PT10", "Isotretinoin_and_TAZvsGSK_PT3", "Isotretinoin_and_TAZvsGSK_PT7", "Isotretinoin_and_TAZvsGSK_PT10")

setwd("~/Neuroblastoma RNAseq analysis/Cluster profiler")


dge_files <- c(TAZvsDMSO_PT3, TAZvsDMSO_PT7, TAZvsDMSO_PT10)
filenames <- c("TAZvsDMSO_PT3", "TAZvsDMSO_PT7", "TAZvsDMSO_PT10")
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")

}

```



# GSEA REACTOME with geneid
```{r}
################################ REACTOME #################################

#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.
setwd("~/Neuroblastoma RNAseq analysis/GSEA REACTOME")
#treatment vs control
filenames <- c("Isotretinoin_and_TAZvsDMSO_PT3", "Isotretinoin_and_TAZvsDMSO_PT7", "Isotretinoin_and_TAZvsDMSO_PT10", "Isotretinoin_and_GSKvsDMSO_PT3", "Isotretinoin_and_GSKvsDMSO_PT7", "Isotretinoin_and_GSKvsDMSO_PT10", "GSKvsDMSO_PT3", "GSKvsDMSO_PT7", "GSKvsDMSO_PT10", "IsotretinoinvsDMSO_PT3", "IsotretinoinvsDMSO_PT7", "IsotretinoinvsDMSO_PT10", "TAZvsDMSO_PT3", "TAZvsDMSO_PT7", "TAZvsDMSO_PT10")

#treatment vs treatment 
filenames <- c("Isotretinoin_and_TAZvsIsotretinoin_PT3", "Isotretinoin_and_TAZvsIsotretinoin_PT7", "Isotretinoin_and_TAZvsIsotretinoin_PT10", "Isotretinoin_and_GSKvsIsotretinoin_PT3", "Isotretinoin_and_GSKvsIsotretinoin_PT7", "Isotretinoin_and_GSKvsIsotretinoin_PT10", "Isotretinoin_and_TAZvsTAZ_PT3", "Isotretinoin_and_TAZvsTAZ_PT7", "Isotretinoin_and_TAZvsTAZ_PT10", "Isotretinoin_and_GSKvsGSK_PT3", "Isotretinoin_and_GSKvsGSK_PT7", "Isotretinoin_and_GSKvsGSK_PT10")

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$pathway <- sub(" ", "_", gsea_df$pathway)
gsea_df<-gsea_df[order(gsea_df$FDR),]
gsea_df$pathway <- gsub("_", " ", gsea_df$pathway )


gsea_df$pathway = str_wrap(gsea_df$pathway, width = 70)
gsea_df

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

    }



### 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 with geneid

```{r}
################################ Hallmarks #################################

#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("Isotretinoin_and_TAZvsDMSO_PT3", "Isotretinoin_and_TAZvsDMSO_PT7", "Isotretinoin_and_TAZvsDMSO_PT10", "Isotretinoin_and_GSKvsDMSO_PT3", "Isotretinoin_and_GSKvsDMSO_PT7", "Isotretinoin_and_GSKvsDMSO_PT10", "GSKvsDMSO_PT3", "GSKvsDMSO_PT7", "GSKvsDMSO_PT10", "IsotretinoinvsDMSO_PT3", "IsotretinoinvsDMSO_PT7", "IsotretinoinvsDMSO_PT10", "TAZvsDMSO_PT3", "TAZvsDMSO_PT7", "TAZvsDMSO_PT10")

#treatment vs treatment 
filenames <- c("Isotretinoin_and_TAZvsIsotretinoin_PT3", "Isotretinoin_and_TAZvsIsotretinoin_PT7", "Isotretinoin_and_TAZvsIsotretinoin_PT10", "Isotretinoin_and_GSKvsIsotretinoin_PT3", "Isotretinoin_and_GSKvsIsotretinoin_PT7", "Isotretinoin_and_GSKvsIsotretinoin_PT10", "Isotretinoin_and_TAZvsTAZ_PT3", "Isotretinoin_and_TAZvsTAZ_PT7", "Isotretinoin_and_TAZvsTAZ_PT10", "Isotretinoin_and_GSKvsGSK_PT3", "Isotretinoin_and_GSKvsGSK_PT7", "Isotretinoin_and_GSKvsGSK_PT10")


for (i in filenames){
# Get pathways from gmt file
Hallmarks_pathways <- gmtPathways("C:/Users/Christina/Documents/GSEA .gmt/h.all.v2022.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_readable <- as.character(row.names(gsea_df))
gsea_df$pathway_readable <- gsub("_", " ", gsea_df$pathway_readable )

#create barplot of enriched pathways:
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=15))
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}


setwd("~/Neuroblastoma RNAseq analysis/GSEA GO terms")
################################ 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("Isotretinoin_and_TAZvsDMSO_PT3", "Isotretinoin_and_TAZvsDMSO_PT7", "Isotretinoin_and_TAZvsDMSO_PT10", "Isotretinoin_and_GSKvsDMSO_PT3", "Isotretinoin_and_GSKvsDMSO_PT7", "Isotretinoin_and_GSKvsDMSO_PT10", "GSKvsDMSO_PT3", "GSKvsDMSO_PT7", "GSKvsDMSO_PT10", "IsotretinoinvsDMSO_PT3", "IsotretinoinvsDMSO_PT7", "IsotretinoinvsDMSO_PT10", "TAZvsDMSO_PT3", "TAZvsDMSO_PT7", "TAZvsDMSO_PT10")

#treatment vs treatment 
#filenames <- c("Isotretinoin_and_TAZvsIsotretinoin_PT3", "Isotretinoin_and_TAZvsIsotretinoin_PT7", "Isotretinoin_and_TAZvsIsotretinoin_PT10", "Isotretinoin_and_GSKvsIsotretinoin_PT3", "Isotretinoin_and_GSKvsIsotretinoin_PT7", "Isotretinoin_and_GSKvsIsotretinoin_PT10", "Isotretinoin_and_TAZvsTAZ_PT3", "Isotretinoin_and_TAZvsTAZ_PT7", "Isotretinoin_and_TAZvsTAZ_PT10", "Isotretinoin_and_GSKvsGSK_PT3", "Isotretinoin_and_GSKvsGSK_PT7", "Isotretinoin_and_GSKvsGSK_PT10")


for (i in filenames){
# Get pathways from gmt file
GO_pathways <- gmtPathways("C:/Users/Christina/Documents/GSEA .gmt/c5.go.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 )
gsea_df$pathway_readable = str_wrap(gsea_df$pathway_readable, width = 70)

 #if more than 20 significant terms only plot the 20 most significant
    if(nrow(gsea_df)>=20){
      p <- ggplot(head(gsea_df, 20), 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=15))
      try(print(p))
 try(ggsave(p, file=paste0("GO_barplot_", i,".png"), width = 35, height = 23, 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=15))
       try(print(p))
try(ggsave(p, file=paste0("GO_barplot_", i,".png"), width = 35, height = 23, units = "cm"))

    }

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


```

#BARPLOT GSEA
```{r}
#plot immune only related TAZ gene sets

#read in csv
#gsea_df <- read.csv("GO_TAZvsDMSO_PT10_immune_only.csv", header = TRUE)
gsea_df <- read_xlsx("GSEA_COMB1_PT10_DNA_REP.xlsx")

gsea_df$pathway_readable <- gsub("_", " ", gsea_df$pathway )


#gsea_df$pathway_readable = str_wrap(gsea_df$pathway_readable, width = 60)

#create barplot of enriched pathways:
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("Isotretinoin + TAZ post-treatment day 10 DNA replication-related gene sets") + 
theme(axis.text=element_text(size=14))

try(ggsave(p, file=paste0("GO_barplot_COMB1_DNA_rep_only_PT10.png"), width = 35, height = 20, units = "cm"))


```



# GSEA 59 NB DIFF GENE SET CUSTOM
```{r}
########################### 59 NB DIFF GENE SET CUSTOM ############################

library(fgsea)
#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.


for (i in filenames){
  
  library(readxl)
  
  # for reactome
  diff_geneset <- read_excel("~/Neuroblastoma RNAseq analysis/GSEA custom gene sets/diff_gene_set.xlsx")
  diff_geneset <- as.list(diff_geneset)
  diff_geneset.ind <- ids2indices( diff_geneset, rownames(edgeRData$counts))
  # Run analysis:
gst.camera <- camera.DGEList(edgeRData, design= design, index=diff_geneset.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)
  
   diff_geneset.ind <- as.vector(diff_geneset.ind)
   
   write.csv(sig_genesets, file=paste0("NB_diff_sig_", i,".csv") )
# p <-   barcodeplot(dge_Isotretinoin_and_TAZ_DMSO_PT10$table$logFC, index=diff_geneset.ind, main="NEUROBLASTOMA DIFFERENTIATION UP",
# labels=c("Negative logFC", "Positive logFC"),
#  quantile=c(-0.5,0.5))
# print(p)

}

```


```{r}

dge_files <- list(dge_Isotretinoin_and_TAZ_RA_PT10)


data <- dge_Isotretinoin_and_TAZ_DMSO_PT10$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"))

up_dge_Isotretinoin_and_TAZ_DMSO_PT10 <- data[data$Expression == 'Up-regulated',]


data <- dge_TAZ_DMSO_PT10$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"))
up_dge_TAZ_DMSO_PT10 <- data[data$Expression == 'Up-regulated',]

data <- dge_GSK_DMSO_PT10$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"))
up_dge_GSK_DMSO_PT10 <- data[data$Expression == 'Up-regulated',]

data <- dge_RA_DMSO_PT10$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"))
up_dge_RA_DMSO_PT10 <- data[data$Expression == 'Up-regulated',]


data <- dge_Isotretinoin_and_GSK_DMSO_PT10$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"))
up_dge_Isotretinoin_and_GSK_DMSO_PT10 <- data[data$Expression == 'Up-regulated',]

df_list <- list(up_dge_Isotretinoin_and_TAZ_DMSO_PT10, up_dge_Isotretinoin_and_GSK_DMSO_PT10)
                
                
#merge all data frames in list
overlapping_DEGS <- df_list %>% reduce(merge, by="geneID")


i <- dge_Isotretinoin_and_TAZ_RA_PT10

  i$table$geneID <- row.names(i$table)
  
df_new <- i[i$table$geneID %in% diff_geneset$NEURONAL_DIFF,]
  print(df_new$table)
  
  
  
  df_new$table$FDR <- p.adjust(data$PValue, method = "BH")
  x <- df_new$table[df_new$table$FDR < 0.05,]
   x <- x[x$logFC > 1,]


```

# GSEA NEURONAL DIFF GENE SET

```{r}
############################ NEURONAL DIFF GENE SET #################

#use design matrix from GLM
design <- model.matrix(~0_and_treatment_timepoint)

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

##### samples removed for which there are no signicant gene sets enriched  "GSK_DMSO_PT3",
filenames <- c("Isotretinoin_and_TAZ_DMSO_PT3",  "Isotretinoin_and_TAZ_DMSO_PT7" , "Isotretinoin_and_TAZ_DMSO_PT10",
               "Isotretinoin_and_GSK_DMSO_PT3" , "Isotretinoin_and_GSK_DMSO_PT7" , "Isotretinoin_and_GSK_DMSO_PT10"  ,   "GSK_DMSO_PT7" ,   "GSK_DMSO_PT10",   "TAZ_DMSO_PT3"  ,  "TAZ_DMSO_PT7"  ,  "TAZ_DMSO_PT10"  ,
               "RA_DMSO_PT3"   ,  "RA_DMSO_PT7"  ,   "RA_DMSO_PT10", "Isotretinoin_and_TAZ_RA_PT3", "Isotretinoin_and_TAZ_RA_PT7", "Isotretinoin_and_TAZ_RA_PT10", "Isotretinoin_and_GSK_RA_PT3", "Isotretinoin_and_GSK_RA_PT7", "Isotretinoin_and_GSK_RA_PT10")


for (i in filenames){
  
  library(readxl)
  
  # for reactome
  LE_NEURONAL_DIFFERENTIATION_UP <- read_excel("~/Neuroblastoma RNAseq analysis/LE_NEURONAL_DIFFERENTIATION_UP.xlsx")
  LE_NEURONAL_DIFFERENTIATION_UP <- as.list(LE_NEURONAL_DIFFERENTIATION_UP)
  LE_NEURONAL_DIFFERENTIATION_UP.ind <- ids2indices( LE_NEURONAL_DIFFERENTIATION_UP, rownames(edgeRData$counts))
  # Run analysis:
  gst.camera <- camera.DGEList(edgeRData, design= design, index=LE_NEURONAL_DIFFERENTIATION_UP.ind,contrast=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)
  

  
}

```




# NB differentation gene sig

```{r}

#test for 59 NB differentation sig -------------------------------------------------------

#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("Isotretinoin_and_TAZvsDMSO_PT3", "Isotretinoin_and_TAZvsDMSO_PT7", "Isotretinoin_and_TAZvsDMSO_PT10", "Isotretinoin_and_GSKvsDMSO_PT3", "Isotretinoin_and_GSKvsDMSO_PT7", "Isotretinoin_and_GSKvsDMSO_PT10", "GSKvsDMSO_PT3", "GSKvsDMSO_PT7", "GSKvsDMSO_PT10", "IsotretinoinvsDMSO_PT3", "IsotretinoinvsDMSO_PT7", "IsotretinoinvsDMSO_PT10", "TAZvsDMSO_PT3", "TAZvsDMSO_PT7", "TAZvsDMSO_PT10", "Isotretinoin_and_TAZvsIsotretinoin_PT3", "Isotretinoin_and_TAZvsIsotretinoin_PT7", "Isotretinoin_and_TAZvsIsotretinoin_PT10", "Isotretinoin_and_GSKvsIsotretinoin_PT3", "Isotretinoin_and_GSKvsIsotretinoin_PT7", "Isotretinoin_and_GSKvsIsotretinoin_PT10", "TAZvsIsotretinoin_PT3", "TAZvsIsotretinoin_PT7", "TAZvsIsotretinoin_PT10", "GSKvsIsotretinoin_PT3", "GSKvsIsotretinoin_PT7", "GSKvsIsotretinoin_PT10", "Isotretinoin_and_TAZvsTAZ_PT3", "Isotretinoin_and_TAZvsTAZ_PT7", "Isotretinoin_and_TAZvsTAZ_PT10", "Isotretinoin_and_GSKvsTAZ_PT3", "Isotretinoin_and_GSKvsTAZ_PT7", "Isotretinoin_and_GSKvsTAZ_PT10", "GSKvsTAZ_PT3", "GSKvsTAZ_PT7", "GSKvsTAZ_PT10", "Isotretinoin_and_TAZvsGSK_PT3", "Isotretinoin_and_TAZvsGSK_PT7", "Isotretinoin_and_TAZvsGSK_PT10", "Isotretinoin_and_TAZvsGSK_PT3", "Isotretinoin_and_TAZvsGSK_PT7", "Isotretinoin_and_TAZvsGSK_PT10")
for (i in filenames){
  
  diff_geneset <- read_excel("~/Neuroblastoma RNAseq analysis/GSEA custom gene sets/diff_gene_set.xlsx")
diff_geneset <- diff_geneset$NEURONAL_DIFF
   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, "NB_gene_sig_results.csv", append = TRUE,sep = ",",col.names = FALSE)

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


  
  
}



```


# random gene signature
Testing a random gene signature to see if it is significantly enriched.
```{r}

filenames <- c("Isotretinoin_and_TAZvsDMSO_PT3", "Isotretinoin_and_TAZvsDMSO_PT7", "Isotretinoin_and_TAZvsDMSO_PT10", "Isotretinoin_and_GSKvsDMSO_PT3", "Isotretinoin_and_GSKvsDMSO_PT7", "Isotretinoin_and_GSKvsDMSO_PT10", "GSKvsDMSO_PT3", "GSKvsDMSO_PT7", "GSKvsDMSO_PT10", "IsotretinoinvsDMSO_PT3", "IsotretinoinvsDMSO_PT7", "IsotretinoinvsDMSO_PT10", "TAZvsDMSO_PT3", "TAZvsDMSO_PT7", "TAZvsDMSO_PT10")

#get list of all genes after filtering 
all_genes <- rownames(counts.keep)
random_genes <- sample(all_genes, 59)

  #random_genes <- as.list(random_genes)
  random_genes.ind <- ids2indices(random_genes, rownames(edgeRData$counts))

for (i in filenames){

#for 59 random genes --------------------------------------------------------------------------
  
  ##select 59 random genes and test for enrichment of these as a signature

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

    gst.camera <- camera(edgeRData, design= design, index=random_genes.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
   try(rownames(sig_genesets) <- paste0(i))

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

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



#Bootstap attempt

```{r}

filenames <- c("Isotretinoin_and_TAZvsDMSO_PT3")

#get list of all genes after filtering 
all_genes <- rownames(counts.keep)
random_genes <- sample(all_genes, 59)

  #random_genes <- as.list(random_genes)
  random_genes.ind <- ids2indices(random_genes, rownames(edgeRData$counts))

for (i in filenames){

#for 59 random genes --------------------------------------------------------------------------
  
  ##select 59 random genes and test for enrichment of these as a signature

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

    gst.camera <- camera(edgeRData, design= design, index=random_genes.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
   try(rownames(sig_genesets) <- paste0(i))

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

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

# Zage NB differentation signature

```{r}

#test fo NB differentation sig by Zage et al. -------------------------------------------------------

#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("Isotretinoin_and_TAZvsDMSO_PT3", "Isotretinoin_and_TAZvsDMSO_PT7", "Isotretinoin_and_TAZvsDMSO_PT10", "Isotretinoin_and_GSKvsDMSO_PT3", "Isotretinoin_and_GSKvsDMSO_PT7", "Isotretinoin_and_GSKvsDMSO_PT10", "GSKvsDMSO_PT3", "GSKvsDMSO_PT7", "GSKvsDMSO_PT10", "IsotretinoinvsDMSO_PT3", "IsotretinoinvsDMSO_PT7", "IsotretinoinvsDMSO_PT10", "TAZvsDMSO_PT3", "TAZvsDMSO_PT7", "TAZvsDMSO_PT10", "Isotretinoin_and_TAZvsIsotretinoin_PT3", "Isotretinoin_and_TAZvsIsotretinoin_PT7", "Isotretinoin_and_TAZvsIsotretinoin_PT10", "Isotretinoin_and_GSKvsIsotretinoin_PT3", "Isotretinoin_and_GSKvsIsotretinoin_PT7", "Isotretinoin_and_GSKvsIsotretinoin_PT10", "TAZvsIsotretinoin_PT3", "TAZvsIsotretinoin_PT7", "TAZvsIsotretinoin_PT10", "GSKvsIsotretinoin_PT3", "GSKvsIsotretinoin_PT7", "GSKvsIsotretinoin_PT10", "Isotretinoin_and_TAZvsTAZ_PT3", "Isotretinoin_and_TAZvsTAZ_PT7", "Isotretinoin_and_TAZvsTAZ_PT10", "Isotretinoin_and_GSKvsTAZ_PT3", "Isotretinoin_and_GSKvsTAZ_PT7", "Isotretinoin_and_GSKvsTAZ_PT10", "GSKvsTAZ_PT3", "GSKvsTAZ_PT7", "GSKvsTAZ_PT10", "Isotretinoin_and_TAZvsGSK_PT3", "Isotretinoin_and_TAZvsGSK_PT7", "Isotretinoin_and_TAZvsGSK_PT10", "Isotretinoin_and_TAZvsGSK_PT3", "Isotretinoin_and_TAZvsGSK_PT7", "Isotretinoin_and_TAZvsGSK_PT10")
for (i in filenames){
  
  diff_geneset <- read_excel("~/Neuroblastoma RNAseq analysis/Zage NB differentiation signature/Zage_NB_diff_sig.xlsx") %>% na.omit() %>% unique()
  

diff_geneset <- diff_geneset$Genes
   diff_geneset.ind <- ids2indices( diff_geneset, rownames(edgeRData$counts))
  # Run analysis:
  gst.camera <- camera.DGEList(edgeRData, design= design, index=diff_geneset.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, "NB_gene_sig_results.csv", append = TRUE,sep = ",",col.names = FALSE))

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


  
  
}

  
```




# Unique GO terms
```{r}


####################### for GO terms#######################

#use design matrix from GLM
design <- model.matrix(~0_and_treatment_timepoint)

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

##### samples removed for which there are no signicant gene sets enriched 
filenames <- c("Isotretinoin_and_TAZ_DMSO_PT3"
              ,  "Isotretinoin_and_TAZ_DMSO_PT7" , "Isotretinoin_and_TAZ_DMSO_PT10",
               "Isotretinoin_and_GSK_DMSO_PT3" , "Isotretinoin_and_GSK_DMSO_PT7" , "Isotretinoin_and_GSK_DMSO_PT10"  ,"GSK_DMSO_PT3", "GSK_DMSO_PT7" ,   "GSK_DMSO_PT10",   "TAZ_DMSO_PT3"  ,  "TAZ_DMSO_PT7"  ,  "TAZ_DMSO_PT10"  ,
               "RA_DMSO_PT3"   ,  "RA_DMSO_PT7"  ,   "RA_DMSO_PT10")

#filenames <-  c("Isotretinoin_and_TAZ_TAZ_PT3", "Isotretinoin_and_TAZ_TAZ_PT7", "Isotretinoin_and_TAZ_TAZ_PT10")

filenames <- c("Isotretinoin_and_TAZ_DMSO_PT10")

for (i in filenames){
  
  
  
# 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[, 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(i)
print(sig_genesets)

#sig_genesets <- sig_genesets[rownames(sig_genesets) %like% "DIFFERENTIATION", ]
 
TAZ_DMSO_PT10_gsea_df <- sig_genesets
gsea_df$padj <- as.numeric(gsea_df$FDR)
gsea_df$pathway <- as.character(row.names(gsea_df))

  #barplot
 p <- ggplot(gsea_df, aes(x = fct_reorder(pathway, - FDR), y = NGenes))_and_
  geom_bar(aes(fill = gsea_df$padj), stat = "identity", width = 0.8) _and_
  scale_fill_continuous(low = "red", high = "blue", na.value = NA)_and_
  scale_x_discrete(expand = expansion(add = .5)) _and_
  coord_flip()_and_ try(facet_grid(~Direction))  _and_ theme_bw() _and_ labs(y="Count", x="")  _and_ggtitle(i)
  try(print(p))
 # try(ggsave(p, file=paste0("GO_terms_barplot_", i,".png"), width = 50, height = 40, units = "cm"))

}




  Isotretinoin_and_TAZ_DMSO_PT3_gsea_df$pathway <- row.names(Isotretinoin_and_TAZ_DMSO_PT3_gsea_df)
    TAZ_DMSO_PT3_gsea_df$pathway <- row.names(TAZ_DMSO_PT3_gsea_df)
  

 x <- unique(TAZ_DMSO_PT3_gsea_df$pathway, Isotretinoin_and_TAZ_DMSO_PT3_gsea_df$pathway)
# write.csv(x, "~/Neuroblastoma RNAseq analysis/GSEA plots saved/TAZ_not_Isotretinoin_and_TAZ_PT3.csv")
 
 
 ############################# PT7 #############################
 Isotretinoin_and_TAZ_DMSO_PT7_gsea_df$pathway <- row.names(Isotretinoin_and_TAZ_DMSO_PT7_gsea_df)
    TAZ_DMSO_PT7_gsea_df$pathway <- row.names(TAZ_DMSO_PT7_gsea_df)
  

 x <- unique(TAZ_DMSO_PT7_gsea_df$pathway, Isotretinoin_and_TAZ_DMSO_PT7_gsea_df$pathway)
 #write.csv(x, "~/Neuroblastoma RNAseq analysis/GSEA plots saved/TAZ_not_Isotretinoin_and_TAZ_PT7.csv")
 
 
 
 ############################# PT10 #############################
 Isotretinoin_and_TAZ_DMSO_PT10_gsea_df$pathway <- row.names(Isotretinoin_and_TAZ_DMSO_PT10_gsea_df)
    TAZ_DMSO_PT10_gsea_df$pathway <- row.names(TAZ_DMSO_PT10_gsea_df)
  

 x <- unique(TAZ_DMSO_PT10_gsea_df$pathway, Isotretinoin_and_TAZ_DMSO_PT10_gsea_df$pathway)
 write.csv(x, "~/Neuroblastoma RNAseq analysis/GSEA plots saved/TAZ_not_Isotretinoin_and_TAZ_PT10.csv")
 
 
```


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

#use design matrix from GLM
design <- model.matrix(~0_and_treatment_timepoint)

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

##### samples removed for which there are no signicant gene sets enriched 
filenames <- c("Isotretinoin_and_TAZ_DMSO_PT3"
              ,  "Isotretinoin_and_TAZ_DMSO_PT7" , "Isotretinoin_and_TAZ_DMSO_PT10",
               "Isotretinoin_and_GSK_DMSO_PT3" , "Isotretinoin_and_GSK_DMSO_PT7" , "Isotretinoin_and_GSK_DMSO_PT10"  ,"GSK_DMSO_PT3", "GSK_DMSO_PT7" ,   "GSK_DMSO_PT10",   "TAZ_DMSO_PT3"  ,  "TAZ_DMSO_PT7"  ,  "TAZ_DMSO_PT10"  ,
               "RA_DMSO_PT3"   ,  "RA_DMSO_PT7"  ,   "RA_DMSO_PT10")

#filenames <-  c("Isotretinoin_and_TAZ_TAZ_PT3", "Isotretinoin_and_TAZ_TAZ_PT7", "Isotretinoin_and_TAZ_TAZ_PT10")

filenames <- c("GSK_DMSO_PT10")

for (i in filenames){
  
  
  
# 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[, 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(i)
print(sig_genesets)

#sig_genesets <- sig_genesets[rownames(sig_genesets) %like% "DIFFERENTIATION", ]
 
GSK_DMSO_PT10_gsea_df <- sig_genesets
gsea_df$padj <- as.numeric(gsea_df$FDR)
gsea_df$pathway <- as.character(row.names(gsea_df))

  #barplot
 p <- ggplot(gsea_df, aes(x = fct_reorder(pathway, - FDR), y = NGenes))_and_
  geom_bar(aes(fill = gsea_df$padj), stat = "identity", width = 0.8) _and_
  scale_fill_continuous(low = "red", high = "blue", na.value = NA)_and_
  scale_x_discrete(expand = expansion(add = .5)) _and_
  coord_flip()_and_ try(facet_grid(~Direction))  _and_ theme_bw() _and_ labs(y="Count", x="")  _and_ggtitle(i)
  try(print(p))
 # try(ggsave(p, file=paste0("GO_terms_barplot_", i,".png"), width = 50, height = 40, units = "cm"))

}




  Isotretinoin_and_GSK_DMSO_PT3_gsea_df$pathway <- row.names(Isotretinoin_and_GSK_DMSO_PT3_gsea_df)
    GSK_DMSO_PT3_gsea_df$pathway <- row.names(GSK_DMSO_PT3_gsea_df)
  

 x <- unique(GSK_DMSO_PT3_gsea_df$pathway, Isotretinoin_and_GSK_DMSO_PT3_gsea_df$pathway)
# write.csv(x, "~/Neuroblastoma RNAseq analysis/GSEA plots saved/GSK_not_Isotretinoin_and_GSK_PT3.csv")
 
 
 ############################# PT7 #############################
 Isotretinoin_and_GSK_DMSO_PT7_gsea_df$pathway <- row.names(Isotretinoin_and_GSK_DMSO_PT7_gsea_df)
    GSK_DMSO_PT7_gsea_df$pathway <- row.names(GSK_DMSO_PT7_gsea_df)
  

 x <- unique(GSK_DMSO_PT7_gsea_df$pathway, Isotretinoin_and_GSK_DMSO_PT7_gsea_df$pathway)
 write.csv(x, "~/Neuroblastoma RNAseq analysis/GSEA plots saved/GSK_not_Isotretinoin_and_GSK_PT7.csv")
 
 
 
 ############################# PT20 #############################
 Isotretinoin_and_GSK_DMSO_PT10_gsea_df$pathway <- row.names(Isotretinoin_and_GSK_DMSO_PT10_gsea_df)
    GSK_DMSO_PT10_gsea_df$pathway <- row.names(GSK_DMSO_PT10_gsea_df)
  

 x <- unique(GSK_DMSO_PT10_gsea_df$pathway, Isotretinoin_and_GSK_DMSO_PT10_gsea_df$pathway)
 write.csv(x, "~/Neuroblastoma RNAseq analysis/GSEA plots saved/GSK_not_Isotretinoin_and_GSK_PT10.csv")
```



```{r}
########### FOR GO DO TOP 15  ###################
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))) _and_ theme_bw(base_size = 5) _and_
  geom_point(aes(size = top_10_GO$NGenes, color = top_10_GO$FDR)) _and_
 theme(axis.text.y = element_text(size =6)) _and_
  
  scale_colour_gradient( low="red", high = "blue", name = "p.adjust")_and_
  ylab(NULL) _and_ facet_grid(.~Direction) _and_  labs(size="Gene set size") _and_  xlab("Count") _and_ggtitle(i)

# print(p)
 
gsea_df <- top_10_GO
gsea_df$padj <- as.numeric(gsea_df$FDR)
gsea_df$pathway <- as.character(row.names(gsea_df)) 

 #barplot
 p <- ggplot(gsea_df, aes(x = fct_reorder(pathway, - FDR), y = NGenes))_and_
  geom_bar(aes(fill = gsea_df$padj), stat = "identity", width = 0.8) _and_
  scale_fill_continuous(low = "red", high = "blue", na.value = NA, name = "p.adjust")_and_
  scale_x_discrete(expand = expansion(add = .5)) _and_
  coord_flip()_and_ facet_grid(~Direction)  _and_ theme_bw() _and_ labs(y="Count", x="")  _and_ggtitle("Isotretinoin _and_ TAZ PT10")
  print(p)

}

```






#Using fgsea for GSEA





```{r}

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


dge_files <- list(dge_Isotretinoin_and_TAZ_DMSO_PT3, dge_Isotretinoin_and_TAZ_DMSO_PT7, dge_Isotretinoin_and_TAZ_DMSO_PT10, dge_Isotretinoin_and_GSK_DMSO_PT3, dge_Isotretinoin_and_GSK_DMSO_PT7, dge_Isotretinoin_and_GSK_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_Isotretinoin_and_TAZ_DMSO_PT3", "dge_Isotretinoin_and_TAZ_DMSO_PT7", "dge_Isotretinoin_and_TAZ_DMSO_PT10", "dge_Isotretinoin_and_GSK_DMSO_PT3", "dge_Isotretinoin_and_GSK_DMSO_PT7", "dge_Isotretinoin_and_GSK_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))) _and_ 
  geom_point(aes(size = gsea_df$size, color = gsea_df$padj)) _and_
  theme_bw(base_size = 5) _and_  theme(axis.text.y = element_text(size =6))_and_
  scale_colour_gradient( low="red", high = "blue", name = "p.adjust")_and_
  ylab(NULL) _and_ facet_grid(.~type) _and_ labs(size="Gene set size") _and_  xlab("Count") _and_ ggtitle(a)#_and_ ggtitle(filenames[[i]])

 print(p)
 
 }


```


# Plot GSEA barplot for just immune-related terms
```{r}


Reactome_pathways <- gmtPathways("TAZ_PT10_immune_only.xlsx")
sig_genesets <- read_xlsx("TAZ_PT10_immune_only.xlsx")
sig_genesets <- as.data.frame(sig_genesets)
sig_genesets$pathway <- sig_genesets$...1

gsea_df <- sig_genesets
gsea_df$padj <- as.numeric(gsea_df$FDR)
#gsea_df$pathway <- as.character(row.names(gsea_df))

p <- ggplot(gsea_df, aes(x = fct_reorder(pathway, - FDR), y = NGenes))_and_
  geom_bar(aes(fill = FDR), stat = "identity", width = 0.8) _and_
  scale_fill_continuous(low = "red", high = "blue", na.value = NA, name = "p.adjust")_and_
  scale_x_discrete(expand = expansion(add = .5)) _and_
  coord_flip()_and_ try(facet_grid(~Direction))  _and_ theme_bw() _and_ labs(y="Count", x="")  _and_ggtitle("TAZ PT10 immune-related terms") 

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



Reactome_pathways <- gmtPathways("GSK_PT10_immune_only.xlsx")
sig_genesets <- read_xlsx("GSK_PT10_immune_only.xlsx")
sig_genesets <- as.data.frame(sig_genesets)
sig_genesets$pathway <- sig_genesets$...1

gsea_df <- sig_genesets
gsea_df$padj <- as.numeric(gsea_df$FDR)
#gsea_df$pathway <- as.character(row.names(gsea_df))

p <- ggplot(gsea_df, aes(x = fct_reorder(pathway, - FDR), y = NGenes))_and_
  geom_bar(aes(fill = FDR), stat = "identity", width = 0.8) _and_
  scale_fill_continuous(low = "red", high = "blue", na.value = NA, name = "p.adjust")_and_
  scale_x_discrete(expand = expansion(add = .5)) _and_
  coord_flip()_and_ try(facet_grid(~Direction))  _and_ theme_bw() _and_ labs(y="Count", x="")  _and_ggtitle("GSK PT10 immune-related terms") 

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


```







#maSigPro

```{r}
library(maSigPro) # load maSigPro library

#get design array
edesign <- read_xlsx("C:/Users/Christina/Documents/Neuroblastoma RNAseq analysis/edesign.xlsx")


 edesign <- column_to_rownames(edesign, var="...1")

#get matrix with normalized gene expression data. Genes must be in rows and arrays in columns. 
 
#get normalised filtered count matrix
rownames(edgeRData) <- paste0("X", seq_len(nrow(edgeRData)))
#getting for all counts only (rawdata)
edgeRData <- DGEList(counts=rawData, genes=rownames(rawData), group = treatment)
edgeRData <- calcNormFactors(edgeRData, method="TMM")
CPM <- cpm(edgeRData, log = FALSE)
data.NB <- CPM
 
#check if colnames in CPM are the same as rownames in edesign
identical(colnames(data.NB), rownames(edesign))



#Create a regression matrix for the full regression model:
design <- make.design.matrix(edesign, degree = 2)

#compute a regression fit for each gene. This is done by the function p.vector().This function also computes the p-value associated to the F-Statistic of the model, which is used to select significant genes. By default maSigPro corrects this p-value for multiple comparisons by applying the linear step-up (B-H) false discovery rate (FDR) procedure (Benjamini and Hochberg, 1995)

#set counts =T for RNA seq data adn theta to 10. Use's a negative binomial model 
fit <- p.vector(data.NB, design, Q = 0.05, MT.adjust = "BH", min.obs = 20, counts= TRUE,  theta=10)

fit$i # returns the number of significant genes
fit$alfa # gives p-value at the Q false discovery control level
fit$SELEC # is a matrix with the significant genes and their expression values

#Once significant genes have been found, maSigPro applies a variable selection procedure to find significant variables for each gene. This will ultimatelly be used to find which are the profile differences between experimental groups. This step is done by the T.fit() function.

tstep <- T.fit(data = fit, design = design$dis, step.method = "forward", alfa = 0.05)

#get sig genes
 sigs <- get.siggenes(tstep, rsq = 0.6, vars = "groups")
 names(sigs$sig.genes)

 #venn diagrams
 suma2Venn(sigs$summary[, c(2:4)])
  suma2Venn(sigs$summary[, c(2:6)])
 
 #see.genes() performs a cluster analysis to group genes by similar profiles
 see.genes(sigs$sig.genes$Isotretinoin_and_TAZvsControl, show.fit = F, dis =design$dis,
cluster.method="hclust" ,cluster.data = 1, k.mclust=TRUE)

 PlotProfiles()
 
 
#as experiment-wide expression profiles
and as by-groups profiles. The first plot (Figure 3) will help to evaluate the consistency of
the clusters while the second plot shows clearly the differences between groups (Figure 4).

```

#Patterns

```{r}
#temporal analsyis using the patterns package
#load package
library("Patterns")

#get normalised filtered count matrix
rownames(edgeRData) <- paste0("X", seq_len(nrow(edgeRData)))
#getting for all counts only (rawdata)
edgeRData <- DGEList(counts=rawData, genes=rownames(rawData), group = treatment)
edgeRData <- calcNormFactors(edgeRData, method="TMM")
CPM <- cpm(edgeRData, log = FALSE)

#turn data into omics array object
micro_US<-as.omics_array(CPM, time=c(3,3, 3, 7, 7, 7,10, 10,10),subject=6)
summary(micro_US)
#plot(micro_US)

#summary(M)

#select genes

#S4 method for signature 'omics_array,omics_array,numeric'
selection1<-geneSelection(
micro_US,
micro_US,
tot.number =1,
data_log = TRUE,
wanted.patterns = NULL,
forbidden.patterns = NULL,
peak = NULL,
alpha = 0.05,
Design = NULL,
lfc = 0
)

geneSelection(
micro_US,
micro_US,
tot.number = 1,
data_log = TRUE,
wanted.patterns = NULL,
forbidden.patterns = NULL,
peak = NULL,
alpha = 0.05,
Design = NULL,
lfc = 0
)


We then merge these selections:
Selection<-unionOmics(list(Selection1,Selection2,Selection3,Selection4,Selection5))
print(Selection)


```

#K-means clustering
```{r}

#Genes significantly expressed in at least one of these comparisons were retained for clustering analysis.

#get list of genes DE in at least one condition
files <- list.files(path = "~/NB/SIGNIFICANT_DEGs", pattern = "vsDMSO" )


#create empty list
my_list <- list()

for (i in files){ 
  #read file 
  DEGs <- read.csv(i)
  
  #get gene names
  genes <- DEGs$X

  #append into a list
  my_list[[length(my_list) + 1]] <- genes   # Append new list element
  }

#merge list
merged.data.frame = Reduce(function(...) merge(..., all=T), my_list)
#remove duplicated gene names
all_DEGs <- merged.data.frame[!duplicated(merged.data.frame$x), ]
all_DEGs <- all_DEGs$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 = TRUE)

 #subset raw data for the list
 counts_DEGs <- CPM[row.names(CPM) %in% all_DEGs,]
 
#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(counts_DEGs)
scaledata <- t(scale(t(DEgenes)))

#hcluster ------------------------------------------------------------------------

#cluster samples
hc <- hclust(as.dist(1-cor(scaledata, method="spearman")), method="complete") # Clusters columns by Spearman correlation.

sampleTree = as.dendrogram(hc, method="average")
plot(sampleTree,
     main = "Sample Clustering",
     ylab = "Height")

#cluster genes
hr <- hclust(as.dist(1-cor(t(scaledata), method="pearson")), method="complete") # Cluster rows by Pearson correlation.

#heatmap ----------------------------------------------------------------------
heatmap.2(DEgenes,
          Rowv=as.dendrogram(hr), 
          Colv=as.dendrogram(hc),
          col=redgreen(100),
          scale="row",
          margins = c(7, 7),
          cexCol = 0.7,
          labRow = F,
          main = "Heatmap.2",
          trace = "none")


#kmeans -----------------------------------------------------------------
# determine optimal number of clusters
#NOTE: with the selection of significantly differentially expressed either up or down regulated genes here, two clusters will dominate first


png('kmeans_clusters_elbowmethod_plot.png',width=400,height=250)
fviz_nbclust(scaledata, kmeans, method = "wss", k.max = 10, ) + geom_vline(xintercept = 3, linetype = 2)+ labs(subtitle = "Elbow method")

dev.off()

# visual check of the value of k choice
k<-3
res.k<-hkmeans(scaledata,k,hc.method='complete')
png(paste0('visualisation_of_clusters',k,'_clusters.png'),width=350,height=300)
print(fviz_cluster(res.k, ellipse=FALSE,labelsize=6,ggtheme = theme_classic(),geom=c('point'),pointsize=1,stand=TRUE,shape='circle')+scale_colour_manual(values = c("darkorange", "green4","mediumvioletred",'royalblue3','darkgoldenrod2','red3','grey78','grey28')))
dev.off()
         
# Silhouette method
png('Silhouette_method.png',width=400,height=250)
fviz_nbclust(scaledata, kmeans, method = "silhouette")+
  labs(subtitle = "Silhouette method")
dev.off()


kdf <- scaledata
kdf <- as.data.frame(kdf)

# conduct clustering
res_k<-hkmeans(kdf,k,hc.method='complete')
kdf$k_clusters<-res_k$cluster
kdf$k_clusters<-as.character(kdf$k_clusters)

kdf<-kdf[order(kdf$k_clusters),]

# extract number of genes in each cluster
for(i in 1:k){
  print(nrow(subset(kdf,k_clusters==i)))
}       


# format and write output table
kdf_table<-kdf
kdf_table$gene<-rownames(kdf)
kdf_table<-kdf_table[,c('gene','k_clusters')]
names(kdf_table)<-c('gene','cluster')
write.table(kdf_table,file=paste0('cluster_details_k',k,'_alllines.txt'),row.names=FALSE,quote=FALSE,col.names=TRUE)


# plot clusters as heatmap ------------------------------------------------

 # colnames(kdf)<- sub("Isotretinoin_and_TAZ", "RA_and_TAZ", colnames(kdf))
 # colnames(kdf)<- sub("Isotretinoin_and_GSK", "RA_and_GSK", colnames(kdf))

#reorder columns to time order
 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", "k_clusters")
kdf <- kdf[, col_order]

#check reorder
colnames(kdf)

png(filename="heatmap_3_test.png", width=3600, bg="white")
h1<-Heatmap(as.matrix(kdf[,1:(ncol(kdf)-1)]),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"))
h2<-Heatmap(as.matrix(kdf[,'k_clusters']),cluster_columns=FALSE,cluster_rows=FALSE, heatmap_legend_param = list(title = "Cluster"))


hlist <- h1+h2
tiff(paste0('Heatmap_k_',k,'_clusters.tiff'),width=3600,height=3000,res=300)
draw(hlist, row_split = kdf$hk_clusters, cluster_row_slices = TRUE,cluster_rows=FALSE, row_dend_reorder = TRUE,show_row_dend = FALSE)
dev.off()





#with just cluster 2
#subset data for just cluster 2
test <- kdf[kdf$k_clusters == '2',]

h1<-Heatmap(as.matrix(test[,1:(ncol(test)-1)]),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('RA',9),rep('TAZ',9), rep("GSK", 9), rep("RA_and_TAZ", 9), rep("RA_and_GSK",9)),levels=c('DMSO','RA','TAZ', 'GSK','RA_and_TAZ','RA_and_GSK')),show_row_names=FALSE,show_column_names = FALSE, heatmap_legend_param = list(title = "Z-score"))
h2<-Heatmap(as.matrix(test[,'k_clusters']),cluster_columns=FALSE,cluster_rows=FALSE, heatmap_legend_param = list(title = "Cluster"))


hlist <- h1+h2
tiff(paste0('cluster_2_only_Hmap_kmeans_k',k,'_clusters.tiff'),width=2000,height=3000,res=300)
draw(hlist, row_split = kdf$hk_clusters, cluster_row_slices = TRUE,cluster_rows=FALSE, row_dend_reorder = TRUE,show_row_dend = FALSE)
dev.off()

#enrichment analysis on clusters---------------------------------------------------

library(clusterProfiler)

bg_universe<-DEgenes
length(bg_universe)

GO_pathways <- gmtPathways("C:/Users/Christina/Documents/GSEA .gmt/c5.go.v2024.1.Hs.symbols.gmt")


clusters<-read.table("~/NB/SIGNIFICANT_DEGs/cluster_details_k3_alllines.txt",stringsAsFactors = FALSE,header=TRUE)

GOtypes<-c('ALL')
for(GOtype in GOtypes){
  for(i in 1:k){
    geneList<-subset(clusters,cluster==i)$gene
    ego <- enrichGO(geneList, OrgDb = "org.Hs.eg.db", ont=GOtype, pAdjustMethod='fdr', pvalueCutoff=0.05,universe=bg_universe,keyType='SYMBOL')
    ego_df<-data.frame(ego,stringsAsFactors = FALSE)
    ego_df<-ego_df[order(ego_df$p.adjust),]
    write.csv(ego_df[,c(2,3,4,6,8)],file = paste0('clusters_',i,'_',GOtype,'GO.csv'),quote=FALSE,row.names=FALSE)
    ego_df$neglog10padj<--log(ego_df$p.adjust,10)
    ego_df<-ego_df[order(ego_df$p.adjust,decreasing=TRUE),] #order table by p-value
    print(nrow(ego_df))
    
    ego_df$Description<-factor(ego_df$Description,level=ego_df$Description)
    
    ego_df$gene_ratio<-sapply(ego_df$GeneRatio, function(x) as.numeric(strsplit(x,split='/')[[1]][1])/as.numeric(strsplit(x,split='/')[[1]][2]))
    
    #if more than 10 significant terms only plot the 20 most significant
    if(nrow(ego_df)>=10){
      png(paste0('clusters',k,'_',i,'_',GOtype,'.png'),width=650,height=475) 
      print(ggplot(ego_df[((nrow(ego_df)-19):nrow(ego_df)),],aes(x=neglog10padj,y=Description,fill=gene_ratio))+geom_bar(stat='identity')+theme_classic()+
              scale_x_continuous(expand=c(0,0))+labs(fill = 'Gene Ratio')+xlab('-log10 adjusted p-val')+ylab('')+theme(axis.text=element_text(size=12),axis.title=element_text(size=12)))
      dev.off()
    }else{
      png(paste0('clusters',k,'_',i,'_',GOtype,'.png'),width=650,height=475) 
      print(ggplot(ego_df,aes(x=neglog10padj,y=Description,fill=gene_ratio))+geom_bar(stat='identity')+theme_classic()+
              scale_x_continuous(expand=c(0,0))+labs(fill = 'Gene Ratio')+xlab('-log10 adjusted p-val')+ylab('')+theme(axis.text=element_text(size=12),axis.title=element_text(size=12)))
      dev.off()
    }
    rm(ego,ego_df,geneList)
  }}




{
  for(i in 1:k){
    i <- "3"
    
    geneList<-subset(clusters,cluster==i)$gene
#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)
#write.csv(ego2, "GO_hub_genes_orange.csv")

GO <- try(barplot(ego2,  font.size = 16 ,showCategory=10, x="Count"))+ scale_y_discrete(labels=function(x) str_wrap(x, width=80))

print(GO)
#try(ggsave(GO, file=paste0("GO_hub_genes_", i,".png"), width = 25, height = 30, 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)


em@result$Description<- gsub("_", " ", em@result$Description )

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




    
    
    ego <- enrichGO(geneList, OrgDb = "org.Hs.eg.db", ont=GOtype, pAdjustMethod='fdr', pvalueCutoff=0.05,universe=bg_universe,keyType='SYMBOL')
    ego_df<-data.frame(ego,stringsAsFactors = FALSE)
    ego_df<-ego_df[order(ego_df$p.adjust),]
    write.csv(ego_df[,c(2,3,4,6,8)],file = paste0('clusters_',i,'_',GOtype,'GO.csv'),quote=FALSE,row.names=FALSE)
    ego_df$neglog10padj<--log(ego_df$p.adjust,10)
    ego_df<-ego_df[order(ego_df$p.adjust,decreasing=TRUE),] #order table by p-value
    print(nrow(ego_df))
    
    ego_df$Description<-factor(ego_df$Description,level=ego_df$Description)
    
    ego_df$gene_ratio<-sapply(ego_df$GeneRatio, function(x) as.numeric(strsplit(x,split='/')[[1]][1])/as.numeric(strsplit(x,split='/')[[1]][2]))
    
    #if more than 10 significant terms only plot the 20 most significant
    if(nrow(ego_df)>=10){
      png(paste0('clusters',k,'_',i,'_',GOtype,'.png'),width=650,height=475) 
      print(ggplot(ego_df[((nrow(ego_df)-19):nrow(ego_df)),],aes(x=neglog10padj,y=Description,fill=gene_ratio))+geom_bar(stat='identity')+theme_classic()+
              scale_x_continuous(expand=c(0,0))+labs(fill = 'Gene Ratio')+xlab('-log10 adjusted p-val')+ylab('')+theme(axis.text=element_text(size=12),axis.title=element_text(size=12)))
      dev.off()
    }else{
      png(paste0('clusters',k,'_',i,'_',GOtype,'.png'),width=650,height=475) 
      print(ggplot(ego_df,aes(x=neglog10padj,y=Description,fill=gene_ratio))+geom_bar(stat='identity')+theme_classic()+
              scale_x_continuous(expand=c(0,0))+labs(fill = 'Gene Ratio')+xlab('-log10 adjusted p-val')+ylab('')+theme(axis.text=element_text(size=12),axis.title=element_text(size=12)))
      dev.off()
    }
    rm(ego,ego_df,geneList)
  }}







# Reactome ORA

for (i in 1:k) {  
  
  #get cluster name
  
#select genes just in cluster
geneList<-subset(clusters,cluster==i)$gene 

#convert to symbols
symbols <- mapIds(org.Hs.eg.db, geneList, 'ENTREZID', 'SYMBOL')
symbols <- as.data.frame(symbols)
geneList <- symbols$symbols

#perform reactome ORA on cluster
reactome_ORA <- enrichPathway(gene=geneList, pvalueCutoff = 0.05, readable=TRUE)
head(reactome_ORA)
write.csv(reactome_ORA, file = paste0('clusters_',i,'_Reactome.csv'),quote=FALSE,row.names=FALSE)

#get dotplot
   R <-   print(ggplot(reactome_ORA,aes(x=GeneRatio,y=Description,fill=p.adjust))+geom_bar(stat='identity')+theme_classic()+
              scale_x_continuous(expand=c(0,0))+labs(fill = 'p.adjust')+xlab('Gene Ratio')+ylab('')+theme(axis.text=element_text(size=12),axis.title=element_text(size=12)))

      
try(ggsave(R, file=paste0("Reactome_clusters_", i,".png"), width = 30, height = 20, units = "cm"))
      
# Hallmark ORA 
m_t2g <- msigdbr(species = "Homo sapiens", category = "H") %>% dplyr::select(gs_name, gene_symbol)
head(m_t2g)

#input gene list as geneIDs as some are lost when converting to entrez
geneList<-subset(clusters,cluster==i)$gene 
em <- enricher(geneList, pvalueCutoff = 0.05,  pAdjustMethod = "BH", TERM2GENE=m_t2g)
head(em)
write.csv(em, file = paste0('clusters_',i,'_Hallmark.csv'),quote=FALSE,row.names=FALSE)

 try(H <- print(ggplot(em,aes(x=GeneRatio,y=Description,fill=p.adjust))+geom_bar(stat='identity')+theme_classic()+
              scale_x_continuous(expand=c(0,0))+labs(fill = 'p.adjust')+xlab('Gene Ratio')+ylab('')+theme(axis.text=element_text(size=12),axis.title=element_text(size=12))))
     

try(ggsave(H, file=paste0("Hallmarks_clusters_", i,".png"), width = 30, height = 20, units = "cm"))
 dev.off()
}


```




#Heatmap of inflammatory genes

```{r}

genes_of_interest = c("HAS2","OSMR","NLRP3","AXL","LIF","VIP","ICAM4","GPC3","RASGRP1","NOD2","LAMP3","PTGIR","IL2RB","EBI3","NMUR1","BST2","C5AR1","SGMS2")


genes_of_interest <- c("CD4", "CD160", "IL1R1", "MR1", "IL17RE", "IL5RA", "OSMR", "IL11RA", "LEPR", "CSF3R", "IL17RC", "ACKR4", "CD74", "LIFR", "IL1RAP", "CCR6", "GFRA1", "IL17RA", "IL1RL2", "GFRA2", "IL20RA", "FZD4", "EPOR", "IL6R", "PRLR", "IL13RA2", "IL21R", "CXCR4", "IL12RB2", "IL6ST", "CNTFR", "CD44", "IL17RD", "ACKR3", "GHR", "CR2", "IL20RB", "CMKLR1", "CSF2RB", "IFNGR1", "GFRA3", "IL4R", "IL17REL", "CTSH", "HLA-DRB1", "HLA-DQB1", "IL10RA", "IFNAR1", "IFNGR2", "IL2RB", "EBI3", "IL1RAPL2", "IFNAR2", "IL17RB", "ACKR2", 
"IL13RA1", "HLA-DRB3", "IL27RA", "F3", "IFNLR1", "IL10RB", "GPR75", "KLRC2", "CRLF1", "CCR10", "C5AR1")





#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 = "Immune response genes",
          trace = "none")
```


#additive genes


```{r}

#find gene that have a greater LFC in combo than in single agent. Just doing this for TAZ as I'm intertested in the DNA damage response and also since TAZ appears to be having a greater effect, I think it will be more useful to check TAZ and PT10 than GSK. 

#merge dataframe 

df <- merge(Isotretinoin_and_TAZvsDMSO_PT10$table, IsotretinoinvsDMSO_PT10$table, by = 'row.names', all = TRUE)

#set rownames as rownames
df <- df %>% remove_rownames %>% column_to_rownames(var="Row.names")

#merge with TAZ as well
df <- merge(df, TAZvsDMSO_PT10$table, by = 'row.names', all = TRUE)

#rename cols
names(df)[2] <- "logFC.isotretinoin_and_TAZ"
names(df)[6] <- "logFC.isotretinoin"
names(df)[10] <- "logFC.TAZ"

#first get genes that have LFC greater than 0 logFC.x=iso+TAZ. logFC.y=iso logFC=TAZ
x <- filter(df, logFC.isotretinoin_and_TAZ > 0 & logFC.isotretinoin > 0 & logFC.TAZ > 0)
#then get genes that have a greater lFC in combo than single agent
additive_genes_up <- filter(x, logFC.isotretinoin_and_TAZ > logFC.isotretinoin & logFC.isotretinoin_and_TAZ > logFC.TAZ)
write.csv(additive_genes_up, file = paste0('additive_gene_list_up.csv'),quote=FALSE,row.names=FALSE)

#for each gene, find difference in LFC between iso+TAZ and TAZ
  
# Extracting column 1
LFC_iso_and_TAZ=additive_genes_up$logFC.isotretinoin_and_TAZ
 
# Extracting column 2
LFC_TAZ=additive_genes_up$logFC.TAZ

#extratcing column 3
LFC_iso=additive_genes_up$logFC.isotretinoin

# printing difference among
# two columns
difference_LFC_TAZ <- print(LFC_iso_and_TAZ-LFC_TAZ)

difference_LFC_iso <- print(LFC_iso_and_TAZ-LFC_iso)

difference_LFC_total <- c(difference_LFC_TAZ, difference_LFC_iso)
mean(difference_LFC_total)
range(difference_LFC_total)

#then get genes that have LFC less than 0
x <- filter(df, logFC.isotretinoin_and_TAZ < 0 & logFC.isotretinoin < 0 & logFC.TAZ < 0)
#then get genes that have a greater lFC in combo than single agent
additive_genes_down <- filter(x, logFC.isotretinoin_and_TAZ < logFC.isotretinoin & logFC.isotretinoin_and_TAZ < logFC.TAZ)
write.csv(additive_genes_down, file = paste0('additive_gene_list_down.csv'),quote=FALSE,row.names=FALSE)


#for each gene, find difference in LFC between iso+TAZ and TAZ
  
# Extracting column 1
LFC_iso_and_TAZ=additive_genes_down$logFC.isotretinoin_and_TAZ
 
# Extracting column 2
LFC_TAZ=additive_genes_down$logFC.TAZ

#extratcing column 3
LFC_iso=additive_genes_down$logFC.isotretinoin

# printing difference among
# two columns
difference_LFC_TAZ <- print(LFC_iso_and_TAZ-LFC_TAZ)

difference_LFC_iso <- print(LFC_iso_and_TAZ-LFC_iso)

difference_LFC_total <- c(difference_LFC_TAZ, difference_LFC_iso)
mean(difference_LFC_total)
range(difference_LFC_total)






#then do enrichment analysis on the genes--------------------------------------------------------------


#for up 
additive_genes_up <- additive_genes_up$Row.names


library(clusterProfiler)

GO_pathways <- gmtPathways("C:/Users/Christina/Documents/GSEA .gmt/c5.go.v2022.1.Hs.symbols.gmt")

    geneList<-additive_genes_up
    ego <- enrichGO(geneList, OrgDb = "org.Hs.eg.db",  pAdjustMethod='fdr', pvalueCutoff=0.05,keyType='SYMBOL')
    goplot(ego)
    
    dotplot(ego, x = "GeneRatio",
  color = "p.adjust",
  showCategory = 10,
  size = NULL,
  split = NULL,
  font.size = 12,
  title = "GO terms for additive upregulated genes",
  orderBy = "x",
  label_format = 30)
    
  

# Reactome ORA
#convert to symbols
symbols <- mapIds(org.Hs.eg.db, geneList, 'ENTREZID', 'SYMBOL')
symbols <- as.data.frame(symbols)
geneList <- symbols$symbols

#perform reactome ORA on cluster
reactome_ORA <- enrichPathway(gene=geneList, pvalueCutoff = 0.05, readable=TRUE)
head(reactome_ORA)

#get dotplot
   R <-   print(ggplot(reactome_ORA,aes(x=GeneRatio,y=Description,fill=p.adjust))+geom_bar(stat='identity')+theme_classic()+
              scale_x_continuous(expand=c(0,0))+labs(fill = 'p.adjust')+xlab('Gene Ratio')+ylab('')+theme(axis.text=element_text(size=12),axis.title=element_text(size=12))) + ggtitle("Reactome terms for additive upregulated genes")

    
# Hallmark ORA 
m_t2g <- msigdbr(species = "Homo sapiens", category = "H") %>% dplyr::select(gs_name, gene_symbol)
head(m_t2g)

#input gene list as geneIDs as some are lost when converting to entrez
em <- enricher(geneList, pvalueCutoff = 0.05,  pAdjustMethod = "BH", TERM2GENE=m_t2g)
head(em)

 try(H <- print(ggplot(em,aes(x=GeneRatio,y=Description,fill=p.adjust))+geom_bar(stat='identity')+theme_classic()+
              scale_x_continuous(expand=c(0,0)) +ggtitle("Hallmarks for additive upregulated genes")+labs(fill = 'p.adjust')+xlab('Gene Ratio')+ylab('')+theme(axis.text=element_text(size=12),axis.title=element_text(size=12) )))
     

#for down --------------------------------------------------------

#for up 
additive_genes_down <- additive_genes_down$Row.names


library(clusterProfiler)

GO_pathways <- gmtPathways("C:/Users/Christina/Documents/GSEA .gmt/c5.go.v2022.1.Hs.symbols.gmt")

    geneList<-additive_genes_down
    ego <- enrichGO(geneList, OrgDb = "org.Hs.eg.db",  pAdjustMethod='fdr', pvalueCutoff=0.05,keyType='SYMBOL')
    goplot(ego)
    
    dotplot(ego, x = "GeneRatio",
  color = "p.adjust",
  showCategory = 10,
  size = NULL,
  split = NULL,
  font.size = 12,
  title = "GO terms for additive downregulated genes",
  orderBy = "x",
  label_format = 30)
    
  

# Reactome ORA
#convert to symbols
symbols <- mapIds(org.Hs.eg.db, geneList, 'ENTREZID', 'SYMBOL')
symbols <- as.data.frame(symbols)
geneList <- symbols$symbols

#perform reactome ORA on cluster
reactome_ORA <- enrichPathway(gene=geneList, pvalueCutoff = 0.05, readable=TRUE)
head(reactome_ORA)

#get dotplot
   R <-   print(ggplot(reactome_ORA,aes(x=GeneRatio,y=Description,fill=p.adjust))+geom_bar(stat='identity')+theme_classic()+
              scale_x_continuous(expand=c(0,0))+labs(fill = 'p.adjust')+xlab('Gene Ratio')+ylab('')+theme(axis.text=element_text(size=12),axis.title=element_text(size=12))) + ggtitle("Reactome terms for additive downregulated genes")

    
# Hallmark ORA 
m_t2g <- msigdbr(species = "Homo sapiens", category = "H") %>% dplyr::select(gs_name, gene_symbol)
head(m_t2g)

#input gene list as geneIDs as some are lost when converting to entrez
em <- enricher(geneList, pvalueCutoff = 0.05,  pAdjustMethod = "BH", TERM2GENE=m_t2g)
head(em)

 try(H <- print(ggplot(em,aes(x=GeneRatio,y=Description,fill=p.adjust))+geom_bar(stat='identity')+theme_classic()+
              scale_x_continuous(expand=c(0,0)) +ggtitle("Hallmarks for additive downregulated genes")+labs(fill = 'p.adjust')+xlab('Gene Ratio')+ylab('')+theme(axis.text=element_text(size=12),axis.title=element_text(size=12) )))



```


#PCA
```{r}
# https://bioconductor.org/packages/devel/bioc/vignettes/PCAtools/inst/doc/PCAtools.html

#load
library("factoextra")
library("FactoMineR")
  library(PCAtools)
library(magrittr)


#get normalised logged 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")
logCPM <- cpm(edgeRData, log = TRUE)
CPM <- logCPM
scaledata <- t(scale(t(CPM)))



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

  #get treatmenttimepoint as extra col in df
sampleID <- row.names(metadata)
metadata$Treatment_timepoint <- c()
for (i in 1:nrow(metadata)) {
  samplename <- sampleID[i]
  Treatment_timepoint <- if(grepl("Isotretinoin\\_and_TAZ_PT3", samplename)) {print("Isotretinoin_and_TAZ_PT3")} 
  else if(grepl("Isotretinoin\\_and_TAZ_PT7", samplename)) {print("Isotretinoin_and_TAZ_PT7")}
  else if(grepl("Isotretinoin\\_and_TAZ_PT10", samplename)) {print("Isotretinoin_and_TAZ_PT10")}
   else if(grepl("Isotretinoin\\_and_GSK_PT3", samplename)) {print("Isotretinoin_and_GSK_PT3")}
   else if(grepl("Isotretinoin\\_and_GSK_PT7", samplename)) {print("Isotretinoin_and_GSK_PT7")}
    else if(grepl("Isotretinoin\\_and_GSK_PT10", samplename)) {print("Isotretinoin_and_GSK_PT10")}
      else if(grepl("DMSO_PT3", samplename)) {print("DMSO_PT3")}
  else if(grepl("DMSO_PT7", samplename)) {print("DMSO_PT7")}
  else if(grepl("DMSO_PT10", samplename)) {print("DMSO_PT10")}
  else if(grepl("GSK_PT3", samplename)) {print("GSK_PT3")}
  else if(grepl("GSK_PT7", samplename)) {print("GSK_PT7")}
  else if(grepl("GSK_PT10", samplename)) {print("GSK_PT10")}
  else if(grepl("Isotretinoin_PT3", samplename)) {print("Isotretinoin_PT3")}
  else if(grepl("Isotretinoin_PT7", samplename)) {print("Isotretinoin_PT7")}
  else if(grepl("Isotretinoin_PT10", samplename)) {print("Isotretinoin_PT10")}
  else if(grepl("TAZ_PT3", samplename)) {print("TAZ_PT3")}
  else if(grepl("TAZ_PT7", samplename)) {print("TAZ_PT7")}
  else if(grepl("TAZ_PT10", samplename)) {print("TAZ_PT10")}
  metadata$Treatment_timepoint[i] <- Treatment_timepoint
}

  #get treatment as col
sampleID <- row.names(metadata)
metadata$Treatment <- c()
for (i in 1:nrow(metadata)) {
  samplename <- sampleID[i]
  Treatment <- if(grepl("_and_TAZ", samplename)) {print("Isotretinoin_and_TAZ")
  } else if(grepl("_and_GSK", samplename)) {print("Isotretinoin_and_GSK")
    } else if (grepl("GSK", samplename)) {print("GSK")
  } else if(grepl("Isotretinoin", samplename)) {print("Isotretinoin")}
  else if(grepl("DMSO", samplename)) {print("DMSO")
  } 
  else {print("TAZ")}
  metadata$Treatment[i] <- Treatment
}

  #get timepoint as col
sampleID <- row.names(metadata)
metadata$Time <- c()
for (i in 1:nrow(metadata)) {
  samplename <- sampleID[i]
  Time <- if(grepl("PT3", samplename)) {print("PT3")
  } else if(grepl("PT7", samplename)) {print("PT7")}
  else if(grepl("PT10", samplename)) {print("PT10")}
  metadata$Time[i] <- Time
}


  #get timepoint as col
sampleID <- row.names(metadata)
metadata$Replicate <- c()
for (i in 1:nrow(metadata)) {
  samplename <- sampleID[i]
  Replicate <- if(grepl("1$", samplename)) {print("1")
  } else if(grepl("2$", samplename)) {print("2")}
  else if(grepl("3$", samplename)) {print("3")}
  metadata$Replicate[i] <- Replicate
}

  #reorder levels for plots
 metadata$Treatment<- factor(metadata$Treatment,  c("DMSO", "Isotretinoin", "TAZ", "GSK", "Isotretinoin_and_TAZ", "Isotretinoin_and_GSK"))
 metadata$Time<- factor(metadata$Time,  c("PT3", "PT7", "PT10"))
 

#get pca data
## -- removing the lower 10% of variables based on variance
 p <- pca(scaledata, metadata = metadata, removeVar = 0.1)

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

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

  #shape by timepoint
  tiff(file="biplot_pc1vspc2_time.tiff",
width=8, height=5, units="in", res=800)
  biplot(p,
          lab = NULL,
    colby = 'Treatment', colkey = c("DMSO"= "#000000", "Isotretinoin"= "#E69F00", "TAZ"= "#56B4E9", "GSK"= "#F0E442", "Isotretinoin_and_TAZ"= "#D55E00", "Isotretinoin_and_GSK"= "#CC79A7"),
    hline = 0, vline = 0,
    legendPosition = 'right', 
    ellipse = FALSE,
      ellipseType = 't',
      ellipseLevel = 0.95,
      ellipseFill = TRUE,
      ellipseAlpha = 1/4,
      ellipseLineSize = 1.0,
    shape = 'Time',
    pointSize = 5,
    title = 'PCA bi-plot with shape representing timepoint',
    subtitle = 'PC1 versus PC2',)
  dev.off()
  
  #shape by timepoint
    tiff(file="biplot_pc1vspc3_time.tiff",
width=8, height=5, units="in", res=800)
  biplot(p,
         x= "PC1",
         y="PC3",
          lab = NULL,
    colby = 'Treatment', colkey = c("DMSO"= "#000000", "Isotretinoin"= "#E69F00", "TAZ"= "#56B4E9", "GSK"= "#F0E442", "Isotretinoin_and_TAZ"= "#D55E00", "Isotretinoin_and_GSK"= "#CC79A7"),
    hline = 0, vline = 0,
    legendPosition = 'right', 
    ellipse = FALSE,
      ellipseType = 't',
      ellipseLevel = 0.95,
      ellipseFill = TRUE,
      ellipseAlpha = 1/4,
      ellipseLineSize = 1.0,
    shape = 'Time',
    pointSize = 5,
    title = 'PCA bi-plot with shape representing timepoint',
    subtitle = 'PC1 versus PC3',)
  dev.off()
  
   #shape by timepoint
    tiff(file="biplot_pc2vspc3_time.tiff",
width=8, height=5, units="in", res=800)
  biplot(p,
         x= "PC2",
         y="PC3",
          lab = NULL,
    colby = 'Treatment', colkey = c("DMSO"= "#000000", "Isotretinoin"= "#E69F00", "TAZ"= "#56B4E9", "GSK"= "#F0E442", "Isotretinoin_and_TAZ"= "#D55E00", "Isotretinoin_and_GSK"= "#CC79A7"),
    hline = 0, vline = 0,
    legendPosition = 'right', 
    ellipse = FALSE,
      ellipseType = 't',
      ellipseLevel = 0.95,
      ellipseFill = TRUE,
      ellipseAlpha = 1/4,
      ellipseLineSize = 1.0,
    shape = 'Time',
    pointSize = 5,
    title = 'PCA bi-plot with shape representing timepoint',
    subtitle = 'PC2 versus PC3',)
  dev.off()

 #shape by replicate---------------------------------------------------
     tiff(file="biplot_pc1vspc2_rep.tiff",
width=8, height=5, units="in", res=800)
  biplot(p,
         
         x= "PC1",
         y="PC2",
          lab = NULL,
    colby = 'Treatment', colkey = c("DMSO"= "#000000", "Isotretinoin"= "#E69F00", "TAZ"= "#56B4E9", "GSK"= "#F0E442", "Isotretinoin_and_TAZ"= "#D55E00", "Isotretinoin_and_GSK"= "#CC79A7"),
    hline = 0, vline = 0,
    legendPosition = 'right', 
    ellipse = FALSE,
      ellipseType = 't',
      ellipseLevel = 0.95,
      ellipseFill = TRUE,
      ellipseAlpha = 1/4,
      ellipseLineSize = 1.0,
    pointSize = 5,
    shape = 'Replicate',
    title = 'PCA bi-plot with shape representing replicate',
    subtitle = 'PC1 versus PC2',)
  dev.off()
  
   #shape by replicate
  
     tiff(file="biplot_pc2vspc3_rep.tiff",
width=8, height=5, units="in", res=800)
  biplot(p,
         x= "PC2",
         y="PC3",
          lab = NULL,
    colby = 'Treatment', colkey = c("DMSO"= "#000000", "Isotretinoin"= "#E69F00", "TAZ"= "#56B4E9", "GSK"= "#F0E442", "Isotretinoin_and_TAZ"= "#D55E00", "Isotretinoin_and_GSK"= "#CC79A7"),
    hline = 0, vline = 0,
    legendPosition = 'right', 
    ellipse = FALSE,
      ellipseType = 't',
      ellipseLevel = 0.95,
      ellipseFill = TRUE,
      ellipseAlpha = 1/4,
      ellipseLineSize = 1.0,
    pointSize = 5,
    shape = 'Replicate',
    title = 'PCA bi-plot with shape representing replicate',
    subtitle = 'PC2 versus PC3',)
  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:7)),
      triangle = TRUE, trianglelabSize = 12,
      hline = 0, vline = 0,
      titleLabSize= 20,
      pointSize = 2,
      shape = 'Time',
      gridlines.major = FALSE, gridlines.minor = FALSE,
     colby = 'Treatment', colkey = c("DMSO"= "#000000", "Isotretinoin"= "#E69F00", "TAZ"= "#56B4E9", "GSK"= "#F0E442", "Isotretinoin_and_TAZ"= "#D55E00", "Isotretinoin_and_GSK"= "#CC79A7"),
      legendPosition = "none",
      title = 'Pairs plot', plotaxes = FALSE,
      margingaps = unit(c(-0.01, -0.01, -0.01, -0.01), 'cm'))
  dev.off()
  
```

#rrvgo 
```{r}

go_analysis <- read.delim(system.file("extdata/example.txt", package="rrvgo"))

go_analysis <- read.csv("GO_GSKvsDMSO_PT3.csv", header = TRUE)

go_analysis <- subset(go_analysis, startsWith(as.character(pathway),"GOBP"))






simMatrix <- calculateSimMatrix(go_analysis$pathway,
                                orgdb="org.Hs.eg.db",
                                ont="BP",
                                method="Rel")

scores <- setNames(go_analysis$padj, go_analysis$ID)
reducedTerms <- reduceSimMatrix(simMatrix,
                                scores,
                                threshold=0.7,
                                orgdb="org.Hs.eg.db")

```


