---
title: "Microarray analysis Triche"
author: "Christina Putnam"
date: "2023-12-13"
output: html_document
editor_options: 
  chunk_output_type: console
---

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

# Load librarys
```{r}
library(GEOquery)
library("caret")
library(oligo)
library(fgsea)
library("randomcoloR")
library(tidyr)
library(tidyverse)
library(hgu133plus2.db)
library(dplyr)
library(survival)
library(survminer)
library(hgu133a.db)
library(affycoretools)
library("readxl")
library(labelled)
library(tidyverse)
library("fuzzyjoin")
library(GSVA)
library(ggpubr)
library(jetset)
library(viridis)
library("affyio")
library(PCAtools)
library(dplyr)
library(oligoClasses)
```

# Import data from GEO
```{r}

#set wd
setwd("~/RMS microarray/Triche")

#When getGEO() is run with the named argument GSEMatrix=FALSE, it does not parse the Series Matrix files but does parse the SOFT files. The advantage of this method is that the SOFT files may contain more information than the series matrix files, but a disadvantage is that it is much slower and the resulting object is not a list of ExpressionSets. However, parsing the Series Matrix can also introduce errors; parsing the SOFT file is less error-prone.

gse92689 <- getGEO('GSE92689')
length(gse92689)
gse92689 <- gse92689[[1]]

```

#get phenodata and add GEO ID 
```{r}
#metadata---------------------------------------------------------------
varLabels(gse92689)
gse92689$supplementary_file
metadata <- pData(gse92689)

#make col with cell file info
metadata['cel_file'] <- str_split(metadata$supplementary_file,"/") %>% map_chr(tail,1)

#celdata <- read.celfiles(paste0("C:/Users/Christina/Documents/RMS microarray/Triche/GSE92689_RAW/",metadata$cel_file),phenoData=phenoData(gse92689))

#select only cel file column 
metadata <- metadata %>%
  dplyr::select(cel_file)

#get start of cel file name
sample_id <- sapply(strsplit(metadata$cel_file, "_"), "[", 2)
sample_id <- sapply(strsplit(sample_id, "-"), "[", 1)

#add col containing cell file name
metadata$start_CELL_File_name <- sample_id
#add col containing geo ids
metadata$geo_id <- rownames(metadata)

#select only cel file and geo columns
metadata <- metadata %>%
  dplyr::select(start_CELL_File_name, geo_id)

#pd---------------------------------------------------------------------------

# getting phenotype data
pd <- read_xlsx("~/RMS microarray/Triche/Clinical_annotation.xlsx")

#rename cols
colnames(pd)[which(names(pd) == "CELL_File name")]  <- "CELL_File_name"

#get start of cel file name
sample_id <- sapply(strsplit(pd$CELL_File_name, "-"), "[", 1)

#add to new col in df
pd$start_CELL_File_name <- sample_id

#merge dfs based on cel file---------------------------------------------------
#merge dfs
#using regrex to match near patterns instead of exact matches
clinical_annotation_with_geoID <- regex_inner_join(pd, metadata, by = c("start_CELL_File_name")) %>% as.data.frame()
#test <- inner_join(metadata, pd, by = "start_CELL_File_name")

#change rownames to GEO_id
rownames(clinical_annotation_with_geoID) <- clinical_annotation_with_geoID$geo_id

#write to xlsx
write_excel_csv(clinical_annotation_with_geoID, "clinical_annotation_1.csv")

#check how many didn't match
dim(metadata %>% filter(!geo_id %in% clinical_annotation_with_geoID$geo_id))
#13 didn't match (present in metadata but not fully annotated)

#lets look at their metadata
not_fully_annotated <- rownames(metadata %>% dplyr::filter(!geo_id %in% clinical_annotation_with_geoID$geo_id))


#reload metadata
metadata <- pData(gse92689)
#subset for the ones that haven't been fully annotated
metadata <- metadata[ rownames(metadata) %in% not_fully_annotated, ]

# test <-  clinical_annotation_with_geoID[ rownames(clinical_annotation_with_geoID) %in% not_fully_annotated, ]
# #one cell line and 12 primary tumours

```


# Get raw data

```{r}
# Getting the raw data (CEL files for Affymetrix data) is distinct from getting the processed ata. The CEL files can be pretty large, and sometimes the download can fail, even with a good connection.

#Nonetheless, if you want to do it, getGEOSuppFiles() will download all the supplementary files for the GEO accession. Note that it doesn’t process or even parse the files, since there are many different types of supplementary files on GEO and R doesn’t know the format ahead of time. It does, however, return the file paths.

#celdata <- read.celfiles(list.celfiles("~/RMS microarray/Triche/GSE92689_RAW", pattern = ".CEL.gz",listGzipped=TRUE, full.names = TRUE))
#exprs(celdata)

#reload metadata
metadata <- pData(gse92689)
#make col with cell file info
metadata['cel_file'] <- str_split(metadata$supplementary_file,"/") %>% map_chr(tail,1)

celdata <- read.celfiles(paste0("C:/Users/Christina/Documents/RMS microarray/Triche/GSE92689_RAW/",metadata$cel_file),phenoData=phenoData(gse92689))

exprs(celdata)

# Pseudo-image plots are used to assess the spatial distribution of the data on the chips. Due to the
# magnitude of the readings, pseudo-images using data on the original scale often mask spatial features
# that may be present on the arrays. This is why we recommend the use of the default log2
# -scale of the
# image method. One useful alternative for the log2
# -scale pseudo-image is the use of the ranks of the
# observations. This can be achieved by setting the transfo argument on the image method
image(celdata[,1])
  
  #plots of raw data ---------------------------------------------------------

#log 2 transformed hsitogram of raw data
tiff(file="histogram_raw_data_triche.tiff
     ", width=6, height=4, units="in", res=300)
hist(celdata, which='pm', ylab = "log-2 intensity", 
        cex.axis = 1, xlab = "sample", main = "Histogram of log-2 intensities for the raw data (Triche)")
dev.off()

tiff(file="boxplot_raw_data_triche.tiff", width=19, height=6, units="in", res=300)
#boxplot pf log 2 tranformed raw data. for all probes, to select PM or MM add arguement
boxplot(celdata, which='pm', ylab = "log-2 intensity",   las = 2, 
        cex.axis = 0.4, xlab = "sample", main = "Boxplot of log-2 intensities for the raw data (Triche)")
dev.off()

 #arrayQualityMetrics(expressionset = celdata, outdir = "Report_for_celdata", force = TRUE)

```

#background correct
```{r}
oligo_BC <- backgroundCorrect(celdata)

tiff(file="boxplot_BC_triche.jpeg", width=19, height=6, units="in", res=300)
boxplot(oligo_BC, which='pm', ylab = "log-2 intensity",   las = 2, 
        cex.axis = 0.4, xlab = "sample", main = "Boxplot of log-2 intensities after background correction")
dev.off()


#histogram
tiff(file="histogram_BC_triche.jpeg", width=6, height=4, units="in", res=300)
hist(oligo_BC, which='pm', ylab = "log-2 intensity", 
        cex.axis = 1, xlab = "sample", main = "Histogram of log-2 intensities after background correction")
dev.off()

```


#Normalisation
```{r}

oligo_normalised <- normalize(celdata,method='quantile',which='pm')

#boxplot
tiff(file="boxplot_normalised_triche.jpeg", width=19, height=6, units="in", res=300)
boxplot(oligo_normalised, which='pm', ylab = "log-2 intensity",   las = 2, 
        cex.axis = 0.4, xlab = "sample", main = "Boxplot of log-2 intensities after quantile normalisation (Triche)")
dev.off()

#histogram
tiff(file="histogram_normalised_triche.jpeg", width=6, height=4, units="in", res=300)
hist(oligo_normalised, which='pm', ylab = "log-2 intensity", 
        cex.axis = 1, xlab = "sample", main = "Histogram of log-2 intensities after quantile normalisation (Triche)")
dev.off()

```

#Summarisation
```{r}
#from probe level to feature level 

oligo_summarised <- rma(celdata,background=TRUE,normalize=TRUE) #this is doing the normalisation and background correction as well all in one step instead of above where i have done it individually (to make plots)

#boxplot
tiff(file="boxplot_normalised_and_summarised_triche.jpeg", width=19, height=6, units="in", res=300)
boxplot(oligo_summarised, which='pm', ylab = "log-2 intensity",   las = 2, 
        cex.axis = 0.4, xlab = "sample", main = "Boxplot of log-2 intensities after normalisation and summarisation to the gene level")
dev.off()

#histogram
tiff(file="histogram_normalised_and_summarised_triche.jpeg", width=6, height=4, units="in", res=300)
hist(oligo_summarised, ylab = "log-2 intensity", 
        cex.axis = 1, xlab = "sample", main = "Histogram of log-2 intensities after normalisation and summarisation")
dev.off()

```


#Performing RMA
```{r}
# the following one line of code performs RMA normalisation, and returns an ExpressionSet object containing the background corrected, normalised, and summarised expression data

# You can also try, running RMA without background correction, or without quantile normalisation, and plot the densities.

gse92689_eset <- oligo::rma(celdata)

#plot histogram of background corrected data
hist(gse92689_eset)

#boxplot of background corrected data
boxplot(gse92689_eset)

#. When preprocessing
# the data, oligo stores the summaries in a matrix called exprs, present in the assayData data slot
# of the ExpressionSet object.

exprs(gse92689_eset)[10001:10010, 1:2]
boxplot(exprs(gse92689_eset)[10001:10010, 1:2])

#get normalised expression
expr = exprs(gse92689_eset)

#QC normalised data
# arrayQualityMetrics(expressionset = gse92689_eset, outdir = "Report_for_normalised_GSE92689", force = TRUE)

```


#Annotation
```{r}
#We need to map these IDs to gene symbols, which can then be further analyzed downstream.To do so, we will use the annotation package hgu133plus2.db, where hgu133plus2 is the array name. Intuitively, different arrays will have different annotation packages, but they will all end with .db.

#What we see are a bunch of maps from the hgu133plus2 probeset identifiers (which we have) to other identifiers (which we don’t). We need to provide the probe identifiers of interest and then retrieve the other identifiers of interest.
ls('package:hgu133plus2.db')

#If we want to extract the symbols, gene identifiers, and gene names, it’s as simple as using the select() function from AnnotationDbi with the probesets as our identifiers:
expression_data <- exprs(gse92689_eset)
#probes <- row.names(expression_data)

#annotate data
    gset <- annotateEset(gse92689_eset, x = "hgu133a.db" )
    gse92689_eset <- annotateEset(gse92689_eset, x = "hgu133a.db" )
    #get symbols
symbol = gset@featureData@data[["SYMBOL"]]
# remove NA
symbol <- symbol[!is.na(symbol)]


top_probes_from_scores <- function(chip, symbol) {        
  require(jetset)
  require(dplyr)
  top_probes <- function(data) { 
    require(jetset)
    require(dplyr)
    data$probeID = rownames(data)
    t <- as.data.frame(matrix(nrow = length(unique(data$symbol)), ncol = 1 ))
    for (i in 1:nrow(t)) { 
      symbs <- unique(data[ ,"symbol"])
      sym <- symbs[i]
      dat <- data %>% filter(., symbol %in% sym) %>% slice_max(., order_by = overall, n = 1, with_ties = FALSE) %>% pull(probeID)
      t[i , 1] <- dat
    }
    return(t)
  }
  probescores <- jscores(chip = chip, symbol = symbol) 
  t1 <- top_probes(data = probescores)
  colnames(t1) <- "ProbeIDs"
  return(t1)
  ## performing master function
gset <- annotateEset(gset, x = chip.db )
symbol = gset@featureData@data[["SYMBOL"]]
chip = sapply(stringi::stri_split_fixed(chip.db,"."), "[[", 1)
top <- top_probes_from_scores(chip = chip, symbol = symbol )
return(top)
 
}


#run function
top_probes <- top_probes_from_scores("hgu133a", symbol)


#filter gset to select rows with top probes 
 gset_genes <- gset[rownames(gset) %in% top_probes$ProbeIDs,]
 
#get gene symbols for probes
   gset_genes <- annotateEset(gset_genes, x = "hgu133a.db" )
   symbol = gset_genes@featureData@data[["SYMBOL"]]
 
   #replace probes as symbols
 expression_data <- exprs(gset_genes)
 rownames(expression_data) <- symbol
 
```


#Remove cell line/subtypes samples 
```{r}
pd <- clinical_annotation_with_geoID

#rename columns
colnames(pd)[which(names(pd) == "Histology...13")] <- "Histology"
colnames(pd)[which(names(pd) == "Sex...6")] <- "Sex"
colnames(pd)[which(names(pd) == "Age...7")] <- "Age"
colnames(pd)[which(names(pd) == "Invas#")] <- "Invasion"
colnames(pd)[which(names(pd) == "Sample Name")]  <- "sampleID"
colnames(pd)[which(names(pd) == "Cell Type")]  <- "Cell_Type"

#remove cell line data from expression data
 #find samples that are cell lines
 cell_line_samples <- subset(pd, Cell_Type == "tumor cell line") %>% rownames()
expression_data <- expression_data[, !colnames(expression_data) %in% cell_line_samples]

 #remove cell line samples from phenotype
pd <- pd[!rownames(pd) %in% cell_line_samples,]

#check dim
dim(expression_data)
dim(pd)


#remove non-alv/embryonal samples from expression data
 #find samples that are alv/embryonal/mixed
 alv_and_embryonal_samples <- subset(pd, Histology == "Embryonal" |  Histology == "Alveolar"|  Histology == "Mixed Alv/Emb") %>% rownames()
expression_data <- expression_data[, colnames(expression_data) %in% alv_and_embryonal_samples]

#select for only alveolar and embryonal samples in phenotype data
alveolar_and_embryonal_samples <- subset(pd, Histology == "Embryonal" |  Histology == "Alveolar"|  Histology == "Mixed Alv/Emb") %>% rownames()
pd <- pd[ rownames(pd) %in% alveolar_and_embryonal_samples, ]
dim(pd) #125 primary tumour samples that are alv/embry/mixed 
dim(expression_data)


#clean phenotype data
pd[] <- lapply(pd, function(x) gsub("NA", NA, x))

# #NAs
# pd$`diseaselocation:ch1` <- gsub("\\?", NA, pd$`diseaselocation:ch1`)
# pd$`diseaselocation:ch1` <- gsub("--", NA, pd$`diseaselocation:ch1`)
# pd$`diseasestage:ch1` <- gsub("--", NA, pd$`diseasestage:ch1`)
# pd$`diseasestate:ch1` <- gsub("--", NA, pd$`diseasestate:ch1`)
# pd$`Sex:ch1` <- gsub("UNKNOWN", NA, pd$`Sex:ch1`)
# pd$`tumorgrading:ch1` <- gsub("--", NA, pd$`tumorgrading:ch1`)

#change to numeric
pd$Age <- as.numeric(pd$Age)

#change fail to event
pd$Fail <- gsub("\\.", "No event", pd$Fail )
pd$Fail <- gsub("Fail", "Event", pd$Fail )

#change translocation NA to neg
pd$Translocation[is.na(pd$Translocation )] <- "NEG"

```


#Plot clinical charcteristics
```{r}

setwd("~/RMS microarray/Triche/Clinical info")
pd

#number of samples
dim(pd) #125

#age--------------------------------------------------------------

df <- data.frame(
value = as.factor(pd$Age))


  #find length of variables
  n = length(unique(df$value))
#set colours
  #myCols <- setNames(distinctColorPalette(n), unique(df$value))
  
   myCols <- viridis(n, option = "G", direction = -1, end = 0.9)

  p <- ggplot(df, aes(  fill=value, x=value)) + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("Age") + theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18), title = element_text(size=18))+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+xlab("Age") + ylab("Number of samples") + guides(fill="none")
 ggsave(p, file=paste0("age.png"), width = 14, height = 9, units = "cm")

#disease stage-----------------------------------------------------------------
 
 pd$`Disease Stage`

df <- data.frame(
value = as.factor(pd$`Disease Stage`))

  #find length of variables
  n = length(unique(df$value))
#set colours
   myCols <- viridis(n, option = "G", direction = -1, end = 0.9)

  p <- ggplot(df, aes(  fill=value, x=value)) + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("Disease stage",) + theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18) ,title = element_text(size=18))+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+xlab("Disease stage") + ylab("Number of samples") + guides(fill="none")
 ggsave(p, file=paste0("stage.png"), width = 9, height = 10, units = "cm")

 #survival-----------------------------------------------------------------
 
 pd$Survival_rounded <- round(as.numeric(pd$Survival))

df <- data.frame(
value = as.factor(pd$Survival_rounded))

  #find length of variables
  n = length(unique(df$value))
#set colours
 myCols <- viridis(n, option = "G", direction = -1, end = 0.9)

  p <- ggplot(df, aes(  fill=value, x=value)) + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("Survival",) + theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18),title = element_text(size=18))+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+xlab("Survival (years)") + ylab("Number of samples") + guides(fill="none")
 ggsave(p, file=paste0("Survival.png"), width = 14, height = 10, units = "cm")


  #FFS-----------------------------------------------------------------
 
 pd$FFS_rounded <- round(as.numeric(pd$FFS))

df <- data.frame(
value = as.factor(pd$FFS_rounded))

  #find length of variables
  n = length(unique(df$value))
#set colours
 myCols <- viridis(n, option = "G", direction = -1, end = 0.9)

  p <- ggplot(df, aes(  fill=value, x=value)) + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("FFS",) + theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18),title = element_text(size=18))+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+xlab("FFS (years)") + ylab("Number of samples") + guides(fill="none")
 ggsave(p, file=paste0("FFS.png"), width = 14, height = 10, units = "cm")

 

 #histology-------------------------------------------------------------------------------

 pd$Histology

df <- data.frame(
value = pd$Histology)

  #find length of variables
  n = length(unique(df$value))
#set colours
 myCols <- viridis(n, option = "G", direction = -1, end = 0.9)

  p <- ggplot(df, aes(  fill=value, x=value)) + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("Histology",) + theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18),title = element_text(size=18))+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+xlab("Histology") + ylab("Number of samples") + guides(fill="none")
 ggsave(p, file=paste0("histology.png"), width = 8, height = 12, units = "cm")

 #translocation-------------------------------------------------------------------------------

 pd$Translocation

df <- data.frame(
value = pd$Translocation)

  #find length of variables
  n = length(unique(df$value))
#set colours
 myCols <- viridis(n, option = "G", direction = -1, end = 0.9)

  p <- ggplot(df, aes(  fill=value, x=value)) + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("Translocation",) + theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18),title = element_text(size=18))+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+xlab("Translocation") + ylab("Number of samples") + guides(fill="none")
 ggsave(p, file=paste0("Translocation.png"), width = 8, height = 12, units = "cm")
 
 #sex------------------------------------------------------------------------
 pd$Sex

df <- data.frame(
value = pd$Sex)

  #find length of variables
  n = length(unique(df$value))
#set colours
 myCols <- viridis(n, option = "G", direction = -1, end = 0.9)

  p <- ggplot(df, aes(  fill=value, x=value)) + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("Sex",) + theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18),title = element_text(size=18))+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+xlab("Sex") + ylab("Number of samples") + guides(fill="none")
 ggsave(p, file=paste0("sex.png"), width = 8, height = 12, units = "cm")
 
 #tumour grade
 
  pd$`Disease Grade`

df <- data.frame(
value = pd$`Disease Grade`)

  #find length of variables
  n = length(unique(df$value))
#set colours
 myCols <- viridis(n, option = "G", direction = -1, end = 0.9)

  p <- ggplot(df, aes(  fill=value, x=value)) + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("Tumour grade",) + theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18),title = element_text(size=18))+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+xlab("Tumour grade") + ylab("Number of samples") + guides(fill="none")
 ggsave(p, file=paste0("grade.png"), width = 10, height = 10, units = "cm")
 
 
  #location---------------------------------------
 
  pd$`Tumor Location`

df <- data.frame(
value = pd$`Tumor Location`)

  #find length of variables
  n = length(unique(df$value))
#set colours
 myCols <- viridis(n, option = "G", direction = -1, end = 0.9)

  p <- ggplot(df, aes(  fill=value, x=value)) + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("Tumour location",) + theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18),title = element_text(size=18))+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+xlab("Tumour location") + ylab("Number of samples") + guides(fill="none")
 ggsave(p, file=paste0("location.png"), width = 14, height = 13, units = "cm")
 
 
   #risk group---------------------------------------
 
  pd$RG

df <- data.frame(
value = pd$RG)

  #find length of variables
  n = length(unique(df$value))
#set colours
 myCols <- viridis(n, option = "G", direction = -1, end = 0.9)

  p <- ggplot(df, aes(  fill=value, x=value)) + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("Risk group",) + theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18),title = element_text(size=18))+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+xlab("Risk group") + ylab("Number of samples") + guides(fill="none")
 ggsave(p, file=paste0("Riskgroup.png"), width = 8, height = 14, units = "cm")
 
  
   #size---------------------------------------
 
  pd$Size

df <- data.frame(
value = pd$Size)

  #find length of variables
  n = length(unique(df$value))
#set colours
 myCols <- viridis(n, option = "G", direction = -1, end = 0.9)

  p <- ggplot(df, aes(  fill=value, x=value)) + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("Tumour size",) + theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18),title = element_text(size=18))+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+xlab("Tumour size") + ylab("Number of samples") + guides(fill="none")
 ggsave(p, file=paste0("tumoursize.png"), width = 7, height = 12, units = "cm")
 
    #node---------------------------------------
 
  pd$Node

df <- data.frame(
value = pd$Node)

  #find length of variables
  n = length(unique(df$value))
#set colours
 myCols <- viridis(n, option = "G", direction = -1, end = 0.9)

  p <- ggplot(df, aes(  fill=value, x=value)) + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("Node status",) + theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18),title = element_text(size=18))+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+xlab("Node status") + ylab("Number of samples") + guides(fill="none")
 ggsave(p, file=paste0("node.png"), width = 8, height = 12, units = "cm")
 
  #invasion---------------------------------------
 
  pd$Invasion

df <- data.frame(
value = pd$Invasion)

  #find length of variables
  n = length(unique(df$value))
#set colours
 myCols <- viridis(n, option = "G", direction = -1, end = 0.9)

  p <- ggplot(df, aes(  fill=value, x=value)) + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("Invasion",) + theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18),title = element_text(size=18))+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+xlab("Invasion") + ylab("Number of samples") + guides(fill="none")
 ggsave(p, file=paste0("invasion.png"), width = 8, height = 12, units = "cm")
 
 
 
  #Fail---------------------------------------
 
  pd$Fail


df <- data.frame(
value = pd$Fail)

  #find length of variables
  n = length(unique(df$value))
#set colours
 myCols <- viridis(n, option = "G", direction = -1, end = 0.9)

  p <- ggplot(df, aes(  fill=value, x=value)) + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("Event",) + theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18),title = element_text(size=18))+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+xlab("Event") + ylab("Number of samples") + guides(fill="none")
 ggsave(p, file=paste0("event.png"), width = 8, height = 12, units = "cm")
```



#PCA
```{r}

setwd("~/RMS microarray/Triche/PCA")

#get normalised counts and remove dups/NAs
#remove na
expression_data <- subset(expression_data, !is.na(rownames(expression_data)))
expression_data <- expression_data[!duplicated(rownames(expression_data)), ]
                                                           
col_order <- rownames(pd)
expression_data <- expression_data[, col_order]

#run pca
p <- pca(expression_data, metadata = pd, removeVar = 0.1)


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

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

    geom_label(aes(x = horn$n + 0, y = 50,
      label = 'Horn\'s', vjust = -2, size = 8)) +
    geom_label(aes(x = elbow + 1, y = 50,
      label = 'Elbow method', vjust = -1, size = 8))
  
  
#biplots------------------------------------------------------------------
 
make_pca_biplot <- function(clinical_variable, readable_variable){
  
  #find length of variables
  n = length(unique(p$metadata[[paste0(clinical_variable)]]))
 p$metadata[[paste0(clinical_variable)]] <-  as.factor(p$metadata[[paste0(clinical_variable)]])
#set colours
  #myCols <- setNames(distinctColorPalette(n), unique(p$metadata[[i]]))
  myCols <- viridis(n, option = "G", direction = -1, end = 0.9)
  legend_title <- readable_variable
  #make filename title
   filename <- paste0("PCA_biplot_",readable_variable,".png")
   #label <- subset(p$metadata, `histology:ch1` == "Mixed Alv/Emb") %>% rownames() #label only mixed samples
 g <- biplot(p,
         lab = NULL,
         #max.overlaps = 200,
    colby = paste0(clinical_variable),
    colkey = myCols,
    hline = 0, vline = 0,
    legendPosition = 'right', 
    ellipse = FALSE,
      ellipseType = 't',
      ellipseLevel = 0.95,
      ellipseFill = TRUE,
      ellipseAlpha = 1/4,
    colLegendTitle = readable_variable,
      ellipseLineSize = 1.0,
   shape = "Translocation",
   shapeLegendTitle = "Chromosomal translocation",
    pointSize = 5,
    title = readable_variable,
    subtitle = 'PC1 versus PC2',)
 print(g)
ggsave(g, file = filename, width = 10, height = 6, dpi = 300)
}

   variables <- c("histology:ch1", "age:ch1", "celltype:ch1", "diseaselocation:ch1", "diseasestage:ch1", "diseasestate:ch1", "Sex:ch1", "tumorgrading:ch1" , "Survival_rounded")
make_pca_biplot(clinical_variable = "Histology", readable_variable = "Histology")
make_pca_biplot(clinical_variable = "Disease Stage", readable_variable = "Disease Stage")
make_pca_biplot(clinical_variable = "Fail", readable_variable = "Event")

```


#check expression of genes in gene signature
```{r}
head(expression_data)

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

expression_signature <- expression_data[row.names(expression_data) %in% hub_genes,]

```


#plot GSVA score function
To create GSVA score plots vs clinical variable with statistics. 
```{r}

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

#check dim
dim(expression_data)

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

genes <- list(genes = genes)

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

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


ARMS_samples <- subset(pd, Translocation == "P7F" | Translocation  == "P3F") %>% rownames()
ERMS_samples <- subset(pd, Translocation == "NEG" ) %>% rownames()


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

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

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

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

  # #select correct colour for histology
  #    if(all(str_detect(gsva_results$Histology, 'Alveolar|Embryonal|Mixed Alv/Emb'))){histology_colour <- c("violetred1","slateblue1", "navy")}
  #  if(all(!str_detect(gsva_results$Histology, "Alveolar"))){histology_colour <- c("slateblue1")}
  #     if(all(!str_detect(gsva_results$Histology, "Embryonal"))){histology_colour <-  c("violetred1", "navy")}

#then selects options for clinical variable picked
  print(clinical_variable)
  
  #clinical_variable <- c("Stage")
    if (clinical_variable == "Stage") {
         variable_levels <- c("Stage 1", "Stage 2", "Stage 3", "Stage 4", "NA")
    readable_clinical_variable <- c("stage")
          #remove comparisons where no enoguh observations
       #in ARMS samples there is only one stage 1
    if(subtype=="ARMS"){
      my_comparisons <- list(c("Stage 2", "Stage 3"), c("Stage 2", "Stage 4"), c("Stage 3", "Stage 4"))
    }
       if(subtype=="ERMS"){
      my_comparisons <- list( c("Stage 3", "Stage 4"))
    }
       
   } else if(clinical_variable == "Grade"){
      variable_levels<-c("IA", "IB", "IIA", "IIB", "III", "IV", "NA")
  readable_clinical_variable <- c("grade")
   my_comparisons <- list(c("IA", "IB"), c("IA", "IIA"), c("IA", "IIB"),c("IA", "III"), c("IA", "IV"), c("IB", "IIA"), c("IB", "IIB"), c("IB", "III"), c("IB", "IV"), c("IIA", "IIB"), c("IIA", "III"), c("IIA", "IV"), c("IIB", "III"), c("IIB", "IV"), c("III", "IV"))
  
    #remove comparisons where no enoguh observations
  #in ARMS there is only 1 IB and IIB so can't compare that
    if(subtype=="ARMS"){
      my_comparisons <- list( c("IA", "IIA"),c("IA", "III"), c("IA", "IV"), c("IIA", "III"), c("IIA", "IV"), c("III", "IV"))
    }
   if(subtype == "ERMS"){
      my_comparisons <- list( c("IA", "IIA"), c("IA", "IIB"), c("IA", "III"), c("IA", "IV"),c("IIA", "IIB"),c("IIA", "III"),c("IIA", "IV"), c("IIB", "III", c("IIB", "IV"), c("III", "IV")))
  }
   
} else if(clinical_variable == "disease_state"){
  
    my_comparisons <- list( c("Alive", "Dead"))
  variable_levels <- c("Alive", "Dead", "NA")
    readable_clinical_variable <- c("disease state")
    
}  else if(clinical_variable == "RG"){
  
    my_comparisons <- list( c("LOW RISK", "INTERMEDIATE RISK"), c("LOW RISK", "HIGH RISK"), c("INTERMEDIATE RISK", "HIGH RISK"))
  variable_levels <- c("LOW RISK", "INTERMEDIATE RISK", "HIGH RISK", "NA")
    readable_clinical_variable <- c("risk group")
}else if(clinical_variable == "Size"){
  
    my_comparisons <- list( c("<= 5 cm", "> 5 cm"))
  variable_levels <-  c("<= 5 cm", "> 5 cm", "NA")
    readable_clinical_variable <- c("tumour size")
}else if(clinical_variable == "Node"){
  
    my_comparisons <- list( c("N-0", "N-1"))
  variable_levels <-   c("N-0", "N-1", "NA")
    readable_clinical_variable <- c("node")
}else if(clinical_variable == "Invasion"){
  
    my_comparisons <- list( c("T-1", "T-2"))
  variable_levels <-  c("T-1", "T-2", "NA")
    readable_clinical_variable <- c("invasion")
} else if(clinical_variable == "Fail"){
  
    my_comparisons <- list( c("Event", "No event"))
  variable_levels <-   c("Event", "N-1", "No event")
    readable_clinical_variable <- c("event")
}
 
#reorder plot
try(gsva_results[[clinical_variable]] <- factor(gsva_results[[clinical_variable]], variable_levels))
 try(fct_relevel(gsva_results[[clinical_variable]]))
 
    #boxplot
    p <- try(ggboxplot(gsva_results  , x =  paste0(clinical_variable), y = "gsva_results", fill = paste0(clinical_variable), outlier.shape = NA)+labs(title= paste0(resistance_type, " resistance score and ", readable_clinical_variable),
        x =paste0(toupper(substr(readable_clinical_variable, 1, 1)), substr(readable_clinical_variable, 2, nchar(readable_clinical_variable)), sep=""), y = "GSVA score")+ theme( axis.text.x = element_text(angle = 90,  size = 12))+ scale_fill_brewer(palette="BuPu", guide = "none")+geom_jitter(aes(col=Translocation))+
  scale_color_manual(values = histology_colour)+ theme(legend.position = "right")+ stat_compare_means(comparisons = my_comparisons,label = "p.signif", method = paste0(correct_method)))
#stat_pwc(method = paste0(correct_method), label = "p.adj.signif", p.adjust.method = "BH",na.rm = TRUE,) 
    print(p)
ggsave(file= paste0(resistance_type, "_", clinical_variable, "_", filename, ".png"), width = 5, height = 5)


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

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

}

```


#GSVA score plots of gene signatures
```{r}

#using GSVA to calculate a chemo-resistant score
#intrinsic all samples -----------------------------------------------------------
setwd("C:/Users/Christina/Documents/RMS microarray/Triche/GSVA_plots_full_annotation")

#get genes 
genes <- read_xlsx("C:/Users/Christina/Documents/RMS RNA seq/Zoes data/CLONES ONLY ANALYSIS/intrinsic_resistance_signature.xlsx") #intrinsic sig
genes <- genes$Gene
#genes <- read_xlsx("C:/Users/Christina/Documents/RMS RNA seq/Luciles data/RH4 WGCNA/RH4_acquired_res_signature.xlsx") #acquired rh4 sig
#genes <- read_xlsx("C:/Users/Christina/Documents/RMS RNA seq/Luciles data/RMSYM WGCNA/RMSYM_acquired_resistance_gene_signature.xlsx") #RMSYM acquired sig

variables <- c( "Grade", "Stage", "Node", "Invasion", "Survival", "FFS", "Fail")

#for intrinsic ARMS only 

for (i in variables){
  plot_resistance_score_with_statistics(genes = genes, subtype = "ARMS",  clinical_variable = i, resistance_type = "Intrinsic")}

#acquired RH4-----------------------------------------------------------------

#get genes 
#genes <- read_xlsx("C:/Users/Christina/Documents/RMS RNA seq/Zoes data/CLONES ONLY ANALYSIS/intrinsic_resistance_signature.xlsx") #intrinsic sig
genes <- read_xlsx("C:/Users/Christina/Documents/RMS RNA seq/Luciles data/RH4 WGCNA/RH4_acquired_res_signature.xlsx") #acquired rh4 sig
genes <- genes$Gene
#genes <- read_xlsx("C:/Users/Christina/Documents/RMS RNA seq/Luciles data/RMSYM WGCNA/RMSYM_acquired_resistance_gene_signature.xlsx") #RMSYM acquired sig

for (i in variables){
  plot_resistance_score_with_statistics(genes = genes, subtype = "ARMS",  clinical_variable = i, resistance_type = "Acquired")}

#acquired RMSYM -----------------------------------------------------------------

#get genes 
#genes <- read_xlsx("C:/Users/Christina/Documents/RMS RNA seq/Zoes data/CLONES ONLY ANALYSIS/intrinsic_resistance_signature.xlsx") #intrinsic sig
#genes <- read_xlsx("C:/Users/Christina/Documents/RMS RNA seq/Luciles data/RH4 WGCNA/RH4_acquired_res_signature.xlsx") #acquired rh4 sig
genes <- read_xlsx("C:/Users/Christina/Documents/RMS RNA seq/Luciles data/RMSYM WGCNA/RMSYM_acquired_resistance_gene_signature.xlsx") #RMSYM acquired sig
genes <- genes$Gene

for (i in variables){
  plot_resistance_score_with_statistics(genes = genes, subtype = "ERMS",  clinical_variable = i, resistance_type = "Acquired")}
```


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

setwd("~/RMS microarray/Triche/GSVA DEGs")

variables <- c( "Grade", "Stage", "Node", "Invasion", "Survival", "FFS", "Fail")

#for intrinsic DEGs genes all samples #-----------------------------------------------------------

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

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

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

for (i in variables){
  plot_resistance_score_with_statistics(clinical_variable = i, subtype = "ARMS", resistance_type = "Intrinsic",genes = genes,)}

#using overlapping DEGs upregulated in RH4 IFO and VCR resistant----------------------------------------------------

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

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

genes <- intersect(DEGS_RH4_IFO_UP$X,DEGS_RH4_VCR_UP$X)


for (i in variables){
  plot_resistance_score_with_statistics(clinical_variable = i, subtype = "ARMS", resistance_type = "Acquired", genes = genes)}
  
  
#using overlapping DEGs upregulated in RMSYM IFO and VCR resistant-----------------------------------------------------------------------------------

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

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

genes <- intersect(DEGS_RMSYM_IFO_UP$X,DEGS_RMSYM_VCR_UP$X)

for (i in variables){
  plot_resistance_score_with_statistics(clinical_variable = i, subtype = "ERMS", resistance_type = "Acquired", genes = genes)}

```


#Scatterplots
```{r}

expression_data

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



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

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

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

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

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

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



```

#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
input_mat = t(expression_data)
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)



# get trait data 
traitData=pd
dim(traitData)
names(traitData)

#remove columns not needed
allTraits <- traitData %>%
  select("age:ch1", "diseaselocation:ch1", "diseasestage:ch1", "diseasestate:ch1", "histology:ch1", "Sex:ch1", "tumorgrading:ch1")

#write.csv(allTraits, "allTraits.csv")

femaleSamples =rownames(datET);
traitRows =match(femaleSamples, rownames(allTraits));
datTraits = allTraits[traitRows, -1];
rownames(datTraits) = rownames(allTraits)


# datTraits<- as.matrix(datTraits)
allTraits$`age:ch1` <- allTraits$`age:ch1`%>% as.factor() %>% as.numeric()
allTraits$`diseaselocation:ch1` <- allTraits$`diseaselocation:ch1`%>% as.factor() %>% as.numeric()
allTraits$`diseasestage:ch1` <- allTraits$`diseasestage:ch1` %>% as.factor() %>% as.numeric()
allTraits$`diseasestate:ch1` <- allTraits$`diseasestate:ch1` %>% as.factor() %>% as.numeric()
allTraits$Histology <- allTraits$Histology %>% as.factor() %>% as.numeric()
allTraits$`Sex:ch1` <- allTraits$`Sex:ch1` %>% as.factor() %>% as.numeric()
allTraits$`tumorgrading:ch1` <- allTraits$`tumorgrading:ch1` %>% as.factor() %>% as.numeric()

# 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(allTraits, signed = FALSE);
# Plot the sample dendrogram and the colors underneath.# Plot therainbow() sample dendrogram and the colors underneath.
pdf(file="trait_denderogram.pdf", width = 21, height = 12)
plotDendroAndColors(sampleTree2, traitColors,
                    groupLabels = names(allTraits),
                    main = "Sample dendrogram and trait heatmap")
dev.off()

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


#see what genes are contributing to the modules on a heatmap
t(datME) %>% heatmap(col = brewer.pal(11,"RdBu"))

#plot heatmap of genes of interest
h1<-Heatmap(as.matrix(t(datME)),cluster_rows=TRUE,cluster_columns=FALSE,use_raster=FALSE,row_dend_reorder = TRUE,show_row_dend = TRUE,show_row_names=TRUE,show_column_names = TRUE, heatmap_legend_param = list(title = "Z-score"))

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




vis <- exportNetworkToVisANT(adjacency, weighted = TRUE, threshold = 0.55)


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(moduleTraitPvalue, 1),signif(moduleTraitCor, 1) ,sep= "")
dim(textMatrix)= dim(moduleTraitCor)
par(mar= c(7,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(allTraits),
            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"))
dev.off()


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


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 darkmagenta

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

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

 genes <- gene_info$genes
 
#select genes with GS. >0.7 and MM > 0.6 as hub genes
hub_genes <-  subset(gene_info, GS. <-0.6 & MM.darkmagenta >0.6)
 
write.csv(hub_genes, file = "hub_genes_darkmagenta_0.6.csv")
#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.6 as hub genes
hub_genes <-  subset(gene_info, GS. <-0.6 & MM.darkolivegreen >0.6)
 
write.csv(hub_genes, file = "hub_genes_darkolivegreen_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")
```


#WGCNA ARMS
```{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)

#remove ERMS samples (keep ARMS and mixed histology)
embryonal_samples <- subset(pd, Histology == "Embryonal") %>% rownames() #68 ERMS samples
#from expression data
ARMS_expression_data <- expression_data[, !colnames(expression_data) %in% embryonal_samples]
#from pheno data
ARMS_pd <- pd[ !rownames(pd) %in% embryonal_samples,]

#check dim of expressiona and pheno data
dim(ARMS_expression_data) #70 ARMS/mixed samples
dim(ARMS_pd)

#need normalised counts in matrix with rows as samples and columns as genes
input_mat = t(ARMS_expression_data)
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)



# get trait data 
traitData=ARMS_pd
dim(traitData)
names(traitData)

#remove columns not needed
allTraits <- traitData %>%
  select("age:ch1", "diseaselocation:ch1", "diseasestage:ch1", "diseasestate:ch1", "histology:ch1", "Sex:ch1", "tumorgrading:ch1")

#write.csv(allTraits, "allTraits.csv")

femaleSamples =rownames(datET);
traitRows =match(femaleSamples, rownames(allTraits));
datTraits = allTraits[traitRows, -1];
rownames(datTraits) = rownames(allTraits)


# datTraits<- as.matrix(datTraits)
allTraits$`age:ch1` <- allTraits$`age:ch1`%>% as.factor() %>% as.numeric()
allTraits$`diseaselocation:ch1` <- allTraits$`diseaselocation:ch1`%>% as.factor() %>% as.numeric()
allTraits$`diseasestage:ch1` <- allTraits$`diseasestage:ch1` %>% as.factor() %>% as.numeric()
allTraits$`diseasestate:ch1` <- allTraits$`diseasestate:ch1` %>% as.factor() %>% as.numeric()
allTraits$Histology <- allTraits$Histology %>% as.factor() %>% as.numeric()
allTraits$`Sex:ch1` <- allTraits$`Sex:ch1` %>% as.factor() %>% as.numeric()
allTraits$`tumorgrading:ch1` <- allTraits$`tumorgrading:ch1` %>% as.factor() %>% as.numeric()

# 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(allTraits, signed = FALSE);
# Plot the sample dendrogram and the colors underneath.# Plot therainbow() sample dendrogram and the colors underneath.
#get names of traits for the dendrogram
trait_names <- c("Age", "Disease location", "Disease stage", "Disease state", "Histology", "Sex", "Tumour grade")
pdf(file="trait_denderogram.pdf", width = 21, height = 12)
plotDendroAndColors(sampleTree2, traitColors,
                    groupLabels = trait_names,
                    main = "Sample dendrogram and trait heatmap")
dev.off()



grp <- c(1,1,1, 2,2)
k_3 <- cutree(dend,k = 3, order_clusters_as_data = FALSE) 
# The FALSE above makes sure we get the clusters in the order of the
# dendrogram, and not in that of the original data. It is like:
# cutree(dend, k = 3)[order.dendrogram(dend)]
the_bars <- cbind(grp, k_3)
dend %>% set("labels", "") %>% plot
colored_bars(colors = the_bars, dend = dend)




#take a look at the samples that are most disimilar on a PCA
potential_outliers <- c("GSM2435289", "GSM2435368")

    variables <- c("histology:ch1")
  
  for (i in variables)
  {   
    #find length of variables
  n = length(unique(p$metadata[[i]]))
  
 p$metadata[[i]] <-  as.factor(p$metadata[[i]])
  
#set colours
  #myCols <- setNames(distinctColorPalette(n), unique(p$metadata[[i]]))
  myCols <- viridis(n, option = "G", direction = -1, end = 0.9)
  
  legend_title <- sapply(strsplit(i,":"), `[`, 1)
  
  #make filename title
   j <- sapply(strsplit(i,":"), `[`, 1)
   filename <- paste0("PCA_biplot_all_labelled",j,".png")
   
   #label <- subset(p$metadata, Histology == "Mixed Alv/Emb") %>% rownames() #label only mixed samples

 g <- biplot(p,
          selectLab = potential_outliers,
         #lab = rownames(p$metadata),
         labSize = 5,
         max.overlaps = 200,
    colby = paste0(i),
    colkey = myCols,
    hline = 0, vline = 0,
    legendPosition = 'right', 
    ellipse = FALSE,
      ellipseType = 't',
      ellipseLevel = 0.95,
      ellipseFill = TRUE,
      ellipseAlpha = 1/4,
    colLegendTitle = "Histology",
      ellipseLineSize = 1.0,
 #  shape = 'histology:ch1',
 #  shapeLegendTitle = "Histology",
    pointSize = 5,
    title = 'PCA bi-plot',
    subtitle = 'PC1 versus PC2',)
ggsave(g, file = filename, width = 8, height = 8, dpi = 300)
  }




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

powers = c(1:10)  # in practice this should include powers up to 20.
# choose power based on SFT criterion

# Plot the results:
par(mfrow = c(1, 2))
# SFT index as a function of different powers
plot(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2], 
    xlab = "Soft Threshold (power)", ylab = "SFT, signed R^2", type = "n", main = paste("Scale independence"))
text(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2], 
    labels = powers, col = "red")
# this line corresponds to using an R^2 cut-off of h
abline(h = 0.9, col = "red")
# Mean connectivity as a function of different powers
plot(sft$fitIndices[, 1], sft$fitIndices[, 5], type = "n", xlab = "Soft Threshold (power)", 
    ylab = "Mean Connectivity", main = paste("Mean connectivity"))
text(sft$fitIndices[, 1], sft$fitIndices[, 5], labels = powers, col = "red")
     

#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 = 9

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

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;


#see what samples are contributing to the modules on a heatmap
t(datME) %>% heatmap(col = brewer.pal(11,"RdBu"))

#plot heatmap of genes of interest
h1<-Heatmap(as.matrix(t(datME)),cluster_rows=TRUE,cluster_columns=FALSE,use_raster=FALSE,row_dend_reorder = TRUE,show_row_dend = TRUE,show_row_names=TRUE,show_column_names = TRUE, heatmap_legend_param = list(title = "Z-score"))

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




vis <- exportNetworkToVisANT(adjacency, weighted = TRUE, threshold = 0.55)


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, allTraits, 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(7,12, 6, 5))
#display the corelation values with a heatmap plot
#INCLUE THE NEXT LINE TO SAVE TO FILE
pdf("module_trait.pdf", pointsize = 16, width = 10, height = 10)
labeledHeatmap(Matrix= moduleTraitCor,
            xLabels= colnames(allTraits),
            yLabels= names(MEs),
            ySymbols= names(MEs),
            colorLabels= FALSE,
            colors= blueWhiteRed(50),
            textMatrix= textMatrix,
            setStdMargins= TRUE,
            cex.text= 0.5,
            zlim= c(-1,1),
            main= paste("Module-trait relationships"))
dev.off()

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


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 darkmagenta

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

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

 genes <- gene_info$genes
 
#select genes with GS. >0.7 and MM > 0.6 as hub genes
hub_genes <-  subset(gene_info, GS. <-0.6 & MM.darkmagenta >0.6)
 
write.csv(hub_genes, file = "hub_genes_darkmagenta_0.6.csv")
#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.6 as hub genes
hub_genes <-  subset(gene_info, GS. <-0.6 & MM.darkolivegreen >0.6)
 
write.csv(hub_genes, file = "hub_genes_darkolivegreen_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")
```
