---
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(oligo)
library(tidyverse)
library("fuzzyjoin")
library(ggpubr)
library(hgu133plus2.db)
library(affycoretools)
library(viridis)
library("ArrayExpress")
library(tidyverse)
library(PCAtools)
library(GSVA)
#library("arrayQualityMetrics")
library("readxl")
library(jetset)
library("affyio")
library(dplyr)
library(oligoClasses)
```

# Import data from arrayexpress

```{r}

setwd("~/RMS microarray/Williamson")

#The output object is an AffyBatch if the ArrayExpress identifier corresponds to an Affymetrix experiment, it is an ExpressionSet if the identifier corresponds to a one colour non Affymetrix experiment and it is a NChannelSet if the identifier corresponds to a two colours experiment. In this example, ’E-MEXP-21’ being an Affymetrix experiment, ’rawset’ is an AffyBatch. The expression values of ’rawset’ are the content from the CEL files. The phenoData of ’rawset’ contains the whole sample annotation file content from the MAGE-TAB sdrf file. The featureData of ’rawset’ contains the whole feature annotation file content from the MAGE-TAB adf file. The experimentData of ’rawset’ contains the experiment annotation file content from the MAGE-TAB idf file.

#load from Array Express
# rawset = ArrayExpress("E-TABM-1202")
# AEset = rawset
# getGEOSuppFiles("E-TABM-1202")

# Read all CEL files in a directory
data_affy <- ReadAffy(celfile.path = "~/RMS microarray/Williamson/raw_data")
AEset = data_affy

#plots of raw data ---------------------------------------------------------

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

tiff(file="boxplot_raw_data_Williamson.jpeg", 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(AEset, which='pm', ylab = "log-2 intensity",   las = 2, 
        cex.axis = 0.6, main = "Boxplot of log-2 intensities for the raw data (Williamson)")
dev.off()

#QC

 # arrayQualityMetrics(expressionset =AEset,
 # outdir = "Report_forAEset",
 # force = TRUE,
 # do.logtransform = TRUE)

```


#Normalisation
```{r}
#perform normalisation with RMA
oligo_normalised <- normalize(AEset,method='quantile',which='pm')

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

#histogram
tiff(file="histogram_normalised_Williamson.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 (Williamson)")
dev.off()

```

#Summarisation
```{r}
#from probe level to feature level 
oligo_summarised <- rma(AEset,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_Williamson.jpeg", width=19, height=6, units="in", res=300)
boxplot(oligo_summarised, which='pm', ylab = "log-2 intensity",   las = 2, 
        cex.axis = 0.6,  main = "Boxplot of log-2 intensities after quantile normalisation and summarisation")
dev.off()

#histogram
tiff(file="histogram_normalised_and_summarised_Williamson.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 quantile normalisation and summarisation)")
dev.off()

```

#Performing RMA

```{r}
# returns an ExpressionSet object containing the background corrected, normalised, and summarised expression data
AEsetnorm <- rma(AEset)

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

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

#get normalised expression
expr = exprs(AEsetnorm)

#QC on processed data
 # arrayQualityMetrics(expressionset = AEsetnorm,
 # outdir = "Report_for_AEsetnorm",
 # intgroup = c("condition", "batch"), #change these to relpase/met
 # 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(AEsetnorm)
#probes <- row.names(expression_data)

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

#use Jetset for one to one gene-probe mapping
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("hgu133plus2", 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 = "hgu133plus2.db" )
   symbol = gset_genes@featureData@data[["SYMBOL"]]
 
   #replace probes as symbols
 expression_data <- exprs(gset_genes)
 rownames(expression_data) <- symbol
 
dim(expression_data)

```


#get clinical data
```{r}
varLabels(AEset)
pd <- pData(AEset)

pd

#number of samples
dim(pd) #101

#select only cols want
pd <- pd %>%
  dplyr::select(Source.Name, Factor.Value..survival_time., Factor.Value..Status., Factor.Value..LymphNode.)

#add column that has cel file names
pd$cel_file_name <- rownames(pd)

#rename source.name col to match patient data 
colnames(pd)[which(names(pd) == "Source.Name")]  <- "Patient"

#get metadata (full clinical annotation)
metadata <- read_xls("~/RMS microarray/Williamson/ITCC patient clinical annotation.xls")

#merge pd with clinical annotation
#using regrex to match near patterns instead of exact matches
clinical_annotation <- regex_full_join(pd, metadata, by = c("Patient")) %>% as.data.frame()

#check they match
clinical_annotation$Patient.x == clinical_annotation$Patient.y

#remove rows where cel file = NA
clinical_annotation <- clinical_annotation[!is.na(clinical_annotation$Patient.x), ]

#check they match
clinical_annotation$Patient.x == clinical_annotation$Patient.y

pd <- clinical_annotation
rownames(pd) <- pd$cel_file_name
```

# clean data
```{r}
#for clinical data
pd_williamson <- read.csv("~/RMS microarray/Williamson/metadata_williamson.csv")

pd_williamson$Fusion.transcript...gene.name <- gsub("Pax3", "P3F", pd_williamson$Fusion.transcript...gene.name)
pd_williamson$Fusion.transcript...gene.name <- gsub("Pax7", "P7F", pd_williamson$Fusion.transcript...gene.name)
pd_williamson$Fusion.transcript...gene.name[is.na(pd_williamson$Fusion.transcript...gene.name)] <- "NEG"


#make sure rownames pd = colnames EM
rownames(pd_williamson) <- paste0(pd_williamson$X.Oligo.cDNA.Array, ".CEL")
pd <- pd_williamson

#clean 'EM'
pd$Stage...IRS.group <- gsub("EM", NA, pd$Stage...IRS.group)
pd$Metastasis...at.diagnosis..recoded. <- gsub("EM", NA, pd$Metastasis...at.diagnosis..recoded.)
pd$Risk..Group <- gsub("EM", NA, pd$Risk..Group)
pd$Tumor.size..recoded. <-  gsub("EM", NA, pd$Tumor.size..recoded.)
pd$Stage...SIOP <-  gsub("EM", NA, pd$Stage...SIOP)

#replace with NA
pd$ID.Histology <- gsub("NA", NA, pd$ID.Histology)
pd$Organism.Part <- gsub("NA", NA, pd$Organism.Part)
pd$Metastasis...Bone_and_Bone.marrow..recoded.<- gsub("NA", NA, pd$Metastasis...Bone_and_Bone.marrow..recoded.)
pd$Factor.Value..LymphNode. <- gsub(" ", NA, pd$Factor.Value..LymphNode.)
pd$Factor.Value..LymphNode. <- gsub("NA", NA, pd$Factor.Value..LymphNode.)
#rename
pd$Factor.Value..Status. <- gsub("dead from disease", "dead", pd$Factor.Value..Status.)

#change to numeric
pd$Age.at.diagnosis..year. <- as.numeric(pd$Age.at.diagnosis..year.)
pd$Factor.Value..survival_time. <- as.numeric(pd$Factor.Value..survival_time.)


pd$survival_rounded <- round(as.numeric(pd$Factor.Value..survival_time.))
pd$Age.at.diagnosis..year. <- round(as.numeric(pd$Age.at.diagnosis..year.))

pd$Histology...Super.class <- gsub("RMS alveolar", "Alveolar", pd$Histology...Super.class)
pd$Histology...Super.class <- gsub("RMS embryonal", "Embryonal", pd$Histology...Super.class)

```



#plot clinical characteristics
```{r}

colnames(pd)

setwd("~/RMS microarray/Williamson/clinical_info")

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

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

  #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))+ theme_bw()+ theme(axis.text.x=element_text(size=16),axis.text.y=element_text(size=16),
        plot.title=element_text(size=16),axis.title.x=element_text(size=16),axis.title.y=element_text(size=16)) + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("Age") +xlab("Age (years)") + ylab("Number of samples") + guides(fill="none")
  
  
 ggsave(p, file=paste0("age.png"), width = 30, height = 15, units = "cm")

   #tumour location---------------------------------------------
    pd$Location..primary.tumor.
  
  df <- data.frame(
  value = pd$Location..primary.tumor.)
  
    #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))+ theme_bw() + 
    geom_bar(stat = "count")+
    scale_fill_manual(values=myCols) + ggtitle("Tumour location",) + theme(axis.text.x=element_text(size=16,angle = 90, vjust = 0.5, hjust=1),axis.text.y=element_text(size=16),
        plot.title=element_text(size=16),axis.title.x=element_text(size=16),axis.title.y=element_text(size=16)) +xlab("Tumour location") + ylab("Number of samples") + guides(fill="none")
   ggsave(p, file=paste0("location.png"), width = 10, height = 12, units = "cm")

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


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

  #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))+theme_bw() + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("Histology",) + theme(axis.text=element_text(size=16),
        axis.title=element_text(size=16)) + theme(axis.text.x=element_text(size=16,angle = 90, vjust = 0.5, hjust=1),axis.text.y=element_text(size=16),
        plot.title=element_text(size=16),axis.title.x=element_text(size=16),axis.title.y=element_text(size=16))+xlab("Histology") + ylab("Number of samples") + guides(fill="none")
 ggsave(p, file=paste0("histology.png"), width = 8, height = 12, units = "cm")

 #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))+ theme_bw() + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("Sex",) + theme(axis.text=element_text(size=16),
        axis.title=element_text(size=16)) + theme(axis.text=element_text(size=16),
        axis.title=element_text(size=16)) + theme(axis.text.x=element_text(size=16),axis.text.y=element_text(size=16),
        plot.title=element_text(size=16),axis.title.x=element_text(size=16),axis.title.y=element_text(size=16))+xlab("Sex") + ylab("Number of samples") + guides(fill="none")
 ggsave(p, file=paste0("sex.png"), width = 8, height = 10, units = "cm")
 
 
 
#metastasis---------------------------------------------------------------------
df <- data.frame(
value = pd$Metastasis...at.diagnosis..recoded.)

  #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))+theme_bw() + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("Metastasis",)+ theme(axis.text=element_text(size=16),
        axis.title=element_text(size=16)) + theme(axis.text.x=element_text(size=16),axis.text.y=element_text(size=16),
        plot.title=element_text(size=16),axis.title.x=element_text(size=16),axis.title.y=element_text(size=16))+xlab("Metastasis") + ylab("Number of samples") + guides(fill="none")
 ggsave(p, file=paste0("Metastasis.png"), width = 7, height = 10, units = "cm")
 
#IRS stage---------------------------------------------------------------------
df <- data.frame(
value = pd$Stage...IRS.group)

  #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))+theme_bw() + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("IRS group",)+ theme(axis.text=element_text(size=16),
        axis.title=element_text(size=16)) + theme(axis.text.x=element_text(size=16),axis.text.y=element_text(size=16),
        plot.title=element_text(size=16),axis.title.x=element_text(size=16),axis.title.y=element_text(size=16))+xlab("IRS group") + ylab("Number of samples") + guides(fill="none")
 ggsave(p, file=paste0("IRS_group.png"), width = 8, height = 10, units = "cm")
 
 
 #SIOP stage---------------------------------------------------------------------
df <- data.frame(
value = pd$Stage...SIOP)

  #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))+theme_bw() + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("SIOP stage",)+ theme(axis.text=element_text(size=16),
        axis.title=element_text(size=16)) + theme(axis.text.x=element_text(size=16),axis.text.y=element_text(size=16),
        plot.title=element_text(size=16),axis.title.x=element_text(size=16),axis.title.y=element_text(size=16))+xlab("SIOP stage") + ylab("Number of samples") + guides(fill="none")
 ggsave(p, file=paste0("SIOP_stage.png"), width = 8, height = 10, units = "cm")
 
 
 #tumour size-------------------------------------
df <- data.frame(
value = pd$Tumor.size..recoded.)

  #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))+theme_bw() + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("Tumour size",)+ theme(axis.text=element_text(size=16),
        axis.title=element_text(size=16)) + theme(axis.text.x=element_text(size=16,angle = 90, vjust = 0.5, hjust=1),axis.text.y=element_text(size=16),
        plot.title=element_text(size=16),axis.title.x=element_text(size=16),axis.title.y=element_text(size=16))+xlab("Tumour size") + ylab("Number of samples") + guides(fill="none")
 ggsave(p, file=paste0("tumour_size.png"), width = 16, height = 12, units = "cm")
 
 #fusion gene status---------------------------------------------------------------
 df <- data.frame(
value = pd$Fusion.transcript...Presence)

  #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))+theme_bw() + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("Chromosomal\ntranslocation\nstatus",)+ theme(axis.text=element_text(size=16),
        axis.title=element_text(size=16)) + theme(axis.text.x=element_text(size=16),axis.text.y=element_text(size=16),
        plot.title=element_text(size=16),axis.title.x=element_text(size=16),axis.title.y=element_text(size=16))+xlab("Chromosomal\ntranslocation\nstatus") + ylab("Number of samples") + guides(fill="none")
 
  #risk group-----------------------------------------------------------------
 df <- data.frame(
value = pd$Risk..Group)

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

    p <- df %>%
  mutate(value = fct_relevel(value, 
            "standard", "high", "very high")) %>%
ggplot(aes(  fill=value, x=value))+theme_bw() + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("Risk group",)+ theme(axis.text=element_text(size=16),
        axis.title=element_text(size=16)) + theme(axis.text.x=element_text(size=16,angle = 90, vjust = 0.5, hjust=1),axis.text.y=element_text(size=16),
        plot.title=element_text(size=16),axis.title.x=element_text(size=16),axis.title.y=element_text(size=16))+xlab("Risk group") + ylab("Number of samples") + guides(fill="none")
    
 ggsave(p, file=paste0("Risk_group.png"), width = 11, height = 12, units = "cm")
 
  #complete remission acheievd -----------------------------------------------------------------
 df <- data.frame(
value = pd$Event...First.Complete.remission...Achieved)

  #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))+theme_bw() + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("Complete remission\nachieved",)+ theme(axis.text=element_text(size=16),
        axis.title=element_text(size=16)) + theme(axis.text.x=element_text(size=16),axis.text.y=element_text(size=16),
        plot.title=element_text(size=16),axis.title.x=element_text(size=16),axis.title.y=element_text(size=16))+xlab("Complete remission\nachieved") + ylab("Number of samples") + guides(fill="none")
 
   #first event occurred -----------------------------------------------------------------
 df <- data.frame(
value = pd$First.event...occurred)

  #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))+theme_bw() + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("First event\noccurred",)+ theme(axis.text=element_text(size=16),
        axis.title=element_text(size=16)) + theme(axis.text.x=element_text(size=16),axis.text.y=element_text(size=16),
        plot.title=element_text(size=16),axis.title.x=element_text(size=16),axis.title.y=element_text(size=16))+xlab("First event\noccurred") + ylab("Number of samples") + guides(fill="none")
   #first event type -----------------------------------------------------------------
 df <- data.frame(
value = pd$First.event...Type)
  
  df$value = str_wrap(df$value, width = 25)

  #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))+theme_bw() + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("First event type",)+ theme(axis.text=element_text(size=16),
        axis.title=element_text(size=16)) + theme(axis.text.x=element_text(size=16,angle = 90, vjust = 0.5, hjust=1),axis.text.y=element_text(size=16),
        plot.title=element_text(size=16),axis.title.x=element_text(size=16),axis.title.y=element_text(size=16))+xlab("First event type") + ylab("Number of samples") + guides(fill="none")
  
 ggsave(p, file=paste0("first_event_type.png"), width = 22, height = 20, units = "cm")
 
 
 #survival time--------------------------------------------------------------

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

  #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))+theme_bw() + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("Survival time",)+ theme(axis.text=element_text(size=16),
        axis.title=element_text(size=16)) + theme(axis.text.x=element_text(size=16),axis.text.y=element_text(size=16),
        plot.title=element_text(size=16),axis.title.x=element_text(size=16),axis.title.y=element_text(size=16))+xlab("Survival time (years)") + ylab("Number of samples") + guides(fill="none")
  
 ggsave(p, file=paste0("survival.png"), width = 30, height = 15, units = "cm")
 
 
  #disease state --------------------------------------------------------------

df <- data.frame(
value = as.factor(pd$Factor.Value..Status.))

  #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("Disease state",) + theme(axis.text=element_text(size=16),
        axis.title=element_text(size=16))+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+xlab("Disease state") + ylab("Number of samples") + guides(fill="none")+ theme_bw(base_size = 20)
 ggsave(p, file=paste0("disease_state.png"), width = 30, height = 15, units = "cm")
 
   #lymoh node --------------------------------------------------------------

df <- data.frame(
value = as.factor(pd$Factor.Value..LymphNode.))

  #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))+theme_bw() + 
  geom_bar(stat = "count")+
  scale_fill_manual(values=myCols) + ggtitle("Lymph node\ninvolvement",)+ theme(axis.text=element_text(size=16),
        axis.title=element_text(size=16)) + theme(axis.text.x=element_text(size=16),axis.text.y=element_text(size=16),
        plot.title=element_text(size=16),axis.title.x=element_text(size=16),axis.title.y=element_text(size=16))+xlab("Lymph node\ninvolvement") + ylab("Number of samples") + guides(fill="none")
  
```



#PCA
```{r}
setwd("~/RMS microarray/Williamson/PCA")

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

colnames(expression_data)
rownames(pd)

#perform 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))
  
  
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 = "Fusion.transcript...Presence",
   shapeLegendTitle = "Fusion gene presence",
    pointSize = 5,
    title = readable_variable,
    subtitle = 'PC1 versus PC2',)
 print(g)
ggsave(g, file = filename, width = 7, height = 4, dpi = 200)
  }
  
   
   make_pca_biplot("First.event...occurred", "First event occured" )
    make_pca_biplot("Location..primary.tumor.", "Tumour location" )
  make_pca_biplot("Histology...Super.class", "Histology" )
  make_pca_biplot("Sex", "Sex" )
  make_pca_biplot("Tumor.size..recoded.", "Tumour size" )
   make_pca_biplot("Risk..Group", "Risk group" )
  make_pca_biplot("Stage...IRS.group", "IRS group" )
  make_pca_biplot("Metastasis...at.diagnosis..recoded.", "Metastasis at diagnosis" )
  make_pca_biplot("Event...First.Complete.remission...Achieved", "Complete remission achieved" )
  make_pca_biplot("First.event...occurred", "First event occurred" )
  make_pca_biplot("First.event...Type", "First event type" )
  make_pca_biplot("Age.at.diagnosis..year.", "Age" )
   make_pca_biplot("survival_rounded", "Survival" )
  make_pca_biplot("Factor.Value..LymphNode.", "Lymph node involvement" )
  make_pca_biplot("Stage...SIOP", "SIOP stage" )
  
  
 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)

  variables <- c("Factor.Value..Age.", "Factor.Value..Histology.",  "Factor.Value..IRS_TNM_STAGE.", "Factor.Value..primary_site_status.", "Factor.Value..LymphNode.", "Factor.Value..tumor_size.", "Factor.Value..survival_time.", "Factor.Value..chromosomal_translocation.",  "Factor.Value..ClinicalGroup.", "Factor.Value..primary_site." ,"Factor.Value..Metastasis." , "Factor.Value..Sex." ,"Factor.Value..Status." )
  
      variables <- c("Factor.Value..Status.")
    
    for (i in variables)
    {   
      #find length of variables
    n = length(unique(p$metadata[[i]]))
  #set colours
      myCols <- viridis(n, option = "G", direction = -1, end = 0.9)
    
    #make filename title
     j <- sapply(strsplit(i,":"), `[`, 1)
     filename <- paste0("PCA_pairs_plot_",j,".png")
   
    g <- pairsplot(p,
        components = getComponents(p, c(1:10)),
        triangle = TRUE, trianglelabSize = 12,
        hline = 0, vline = 0,
        titleLabSize= 20,
        pointSize = 2,
        gridlines.major = FALSE, gridlines.minor = FALSE,
       colby = paste0(i), 
       colkey = myCols,
        legendPosition = "none",
        title = 'Pairs plot', plotaxes = FALSE,
        margingaps = unit(c(-0.01, -0.01, -0.01, -0.01), 'cm'))
   ggsave(g, file = filename, width = 8.5, height = 5, dpi = 300)
    }
  
```


#function to plot GSVA scores
```{r}

 # For the Williamson dataset, this included metastasis at diagnosis, IRS stage, tumour size, risk group, complete remission achieved, first event occurred, survival time and lymph node involvement.

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

#check dim
dim(expression_data)

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

#change genes to list
genes <- list(genes = genes)

#do GSVA to score samples
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$TNM_stage<-  pd[,'Factor.Value..IRS_TNM_STAGE.']
# pheno_data$metastasis <-  pd[,'Factor.Value..Metastasis.']
 pheno_data$Histology <- pd[,'Histology...Super.class']
# pheno_data$clinical_group <- pd[,'Factor.Value..ClinicalGroup.']
# pheno_data$lymphnode <- pd[,'Factor.Value..LymphNode.']
# pheno_data$tumour_size <- pd[,'Factor.Value..tumor_size.']
# pheno_data$disease_state <- pd[,'Factor.Value..Status.']
# pheno_data$age_status <- pd[,'Factor.Value..AgeStatus.']
# pheno_data$primary_site <- pd[,'Factor.Value..primary_site.']
# pheno_data$primary_site_status <- pd[,'Factor.Value..primary_site_status.']
# pheno_data$COG_classifier <- pd[,'Factor.Value..COGClassifier.']
# pheno_data$survival_time <- pd[,'Factor.Value..survival_time.']
pheno_data$first_event_occurred <- pd[,'First.event...occurred']
pheno_data$metastasis <- pd[,'Metastasis...at.diagnosis..recoded.']
pheno_data$IRS_group <- pd[,'Stage...IRS.group']
pheno_data$tumour_size <- pd[,'Tumor.size..recoded.']
pheno_data$risk_group <- pd[,'Risk..Group']
pheno_data$complete_remission_achieved <- pd[,'Event...First.Complete.remission...Achieved']
pheno_data$lymphnode <- pd[,'Factor.Value..LymphNode.']
pheno_data$Survival <- pd[,'Factor.Value..survival_time.']
pheno_data$Translocation <- pd[,'Fusion.transcript...gene.name']
pheno_data$SIOP_stage <- pd[,'Stage...SIOP']

pheno_data <- pheno_data %>% remove_rownames %>% column_to_rownames(var=".")

ARMS_samples <- subset(pd, Fusion.transcript...gene.name == "P7F" | Fusion.transcript...gene.name  == "P3F"| Fusion.transcript...gene.name  == "Ncoa1P3F") %>% rownames()
ERMS_samples <- subset(pd, Fusion.transcript...gene.name == "NEG") %>% rownames()


if(subtype=="ARMS"){
  #keep ARMS samples
pheno_data <- pheno_data[rownames(pheno_data) %in% ARMS_samples, ]
histology_colour <- c("#009E73","#CC79A7", "#0072B2")
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}"
  }

#then selects options for clinical variable picked
  
  print(clinical_variable)

  if (clinical_variable == "clinical_group") {
    #specify variables
    my_comparisons <- list( c("I", "II,III"))
   variable_levels <- c("I", "II,III", "NA")
    readable_clinical_variable <- c("clinical group")
  } else if(clinical_variable == "TNM_stage"){
  my_comparisons <- list( c("1", "2"), c("1", "3"), c("1", "4"), c("2", "3"), c("2", "4"), c("3", "4"))
  variable_levels<-c("1", "2", "3", "4", "NA")
  readable_clinical_variable <- c("TNM stage")
} else if(clinical_variable == "metastasis"){
    readable_clinical_variable <- c("metastasis")
    my_comparisons <- list( c("0", "1"))
    variable_levels <- c("0", "1", "EM")
  } else if(clinical_variable == "lymphnode"){
   my_comparisons <- list(c("N0", "N1"))
    variable_levels <- c("N0", "N1","NA")
    readable_clinical_variable <- c("lymph node")
  } else if(clinical_variable == "tumour_size"){
    my_comparisons <- list( c("<= 5 cm", "> 5 cm"))
   variable_levels <- c("<= 5 cm", "> 5 cm", "EM")
    readable_clinical_variable <- c("tumour size")
  } 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 == "COG_classifier"){
    readable_clinical_variable <- c("COG classifier")
    #ARMS sample have no low risk so need to remove that comparison and ERMS have no intermediate so need to remove that comparison
      if(all(str_detect(gsva_results$Histology, 'ARMS|ERMS'))){    
my_comparisons <- list( c("LOW RISK", "INTERMEDIATE RISK"), c("LOW RISK", "HIGH RISK") ,c("HIGH RISK", "INTERMEDIATE RISK"))
variable_levels <- c("LOW RISK", "INTERMEDIATE RISK", "HIGH RISK", "NA")}
   if(all(!str_detect(gsva_results$Histology, "ARMS"))){    
my_comparisons <- list( c("LOW RISK", "HIGH RISK"))
variable_levels <- c("LOW RISK","HIGH RISK", "NA")}
      if(all(!str_detect(gsva_results$Histology, "ERMS"))){
     my_comparisons <- list(c("HIGH RISK", "INTERMEDIATE RISK"))
variable_levels <- c("INTERMEDIATE RISK", "HIGH RISK", "NA")}

} else if(clinical_variable == "primary_site_status"){
my_comparisons <- list( c("Unfavourable", "Favourable"))
variable_levels <- c("Unfavourable", "Favourable", "NA")
readable_clinical_variable <- c("primary site status")
  } else if(clinical_variable == "first_event_occurred"){
my_comparisons <- list( c("Y", "N"))
variable_levels <- c("Y", "N")
readable_clinical_variable <- c("first event occurred")
  }else if(clinical_variable == "IRS_group"){
my_comparisons <- list( c("1", "2"),c("1", "3"),c("1", "4"),c("2", "3"),c("2", "4"),c("3", "4"))
variable_levels <- c("1", "2", "3", "4", "EM")
readable_clinical_variable <- c("IRS group")
  }else if(clinical_variable == "risk_group"){
    if(subtype=="ARMS"){
      my_comparisons <- list( c("high", "very high"))
      variable_levels <- c( "high", "very high", "EM")
    }else{
my_comparisons <- list( c("standard", "high"),c("high", "very high"),c("standard", "very high"))
variable_levels <- c("standard", "high", "very high", "EM")}
readable_clinical_variable <- c("risk group")

  }else if(clinical_variable == "complete_remission_achieved"){
    
my_comparisons <- list( c("Y", "N"))
variable_levels <- c("Y", "N", "NA")
readable_clinical_variable <- c("complete remission achieved")

  }else if(clinical_variable == "SIOP_stage"){
    
my_comparisons <- list( c("1", "2"),  c("1", "3"),  c("1", "4"),  c("2", "3"), c("2", "4"), c("3", "4"))
variable_levels <- c("1", "2","3", "4", "NA")
readable_clinical_variable <- c("SIOP stage")

  }
  
  
#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\n", readable_clinical_variable),
        x =paste0(toupper(substr(readable_clinical_variable, 1, 1)), substr(readable_clinical_variable, 2, nchar(readable_clinical_variable)), sep=""), y = "GSVA score")+ theme( axis.text.x = element_text(angle = 90,  size = 12))+ try(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,) 
    try(print(p))
try(ggsave(file= paste0(resistance_type, "_", clinical_variable,"_", filename, ".png"), width = 5, height = 5))
  
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, name="Translocation") + 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)"))+ scale_fill_discrete(name = "New Legend Title")
 ggsave(s, file= paste0(resistance_type, "_scatterplot_", clinical_variable,"_", subtype, ".png")) 
  #+ scale_color_discrete(name="Histology")
} 
}

```


#GSVA score plots of signatures
Get GSVA scores of gene signatures as plots using the function
```{r}

setwd("~/RMS microarray/Williamson/GSVA with full clinical annotation")

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

variables <- c("first_event_occurred", "metastasis", "IRS_group", "tumour_size", "risk_group", "complete_remission_achieved", "lymphnode", "Survival", "SIOP_stage")

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


#for ARMS acquired samples ---------------------------------------------
hub_genes <- read_xlsx("C:/Users/Christina/Documents/RMS RNA seq/Luciles data/RH4 WGCNA/RH4_acquired_res_signature.xlsx") #acquired rh4 sig
genes <- hub_genes$Gene

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


#for ERMS acquired samples -------------------------------
hub_genes <- read_xlsx("C:/Users/Christina/Documents/RMS RNA seq/Luciles data/RMSYM WGCNA/RMSYM_acquired_resistance_gene_signature.xlsx") #RMSYM acquired sig
genes <- hub_genes$Gene

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


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

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

#for intrinsic DEGs genes all sample ---------------------------------------

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

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

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

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

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

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

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

```
