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

include=FALSE:

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

#load libraries
```{r}
library(dplyr)
library(readr)
library(tidyverse)
library(RColorBrewer)
library(GSVA)

library(fgsea)
library(edgeR)
library("gplots")
library("ggrepel")
library("ggplot2")
library('knitr')
```

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

setwd("~/RMS RNA seq/Luciles data/All_Raw_Reads")
# #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")
# 
#read in count matrix
counts <- read.csv('count_matrix_Luciles.csv')

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

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

#write.csv(rawData, "count_matrix_Luciles.csv")



```



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

# DESCRIBE EXPERIMENTAL GROUPS
```{r}
########## DESCRIBE EXPERIMENTAL GROUPS #########################

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

#for sample treatment

sampleID <- colnames(rawData)
treatment <- c()
for (i in 1:ncol(rawData)) {
  samplename <- sampleID[i]
  group <- if(grepl("RH4ctrl", samplename)) {print("RH4ctrl")
  } else if(grepl("RH4ifo", samplename)) {print("RH4ifo")}
  else if(grepl("RH4vcr", samplename)) {print("RH4vcr")
  } else if(grepl("RMSYMctrl", samplename)) {print("RMSYMctrl")
  } else if(grepl("RMSYMifo", samplename)) {print("RMSYMifo")
  } else {print("RMSYMvcr")}
  treatment[i] <- group
}



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


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



```

#Filtering to remove genes with low expression

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


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

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

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

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


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

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

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


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

```


#Create DGEList object
```{r}

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



#Density plots

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

```{r}

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

print(p)


#first peak has been greatly reduced
```


QC Plots before normalisation

Barplot of raw reads

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


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

#make a bar plot of the number of reads, i.e. library size, in each sample
png(filename="lib_size.png", width=900, bg="white")
par(mar=c(11,6,4,1)+.1)
barplot(data$lib.size,names=data$rowname, ylim=c(0,60000000),las=2,   cex.names = 0.8, ylab="Library size", xlab = "", cex.axis = 0.8);
title("Barplot of library sizes")
title(xlab="Sample", line=5, cex.lab=1)
dev.off()



```

#Boxplot of logCPM before normalisation

```{r}

#Sample outliers can be detected based on the distribution of read counts between samples which are visualised using a box and whisker plot of log2 counts per milliom (CPM) for each sample. CPM is used to correct for differences in library size (ie depth of sequencing per sample) which are logged because the count data are not normally distributed. We will use the cpm function in edgeR to get log2(CPM) which also adds a small offset to avoid taking a log of zero.
# Get log2 counts per million
logcounts <- cpm(edgeRData,log=TRUE)
# Check distributions of samples using boxplots
boxplot(logcounts, xlab="", ylab="Log2 counts per million",las=2)
# Let's add a blue horizontal line that corresponds to the median logCPM
abline(h=median(logcounts),col="blue")
title(xlab="Sample", line=6.5, cex.lab=1.1)
title("Boxplots of log-2 CPMs before normalisation")
```

#MDS plot

```{r}


treatment_names <- c("RH4ctrl", "RH4ifo", "RH4vcr", "RMSYMctrl", "RMSYMifo", "RMSYMvcr")


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

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

```

#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("orange",  "red", "green",  "cyan",  "#FF69B4",  "purple")[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}
design <- model.matrix(~0+group, data=edgeRData$samples)
colnames(design) <- levels(edgeRData$samples$group)
design

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


# note values of common dispersion and BCV
```


#Boxplot of logCPM

```{r}
#Sample outliers can be detected based on the distribution of read counts between samples which are visualised using a box and whisker plot of log2 counts per milliom (CPM) for each sample. CPM is used to correct for differences in library size (ie depth of sequencing per sample) which are logged because the count data are not normally distributed. We will use the cpm function in edgeR to get log2(CPM) which also adds a small offset to avoid taking a log of zero.
# Get log2 counts per million
logcounts <- cpm(edgeRData,log=TRUE)
# Check distributions of samples using boxplots
boxplot(logcounts, xlab="", ylab="Log2 counts per million",las=2)
# Let's add a blue horizontal line that corresponds to the median logCPM
abline(h=median(logcounts),col="blue")
title(xlab="Sample", line=6, cex.lab=1.1)
title("Boxplots of log-2 CPMs after normalisation")
```

#BCV plot

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

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

 
  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}

design <- model.matrix(~0+group, data=edgeRData$samples)
colnames(design) <- levels(edgeRData$samples$group)
design


# create contrast matrix
my.contrasts <- makeContrasts(RH4ctrl_RH4ifo=RH4ifo-RH4ctrl,
  RH4ctrl_RH4vcr=RH4vcr-RH4ctrl,RMSYMctrl_RMSYMifo=RMSYMifo-RMSYMctrl,
  RMSYMctrl_RMSYMvcr=RMSYMvcr-RMSYMctrl, 
  RH4ctrl_RMSYMctrl=RH4ctrl-RMSYMctrl,
  levels=design)



fit <- glmQLFit(edgeRData, design)
RH4ctrl_RH4ifo <- glmQLFTest(fit, contrast=my.contrasts[,"RH4ctrl_RH4ifo"])
RH4ctrl_RH4vcr <- glmQLFTest(fit, contrast=my.contrasts[,"RH4ctrl_RH4vcr"])
RMSYMctrl_RMSYMifo <- glmQLFTest(fit, contrast=my.contrasts[,"RMSYMctrl_RMSYMifo"])
RMSYMctrl_RMSYMvcr <- glmQLFTest(fit, contrast=my.contrasts[,"RMSYMctrl_RMSYMvcr"])
RH4ctrl_RMSYMctrl <- glmQLFTest(fit, contrast=my.contrasts[,"RH4ctrl_RMSYMctrl"])

```


```{r}
# 
# ###################### DGE- exact tests ###############################################
# 
# #RH4ifo vs control
# RH4ctrl_RH4ifo<- exactTest(edgeRData, pair=c("RH4ctrl", "RH4ifo"))
# 
# #RH4vcr vs control
# RH4ctrl_RH4vcr<- exactTest(edgeRData, pair=c("RH4ctrl", "RH4vcr"))
# 
# #RMSYMifo vs control
# RMSYMctrl_RMSYMifo<- exactTest(edgeRData, pair=c("RMSYMctrl", "RMSYMifo"))
# 
# #RMSYMvcr vs control
# RMSYMctrl_RMSYMvcr<- exactTest(edgeRData, pair=c("RMSYMctrl", "RMSYMvcr"))

```


#DGE visualisation


```{r}

dge_files <- list(RH4ctrl_RH4ifo, RH4ctrl_RH4vcr, RMSYMctrl_RMSYMifo,RMSYMctrl_RMSYMvcr)



filenames <- c("RH4ctrl_RH4ifo", "RH4ctrl_RH4vcr", "RMSYMctrl_RMSYMifo" , "RMSYMctrl_RMSYMvcr", "RH4ctrl_RMSYMctrl")


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



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

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

  # Volcano Plots-----------------------------------------------
  

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

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

  )
  top_genes %>%
    kable()

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


  #vol plot with labels of specific genes---------------------------------------------------------------

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


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

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

 #plot vol plot highlighting genelist and adding labels
 p3 <- ggplot(data = vol_plot_data, # Original data  
       aes(x = logFC, y = -log(FDR,10))) + 
  geom_hline(yintercept = -log10(0.05),
             linetype = "dashed") + 
  geom_vline(xintercept = c(-1.5, 1.5),
             linetype = "dashed") +  
 scale_colour_manual(values = cols)+
    ggtitle(filenames[[i]]) +
    geom_point(aes(colour = Expression), 
             alpha = 0.8, 
             shape = 16,
             size = 1)  + geom_point(data = ils, # New layer containing data subset il_genes
             size = 2,
             shape = 21,
             fill = "yellow",
             colour = "black") +  geom_label_repel(data = ils, label.size = 0.5, force = 0.5, max.overlaps = 15, 
                     mapping = aes(logFC, -log(FDR,10), label = rownames(ils)),
                     size = 3.8) # Add labels last to appear as the toplayer
                  
 print(p3)
  
 ggsave(p3, file=paste0("Volcano_plot_MDR_genes", j,".png"), width = 25, height = 25, units = "cm")
 # dev.off()
  
}

```


#Number of DEGs

```{r}

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


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

```{r}

######### for RH4 VCR and RH4 IFO #######################
RH4ctrl_RH4ifo$table$FDR <- p.adjust(RH4ctrl_RH4ifo$table$PValue, method = "BH")
 RH4ctrl_RH4ifo$table <- RH4ctrl_RH4ifo$table %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "significant",
                           logFC < -1 & FDR < 0.05 ~ "significant",
                           TRUE ~ "None") )
 
RH4ctrl_RH4vcr$table$FDR <- p.adjust(RH4ctrl_RH4vcr$table$PValue, method = "BH")
 RH4ctrl_RH4vcr$table <- RH4ctrl_RH4vcr$table %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "significant",
                           logFC < -1 & FDR < 0.05 ~ "significant",
                           TRUE ~ "None") )
 
df <- merge(RH4ctrl_RH4vcr$table, RH4ctrl_RH4ifo$table, by = 'row.names', all = TRUE)

df <- df %>%
   mutate(sig = case_when(Expression.x == "significant" & Expression.y == "significant" ~ "Both",
                          Expression.x == "significant" ~ "RH4 VCR-resistant",
                           Expression.y == "significant" ~ "RH4 IFO-resistant",
                          
                           TRUE ~ "None") )
  


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("orange", "grey", "red", "green" )) + ggtitle("Differentially expressed genes for RH4 VCR-resistant vs control and RH4 IFO-resistant vs control") +
  xlab("Log-2 fold change RH4 IFO-resistant vs control") + ylab("Log-2 fold change RH4 VCR-resistant vs control") +theme(axis.text=element_text(size=12),
        axis.title=element_text(size=12)) + guides(color=guide_legend(title="Condition in which the gene is differentially expressed")) 

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

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




# for RMSYM VCR and RMSYM IFO --------------------------------------------------------
RMSYMctrl_RMSYMifo$table$FDR <- p.adjust(RMSYMctrl_RMSYMifo$table$PValue, method = "BH")
 RMSYMctrl_RMSYMifo$table <- RMSYMctrl_RMSYMifo$table %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "significant",
                           logFC < -1 & FDR < 0.05 ~ "significant",
                           TRUE ~ "None") )
 
RMSYMctrl_RMSYMvcr$table$FDR <- p.adjust(RMSYMctrl_RMSYMvcr$table$PValue, method = "BH")
 RMSYMctrl_RMSYMvcr$table <- RMSYMctrl_RMSYMvcr$table %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "significant",
                           logFC < -1 & FDR < 0.05 ~ "significant",
                           TRUE ~ "None") )
 
df <- merge(RMSYMctrl_RMSYMvcr$table, RMSYMctrl_RMSYMifo$table, by = 'row.names', all = TRUE)

df <- df %>%
   mutate(sig = case_when(Expression.x == "significant" & Expression.y == "significant" ~ "Both",
                          Expression.x == "significant" ~ "RMSYM VCR-resistant",
                           Expression.y == "significant" ~ "RMSYM IFO-resistant",
                          
                           TRUE ~ "None") )
  


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("violet", "grey", "deeppink1", "purple" )) + ggtitle("Differentially expressed genes for RMSYM VCR-resistant vs control and RMSYM IFO-resistant vs control") +
  xlab("Log-2 fold change RMSYM IFO-resistant vs control") + ylab("Log-2 fold change RMSYM VCR-resistant vs control") +theme(axis.text=element_text(size=12),
        axis.title=element_text(size=12)) + guides(color=guide_legend(title="Condition in which the gene is differentially expressed")) 

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)





# for RH4 IFO and RMSYM IFO --------------------------------------------------------
RMSYMctrl_RMSYMifo$table$FDR <- p.adjust(RMSYMctrl_RMSYMifo$table$PValue, method = "BH")
 RMSYMctrl_RMSYMifo$table <- RMSYMctrl_RMSYMifo$table %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "significant",
                           logFC < -1 & FDR < 0.05 ~ "significant",
                           TRUE ~ "None") )
 
RH4ctrl_RH4ifo$table$FDR <- p.adjust(RH4ctrl_RH4ifo$table$PValue, method = "BH")
 RH4ctrl_RH4ifo$table <- RH4ctrl_RH4ifo$table %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "significant",
                           logFC < -1 & FDR < 0.05 ~ "significant",
                           TRUE ~ "None") )
 
df <- merge(RH4ctrl_RH4ifo$table, RMSYMctrl_RMSYMifo$table, by = 'row.names', all = TRUE)

df <- df %>%
   mutate(sig = case_when(Expression.x == "significant" & Expression.y == "significant" ~ "Both",
                          Expression.x == "significant" ~ "RH4 IFO-resistant",
                           Expression.y == "significant" ~ "RMSYM IFO-resistant",
                          
                           TRUE ~ "None") )
  
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("orange", "grey", "deeppink1", "red" )) + ggtitle("Differentially expressed genes for RH4 IFO-resistant vs control and RMSYM IFO-resistant vs control") +
  xlab("Log-2 fold change RMSYM IFO-resistant vs control") + ylab("Log-2 fold change RH4 IFO-resistant vs control") +theme(axis.text=element_text(size=12),
        axis.title=element_text(size=12)) + guides(color=guide_legend(title="Condition in which the gene is differentially expressed")) 

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)


# for RH4 VCR and RMSYM VCR --------------------------------------------------------
RMSYMctrl_RMSYMvcr$table$FDR <- p.adjust(RMSYMctrl_RMSYMvcr$table$PValue, method = "BH")
 RMSYMctrl_RMSYMvcr$table <- RMSYMctrl_RMSYMvcr$table %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "significant",
                           logFC < -1 & FDR < 0.05 ~ "significant",
                           TRUE ~ "None") )
 
RH4ctrl_RH4vcr$table$FDR <- p.adjust(RH4ctrl_RH4vcr$table$PValue, method = "BH")
 RH4ctrl_RH4vcr$table <- RH4ctrl_RH4vcr$table %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "significant",
                           logFC < -1 & FDR < 0.05 ~ "significant",
                           TRUE ~ "None") )
 
df <- merge(RH4ctrl_RH4vcr$table, RMSYMctrl_RMSYMvcr$table, by = 'row.names', all = TRUE)

df <- df %>%
   mutate(sig = case_when(Expression.x == "significant" & Expression.y == "significant" ~ "Both",
                          Expression.x == "significant" ~ "RH4 VCR-resistant",
                           Expression.y == "significant" ~ "RMSYM VCR-resistant",
                          
                           TRUE ~ "None") )
  
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("deepskyblue", "grey", "purple", "green" )) + ggtitle("Differentially expressed genes for RH4 VCR-resistant vs control and RMSYM VCR-resistant vs control") +
  xlab("Log-2 fold change RMSYM VCR-resistant vs control") + ylab("Log-2 fold change RH4 VCR-resistant vs control") +theme(axis.text=element_text(size=12),
        axis.title=element_text(size=12)) + guides(color=guide_legend(title="Condition in which the gene is differentially expressed")) 

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)


```

  
  

#Compare lists of DEGs to find overlapping DEGs in ALL samples


```{r}



#get list of DEGS

#for RH4 control vs vcr
data <- RH4ctrl_RH4vcr$table
   data$FDR <- p.adjust(data$PValue, method = "BH")
   data$geneID <- row.names(data)
  data <- data %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "Up-regulated",
                           logFC < -1 & FDR < 0.05 ~ "Down-regulated",
                           TRUE ~ "Unchanged"))
DEG_list_RH4ctrl_RH4vcr <- data[data$Expression == 'Up-regulated' | data$Expression == 'Down-regulated',]
UP_DEG_list_RH4ctrl_RH4vcr <- data[data$Expression == 'Up-regulated',]
DOWN_DEG_list_RH4ctrl_RH4vcr <- data[data$Expression == 'Down-regulated',]

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

# for RH4 control vs ifo
data <- RH4ctrl_RH4ifo$table
   data$FDR <- p.adjust(data$PValue, method = "BH")
   data$geneID <- row.names(data)
  data <- data %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "Up-regulated",
                           logFC < -1 & FDR < 0.05 ~ "Down-regulated",
                           TRUE ~ "Unchanged"))
DEG_list_RH4ctrl_RH4ifo <- data[data$Expression == 'Up-regulated' | data$Expression == 'Down-regulated',]
UP_DEG_list_RH4ctrl_RH4ifo <- data[data$Expression == 'Up-regulated',]
DOWN_DEG_list_RH4ctrl_RH4ifo <- data[data$Expression == 'Down-regulated',]

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

#for RMSYM control vs vcr
data <- RMSYMctrl_RMSYMvcr$table
   data$FDR <- p.adjust(data$PValue, method = "BH")
   data$geneID <- row.names(data)
  data <- data %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "Up-regulated",
                           logFC < -1 & FDR < 0.05 ~ "Down-regulated",
                           TRUE ~ "Unchanged"))
DEG_list_RMSYMctrl_RMSYMvcr <- data[data$Expression == 'Up-regulated' | data$Expression == 'Down-regulated',]
UP_DEG_list_RMSYMctrl_RMSYMvcr <- data[data$Expression == 'Up-regulated',]
DOWN_DEG_list_RMSYMctrl_RMSYMvcr <- data[data$Expression == 'Down-regulated',]

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

#for RMSYM control vs ifo
data <- RMSYMctrl_RMSYMifo$table
   data$FDR <- p.adjust(data$PValue, method = "BH")
   data$geneID <- row.names(data)
  data <- data %>% 
  mutate(
    Expression = case_when(logFC > 1 & FDR < 0.05 ~ "Up-regulated",
                           logFC < -1 & FDR < 0.05 ~ "Down-regulated",
                           TRUE ~ "Unchanged"))
DEG_list_RMSYMctrl_RMSYMifo <- data[data$Expression == 'Up-regulated' | data$Expression == 'Down-regulated',]
UP_DEG_list_RMSYMctrl_RMSYMifo <- data[data$Expression == 'Up-regulated',]
DOWN_DEG_list_RMSYMctrl_RMSYMifo <- data[data$Expression == 'Down-regulated',]

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



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

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



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


```


#Upset plots

```{r upset plots}
#Upset plot 
# load libraries
library(plyr)
library(reshape2)
library(UpSetR)
library(ComplexHeatmap)
library("viridis")

#make sure directory is set to where list of signif gene lists

  # make a list of input files to be read that contain list of sig genes
  filelist = list.files(pattern = "*.csv")
  
  # make a 3 column table of listname,gene,1
  res<-lapply(filelist, function(x){
    data.frame(
      set=x,
      geneID=as.character(read.table(x)[,1]),
      val=1)
  })
  
  res<-ldply(res)
  
  # turn the 3 column long table into a wide
  res1<-acast(res,geneID~set,value.var="val",fill=0) 
  
  
  # force as dataframe
  res1<-as.data.frame(res1)
  
  # 1st column must be a name
  res1$name=rownames(res1)
  
  # rearrange columns
  res2=res1[,c(ncol(res1),1:(ncol(res1)-1))]
  
  #change column name
  old_colnames <- colnames(res2)
  x <- gsub('.csv','',old_colnames)
   y <-  gsub('DEG_list_','',x)
   
  colnames(res2) <- c(y)
  
 
  #make combination matrix
m = make_comb_mat(res2, mode = "intersect", universal_set = )
  m1 = make_comb_mat(res2, mode = "distinct",  universal_set =)
  

#add colour palette
  y <- viridis(6)
  #plot upset
  Intersect_upset_plot <- UpSet(m, column_title = "intersect mode", pt_size = unit(5, "mm"), lwd = 3, bg_col = "#F0F0FF", bg_pt_col = "#CCCCFF" , 
        comb_col = y[comb_degree(m)])
  
 # grid.text(,x = 0.65, y=0.95, gp=gpar(fontsize=20))
   print(Intersect_upset_plot)
   
    Distinct_upset_plot <- UpSet(m1, column_title = "distinct mode", pt_size = unit(5, "mm"), lwd = 3, bg_col = "#F0F0FF", bg_pt_col = "#CCCCFF" , 
        comb_col = y[comb_degree(m)])
  
 # grid.text(,x = 0.65, y=0.95, gp=gpar(fontsize=20))
   print(Distinct_upset_plot)


```


```{r}

###################### For VCR UP ############################
#put all data frames into list
df_list <- list(UP_DEG_list_RH4ctrl_RH4vcr, UP_DEG_list_RMSYMctrl_RMSYMvcr)
#merge all data frames in list
overlapping_DEGS <- df_list %>% reduce(merge, by="geneID")
dim(overlapping_DEGS)


###################### For IFO UP ############################
#put all data frames into list
df_list <- list(UP_DEG_list_RH4ctrl_RH4ifo, UP_DEG_list_RMSYMctrl_RMSYMifo)
#merge all data frames in list
overlapping_DEGS <- df_list %>% reduce(merge, by="geneID")
dim(overlapping_DEGS)



###################### For RH4 UP ############################
#put all data frames into list
df_list <- list(UP_DEG_list_RH4ctrl_RH4vcr, UP_DEG_list_RH4ctrl_RH4ifo)
#merge all data frames in list
overlapping_DEGS <- df_list %>% reduce(merge, by="geneID")
dim(overlapping_DEGS)

###################### For RMSY UP ############################
#put all data frames into list
df_list <- list(UP_DEG_list_RMSYMctrl_RMSYMifo, UP_DEG_list_RMSYMctrl_RMSYMvcr)
#merge all data frames in list
overlapping_DEGS <- df_list %>% reduce(merge, by="geneID")
dim(overlapping_DEGS)


```


```{r}
###################### For VCR DOWN ############################
#put all data frames into list
df_list <- list(DOWN_DEG_list_RH4ctrl_RH4vcr, DOWN_DEG_list_RMSYMctrl_RMSYMvcr)
#merge all data frames in list
overlapping_DEGS <- df_list %>% reduce(merge, by="geneID")
dim(overlapping_DEGS)


###################### For IFO DOWN ############################
#put all data frames into list
df_list <- list(DOWN_DEG_list_RH4ctrl_RH4ifo, DOWN_DEG_list_RMSYMctrl_RMSYMifo)
#merge all data frames in list
overlapping_DEGS <- df_list %>% reduce(merge, by="geneID")
dim(overlapping_DEGS)



###################### For RH4 DOWN ############################
#put all data frames into list
df_list <- list(DOWN_DEG_list_RH4ctrl_RH4vcr, DOWN_DEG_list_RH4ctrl_RH4ifo)
#merge all data frames in list
overlapping_DEGS <- df_list %>% reduce(merge, by="geneID")
dim(overlapping_DEGS)

###################### For RMSY DOWN ############################
#put all data frames into list
df_list <- list(DOWN_DEG_list_RMSYMctrl_RMSYMifo, DOWN_DEG_list_RMSYMctrl_RMSYMvcr)
#merge all data frames in list
overlapping_DEGS <- df_list %>% reduce(merge, by="geneID")
dim(overlapping_DEGS)

```


```{r}
###################### For VCR ALL ############################
#put all data frames into list
df_list <- list(DEG_list_RH4ctrl_RH4vcr, DEG_list_RMSYMctrl_RMSYMvcr)
#merge all data frames in list
overlapping_DEGS <- df_list %>% reduce(merge, by="geneID")
dim(overlapping_DEGS)


###################### For IFO DOWN ############################
#put all data frames into list
df_list <- list(DEG_list_RH4ctrl_RH4ifo, DEG_list_RMSYMctrl_RMSYMifo)
#merge all data frames in list
overlapping_DEGS <- df_list %>% reduce(merge, by="geneID")
dim(overlapping_DEGS)



###################### For RH4 DOWN ############################
#put all data frames into list
df_list <- list(DEG_list_RH4ctrl_RH4vcr, DEG_list_RH4ctrl_RH4ifo)
#merge all data frames in list
overlapping_DEGS <- df_list %>% reduce(merge, by="geneID")
dim(overlapping_DEGS)

###################### For RMSY DOWN ############################
#put all data frames into list
df_list <- list(DEG_list_RMSYMctrl_RMSYMifo, DEG_list_RMSYMctrl_RMSYMvcr)
#merge all data frames in list
overlapping_DEGS <- df_list %>% reduce(merge, by="geneID")
dim(overlapping_DEGS)

```

#Sense check genes


```{r}

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

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

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

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

  
 


#Number of DEGs

```{r}

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

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

p


```





#GSEA REACTOME
```{r}

################################ REACTOME #################################


# Check contrasts:
print(colnames(design))

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

filenames <- c( "RH4ctrl_RH4vcr", "RH4ctrl_RH4ifo", "RMSYMctrl_RMSYMifo" , "RMSYMctrl_RMSYMvcr")



for (i in filenames){
# Get pathways from gmt file
Reactome_pathways <- gmtPathways("C:/Users/Christina/Documents/GSEA .gmt/c2.cp.reactome.v2022.1.Hs.symbols.gmt")
#create index 
reactome.ind <- ids2indices( Reactome_pathways, rownames(edgeRData$counts))
# Run analysis:
gst.camera <- camera.DGEList(edgeRData, design= design, index=reactome.ind,contrast=my.contrasts[, i], inter.gene.cor=0.05)

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

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

#  #if more than 20 significant terms only plot the 20 most significant
#     if(nrow(gsea_df)>=20){
#       p <- ggplot(gsea_df[((nrow(gsea_df)-19):nrow(gsea_df)),], aes(x = fct_reorder(pathway, - 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("Reactome_barplot_", i,".png"), width = 35, height = 20, units = "cm")) 
#     }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=6))
# try(ggsave(p, file=paste0("Reactome_barplot_", i,".png"), width = 35, height = 20, units = "cm"))
# 
#     }

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

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

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

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

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


#GSEA HALLMARKS

```{r}

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


filenames <- c( "RH4ctrl_RH4vcr", "RH4ctrl_RH4ifo", "RMSYMctrl_RMSYMifo" , "RMSYMctrl_RMSYMvcr")

#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){
# 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 <- as.character(row.names(gsea_df))

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

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

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

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

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


```



#GSEA GO

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

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


for (i in filenames){
# Get pathways from gmt file
GO_pathways <- gmtPathways("C:/Users/Christina/Documents/GSEA .gmt/c5.go.v2022.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),]

 #if more than 20 significant terms only plot the 20 most significant
    if(nrow(gsea_df)>=20){
      p <- ggplot(gsea_df[((nrow(gsea_df)-19):nrow(gsea_df)),], aes(x = fct_reorder(pathway, - 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 = 35, height = 20, units = "cm")) 
    }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=15))
try(ggsave(p, file=paste0("GO_barplot_", i,".png"), width = 45, height = 20, units = "cm"))

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

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

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

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

```


Find overlapping GSEA GO terms


```{r}

############################### for "RH4ctrl_RH4vcr"###################################
# 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=contrasts[, "RH4ctrl_RH4vcr"], inter.gene.cor=0.05)

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

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




###################################### for "RMSYMctrl_RMSYMifo"############################
# 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=contrasts[, "RMSYMctrl_RMSYMifo"], inter.gene.cor=0.05)

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

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

############################ for "RMSYMctrl_RMSYMvcr" #############################################

# 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=contrasts[, "RMSYMctrl_RMSYMvcr"], inter.gene.cor=0.05)

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

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






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







```


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


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

colourBy <- top_10_GO$FDR

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

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

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

}
```






#fgsea for GSEA


```{r}


dge_files <- list(RH4ctrl_RH4ifo, RH4ctrl_RH4vcr, RMSYMctrl_RMSYMifo, RMSYMctrl_RMSYMvcr)

for (i in dge_files) {

  set.seed(42)
  fgseaRes <- fgsea(examplePathways, exampleRanks,  maxSize=500, nproc=1)
  
  
  # rank by LFC
  rankData <- i$table$logFC
  names(rankData) <- i$genes$genes
  head(rankData)
  
  # for hallmarks
#load pathways

#GO_pathways <- gmtPathways("C:/Users/Christina/Documents/Neuroblastoma RNAseq analysis/GSEA .gmt/c5.go.v2022.1.Hs.symbols.gmt")
reactome_pathways <- gmtPathways("C:/Users/Christina/Documents/Neuroblastoma RNAseq analysis/GSEA .gmt/c2.cp.reactome.v2022.1.Hs.symbols.gmt")

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



  # look at top 10 results
x <- fgseaRes[order(fgseaRes$padj),]
x <- x[1:10, ]
#print(i$comparison)
#print(x)



gsea_df <- x[padj < 0.05,]
gsea_df$padj <- as.numeric(gsea_df$padj)

print(gsea_df)


gsea_df <- gsea_df %>%
  dplyr::select(pathway, NES, padj) %>%
  dplyr::mutate(
    padj.group = cut(.$padj, breaks = c(-Inf, 0.25, Inf), labels = c(1, 0)),
    nes.group = cut(.$NES, breaks = c(-Inf, 0, Inf), labels = c(0, 1)),
    comb.group = paste0(padj.group, nes.group)
  ) %>%
  dplyr::mutate(group = dplyr::case_when(
    comb.group == "10" ~ "A", # A is nes<0 & p<0.25
    comb.group == "11" ~ "B", # B is nes>0 & p<0.25
    TRUE ~ "C"
  )) %>% # C is grey
  dplyr::arrange(NES) %>%
  dplyr::mutate(index = 1:dplyr::n())


p <- ggplot(gsea_df, aes(x = index, y = NES)) +
  geom_bar(aes(fill = gsea_df$padj), stat = "identity", width = 0.8) +
  scale_fill_continuous(low = "red", high = "blue", na.value = NA)+
  scale_x_discrete(expand = expansion(add = .5)) +
  scale_y_continuous(breaks = seq(
    floor(min(gsea_df$NES)), ceiling(max(gsea_df$NES)),
    ceiling((ceiling(max(gsea_df$NES)) - floor(min(gsea_df$NES))) / 6)
  )) +
  coord_flip()


pos_new <- sum(gsea_df$NES>0); neg_nes <- sum(gsea_df$NES<0)
if(neg_nes == 0 & pos_new != 0 ){
  p <- p +
    geom_text(
      data = subset(gsea_df, NES > 0),
      aes(x = index, y = 0, label = paste0("  ", pathway), color = padj.group),
      size = 2.5, hjust = "outward"
    )
}else if(neg_nes != 0 & pos_new == 0 ){
  p <- p +
    geom_text(
      data = subset(gsea_df, NES < 0),
      aes(x = index, y = 0, label = paste0("  ", pathway), color = padj.group),
      size = 2.5, hjust = "inward"
    )
}else{
  
  if(pos_new<=neg_nes){
    p <- p +
      geom_text(
        data = subset(gsea_df, NES > 0),
        aes(x = index, y = 0, label = paste0(pathway, "  "), color = padj.group),
        size = 2.5,
        hjust = "inward"
      ) +
      geom_text(
        data = subset(gsea_df, NES < 0),
        aes(x = index, y = 0, label = paste0("  ", pathway), color = padj.group),
        size = 2.5, hjust = "outward"
      )
  }else if (pos_new > neg_nes){
    p <- p +
      geom_text(
        data = subset(gsea_df, NES > 0),
        aes(x = index, y = 0, label = paste0(pathway, "  "), color = padj.group),
        size = 2.5,
        hjust = "outward"
      ) +
      geom_text(
        data = subset(gsea_df, NES < 0),
        aes(x = index, y = 0, label = paste0("  ", pathway), color = padj.group),
        size =  2.5, hjust = "inward"
      )
  }}

p <- p +
  scale_colour_manual(values = c("black")) +
  labs(x = "", y = "Normalized enrichment score") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
                                                          panel.background = element_blank())  




print(p)
}

```

# WGCNA

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

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

#need to remove control samples (??? don't think I do)
#logCPM <- logCPM[, !colnames(logCPM) %in% c("RH30_P_1", "RH30_P_2", "RH30_P_3", "RH30_P_4", "RH30_P_5", "RH30_P_6", "RH30_P_6")]

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


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

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


# create trait data 
datTraits =read_xlsx("Luciles_trait_data.xlsx")
 datTraits <- as.matrix(datTraits)
row.names(datTraits) = row.names(input_mat)
#check trait data
datTraits


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





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

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

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

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

#Select or set soft threshold power 
#beta = 14

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

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



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


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




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


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


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


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

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

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



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



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



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


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



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



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



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

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





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



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


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



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


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


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

#summary output 

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

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

list_module_sizes[i] <- module_size
}

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



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



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

#find average module size
average_module_size/number_of_modules

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

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

#view geneinfo file
 fix(geneInfo)

 
 
 #getting genes from interesting modules

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

  
 #get gene info for module selected
 gene_info_darkslateblue <- geneInfo[geneInfo$genes %in% module_genes,] 
  write.csv(gene_info_darkslateblue, "gene_info_darkslateblue.csv")

 
 
 genes_darkslateblue <- row.names(gene_info_darkslateblue)
 
#select genes with GS. >0.7 and MM > 0.7 as hub genes
hub_genes <-  subset(gene_info_darkslateblue, GS. >= 0.7 & MM.darkslateblue >= 0.5)
 


```


#WGCNA RH4 split on cell line
```{r}

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

# #need normalised counts in matrix with rows as samples and columns as genes
# #get normalised counts
# rownames(edgeRData) <- paste0("X", seq_len(nrow(edgeRData)))
# edgeRData <- DGEList(counts=counts.keep, genes=rownames(counts.keep), group = treatment)
# edgeRData <- calcNormFactors(edgeRData, method="TMM")
# CPM <- cpm(edgeRData, log = TRUE)
# 
# #remove embryonal cell line RMSYM
# CPM <- CPM[,!grepl("RMSYM", colnames(CPM))]
# #check colnmaes
# colnames(CPM)



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

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

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

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

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

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


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

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

sampleID <- colnames(rawData)
treatment <- c()
for (i in 1:ncol(rawData)) {
  samplename <- sampleID[i]
  group <- if(grepl("RH4ctrl", samplename)) {print("RH4ctrl")
  } else if(grepl("RH4ifo", samplename)) {print("RH4ifo")}
  else if(grepl("RH4vcr", samplename)) {print("RH4vcr")
  } else if(grepl("RMSYMctrl", samplename)) {print("RMSYMctrl")
  } else if(grepl("RMSYMifo", samplename)) {print("RMSYMifo")
  } else {print("RMSYMvcr")}
  treatment[i] <- group
}

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


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


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


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

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


# create trait data 
datTraits =read_xlsx("Luciles_RH4_trait_data.xlsx")
 datTraits <- as.matrix(datTraits)
row.names(datTraits) = row.names(input_mat)
#check trait data
datTraits

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





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

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

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

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

#Select or set soft threshold power 
#beta = 14

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

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


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

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

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


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

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


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

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

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



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



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



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



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



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



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

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



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


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



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


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


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


module_size <- colnames(datME)[moduleColors=="steelblue"]
length(module_size)



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

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

list_module_sizes[i] <- module_size
}

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


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

#find average module size
average_module_size/number_of_modules




# Create the data for the chart 
name   <- module_list
value=list_module_sizes$list_module_sizes 
GFG_Name <- c("lightgreen"  ,   "steelblue"     , "orange"      ,   "sienna3"       , "darkmagenta" ,   "darkolivegreen" ,"darkorange")
X <- data.frame(name,GFG_Name,value) 
 
  
  
# Plot the bar chart  

par(mar=c(11,6,4,1)+.1)
p <- barplot(value,las=2,names.arg=GFG_Name,ylab="Size", ylim = c(0,5000)
        ,col=c("lightgreen"   ,  "steelblue"    ,  "orange"   ,      "sienna3"     ,   "darkmagenta" ,   "darkolivegreen", "darkorange"), 
main="Module sizes of the RH4 acquired resistance data") 
title(xlab="Module name", line=6, cex.lab=1)
jpeg(file="module_size_RH4_luciles.jpeg")



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

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

#view geneinfo file
 fix(geneInfo)
 
geneInfo <-  read.csv("RH4_geneInfo.csv")
 

 #getting genes from interesting modules--------------------------------------------------
 #for orange

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

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

 genes <- row.names(gene_info)
 
#select genes with GS. >0.7 and MM > 0.6 as hub genes
hub_genes <-  subset(gene_info, GS. >0.6 & MM.orange >0.6)
 
write.csv(hub_genes, file = "hub_genes_orange_0.6.csv")


 #getting genes from interesting modules--------------------------------------------------
 #for lightgreen

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

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

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

#getting genes from interesting modules--------------------------------------------------
 #for steelblue

  module_genes <- colnames(datET)[moduleColors=="steelblue"] 
 


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

 genes <- row.names(gene_info)
 
#select genes with GS. >0.7 and MM > 0.6 as hub genes
hub_genes <-  subset(gene_info, GS. >0.6 & MM.steelblue >0.6)
 
write.csv(hub_genes, file = "hub_genes_steelblue_0.6.csv")
```


#GSEA on hub genes
```{r}

#use clusterprofiler to do GSEA on module associated with resistance

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

gene_info <- read.csv("C:/Users/Christina/Documents/RMS RNA seq/Luciles data/RH4 WGCNA/RH4_geneInfo.csv")

options(enrichplot.colours = c("red","blue"))
modules <- c("orange", "lightgreen")


for (i in modules){
  #get hub genes from module
hub_genes <- read.csv(file=paste0("C:/Users/Christina/Documents/RMS RNA seq/Luciles data/RH4 WGCNA/hub_genes_", i,"_0.6.csv") )
#get just genes 
geneList <- hub_genes$genes


#Enrichment analysis on hub genes 

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

GO <- try(barplot(ego2,  font.size = 16 ,showCategory=20, x="Count"))
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)
write.csv(em, file=paste0("Hallmark_hub_genes_", i,".csv"))

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



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

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

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

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

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

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

}


```


#Get gene sig score

```{r}

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

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

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


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

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

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


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



```



# Score
```{r}

#using GSVA to calculate a chemo-resistant score


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

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

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

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


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

#get sample group as extra col in df
sampleID <- row.names(gsva_results)
gsva_results$sample <- c()
for (i in 1:nrow(gsva_results)) {
  samplename <- sampleID[i]
  group <- if(grepl("RH4ctrl", samplename)) {print("RH4 control")
  } else if(grepl("RH4ifo", samplename)) {print("RH4 ifosfamide-resistant")}
  else if(grepl("RH4vcr", samplename)) {print("RH4 vincristine-resistant")
  } else if(grepl("RMSYMctrl", samplename)) {print("RMSYM control")
  } else if(grepl("RMSYMifo", samplename)) {print("RMSYM ifosfamide-resistant")
  } else {print("RMSYM vincristine-resistant")}
  gsva_results$sample[i] <- group
}

#get treatment  as extra col in df
sampleID <- row.names(gsva_results)
gsva_results$treatment <- c()
for (i in 1:nrow(gsva_results)) {
  samplename <- sampleID[i]
  group <- if(grepl("RH4ctrl", samplename)) {print("Control")
  } else if(grepl("RH4ifo", samplename)) {print("Ifosfamide-resistant")}
  else if(grepl("RH4vcr", samplename)) {print("Vincristine-resistant")
  } else if(grepl("RMSYMctrl", samplename)) {print("Control")
  } else if(grepl("RMSYMifo", samplename)) {print("Ifosfamide-resistant")
  } else {print("Vincristine-resistant")}
  gsva_results$treatment[i] <- group
}


#quick plot to look at scores
Rh4_acquired_resistance_score_plot=ggplot(data=gsva_results, mapping=aes(x=sample, y=gsva_results, fill=treatment))+geom_boxplot() +
  stat_summary(fun = "mean", geom = "point", shape = 8,
               size = 2, color = "white") +labs(title="RH4 acquired resistance score: Lucile's data",
        x ="Sample", y = "RH4 acquired resistance score GSVA score")

Rh4_acquired_resistance_score_plot



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

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

shapiro.test(gsva_results$gsva_results)

qqPlot(gsva_results$gsva_results)

#data not normally distributed 

#we do non-parametric anova or kruskal test

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

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


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


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


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

# First do nromality test then a global test

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



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


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


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

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


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






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

#read in GSVA scores for Zoes data
#get hub genes 
gsva_results_zoes_data <- read.csv("C:/Users/Christina/Documents/RMS RNA seq/Zoes data/Zoes_RH4_acquired_score") 

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


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


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


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

#frist normality test 
shapiro.test(gsva_scores$gsva_results)

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

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


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


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



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

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


```



#Overlap DEGs RH4
```{r}
#getting hub genes from overalapping DEGs and WGCNA hub genes 


#get all hub genes
  x <- read.csv("hub_geneslightgreen.csv")
  y <- read.csv("hub_genesorange.csv")
z <- read.csv("hub_genessteelblue.csv")

all_clone_hub_genes <- c(x$x, y$x, z$x)


#get list of DEGs for clones C7 and C13 
#get list of genes DE in at least one condition
files <- list.files( pattern = "RH4ctrl" )



#create empty list
my_list <- list()

for (i in files){ 
  #read file 
  DEGs <- read.csv(i)
  
  #filter for up/down degs
  DEGs <- DEGs[DEGs$Expression != "Unchanged",]
  
  #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 overlap in DEG genes and genes from signif resistant modules from WGCNA
gene_signature <- Reduce(intersect,list(all_DEGs,all_clone_hub_genes))

write.csv(gene_signature, file = "RH4_gene_signature_by_overlap_with_DEGS.csv")



```


#PPI
```{r}

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

"orange"
modules <- c("lightgreen" )
# geneList <- all_clone_hub_genes
# i <- "all_clone_hub_genes"

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

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



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


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

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

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

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

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

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

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

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

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


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

nodes_vis <- left_join(nodes_vis, colour_df, by = "label") #using the left_join() function to merge the nodes_vis and the colour_df, visnetwork will then recognise the colour value
#vis of network
 projection <- visNetwork(nodes_vis, edges_vis, width = "100%") %>%
    visPhysics(stabilization = FALSE, timestep = .1,barnesHut = list(
      centralGravity = 0.5, 
      gravitationalConstant = -30000, 
      springLength = 1000, 
      springConstant = 0.04, 
      avoidOverlap = 0)) %>% 
    visNodes(shape = "dot",
             color = list(
               background = "#0085AF",
               border = "#013848",
               highlight = "#FF8000"
             ),
             scaling = list(min = 20, max = 100),
             shadow = list(enabled = FALSE)
    ) %>%
    visEdges(
      smooth = T,
      shadow = FALSE,
      color = list(highlight = "#C62F4B")
    ) %>%
    visOptions(highlightNearest = list(enabled = TRUE) ) %>% 
    visLayout(randomSeed = 123)
  projection



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

#select genes with score greater than or equal to 6
top_nodes <- as.data.frame(top_nodes)
   table(top_nodes$top_nodes >= 8) #get number for which this is true
  selected_nodes <- subset(top_nodes, top_nodes >= 8) #subset data to select for score >=8
  write.csv(selected_nodes, file=paste0("selected_nodes_", i,".csv")) #save selected nodes
  
  #find overlap between selected nodes and hub genes
  selected_nodes <- tibble::rownames_to_column(selected_nodes, "genes")#chnage rownames to col so all data is retained when subset
  gene_sig_using_PPI <- selected_nodes[selected_nodes$genes %in% hub_genes,]
  #save results
 write.csv(gene_sig_using_PPI, file=paste0("gene_sig_using_PPI_", i,".csv"))
  
  
}
```



#WGCNA RMSYM
```{r}

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

# #need normalised counts in matrix with rows as samples and columns as genes
# #get normalised counts
# rownames(edgeRData) <- paste0("X", seq_len(nrow(edgeRData)))
# edgeRData <- DGEList(counts=counts.keep, genes=rownames(counts.keep), group = treatment)
# edgeRData <- calcNormFactors(edgeRData, method="TMM")
# CPM <- cpm(edgeRData, log = TRUE)
# 
# #remove embryonal cell line RMSYM
# CPM <- CPM[,!grepl("RMSYM", colnames(CPM))]
# #check colnmaes
# colnames(CPM)



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

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

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

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

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

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


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

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

sampleID <- colnames(rawData)
treatment <- c()
for (i in 1:ncol(rawData)) {
  samplename <- sampleID[i]
  group <- if(grepl("RH4ctrl", samplename)) {print("RH4ctrl")
  } else if(grepl("RH4ifo", samplename)) {print("RH4ifo")}
  else if(grepl("RH4vcr", samplename)) {print("RH4vcr")
  } else if(grepl("RMSYMctrl", samplename)) {print("RMSYMctrl")
  } else if(grepl("RMSYMifo", samplename)) {print("RMSYMifo")
  } else {print("RMSYMvcr")}
  treatment[i] <- group
}

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


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


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


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

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


# create trait data 
datTraits =read_xlsx("Luciles_RMSYM_trait_data.xlsx")
 datTraits <- as.matrix(datTraits)
row.names(datTraits) = row.names(input_mat)
#check trait data
datTraits

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





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

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

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

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

#Select or set soft threshold power 
#beta = 14

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

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



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





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

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


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


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

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

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



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



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



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


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



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



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



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

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





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



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


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



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


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


module_size <- colnames(datET)[moduleColors=="brown4"]
length(module_size)
#1249 genes 

#summary output 

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

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

list_module_sizes[i] <- module_size
}

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



# Create the data for the chart 
name   <- module_list
value=list_module_sizes$list_module_sizes 
GFG_Name <- c("plum1"     ,     "skyblue"     ,   "turquoise"   ,   "honeydew1"    ,  "brown2"    ,     "darkorange2"   ,
 "blue2"    ,      "yellow4"      ,  "maroon"    ,     "navajowhite2" ,  "coral1"     , "darkslateblue" ,
 "plum"       ,    "thistle2" ,      "darkgrey"   ,    "green"     ,     "darkolivegreen" ,"brown4"    ,    
 "lightcoral" ,   "grey")
X <- data.frame(name,GFG_Name,value) 
 
  
  
# Plot the bar chart  
par(mar=c(11,6,4,1)+.1)
barplot(value,las=2,names.arg=GFG_Name,ylab="Size", ylim = c(0,3000)
        ,col=c("plum1"     ,     "skyblue"     ,   "turquoise"   ,   "honeydew1"    ,  "brown2"    ,     "darkorange2"   ,
 "blue2"    ,      "yellow4"      ,  "maroon"    ,     "navajowhite2" ,  "coral1"     , "darkslateblue" ,
 "plum"       ,    "thistle2" ,      "darkgrey"   ,    "green"     ,     "darkolivegreen" ,"brown4"    ,    
 "lightcoral" ,   "grey"), 
main="Module sizes of the RMSYM acquired resistance data") 
title(xlab="Module name", line=6, cex.lab=1)













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

#find average module size
average_module_size/number_of_modules

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

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

#view geneinfo file
 fix(geneInfo)

 
 geneInfo <- read.csv("C:/Users/Christina/Documents/RMS RNA seq/Luciles data/RMSYM WGCNA/RMSYM_geneInfo.csv")
 
 #getting genes from interesting modules--------------------------------------------------
 #for darkolivegreen

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

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

 genes <- row.names(gene_info)
 
#select genes with GS. >0.7 and MM > 0.7 as hub genes
hub_genes <-  subset(gene_info, GS. >0.7 & MM.darkolivegreen >0.7)
 
write.csv(hub_genes, file = "hub_genes_darkolivegreen_0.7.csv")


 #getting genes from interesting modules--------------------------------------------------
 #for green

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

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

 genes <- row.names(gene_info)
 
#select genes with GS. >0.7 and MM > 0.7 as hub genes
hub_genes <-  subset(gene_info, GS. >0.7 & MM.green >0.7)
 
write.csv(hub_genes, file = "hub_genes_green_0.7.csv")


 #getting genes from interesting modules--------------------------------------------------
 #for lightcoral

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

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

 genes <- row.names(gene_info)
 
#select genes with GS. >0.7 and MM > 0.7 as hub genes
hub_genes <-  subset(gene_info, GS. >0.7 & MM.lightcoral >0.7)
 
write.csv(hub_genes, file = "hub_genes_lightcoral_0.7.csv")

#getting genes from interesting modules--------------------------------------------------
 #for brown4

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

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

 genes <- row.names(gene_info)
 
#select genes with GS. >0.7 and MM > 0.7 as hub genes
hub_genes <-  subset(gene_info, GS. >0.7 & MM.brown4 >0.7)
 
write.csv(hub_genes, file = "hub_genes_brown4_0.7.csv")

```



#GSEA on hub genes
```{r}

#use clusterprofiler to do GSEA on module associated with resistance

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

#read in geneInfo
gene_info <- read.csv("C:/Users/Christina/Documents/RMS RNA seq/Luciles data/RMSYM WGCNA/RMSYM_geneInfo.csv")


modules <- c("green", "darkolivegreen", "brown4", "lightcoral")

for (i in modules){
  #get hub genes from module
hub_genes <- read.csv(file=paste0("C:/Users/Christina/Documents/RMS RNA seq/Luciles data/RMSYM WGCNA/hub_genes_", i,"_0.6.csv") )
#get just genes 
geneList <- hub_genes$genes


#Enrichment analysis on hub genes 

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

GO <- try(barplot(ego2,  font.size = 16 ,showCategory=20, x="Count"))
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)
write.csv(em, file=paste0("Hallmark_hub_genes_", i,".csv"))

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



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

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

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

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

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

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

}


```


#PPI
```{r}

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

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

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

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



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


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

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

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

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

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

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

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

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

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


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

nodes_vis <- left_join(nodes_vis, colour_df, by = "label") #using the left_join() function to merge the nodes_vis and the colour_df, visnetwork will then recognise the colour value
#vis of network
 projection <- visNetwork(nodes_vis, edges_vis, width = "100%") %>%
    visPhysics(stabilization = FALSE, timestep = .1,barnesHut = list(
      centralGravity = 0.5, 
      gravitationalConstant = -30000, 
      springLength = 1000, 
      springConstant = 0.04, 
      avoidOverlap = 0)) %>% 
    visNodes(shape = "dot",
             color = list(
               background = "#0085AF",
               border = "#013848",
               highlight = "#FF8000"
             ),
             scaling = list(min = 20, max = 100),
             shadow = list(enabled = FALSE)
    ) %>%
    visEdges(
      smooth = T,
      shadow = FALSE,
      color = list(highlight = "#C62F4B")
    ) %>%
    visOptions(highlightNearest = list(enabled = TRUE) ) %>% 
    visLayout(randomSeed = 123)
  projection



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

#select genes with score greater than or equal to 6
top_nodes <- as.data.frame(top_nodes)
   table(top_nodes$top_nodes >= 8) #get number for which this is true
  selected_nodes <- subset(top_nodes, top_nodes >= 8) #subset data to select for score >=8
  write.csv(selected_nodes, file=paste0("selected_nodes_", i,".csv")) #save selected nodes
  
  #find overlap between selected nodes and hub genes
  selected_nodes <- tibble::rownames_to_column(selected_nodes, "genes")#chnage rownames to col so all data is retained when subset
  gene_sig_using_PPI <- selected_nodes[selected_nodes$genes %in% hub_genes,]
  #save results
 write.csv(gene_sig_using_PPI, file=paste0("gene_sig_using_PPI_", i,".csv"))
  
  
}
```




#Get gene sig 

```{r}

#for darkolivegreen------------------------------------------------------
 #read in hub proteins 

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

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


 gene_sig_RMSYM <- module_genes[module_genes$genes %in% hub_proteins,]
  gene_sig_RMSYM_darkolive_green <- gene_sig_RMSYM$X
  
  
write.csv(gene_sig_RMSYM_darkolive_green, "gene_sig_RMSYM_darkolive_green.csv")
  
#for green-----------------------------------------------------------------
 #read in hub proteins 

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

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


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



```

#score ERMS
```{r}

#using GSVA to calculate a chemo-resistant score


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

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

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

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


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

#get sample group as extra col in df
sampleID <- row.names(gsva_results)
gsva_results$sample <- c()
for (i in 1:nrow(gsva_results)) {
  samplename <- sampleID[i]
  group <- if(grepl("RH4ctrl", samplename)) {print("RH4 control")
  } else if(grepl("RH4ifo", samplename)) {print("RH4 ifosfamide-resistant")}
  else if(grepl("RH4vcr", samplename)) {print("RH4 vincristine-resistant")
  } else if(grepl("RMSYMctrl", samplename)) {print("RMSYM control")
  } else if(grepl("RMSYMifo", samplename)) {print("RMSYM ifosfamide-resistant")
  } else {print("RMSYM vincristine-resistant")}
  gsva_results$sample[i] <- group
}

#get treatment  as extra col in df
sampleID <- row.names(gsva_results)
gsva_results$treatment <- c()
for (i in 1:nrow(gsva_results)) {
  samplename <- sampleID[i]
  group <- if(grepl("RH4ctrl", samplename)) {print("Control")
  } else if(grepl("RH4ifo", samplename)) {print("Ifosfamide-resistant")}
  else if(grepl("RH4vcr", samplename)) {print("Vincristine-resistant")
  } else if(grepl("RMSYMctrl", samplename)) {print("Control")
  } else if(grepl("RMSYMifo", samplename)) {print("Ifosfamide-resistant")
  } else {print("Vincristine-resistant")}
  gsva_results$treatment[i] <- group
}



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

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

shapiro.test(gsva_results$gsva_results)

qqPlot(gsva_results$gsva_results)

#data not normally distributed 

#we do non-parametric anova or kruskal test

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

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


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


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


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

# First do nromality test then a global test

compare_means(gsva_results ~ sample,  data = gsva_results, method = "anova")



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


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


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

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


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






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

#read in GSVA scores for Zoes data
#get hub genes 
gsva_results_zoes_data <- read.csv("C:/Users/Christina/Documents/RMS RNA seq/Luciles data/RMSYM WGCNA/Zoes_RMSYM_acquired_score.csv") 

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


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


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


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

#frist normality test 
shapiro.test(gsva_scores$gsva_results)

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

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


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


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



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


write.csv(gsva_scores, "gsva_scores_ERMS_acquired_all_samples.csv")


```






#Overlap DEGs RMSYM
```{r}
#getting hub genes from overalapping DEGs and WGCNA hub genes 


#get all hub genes
  x <- read.csv("hub_genesbrown4.csv")
  y <- read.csv("hub_genesdarkolivegreen.csv")
z <- read.csv("hub_genesdarkorange2.csv")
a <- read.csv("hub_genesgreen.csv")
b <- read.csv("hub_geneslightcoral.csv")

all_clone_hub_genes <- c(x$x, y$x, z$x, a$x, b$x)


#get list of DEGs for clones C7 and C13 
#get list of genes DE in at least one condition
files <- list.files( pattern = "RMSYMctrl" )



#create empty list
my_list <- list()

for (i in files){ 
  #read file 
  DEGs <- read.csv(i)
  
  #filter for up/down degs
  DEGs <- DEGs[DEGs$Expression != "Unchanged",]
  
  #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 overlap in DEG genes and genes from signif resistant modules from WGCNA
gene_signature <- Reduce(intersect,list(all_DEGs,all_clone_hub_genes))

write.csv(gene_signature, file = "RMSYM_gene_signature_by_overlap_with_DEGS.csv")



```


#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



#remove RMSYM samples
CPM <- CPM[,!grepl("RMSYM", colnames(CPM))]

CPM <- CPM[,!grepl("RH4", colnames(CPM))]

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

  #get sample group as extra col in df
sampleID <- row.names(metadata)
metadata$Treatment_condition <- c()
for (i in 1:nrow(metadata)) {
  samplename <- sampleID[i]
  Treatment_condition <- if(grepl("RH4ctrl", samplename)) {print("RH4ctrl")
  } else if(grepl("RH4ifo", samplename)) {print("RH4ifo")}
  else if(grepl("RH4vcr", samplename)) {print("RH4vcr")
  } else if(grepl("RMSYMctrl", samplename)) {print("RMSYMctrl")
  } else if(grepl("RMSYMifo", samplename)) {print("RMSYMifo")
  } else {print("RMSYMvcr")}
  metadata$Treatment_condition[i] <- Treatment_condition
}


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


  screeplot(p, axisLabSize = 18, titleLabSize = 22)
 
 
   biplot(p, showLoadings = TRUE, lab = NULL)
   
    plotloadings(p, labSize = 3)
 
    
    ########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))
 
 
  ## biplot
   biplot(p,
    lab = paste0(p$metadata$Group),
    colby = 'Group',
    hline = 0, vline = 0,
    legendPosition = 'right')
  
  biplot(p,
          lab = NULL,
    colby = 'Treatment_condition',
    hline = 0, vline = 0,
    legendPosition = 'right', 
    
    ellipse = TRUE,
      ellipseType = 't',
      ellipseLevel = 0.95,
      ellipseFill = TRUE,
      ellipseAlpha = 1/4,
      ellipseLineSize = 1.0,
    title = 'PCA bi-plot for R acquired resistance data',
    subtitle = 'PC1 versus PC2',)
  
  
   
     biplot(p,
               showLoadings = TRUE,
      lengthLoadingsArrowsFactor = 1.5,
      sizeLoadingsNames = 5,
      colLoadingsNames = 'red4',

    lab = NULL,
    colby = 'Group',
    hline = 0, vline = 0,
    legendPosition = 'right', 
    ellipse = TRUE,
      ellipseType = 't',
      ellipseLevel = 0.95,
      ellipseFill = TRUE,
      ellipseAlpha = 1/4,
      ellipseLineSize = 1.0,
    title = 'PCA bi-plot',
      subtitle = 'PC1 versus PC2',)
  
 
 
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)


  pairsplot(p,
      components = getComponents(p, c(1:3)),
      triangle = TRUE, trianglelabSize = 12,
      hline = 0, vline = 0,
      pointSize = 1.5,
      gridlines.major = FALSE, gridlines.minor = FALSE,
      colby = 'Group',
      legendPosition = "none",
      title = 'Pairs plot', plotaxes = FALSE,
      margingaps = unit(c(-0.01, -0.01, -0.01, -0.01), 'cm'))
```




#score intrinsic

```{r}
#using GSVA to calculate a chemo-resistant score

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

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

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

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


#get sample group as extra col in df
sampleID <- row.names(gsva_results)
gsva_results$sample <- c()
for (i in 1:nrow(gsva_results)) {
  samplename <- sampleID[i]
  group <- if(grepl("RH4ctrl", samplename)) {print("RH4 control")
  } else if(grepl("RH4ifo", samplename)) {print("RH4 ifosfamide-resistant")}
  else if(grepl("RH4vcr", samplename)) {print("RH4 vincristine-resistant")
  } else if(grepl("RMSYMctrl", samplename)) {print("RMSYM control")
  } else if(grepl("RMSYMifo", samplename)) {print("RMSYM ifosfamide-resistant")
  } else {print("RMSYM vincristine-resistant")}
  gsva_results$sample[i] <- group
}

#get treatment  as extra col in df
sampleID <- row.names(gsva_results)
gsva_results$treatment <- c()
for (i in 1:nrow(gsva_results)) {
  samplename <- sampleID[i]
  group <- if(grepl("RH4ctrl", samplename)) {print("Control")
  } else if(grepl("RH4ifo", samplename)) {print("Ifosfamide-resistant")}
  else if(grepl("RH4vcr", samplename)) {print("Vincristine-resistant")
  } else if(grepl("RMSYMctrl", samplename)) {print("Control")
  } else if(grepl("RMSYMifo", samplename)) {print("Ifosfamide-resistant")
  } else {print("Vincristine-resistant")}
  gsva_results$treatment[i] <- group
}


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


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


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



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


intrinsic_resistance_score_plot

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


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



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

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

shapiro.test(gsva_results$gsva_results)

qqPlot(gsva_results$gsva_results)

#data not normally distributed 

#we do non-parametric anova or kruskal test

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

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


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


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

```


#info
```{r}
#read in gene sig
gene_sig <- read_xlsx("intrinsic_resistance_signature.xlsx")
genes <- gene_sig$Gene

#read in geneinfo
gene_info <- read.csv("geneInfo_clones.csv")

#read in degs info
C7 <- read.csv("unfiltered_DEGs_RH30_P_vs_RH_C7.csv")


#subset data for genes in col
de_info <- C7[C7$X %in% genes,]
de_info <- de_info[c("X", "logFC","FDR")]

write.csv(de_info, "C7_gene_sig_info.csv")

#read in degs info
C13 <- read.csv("unfiltered_DEGs_RH30_P_vs_RH_C13.csv")


#subset data for genes in col
de_info <- C13[C13$X %in% genes,]
de_info <- de_info[c("X", "logFC","FDR")]

write.csv(de_info, "C13_gene_sig_info.csv")


#subset data for genes in col
info <- gene_info[gene_info$genes %in% genes,]
info <- info[c("genes", "GS.","MM.darkgrey","MM.darkgreen")]

write.csv(info, "gene_sig_genes_info.csv")




```








#MDR genes
```{r}


MDR_genes <- read_delim("MDR_genes.txt")
MDR_genes <- MDR_genes$Gene


```

