---
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(recipes)
library(readr)
#library(tidymodels)
library(caret)
library(clusterProfiler)
library(rsample)
library(msigdb)
library(modeldata)
library(tidyverse)
library(sva)
library("fuzzyjoin")
library(ggpubr)
library(hgu133plus2.db)
library(affycoretools)
library(affy)
library(viridis)
library("ArrayExpress")
library(sva)
library(tidyverse)
library(PCAtools)
library(GSVA)
#library("arrayQualityMetrics")
library("readxl")
library(jetset)
library("affyio")
library(dplyr)
library(oligoClasses)
```


#clean pd
```{r}

pd_triche <- read.csv("~/RMS microarray/Triche/metadata_Triche.csv")
#for clinical data
pd_williamson <- read.csv("~/RMS microarray/Williamson/metadata_williamson.csv")

#rename cols going to merge so they are the same and change so match
#sex--------------------------------------
pd_triche$Sex

pd_triche$Sex <- gsub("female", "F", pd_triche$Sex)
pd_triche$Sex <- gsub("male", "M", pd_triche$Sex)

#age------------------------------
pd_williamson$Age.at.diagnosis..year.
pd_triche$Age
pd_williamson <- dplyr::rename(pd_williamson, Age = Age.at.diagnosis..year.)

#histology-------------------------------------------------
pd_williamson$Histology...Super.class
pd_triche$Histology

pd_williamson <- dplyr::rename(pd_williamson, Histology = Histology...Super.class)

pd_williamson$Histology <- gsub("RMS alveolar", "Alveolar", pd_williamson$Histology)
pd_williamson$Histology <- gsub("RMS embryonal", "Embryonal", pd_williamson$Histology)

pd_triche$Histology <- gsub("Mixed Alv/Emb", "Alveolar", pd_triche$Histology)

#translocation---------------------------------
pd_williamson$Fusion.transcript...gene.name
pd_triche$Translocation

pd_williamson <- dplyr::rename(pd_williamson, Translocation = Fusion.transcript...gene.name)

pd_williamson$Translocation <- gsub("Pax3", "P3F", pd_williamson$Translocation)
pd_williamson$Translocation <- gsub("Pax7", "P7F", pd_williamson$Translocation)
pd_williamson$Translocation <- gsub("Ncoa1P3F", "P3F", pd_williamson$Translocation)
pd_williamson$Translocation[is.na(pd_williamson$Translocation)] <- "NEG"

pd_triche$Translocation[is.na(pd_triche$Translocation )] <- "NEG"

#Fusion status--------------------------------
pd_williamson$Fusion_status <- pd_williamson$Translocation
pd_triche$Fusion_status <- pd_triche$Translocation

pd_williamson$Fusion_status <- gsub("P3F", "FP", pd_williamson$Fusion_status)
pd_williamson$Fusion_status <- gsub("P7F", "FP", pd_williamson$Fusion_status)
pd_williamson$Fusion_status <- gsub("NEG", "FN", pd_williamson$Fusion_status)

pd_triche$Fusion_status <- gsub("P3F", "FP",pd_triche$Fusion_status)
pd_triche$Fusion_status <- gsub("P7F", "FP", pd_triche$Fusion_status)
pd_triche$Fusion_status <- gsub("NEG", "FN", pd_triche$Fusion_status)


#survival-----------------------------------------------
pd_williamson$Factor.Value..survival_time.
pd_triche$Survival

pd_williamson <- dplyr::rename(pd_williamson, Survival = Factor.Value..survival_time.)

#size--------------------------------
pd_williamson$Tumor.size..recoded.
pd_triche$Size

pd_williamson <- dplyr::rename(pd_williamson, Size = Tumor.size..recoded.)

pd_williamson$Size <- gsub("EM", NA, pd_williamson$Size )
pd_williamson$Size <- gsub("EM", NA, pd_williamson$Size )

#node-------------------------------------------------------------
pd_williamson$Factor.Value..LymphNode.
pd_triche$Node

pd_williamson <- dplyr::rename(pd_williamson, Node = Factor.Value..LymphNode.)

pd_triche$Node <- gsub("N-0", "N0", pd_triche$Node )
pd_triche$Node <- gsub("N-1", "N1", pd_triche$Node )

#fail------------------------------------------------------
pd_williamson$First.event...occurred
pd_triche$Fail

pd_williamson <- dplyr::rename(pd_williamson, Fail = First.event...occurred)

pd_triche$Fail <- gsub("\\.", "Not_fail", pd_triche$Fail )
pd_williamson$Fail <- gsub("Y", "Fail", pd_williamson$Fail )
pd_williamson$Fail <- gsub("N", "Not_fail", pd_williamson$Fail )


#locati

```


#impute Triche 
```{r}


#read in triche
EM_triche <- read.csv("~/RMS microarray/Triche/EM_Triche.csv") %>% remove_rownames %>% column_to_rownames(var="X")


#get top 50 genes for imputation to reduce curse of dimentionality
# Get variance per gene
gene_variances <- apply(EM_triche, 1, var)

#Sort
top50_genes <- names(sort(gene_variances, decreasing = TRUE))[1:50]

# Subset expression matrix to top 50 variable genes
EM_triche <- EM_triche[top50_genes, ]

#only use informative clinical variables
pd_triche

# Select columns:
pd_triche <- pd_triche %>% dplyr::select(Sex, Age, Disease.Stage, Disease.Grade, Tumor.Location, Histology, Translocation, Survival, Fail, FFS, RG,  Size, Node, Invasion, geo_id, Fusion_status)

#find % of missing data 
Triche_important_cols <- pd_triche %>% dplyr::select(Age, Sex, Node, Size, Translocation, Fusion_status, Histology )
na_percent <- sum(is.na(Triche_important_cols)) / (nrow(Triche_important_cols) * ncol(Triche_important_cols)) * 100


#make sure rownames pd = colnames EM
rownames(pd_triche) <- pd_triche$geo_id
rownames(pd_triche) == colnames(EM_triche)

#join Triche EM and metadata-------------------------------------------------------
EM_triche <- EM_triche %>% t()
# pd_triche <- pd_triche %>% dplyr::select(Age, Histology, Translocation, Survival, Fail, Size, Node, Sex)
rownames(EM_triche) == rownames(pd_triche)

Triche_EM_metadata_merged <- merge(EM_triche,pd_triche, by = 0, all=TRUE )%>% remove_rownames %>% column_to_rownames(var="Row.names")

#impute Triche-------------------------------------------------------------------
rec_obj <- recipe(Fail ~ ., data = Triche_EM_metadata_merged)
rec_obj

#impute with
imputed <- rec_obj %>%
  step_impute_knn(all_predictors()) 
imputed

trained_rec <- prep(imputed, training = Triche_EM_metadata_merged)

train_data <- bake(trained_rec, new_data = Triche_EM_metadata_merged)

Triche_EM_metadata_imputed <- train_data %>% as.data.frame()

rownames(Triche_EM_metadata_imputed) <- rownames(Triche_EM_metadata_merged)

#remove NA fail
keep <- subset(Triche_EM_metadata_imputed, Fail != "Fail" | Fail != "Not_fail") %>% rownames()
Triche_EM_metadata_imputed <- Triche_EM_metadata_imputed[rownames(Triche_EM_metadata_imputed) %in% keep,]

write.csv(Triche_EM_metadata_imputed, "Triche_EM_metadata_imputed.csv")
```

#imputed Williamson 
```{r}

#read in williamson
EM_williamson <- read.csv("~/RMS microarray/Williamson/expression_data_williamson.csv") %>% remove_rownames %>% column_to_rownames(var="X")


#get top 50 genes for imputation to reduce curse of dimentionality
# Get variance per gene
gene_variances <- apply(EM_williamson, 1, var)

#Sort
top50_genes <- names(sort(gene_variances, decreasing = TRUE))[1:50]

# Subset expression matrix to top 50 variable genes
EM_williamson <- EM_williamson[top50_genes, ]

#only use informative clinical variables
pd_williamson

rownames(pd_williamson) <- paste0(pd_williamson$X.Oligo.cDNA.Array, ".CEL")

# Select columns:
pd_williamson <- pd_williamson %>% dplyr::select(Survival, Node, Sex, Histology, Location..primary.tumor., Metastasis...at.diagnosis..recoded., Stage...IRS.group, Stage...SIOP, Size, Translocation, Age, Fail, First.event...Type, Fusion_status, Risk..Group ,Event...First.Complete.remission...Achieved)


#find % of missing data 
Williamson_important_cols <- pd_williamson %>% dplyr::select(Age, Sex, Node, Size, Translocation, Fusion_status, Histology )
na_percent <- sum(is.na(Williamson_important_cols)) / (nrow(Williamson_important_cols) * ncol(Williamson_important_cols)) * 100

#make sure rownames pd = colnames EM

rownames(pd_williamson) == colnames(EM_williamson)

#join Williamson EM and metadata---------------------------------------------------
 EM_williamson <- EM_williamson %>% t()
# pd_williamson <- pd_williamson %>% dplyr::select(Age, Histology, Translocation, Survival, Fail, Size, Node, Sex)
# rownames(EM_williamson) == rownames(pd_williamson)

Williamson_EM_metadata_merged <- merge(EM_williamson,pd_williamson, by = 0, all=TRUE )%>% remove_rownames %>% column_to_rownames(var="Row.names")

rec_obj <- recipe(Fail ~ ., data = Williamson_EM_metadata_merged)
rec_obj

imputed <- rec_obj %>%
  step_impute_knn(all_predictors()) 
imputed

trained_rec <- prep(imputed, training = Williamson_EM_metadata_merged)

train_data <- bake(trained_rec, new_data = Williamson_EM_metadata_merged)

Williamson_EM_metadata_imputed <- train_data %>% as.data.frame()

rownames(Williamson_EM_metadata_imputed) <- rownames(Williamson_EM_metadata_merged)
  
write.csv(Williamson_EM_metadata_imputed, "Williamson_EM_metadata_imputed.csv")
```

#sva
```{r}

#merge EM----------------
#read in triche
EM_triche <- read.csv("~/RMS microarray/Triche/EM_Triche.csv") %>% remove_rownames %>% column_to_rownames(var="X")
#remove NA fail from triche
EM_triche <- EM_triche[,colnames(EM_triche)%in% keep]

#read in williamson
EM_williamson <- read.csv("~/RMS microarray/Williamson/expression_data_williamson.csv") %>% remove_rownames %>% column_to_rownames(var="X")

#merge data by common genes
# Make sure genes are rownames
common_genes <- intersect(rownames(EM_williamson), rownames(EM_triche))

# Subset
EM_williamson <- EM_williamson[common_genes, ]
EM_triche <- EM_triche[common_genes, ]

all(rownames(EM_triche) == rownames(EM_williamson))
EM_merged <- merge(EM_triche, EM_williamson, by = 0, all = TRUE)%>% remove_rownames %>% column_to_rownames(var="Row.names")


#merge metadata-----------------------------------------
#select importnant col
Williamson_EM_metadata_imputed <- Williamson_EM_metadata_imputed %>% dplyr::select(Age, Histology, Translocation, Survival, Fail, Size, Node, Sex, Fusion_status)

Triche_EM_metadata_imputed <- Triche_EM_metadata_imputed %>% dplyr::select(Age, Histology, Translocation, Survival, Fail, Size, Node, Sex, Fusion_status)

combined_metadata <- bind_rows(Triche_EM_metadata_imputed, Williamson_EM_metadata_imputed)

#check EM and metadata is in same order
combined_metadata <- combined_metadata[colnames(EM_merged), ]
all(rownames(combined_metadata) == colnames(EM_merged))

#convert metadata to fcators for sva
convert_char_to_factor <- function(df) {
  for (colname in names(df)) {
    if (is.character(df[[colname]])) {
      df[[colname]] <- as.factor(df[[colname]])
    }
  }
  return(df)
}

combined_metadata <- convert_char_to_factor(combined_metadata)

#add batch
combined_metadata$batch <- c(rep("Triche", ncol(EM_triche)), rep("Williamson", ncol(EM_williamson)))

#model matrix for variables you don't want to adjsut batch effect for 
mod <- model.matrix(~Fail + Age + Histology + Translocation + Survival + Size + Node + Sex, data = combined_metadata)

#need to t() EM so genes are rows
#EM_merged <- EM_merged %>% t()

# Apply ComBat to remove batch effects
combat_expr <- ComBat(dat = EM_merged,
                      batch = combined_metadata$batch,
                      mod = mod, 
                      par.prior = TRUE, 
                      prior.plots = FALSE)

write.csv(combat_expr, "combat_EM.csv")

```


#PCA
Check sva hasn't removed too much variation
```{r}

combined_metadata$Event <- ifelse(combined_metadata$Fail == "Fail", "Event", "No event" )

EM_merged

#combat_expr <- combat_expr %>% t()
p <- pca(combat_expr, metadata = combined_metadata, removeVar = 0.1)

  #reorder levels for plots
 # p$metadata$`age:ch1` <- factor(p$metadata$`age:ch1`) %>% order(p$metadata$`age:ch1`)
 # 
 # metadata$Time<- factor(metadata$Time,  c("PT3", "PT7", "PT10"))
  
#biplots------------------------------------------------------------------
#   
#   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.", "Factor.Value..COGClassifier." )
#   
#    variables <- c("Factor.Value..COGClassifier.")
#   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",j,".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(i),
#     colkey = myCols,
#     hline = 0, vline = 0,
#     legendPosition = 'right', 
#     ellipse = FALSE,
#       ellipseType = 't',
#       ellipseLevel = 0.95,
#       ellipseFill = TRUE,
#       ellipseAlpha = 1/4,
#     colLegendTitle = "COG classifier",
#       ellipseLineSize = 1.0,
#  #  shape = 'histology:ch1',
#  #  shapeLegendTitle = "Histology",
#     pointSize = 5,
#     title = 'PCA bi-plot',
#     subtitle = 'PC1 versus PC2',)
#  print(g)
# ggsave(g, file = filename, width = 10, height = 6, dpi = 300)
#   }
  
  
   #clinical variable=colname(pd), readable_variable=name of legend title
   
make_pca_biplot <- function(clinical_variable, plot_title, legend_title, filename){
  
  #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 <- legend_title
  #make filename title
   filename <- paste0("PCA_biplot_",filename,".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 = legend_title,
      ellipseLineSize = 1.0,
   shape = "Translocation",
   shapeLegendTitle = "Translocation",
    pointSize = 5,
    title = plot_title,
    subtitle = 'PC1 versus PC2',)
 print(g)
ggsave(g, file = filename, width = 7, height = 4.5, dpi = 100)
}

make_pca_biplot("batch", "Before batch correction" , "Batch", "batch_before_correction")
  
   
   make_pca_biplot("Fail", "Fail" )
   
   make_pca_biplot("Translocation",  "After batch correction" ,"Translocation","translocation_after_correction" )
   
   make_pca_biplot("batch", "After batch correction" , "Batch", "batch_after_correction")
   
   make_pca_biplot("Event",  "After batch correction" ,"Event","Event_after_correction" )
   
    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" )
  
  #look at the samples with mixed hsitology 
      #label <- subset(p$metadata, `histology:ch1` == "Mixed Alv/Emb") %>% rownames() #label only mixed samples
    
  #  x <- p$metadata[row.names(p$metadata) %in% label,]
  # 
  #    ERMS <- c("GSM2435372", "GSM2435277")
  # 
  # x <- p$metadata[rownames(p$metadata)%in% ERMS,]
  # 
  # write.csv(p$metadata, "metadata.csv")
  
 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)
    }
  
```


#merge EM and metadata
```{r}

combat_expr <- combat_expr %>% t()

combined_EM_metadata <- merge(combat_expr, combined_metadata, by = 0 , all=TRUE)%>% remove_rownames %>% column_to_rownames(var="Row.names")

combined_EM_metadata <- combined_EM_metadata %>% dplyr::select(-Event)

write.csv(combined_EM_metadata, "combat_combined_EM_metadata.csv")
```

#get metadata characteristics
```{r}

metadata <- combined_metadata %>% dplyr::select(-Fail, -Survival)

meta_long <- metadata  %>%
  mutate(across(everything(), as.character)) %>%
  pivot_longer(cols = everything(), names_to = "Variable", values_to = "Category") %>%
  group_by(Variable, Category) %>%
  summarise(Count = n(), .groups = "drop") %>%
  group_by(Variable) %>%
  mutate(Percentage = round(100 * Count / sum(Count), 1)) %>%
  ungroup()



# Step 1: Convert Category to numeric and group by age ranges
summary_age_grouped <- meta_long %>%
  filter(Variable == "Age") %>%
  mutate(
    Age = as.numeric(Category),
    Age_Group = cut(
      Age,
      breaks = c(-0.1, 4, 9, 14, 19, 24),  # adjust upper bound to cover up to 21.6
      labels = c("0–4", "5–9", "10–14", "15–19", "20–24"),
      right = TRUE
    )
  ) %>%
  group_by(Variable = "Age", Category = Age_Group) %>%
  summarise(
    Count = sum(Count),
    Percentage = sum(Percentage),
    .groups = "drop"
  )


summary_table_grouped <- meta_long %>%
  filter(Variable != "Age") %>%
  bind_rows(summary_age_grouped)


library(flextable)

table_flex <-library(flextable)

table_flex <- summary_table_grouped %>%
  mutate(
    Percentage = paste0(round(Percentage, 1), "%")  # Format Percentage
  ) %>%
  flextable() %>%
  set_header_labels(
    Variable = "Variable",
    Category = "Category",
    Count = "N",
    Percentage = "%"
  ) %>%
  autofit()


library(officer)

doc <- read_docx() %>%
  body_add_flextable(value = table_flex)
print(doc, target = "age_summary.docx")

write.csv(summary_table_grouped, "summary_table_grouped.csv")

```

#split data
split into train and test. The function createDataPartition can be used to create balanced splits of the data. If the y argument to this function is a factor, the random sampling occurs within each class and should preserve the overall class distribution of the data. For example, to create a single 80/20% split:
```{r}

set.seed(285)
 trainIndex <- createDataPartition(combined_EM_metadata$Fail, p = .8, 
                                  list = TRUE,
                                  times = 100)
head(trainIndex)

for  (name in names(trainIndex)){
  
combined_EM_metadata_Train <- combined_EM_metadata[ trainIndex[[name]],]
combined_EM_metadata_Test  <- combined_EM_metadata[-trainIndex[[name]],]

write.csv(combined_EM_metadata_Train, paste0(name, "_Train.csv"))


write.csv(combined_EM_metadata_Test, paste0(name, "_Test.csv"))
}





for (i in 1:100) {
  # Create partition
  trainIndex <- createDataPartition(combined_EM_metadata$Fail, p = 0.8, list = FALSE)
  
  # Subset data
  train_data <- combined_EM_metadata[trainIndex, ]
  test_data  <- combined_EM_metadata[-trainIndex, ]
  
  # Save to CSV files (change path if needed)
  train_filename <- paste0("combined_EM_metadata_Train_", i, ".csv")
  test_filename  <- paste0("combined_EM_metadata_Test_", i, ".csv")
  
  write.csv(train_data, train_filename, row.names = FALSE)
  write.csv(test_data, test_filename, row.names = FALSE)
}

```

#Logistic regression
```{r}

combined_EM_metadata_Train <- read.csv("combined_EM_metadata_Train.csv")%>% remove_rownames %>% column_to_rownames(var="X")


combined_EM_metadata_Train <- combined_EM_metadata_Train %>% dplyr::select(-Survival)

#add factors to all categorical variables
combined_EM_metadata_Train[] <- lapply(combined_EM_metadata_Train, function(x) if(is.character(x)) factor(x) else x)

levels(combined_EM_metadata_Train$Fail)

#for logistic regression, fail is second as by default, in a binary factor, the first factor level is the reference level (coded as 0), and the second level is the "event" (coded as 1)
combined_EM_metadata_Train$Fail <- fct_relevel(combined_EM_metadata_Train$Fail, "Not_fail", "Fail")

  #run univariate logistic regression for 
predictors <- setdiff(names(combined_EM_metadata_Train), "Fail")

# Initialize a vector to store p-values
pvals <- numeric(length(predictors))
names(pvals) <- predictors

odds_ratio <- numeric(length(predictors))
names(odds_ratio) <- predictors

CI_lower <- numeric(length(predictors))
names(CI_lower) <- predictors

CI_upper <- numeric(length(predictors))
names(CI_upper) <- predictors

# Run logistic regression for each predictor
for (var in predictors) {
  
 formula <- as.formula(paste("Fail ~ `", var, "`", sep = ""))

  # Fit the model
  fit <- glm(formula, data = combined_EM_metadata_Train, family = "binomial")
  
  pvals[var] <- summary(fit)$coefficients[2, 4]  # Extract the p-value
  #coef[var] <- unname(coef(fit)[2])  # Extract the p-value
  #coef <- coef(fit)[2]  # get the coefficient for the feature
  odds_ratio[var] <- exp(unname(coef(fit)[2]))
  
CI_lower[var] <- exp(confint(fit)[2,1 ])
CI_upper[var] <- exp(confint(fit)[2,2 ])

}

# View the p-values
pvals


pvals_adj <- p.adjust(pvals, method = "BH")  # FDR adjustment

results <- data.frame(
  variable = names(pvals),
  p_value = pvals,
   OddsRatio = odds_ratio,
  p_adj = pvals_adj,
  CI_lower = CI_lower, 
  CI_upper= CI_upper
)

# Sort by adjusted p-value
results <- results[order(results$p_adj), ]

#write.csv(results, "LR_results.csv")

# View significant results (e.g., FDR < 0.01)
significant <- subset(results, p_value < 0.01)

#find expression values
combined_EM_metadata_Train_t <- combined_EM_metadata_Train[,significant$variable]

range(unlist(combined_EM_metadata_Train_t[sapply(combined_EM_metadata_Train_t, is.numeric)]), na.rm = TRUE)
mean(unlist(combined_EM_metadata_Train_t[sapply(combined_EM_metadata_Train_t, is.numeric)]), na.rm = TRUE)




# Top 10 with highest ORs
top_highest <- significant[order(-significant$OddsRatio), ][1:10, ]

# Top 10 with lowest ORs
top_lowest <- significant[order(significant$OddsRatio), ][1:10, ]

# Combine
top_combined <- rbind(top_highest, top_lowest)

top_combined$variable <- factor(top_combined$variable, levels = top_combined$variable[order(top_combined$OddsRatio)])

# Log-transform OR
top_combined$log_OR <- log(top_combined$OddsRatio)


significant$variable <- factor(significant$variable, levels = significant$variable[order(significant$OddsRatio)])


# HR plot no padj colour
ggplot(significant, aes(x = variable, y = OddsRatio)) +
  geom_point(size = 3) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "red")+
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2) +
  coord_flip() +
  labs(
    title = "Hazard Ratios from Cox Model",
    x = "Gene",
    y = "Hazard Ratio (95% CI)"
  ) +
  theme_minimal()

#HR plot colour y padjust
OR_plot <- ggplot(top_combined, aes(x = variable, y = log_OR, color = p_value )) +
  geom_point(size = 3) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "red")+
  geom_errorbar(aes(ymin = log(CI_lower), ymax = log(CI_upper)), width = 0.2) +
  coord_flip() +
  labs(
    title = "Odds Ratios from logistic regression model",
    x = "Variable",
    y = "Odds Ratio (95% CI)",
    color = "P.adjust", size=9
  ) +
  theme_minimal()+
  scale_color_gradient(low = "blue", high = "red")+theme(plot.title = element_text(size = 12, face = "bold"),
    legend.title=element_text(size=10), 
    legend.text=element_text(size=10), 
    axis.text.x=element_text(size=10),
    axis.text.y=element_text(size=10),)
ggsave(OR_plot, filename = "OR_plot.png", bg= "white", height = 3.5, width=6)



```

#ORA on logistic regression results
```{r}

#select for just OR > 1
results_df <- significant[significant$OddsRatio > 1,]
geneList <- results_df$variable
filename <- c("ORA_increased")

#select for just OR < 1
results_df <- significant[significant$OddsRatio < 1,]
geneList <- results_df$variable
filename <- c("ORA_decreased")



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

write.csv(ego2, file=paste0("GO_", filename, ".csv"))

GO <- try(barplot(ego2,  font.size = 16 ,showCategory=10, x="Count")) + scale_y_discrete(labels=function(ego2) str_wrap(ego2, width=60))
try(ggsave(GO, file=paste0("GO_", filename, ".png"), width = 20, height = 10, units = "cm"))

#KEGG
# Convert gene symbols to Entrez IDs
gene_df <- bitr(geneList,
                fromType = "SYMBOL",
                toType = "ENTREZID",
                OrgDb = org.Hs.eg.db)

gene_df <- gene_df$ENTREZID

# Check conversion
head(gene_df)
kk <- enrichKEGG(gene         = gene_df,
                 organism     = 'hsa',
                 pvalueCutoff = 0.5, 
                 keyType = "ncbi-geneid")
head(kk)

```

#Filter EM
```{r}

#filter clinical and EM to just include genes significantly correlated with surival to reduce the dimentionality
combined_EM_metadata_Train <- combined_EM_metadata_Train[,colnames(combined_EM_metadata_Train) %in% significant$variable] 

# #read in data
clinical <- read.csv("combined_EM_metadata_Train.csv") %>% remove_rownames %>% column_to_rownames(var="X")

rownames(combined_EM_metadata_Train) == rownames(clinical)

combined_EM_metadata_Train$Fail <- clinical$Fail

#add factors 
combined_EM_metadata_Train$Fail <- factor(combined_EM_metadata_Train$Fail)

#for caret, factors need to be positive first and negative second
combined_EM_metadata_Train$Fail <- fct_relevel(combined_EM_metadata_Train$Fail, "Fail", "Not_fail")

#clinical_and_EM_filtered <- Triche_EM_metadata_merged_filtered
```


#RFE
Using recursive feature elimination to reduce the number of features 
```{r}

set.seed(285)
# Set up control for rfe using random forest function and cross validation 
ctrl <- rfeControl(functions=rfFuncs, method="cv", number=5)  # 5-fold cross-validation

# Run RFE to select the best number of features
n_features <- ncol(combined_EM_metadata_Train) - 1 
set.seed(285)
rfe_results <- rfe(combined_EM_metadata_Train[,1:n_features], combined_EM_metadata_Train$Fail, sizes=seq(1, n_features, by = 1), rfeControl=ctrl)

# View the results of the feature selection
print(rfe_results)
saveRDS(rfe_results, "rfe_results.rds")

rfe_results_df <- rfe_results$results %>% as.data.frame()

write.csv(rfe_results_df, "rfe_results.csv")
write.csv(rfe_results$variables$var, "rfe_variables.csv")

# Plot the results
png(filename=paste0("rfe_plot_new.png"))
plot(rfe_results, type = c("g", "o"))
dev.off()

#get features from optimum model with highest accuracy
predictors(rfe_results)

#get top features
selected_features <- rfe_results$optVariables[1:45]

write.csv(selected_features, "selected_features.csv")

# Subset the dataset to keep only the selected features
subset_data <- combined_EM_metadata_Train[, selected_features]

#add fail back in to df 
subset_data$Fail <- combined_EM_metadata_Train$Fail

# Check the structure of the subset
str(subset_data)

#relevel so first level is event
subset_data$Fail <- fct_relevel(subset_data$Fail, "Fail", "Not_fail")



```

#RF------------------------------------------
#RF standard sampling
```{r}


fitControl <- trainControl(## 5-fold CV
                           method = "repeatedcv",
                           number = 5,
                           ## repeated ten times
                           repeats = 10,
                           savePredictions = T, 
                           classProbs = T,
                           summaryFunction = twoClassSummary 
                          )



gbmGrid <-  expand.grid(mtry = c(5, 10, 15, 20, 25))



#don't need to preprocess as dont need to normalise for decision trees and xgboost
set.seed(825)

rf_standard <- train(Fail ~ ., data = subset_data, 
                 method = 'rf', 
                 trControl = fitControl,
                  tuneGrid = gbmGrid,
                   metric="ROC"
                 )
rf_standard


#plot results of tuning grid
trellis.par.set(caretTheme())

plot(rf_standard)


roc_obj <- roc(response = final_model$pred$obs, predictor = final_model$pred$Fail, levels = c("Fail","Not_fail"))
plot(roc_obj)


```


#RF down sampling
```{r}

fitControl <- trainControl(## 10-fold CV
                           method = "repeatedcv",
                           number = 5,
                           ## repeated ten times
                           repeats = 10,
                           savePredictions = T, 
                           classProbs = T,
                           sampling = "down",
                            summaryFunction = twoClassSummary 
                          )


gbmGrid <-  expand.grid(mtry = c(5, 10, 15, 20, 25))


#don't need to preprocess as dont need to normalise for decision trees and xgboost
set.seed(825)

rf_down <- train(Fail ~ ., data = subset_data, 
                 method = 'rf', 
                 trControl = fitControl,
                  tuneGrid = gbmGrid,
                  metric="ROC")
rf_down


#plot results of tuning grid
trellis.par.set(caretTheme())
plot(rf_down)  


#roc_obj <- roc(response = final_model$pred$obs, predictor = final_model$pred$Fail, levels = c("Fail","Not_fail"))
#plot(roc_obj)


```

#RF up sampling
```{r}

fitControl <- trainControl(## 10-fold CV
                           method = "repeatedcv",
                           number = 5,
                           ## repeated ten times
                           repeats = 10,
                           savePredictions = T, 
                           classProbs = T,
                           sampling = "up",
                            summaryFunction = twoClassSummary 
                          )



gbmGrid <-  expand.grid(mtry = c(5, 10, 15, 20, 25))



#don't need to preprocess as dont need to normalise for decision trees and xgboost
set.seed(825)

rf_up <- train(Fail ~ ., data = subset_data, 
                 method = 'rf', 
                 trControl = fitControl,
                  tuneGrid = gbmGrid,
                metric="ROC")
rf_up

#plot results of tuning grid
trellis.par.set(caretTheme())
plot(rf_up)  


#plot results of tuning grid
trellis.par.set(caretTheme())
plot(rf_up, cex.axis=2)  

```

#compare RF
```{r}

#plot ROC curve of rf for different sampling methods
library(pROC)
# Prevent scientific notation globally for this plotting session
options(scipen = 999)


roc_rose <- plot(roc(response = rf_standard$pred$obs, predictor = rf_standard$pred$Fail, levels = c("Not_fail","Fail")), print.auc = TRUE, col = "#0072B2",legacy.axes = TRUE)

roc_rose <- plot(roc(response = rf_down$pred$obs, predictor = rf_down$pred$Fail, levels = c("Not_fail","Fail")), print.auc = TRUE,legacy.axes = TRUE ,
                 col = "#F0E442", print.auc.y = .4, add = TRUE)#

roc_rose <- plot(roc(response = rf_up$pred$obs, predictor = rf_up$pred$Fail, levels = c("Not_fail","Fail")), print.auc = TRUE, 
                 col = "#009E73", print.auc.y = .6,legacy.axes = TRUE, add = TRUE)
legend(x = "topleft", legend=c("Standard sampling", "Downsampling", "Upsampling"), 
       fill = c("#0072B2","#F0E442" ,"#009E73"))



#plot accurcay and kappa
resamps <- resamples(list(Standard = rf_standard,
                          Downsampling = rf_down, 
                          Upsampling = rf_up))
resamps


theme1 <- trellis.par.get()
theme1$plot.symbol$col = rgb(.2, .2, .2, .4)
theme1$plot.symbol$pch = 16
theme1$plot.line$col = rgb(1, 0, 0, .7)
theme1$plot.line$lwd <- 2
trellis.par.set(theme1)
bwplot(resamps, layout = c(2, 1), metric = c("Sens", "Spec"))


bwplot(resamps, layout = c(3, 1))
```



#Test model
```{r}

combined_EM_metadata_Test <- read.csv("combined_EM_metadata_Test.csv")%>% remove_rownames %>% column_to_rownames(var="X")

final_model <- rf_standard
final_model <- rf_down
final_model <- rf_up

#subset for features in model
combined_EM_metadata_Test_selected <- combined_EM_metadata_Test[,selected_features]

#check they match
colnames(combined_EM_metadata_Test_selected) == colnames(subset_data)

#relevel factor so event is second
combined_EM_metadata_Test$Fail <- fct_relevel(combined_EM_metadata_Test$Fail, "Fail", "Not_fail")

# Make predictions using the trained model on the new dataset
predictions <- predict(final_model, newdata = combined_EM_metadata_Test_selected)

true_labels <- combined_EM_metadata_Test$Fail
conf_matrix <- confusionMatrix(predictions, true_labels,positive = "Fail")
print(conf_matrix)


# Example of adjusting threshold
predictions_prob <- predict(final_model, newdata = combined_EM_metadata_Test_selected, type = "prob")
# 
# library(pROC)
# 
# # Assuming `target` is your actual class and `probs$event` are predicted probabilities
 roc_obj <- roc(response = true_labels, predictor = predictions_prob$Fail, levels = c("Not_fail","Fail"))
# #plot(roc_obj, legacy.axes = TRUE)

# Best threshold based on Youden’s J (sensitivity + specificity - 1)
best_thresh <- coords(roc_obj, "best", ret = "threshold")
best_thresh <- as.numeric(best_thresh)  # ensure it's numeric

# Assume we are predicting the positive class (e.g., the second column in the result)
predictions <- ifelse(predictions_prob$Fail > best_thresh, "Fail", "Not_fail")


predictions<- factor(predictions)
conf_matrix <- confusionMatrix(predictions, true_labels,  positive = "Fail")
print(conf_matrix)




#plot ROC curve of rf for different sampling methods
library(pROC)
roc_rose <- plot(roc(response = true_labels, predictor = predictions_prob$Fail, levels = c("Not_fail","Fail")), print.auc = TRUE, col = "#0072B2",legacy.axes = TRUE)

roc_rose <- plot(roc(response = true_labels, predictor = predictions_prob$Fail, levels = c("Not_fail","Fail")), print.auc = TRUE,legacy.axes = TRUE ,
                 col = "#F0E442", print.auc.y = .4, add = TRUE)#

roc_rose <- plot(roc(response = true_labels, predictor = predictions_prob$Fail, levels = c("Not_fail","Fail")), print.auc = TRUE, 
                 col = "#009E73", print.auc.y = .6,legacy.axes = TRUE, add = TRUE)
legend(x = "topleft", legend=c("Standard sampling", "Downsampling", "Upsampling"), 
       fill = c("#0072B2","#F0E442" ,"#009E73"))
```

#XGBoost------------------------------

#Train model standard
```{r}

preprocess = c("center", "scale")

set.seed(825)


ctrl <- trainControl(
  method = "repeatedcv",
  number = 5,
  repeats = 10,
  classProbs = TRUE,
  summaryFunction = twoClassSummary,
                           savePredictions = T
)


tune_grid <- expand.grid(
  nrounds = c(50, 100),       # number of boosting rounds
  max_depth = c(2, 3, 4),          # shallower trees to prevent overfitting
  eta = c(0.01, 0.05),        # conservative learning rates, default 0.3
  gamma = c(1, 5),                 # regularisation
  colsample_bytree = c(0.3, 0.5),  # % of features used per tree
  min_child_weight = c(5, 10),     # minimum number of instances in leaf node
  subsample = c(0.5, 0.7)          # % of samples used per boosting round
)



xgb <- train(
  Fail ~ .,
  data = subset_data,
  method = "xgbTree",
  trControl = ctrl,
  tuneGrid = tune_grid , 
  preProcess = preprocess, 
  metric = "ROC",
  lambda=1
)


```


#Training upsampling
```{r}

# i <- c("univariate_cox_p.adj_0.05")
#i <- c("univariate_cox_p.val_0.001")
#i <- c("univariate_cox_p.val_0.01")
# i <- c("highly_variable")
# i <- c("univariate_cox_p.val_0.05")


preprocess = c("center", "scale")


set.seed(825)

tune_grid <- expand.grid(
  nrounds = c(50, 100),       # number of boosting rounds
  max_depth = c(2, 3, 4),          # shallower trees to prevent overfitting
  eta = c(0.01, 0.05),        # conservative learning rates
  gamma = c(1, 5),                 # regularization: min loss reduction
  colsample_bytree = c(0.3, 0.5),  # % of genes (features) used per tree
  min_child_weight = c(5, 10),     # minimum number of instances in leaf node
  subsample = c(0.5, 0.7)          # % of samples used per boosting round
)


#control with upsampling
ctrl_up <- trainControl(
  method = "repeatedcv",
  number = 5,
  repeats = 10,
  classProbs = TRUE,
  summaryFunction = twoClassSummary,
  sampling = "up",
  verboseIter = TRUE,
  savePredictions = T
)



# Train with upsampling
xgb_up <- train(
  Fail ~ .,
  data = subset_data,
  method = "xgbTree",
  trControl = ctrl_up,
  metric = "ROC",
  tuneGrid = tune_grid ,
  preProcess = preprocess,
  lambda = 1
)





```

#Training down
```{r}

# i <- c("univariate_cox_p.adj_0.05")
#i <- c("univariate_cox_p.val_0.001")
#i <- c("univariate_cox_p.val_0.01")
# i <- c("highly_variable")
# i <- c("univariate_cox_p.val_0.05")

preprocess = c("center", "scale")



set.seed(825)

tune_grid <- expand.grid(
  nrounds = c(50, 100),       # number of boosting rounds
  max_depth = c(2, 3, 4),          # shallower trees to prevent overfitting
  eta = c(0.01, 0.05),        # conservative learning rates
  gamma = c(1, 5),                 # regularization: min loss reduction
  colsample_bytree = c(0.3, 0.5),  # % of genes (features) used per tree
  min_child_weight = c(5, 10),     # minimum number of instances in leaf node
  subsample = c(0.5, 0.7)          # % of samples used per boosting round
)



#control with upsampling
ctrl_down <- trainControl(
  method = "repeatedcv",
  number = 5,
  repeats = 10,
  classProbs = TRUE,
  summaryFunction = twoClassSummary,
  sampling = "down",
  verboseIter = TRUE,
  savePredictions = T
)



# Train with upsampling
xgb_down <- train(
  Fail ~ .,
  data = subset_data,
  method = "xgbTree",
  trControl = ctrl_down,
  metric = "ROC",
  tuneGrid = tune_grid ,
  preProcess = preprocess,
  lambda = 1
)





```



#compare XGBoost
```{r}

#plot ROC curve of rf for different sampling methods
library(pROC)
roc_rose <- plot(roc(response = xgb$pred$obs, predictor = xgb$pred$Fail, levels = c("Not_fail","Fail")), print.auc = TRUE, col = "#E69F00",legacy.axes = TRUE )

roc_rose <- plot(roc(response = xgb_down$pred$obs, predictor = xgb_down$pred$Fail, levels = c("Not_fail","Fail")), print.auc = TRUE, 
                 col = "#009E73", print.auc.y = .4, add = TRUE,legacy.axes = TRUE )#

roc_rose <- plot(roc(response = xgb_up$pred$obs, predictor = xgb_up$pred$Fail, levels = c("Not_fail","Fail")), print.auc = TRUE, 
                 col = "#CC79A7", print.auc.y = .6, add = TRUE,legacy.axes = TRUE )
legend(x = "topleft", legend=c("Standard", "Downsampling", "Upsampling"), 
       fill = c("#E69F00","#009E73" ,"#CC79A7"))


#plot accurcay and kappa
resamps <- resamples(list(Standard = xgb,
                          Downsampling = xgb_down, 
                          Upsampling = xgb_up))
resamps

theme1 <- trellis.par.get()
theme1$plot.symbol$col = rgb(.2, .2, .2, .4)
theme1$plot.symbol$pch = 16
theme1$plot.line$col = rgb(1, 0, 0, .7)
theme1$plot.line$lwd <- 2
trellis.par.set(theme1)
bwplot(resamps, layout = c(2, 1), metric = c("Sens", "Spec"))

```



#Test model
```{r}

model <- xgb
model <- xgb_down
model <- xgb_up

#subset for features in model
combined_EM_metadata_Test_selected <- combined_EM_metadata_Test[,selected_features]

#check they match
colnames(combined_EM_metadata_Test_selected) == colnames(subset_data)

#relevel factor so event is second
combined_EM_metadata_Test$Fail <- fct_relevel(combined_EM_metadata_Test$Fail, "Fail", "Not_fail")

# Make predictions using the trained model on the new dataset
predictions <- predict(model, newdata = combined_EM_metadata_Test_selected)

true_labels <- combined_EM_metadata_Test$Fail
conf_matrix <- confusionMatrix(predictions, true_labels,positive = "Fail")
print(conf_matrix)


# Example of adjusting threshold
predictions_prob <- predict(model, newdata = combined_EM_metadata_Test_selected, type = "prob")

library(pROC)

# Assuming `target` is your actual class and `probs$event` are predicted probabilities
roc_obj <- roc(response = true_labels, predictor = predictions_prob$Fail, levels = c("Fail", "Not_fail"))
#plot(roc_obj)

# Best threshold based on Youden’s J (sensitivity + specificity - 1)
best_thresh <- coords(roc_obj, "best", ret = "threshold")
best_thresh <- as.numeric(best_thresh)  # ensure it's numeric

# Assume we are predicting the positive class (e.g., the second column in the result)
predictions <- ifelse(predictions_prob$Fail > best_thresh, "Fail", "Not_fail")


predictions<- factor(predictions)
conf_matrix <- confusionMatrix(predictions, true_labels,  positive = "Fail")
print(conf_matrix)



#plot ROC curve of rf for different sampling methods
library(pROC)
roc_rose <- plot(roc(response = true_labels, predictor = predictions_prob$Fail, levels = c("Not_fail","Fail")), print.auc = TRUE, col = "#E69F00",legacy.axes = TRUE )

roc_rose <- plot(roc(response = true_labels, predictor = predictions_prob$Fail, levels = c("Not_fail","Fail")), print.auc = TRUE, 
                 col = "#009E73", print.auc.y = .4, add = TRUE,legacy.axes = TRUE )#

roc_rose <- plot(roc(response = true_labels, predictor = predictions_prob$Fail, levels = c("Not_fail","Fail")), print.auc = TRUE, 
                 col = "#CC79A7", print.auc.y = .6, add = TRUE,legacy.axes = TRUE )
legend(x = "topleft", legend=c("Standard", "Downsampling", "Upsampling"), 
       fill = c("#E69F00","#009E73" ,"#CC79A7"))

```










#variable importance
```{r}

set.seed(825)
#variable importance
#png(filename=paste0("~/RMS microarray/Triche/ML/Variable_importance/", i, ".png"))
gbmImp <- varImp(xgb_up, scale = TRUE)
plot(gbmImp, top = 20)
#dev.off()
#gbmImp <- gbmImp$importance

#write.csv(gbmImp$importance, "important_features.csv")



keep <- subset(gbmImp$importance, Overall != 0) %>% rownames()

# Subset the dataset to keep only the selected features
subset_data <- combined_EM_metadata_Train[, keep]

#add fail back in to df 
subset_data$Fail <- combined_EM_metadata_Train$Fail

# Check the structure of the subset
str(subset_data)

#relevel so first level is event
subset_data$Fail <- fct_relevel(subset_data$Fail, "Fail", "Not_fail")



preprocess = c("center", "scale")

set.seed(825)


ctrl <- trainControl(
  method = "repeatedcv",
  number = 5,
  repeats = 10,
  classProbs = TRUE,
  summaryFunction = twoClassSummary,
                           savePredictions = T
)


tune_grid <- expand.grid(
  nrounds = c(50, 100),       # number of boosting rounds
  max_depth = c(2, 3, 4),          # shallower trees to prevent overfitting
  eta = c(0.01, 0.05),        # conservative learning rates, default 0.3
  gamma = c(1, 5),                 # regularisation
  colsample_bytree = c(0.3, 0.5),  # % of features used per tree
  min_child_weight = c(5, 10),     # minimum number of instances in leaf node
  subsample = c(0.5, 0.7)          # % of samples used per boosting round
)



xgb_important_features <- train(
  Fail ~ .,
  data = subset_data,
  method = "xgbTree",
  trControl = ctrl,
  tuneGrid = tune_grid , 
  preProcess = preprocess, 
  metric = "ROC",
  lambda=1
)


library(pROC)
roc_rose <- plot(roc(response = xgb_important_features$pred$obs, predictor = xgb_important_features$pred$Fail, levels = c("Fail","Not_fail")), print.auc = TRUE, col = "#E69F00",legacy.axes = TRUE )



model <- xgb_important_features

#subset for features in model
combined_EM_metadata_Test_selected <- combined_EM_metadata_Test[,keep]

#check they match
colnames(combined_EM_metadata_Test_selected) == colnames(subset_data)

#relevel factor so event is second
combined_EM_metadata_Test$Fail <- fct_relevel(combined_EM_metadata_Test$Fail, "Fail", "Not_fail")

# Make predictions using the trained model on the new dataset
predictions <- predict(model, newdata = combined_EM_metadata_Test_selected)

true_labels <- combined_EM_metadata_Test$Fail
conf_matrix <- confusionMatrix(predictions, true_labels,positive = "Fail")
print(conf_matrix)


# Example of adjusting threshold
predictions_prob <- predict(model, newdata = combined_EM_metadata_Test_selected, type = "prob")

library(pROC)

# Assuming `target` is your actual class and `probs$event` are predicted probabilities
roc_obj <- roc(response = true_labels, predictor = predictions_prob$Fail, levels = c("Fail", "Not_fail"))
#plot(roc_obj)

# Best threshold based on Youden’s J (sensitivity + specificity - 1)
best_thresh <- coords(roc_obj, "best", ret = "threshold")
print(best_thresh)

#0.4108368, 0.4298175, 0.2695673
# Assume we are predicting the positive class (e.g., the second column in the result)
predictions <- ifelse(predictions_prob[,1] >0.4579119,"Fail", "Not_fail")  # Adjust threshold to 0.6

predictions<- factor(predictions)
conf_matrix <- confusionMatrix(predictions, true_labels,  positive = "Fail")
print(conf_matrix)



#plot ROC curve of rf for different sampling methods
library(pROC)
roc_rose <- plot(roc(response = true_labels, predictor = predictions_prob$Fail, levels = c("Fail","Not_fail")), print.auc = TRUE, col = "#E69F00",legacy.axes = TRUE )

roc_rose <- plot(roc(response = true_labels, predictor = predictions_prob$Fail, levels = c("Fail","Not_fail")), print.auc = TRUE, 
                 col = "#009E73", print.auc.y = .4, add = TRUE,legacy.axes = TRUE )#

roc_rose <- plot(roc(response = true_labels, predictor = predictions_prob$Fail, levels = c("Fail","Not_fail")), print.auc = TRUE, 
                 col = "#CC79A7", print.auc.y = .6, add = TRUE,legacy.axes = TRUE )
legend(x = "topleft", legend=c("Standard", "Downsampling", "Upsampling"), 
       fill = c("#E69F00","#009E73" ,"#CC79A7"))
```

#top features LR plot
```{r}
selected_features <- read.csv("important_features.csv")

#get top 20 features
top_20 <- selected_features$X[1:20]

#read in LR results
survival_results <- read.csv("LR_results.csv")%>% remove_rownames %>% column_to_rownames(var="X")

#subset from lR results
LR_results <- survival_results[survival_results$variable %in% top_20,]

#change to factor
LR_results$variable <- factor(LR_results$variable, levels = LR_results$variable[order(LR_results$OddsRatio)])

#plot to see log regression
#HR plot colour by padjust
OR_plot <- ggplot(LR_results, aes(x = variable, y = OddsRatio, color = p_value )) +
  geom_point(size = 3) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "red")+
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2) +
  coord_flip() +
  labs(
    title = "Odds Ratios from logistic regression model",
    x = "Variable",
    y = "Odds Ratio (95% CI)",
    color = "P.adjust", size=9
  ) +
  theme_minimal()+scale_y_log10()+
  scale_color_gradient(low = "blue", high = "red")+theme(plot.title = element_text(size = 12, face = "bold"),
    legend.title=element_text(size=10), 
    legend.text=element_text(size=10), 
    axis.text.x=element_text(size=10),
    axis.text.y=element_text(size=10),)
ggsave(OR_plot, filename = "OR_plot_log_top_20.png", bg= "white", height = 3.5, width=6)

```

#ORA on important variables
```{r}
#ORA-------------------------------------

#read in features
important_features <- read.csv("important_features.csv")


important_features <- important_features$X

#read in LR results
LR_results <- read.csv("LR_results.csv")

  #subset from lR results
results <- LR_results[LR_results$variable %in% important_features,]


#split by OR > or <1

#select for just OR > 1
results_df <- results[results$OddsRatio > 1,]
geneList <- results_df$variable
filename <- c("ORA_increased_important_variables")

#select for just OR < 1
results_df <- results[results$OddsRatio < 1,]
geneList <- results_df$variable
filename <- c("ORA_decreased_important_variables")



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

write.csv(ego2, file=paste0("GO_", filename, ".csv"))

GO <- try(barplot(ego2,  font.size = 16 , x="Count")) + scale_y_discrete(labels=function(ego2) str_wrap(ego2, width=60))
  try(ggsave(GO, file=paste0("GO_", filename, ".png"), width = 20, height = 5, units = "cm"))


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



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

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



try(ggsave(H, file=paste0("H_", filename, ".png"), width = 20, height = 5, units = "cm"))



```






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


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

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

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

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

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

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

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

```


#Bootstrap
```{r}
#read in data
#combined_EM_metadata <- read.csv("combat_combined_EM_metadata.csv")

#create 100 partitions
set.seed(285)
 trainIndex <- createDataPartition(combined_EM_metadata$Fail, p = .8, 
                                  list = TRUE,
                                  times = 100)
head(trainIndex)

#for each partition, split the data, perform LR, RFE, train and test the model
for  (name in names(trainIndex)){
  
  #split data -------------------------------------------------
combined_EM_metadata_Train <- combined_EM_metadata[ trainIndex[[name]],]
combined_EM_metadata_Test  <- combined_EM_metadata[-trainIndex[[name]],]

#write.csv(combined_EM_metadata_Train, paste0(name, "_Train.csv"))


#write.csv(combined_EM_metadata_Test, paste0(name, "_Test.csv"))

#Logistic regression---------------------------------------------


#combined_EM_metadata_Train <- read.csv("combined_EM_metadata_Train.csv")%>% remove_rownames %>% column_to_rownames(var="X")


combined_EM_metadata_Train <- combined_EM_metadata_Train %>% dplyr::select(-Survival)

#add factors to all categorical variables
combined_EM_metadata_Train[] <- lapply(combined_EM_metadata_Train, function(x) if(is.character(x)) factor(x) else x)

levels(combined_EM_metadata_Train$Fail)

combined_EM_metadata_Train$Fail <- fct_relevel(combined_EM_metadata_Train$Fail, "Not_fail", "Fail")

  #run univariate logistic regression for 
predictors <- setdiff(names(combined_EM_metadata_Train), "Fail")

# Initialize a vector to store p-values
pvals <- numeric(length(predictors))
names(pvals) <- predictors

odds_ratio <- numeric(length(predictors))
names(odds_ratio) <- predictors

CI_lower <- numeric(length(predictors))
names(CI_lower) <- predictors

CI_upper <- numeric(length(predictors))
names(CI_upper) <- predictors

# Run logistic regression for each predictor
for (var in predictors) {
  
 formula <- as.formula(paste("Fail ~ `", var, "`", sep = ""))

  # Fit the model
  fit <- glm(formula, data = combined_EM_metadata_Train, family = "binomial")
  
  pvals[var] <- summary(fit)$coefficients[2, 4]  # Extract the p-value
  #coef[var] <- unname(coef(fit)[2])  # Extract the p-value
  #coef <- coef(fit)[2]  # get the coefficient for the feature
  odds_ratio[var] <- exp(unname(coef(fit)[2]))
  
CI_lower[var] <- exp(confint(fit)[2,1 ])
CI_upper[var] <- exp(confint(fit)[2,2 ])

}

# View the p-values
pvals


pvals_adj <- p.adjust(pvals, method = "BH")  # FDR adjustment

results <- data.frame(
  variable = names(pvals),
  p_value = pvals,
   OddsRatio = odds_ratio,
  p_adj = pvals_adj,
  CI_lower = CI_lower, 
  CI_upper= CI_upper
)

# Sort by adjusted p-value
results <- results[order(results$p_adj), ]

write.csv(results, paste0(name, "_LR_results.csv"))

# View significant results (e.g., FDR < 0.01)
#significant <- subset(results, p_value < 0.01)

# #find expression values
# combined_EM_metadata_Train_t <- combined_EM_metadata_Train[,significant$variable]
# 
# range(unlist(combined_EM_metadata_Train_t[sapply(combined_EM_metadata_Train_t, is.numeric)]), na.rm = TRUE)
# mean(unlist(combined_EM_metadata_Train_t[sapply(combined_EM_metadata_Train_t, is.numeric)]), na.rm = TRUE)
# 
# 
# 
# 
# # Top 10 with highest ORs
# top_highest <- significant[order(-significant$OddsRatio), ][1:10, ]
# 
# # Top 10 with lowest ORs
# top_lowest <- significant[order(significant$OddsRatio), ][1:10, ]
# 
# # Combine
# top_combined <- rbind(top_highest, top_lowest)
# 
# top_combined$variable <- factor(top_combined$variable, levels = top_combined$variable[order(top_combined$OddsRatio)])
# 
# # Log-transform OR
# top_combined$log_OR <- log(top_combined$OddsRatio)
# 
# 
# significant$variable <- factor(significant$variable, levels = significant$variable[order(significant$OddsRatio)])
# 
# 
# # HR plot no padj colour
# ggplot(significant, aes(x = variable, y = OddsRatio)) +
#   geom_point(size = 3) +
#   geom_hline(yintercept = 1, linetype = "dashed", color = "red")+
#   geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.2) +
#   coord_flip() +
#   labs(
#     title = "Hazard Ratios from Cox Model",
#     x = "Gene",
#     y = "Hazard Ratio (95% CI)"
#   ) +
#   theme_minimal()
# 
# #HR plot colour y padjust
# OR_plot <- ggplot(top_combined, aes(x = variable, y = log_OR, color = p_value )) +
#   geom_point(size = 3) +
#   geom_hline(yintercept = 1, linetype = "dashed", color = "red")+
#   geom_errorbar(aes(ymin = log(CI_lower), ymax = log(CI_upper)), width = 0.2) +
#   coord_flip() +
#   labs(
#     title = "Odds Ratios from logistic regression model",
#     x = "Variable",
#     y = "Odds Ratio (95% CI)",
#     color = "P.adjust", size=9
#   ) +
#   theme_minimal()+
#   scale_color_gradient(low = "blue", high = "red")+theme(plot.title = element_text(size = 12, face = "bold"),
#     legend.title=element_text(size=10), 
#     legend.text=element_text(size=10), 
#     axis.text.x=element_text(size=10),
#     axis.text.y=element_text(size=10),)
# ggsave(OR_plot, filename = "OR_plot.png", bg= "white", height = 3.5, width=6)


#get the number of vars with the highest accuracy <50
# Subset the results data frame to those with < 50 variables
subsetted_rfe <- rfe_results$results[rfe_results$results$Variables < 50, ]

# Get number of variables where accuracy is highest
number_of_vars <- subsetted_rfe$Variables[which.max(subsetted_rfe$Accuracy)]

# Get the top 'number_of_vars' selected features
selected_features <- rfe_results$optVariables[1:number_of_vars]

write.csv(selected_features, "selected_features.csv")

# Subset the dataset to keep only the selected features
subset_data <- combined_EM_metadata_Train[, selected_features]

#add fail back in to df 
subset_data$Fail <- combined_EM_metadata_Train$Fail

# Check the structure of the subset
str(subset_data)

#relevel so first level is event
subset_data$Fail <- fct_relevel(subset_data$Fail, "Fail", "Not_fail")

}


```

#Get overlapping LR results
```{r}
# Load required libraries
library(dplyr)

# Set your folder path where all *_LR_results.csv files are stored
lr_results_path <- "/scratch/cmp1g20/ML/log_regression_results/"

# List all LR result CSVs
lr_files <- list.files(lr_results_path, pattern = "_LR_results.csv$", full.names = TRUE)

# Initialize a list to store significant genes per split
signif_gene_lists <- list()

# Loop over each file and extract significant variables
for (file in lr_files) {
  df <- read.csv(file)
  
  # Filter for adjusted p-value < 0.01
  sig_genes <- df %>%
    filter(p_adj < 0.01) %>%
    pull(variable)
  
  # Store in list using file name as identifier
  signif_gene_lists[[basename(file)]] <- sig_genes
}

# Find the intersection (genes common to all splits)
overlapping_genes <- Reduce(intersect, signif_gene_lists)

# Output result
cat("Number of overlapping significant genes across all splits:", length(overlapping_genes), "\n")
print(overlapping_genes)

# Optional: Save to file
#write.csv(overlapping_genes, "/scratch/cmp1g20/ML/overlapping_significant_genes.csv", row.names = FALSE)


```

#Overlapping LR with same OR
```{r}
# Load required libraries
library(dplyr)

# Set path
lr_results_path <- "/scratch/cmp1g20/ML/log_regression_results/"

# List all result files
lr_files <- list.files(lr_results_path, pattern = "_LR_results.csv$", full.names = TRUE)

# To track OR direction per gene
gene_or_direction <- list()

# To track significant genes per split
signif_gene_lists <- list()

# Loop through each file
for (file in lr_files) {
  df <- read.csv(file)
  
  # Filter for adjusted p < 0.01
  sig_df <- df %>%
    filter(p_adj < 0.01)
  
  # Extract gene names
  sig_genes <- sig_df$variable
  signif_gene_lists[[basename(file)]] <- sig_genes
  
  # Track OR direction
  for (i in seq_len(nrow(sig_df))) {
    gene <- sig_df$variable[i]
    or <- sig_df$OddsRatio[i]
    
    # Initialize tracking if not already
    if (is.null(gene_or_direction[[gene]])) {
      gene_or_direction[[gene]] <- c()
    }
    
    # Save OR direction only (not value)
    gene_or_direction[[gene]] <- c(gene_or_direction[[gene]], ifelse(or > 1, "up", "down"))
  }
}

# Step 1: Find genes significant in *all* splits
overlapping_genes <- Reduce(intersect, signif_gene_lists)

# Step 2: Filter for consistent OR direction
# Separate into upregulated (OR > 1) and downregulated (OR < 1)
upregulated_genes <- c()
downregulated_genes <- c()

for (gene in overlapping_genes) {
  directions <- gene_or_direction[[gene]]
  unique_dirs <- unique(directions)
  
  if (length(unique_dirs) == 1) {
    if (unique_dirs == "up") {
      upregulated_genes <- c(upregulated_genes, gene)
    } else if (unique_dirs == "down") {
      downregulated_genes <- c(downregulated_genes, gene)
    }
  }
}

# Output
cat("✅ Overlapping significant genes with OR > 1 (risk):", length(upregulated_genes), "\n")
cat("✅ Overlapping significant genes with OR < 1 (protective):", length(downregulated_genes), "\n")

# Save to separate CSVs
write.csv(upregulated_genes, "overlapping_genes_OR_greater_than_1.csv", row.names = FALSE)
write.csv(downregulated_genes, "overlapping_genes_OR_less_than_1.csv", row.names = FALSE)

```



#GSEA on LR
```{r}

#select for just OR > 1
results_df <- read.csv("overlapping_genes_OR_greater_than_1.csv")
geneList <- results_df$variable
filename <- c("ORA_increased")

#select for just OR < 1
results_df <- read.csv("overlapping_genes_OR_less_than_1.csv")
geneList <- results_df$variable
filename <- c("ORA_decreased")



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

write.csv(ego2, file=paste0("GO_", filename, ".csv"))

GO <- try(barplot(ego2,  font.size = 16 ,showCategory=10, x="Count")) + scale_y_discrete(labels=function(ego2) str_wrap(ego2, width=60))
try(ggsave(GO, file=paste0("GO_", filename, ".png"), width = 20, height = 10, units = "cm"))


```

#RFE feature stability
from the rfe results, check if there are features consistently present in the top 50
```{r}
library(dplyr)
library(ggplot2)

# Path to your selected features files
selected_features_path <- "/mainfs/scratch/cmp1g20/ML/RFE_results/"

# List all selected features files
feature_files <- list.files(selected_features_path, pattern = "_selected_features.csv$", full.names = TRUE)

# Initialize an empty vector to store all features from all splits
all_features <- c()

# Loop over each file and read the selected features, append to vector
for (file in feature_files) {
  features <- read.csv(file, stringsAsFactors = FALSE)[,1]  # assuming features are in first column
  all_features <- c(all_features, features)
}

# Count how often each feature appears
feature_counts <- table(all_features) %>% as.data.frame()
colnames(feature_counts) <- c("Feature", "Count")

# Calculate proportion (frequency / number of splits)
num_splits <- length(feature_files)
feature_counts$Proportion <- feature_counts$Count / num_splits

# Sort by proportion descending
feature_counts <- feature_counts %>% arrange(desc(Proportion))

# View top stable features
head(feature_counts, 20)

# Plot top 20 most stable features
ggplot(feature_counts[1:20, ], aes(x = reorder(Feature, Proportion), y = Proportion)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  coord_flip() +
  labs(title = "Top 20 Most Stable Features Across 100 Splits",
       x = "Feature",
       y = "Proportion of Splits Selected") +
  theme_minimal()

```

#Internal model performance ROC
For the 100 training sets calculated ROC curve
```{r}
library(dplyr)
library(purrr)
library(pROC)
library(ggplot2)

plot_mean_roc_from_csv <- function(roc_csv_folder, auc_csv_folder = NULL) {
  # 1. Read all ROC CSVs
  roc_files <- list.files(roc_csv_folder, pattern = "\\.csv$", full.names = TRUE)
  roc_list <- lapply(roc_files, function(f) {
    # read the CSV as data frame
    df <- read.csv(f, stringsAsFactors = FALSE)
    # Expect columns: fpr, tpr, thresholds, maybe split
    # Make sure numeric
    df <- df %>%
      mutate(across(c(fpr, tpr, thresholds), as.numeric))
    return(df)
  })
  
  # 2. Interpolate to common FPR grid
  fpr_grid <- seq(0, 1, length.out = 100)
  interp_rocs <- lapply(roc_list, function(df) {
    approx(x = df$fpr, y = df$tpr, xout = fpr_grid, ties = mean)$y
  })
  
  tpr_matrix <- do.call(cbind, interp_rocs)
  
  mean_tpr <- rowMeans(tpr_matrix, na.rm = TRUE)
  sd_tpr <- apply(tpr_matrix, 1, sd, na.rm = TRUE)
  n <- ncol(tpr_matrix)
  se_tpr <- sd_tpr / sqrt(n)
  ci_upper <- pmin(mean_tpr + 1.96 * se_tpr, 1)
  ci_lower <- pmax(mean_tpr - 1.96 * se_tpr, 0)
  
  # 3. Prepare summary data frame
  roc_summary_df <- data.frame(
    fpr = fpr_grid,
    mean_tpr = mean_tpr,
    ci_lower = ci_lower,
    ci_upper = ci_upper
  )
  
  # 4. Plot
  p <- ggplot(roc_summary_df, aes(x = fpr, y = mean_tpr)) +
    geom_line(color = "blue", size = 1) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), alpha = 0.2, fill = "lightblue") +
    geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "grey") +
    labs(
      title = "Mean ROC Curve with 95% CI",
      x = "False Positive Rate (FPR)",
      y = "True Positive Rate (TPR)"
    ) +
    theme_minimal()
  
  # 5. Optionally read AUCs if provided and annotate
  if (!is.null(auc_csv_folder)) {
    auc_files <- list.files(auc_csv_folder, pattern = "\\.csv$", full.names = TRUE)
    auc_vals <- sapply(auc_files, function(f) {
      x <- read.csv(f, stringsAsFactors = FALSE)
      # assume AUC in first cell
      as.numeric(x[1, 1])
    })
    mean_auc <- mean(auc_vals, na.rm = TRUE)
    sd_auc <- sd(auc_vals, na.rm = TRUE)
    p <- p + 
      ggtitle(sprintf("Mean ROC Curve\nMean AUC = %.3f ± %.3f", mean_auc, sd_auc))
  }
  
  return(p)
}

# Example usage:
roc_folder <- "/scratch/cmp1g20/ML/test_roc_curves/"
auc_folder <- "/scratch/cmp1g20/ML/roc_cv_auc/"

mean_roc_plot <- plot_mean_roc_from_csv(roc_folder, auc_folder)
print(mean_roc_plot)

```

#Test ROC
```{r}
library(tidyverse)

# Define path to ROC CSVs
roc_dir <- "/scratch/cmp1g20/ML/test_roc_curves"

# List all files ending in .csv
roc_files <- list.files(roc_dir, pattern = "*.csv", full.names = TRUE)

# Read all ROC CSVs into one long data frame
all_roc <- roc_files %>%
  setNames(nm = basename(.)) %>%  # keep filenames as IDs
  map_dfr(read_csv, .id = "file")
# Make sure fpr is numeric and round to avoid floating point mismatch
all_roc <- all_roc %>%
  mutate(fpr = round(fpr, 3))

# Compute mean and CI
roc_summary <- all_roc %>%
  group_by(fpr) %>%
  summarise(
    mean_tpr = mean(tpr, na.rm = TRUE),
    sd_tpr = sd(tpr, na.rm = TRUE),
    n = n(),
    se_tpr = sd_tpr / sqrt(n),
    ci_lower = mean_tpr - 1.96 * se_tpr,
    ci_upper = mean_tpr + 1.96 * se_tpr,
    .groups = "drop"
  )
library(ggplot2)

ggplot(roc_summary, aes(x = fpr, y = mean_tpr)) +
  geom_line(color = "#0072B2", size = 1.2) +
  geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), alpha = 0.3, fill = "#0072B2") +
  labs(
    x = "False Positive Rate",
    y = "True Positive Rate",
    title = "Mean ROC Curve with 95% Confidence Interval (Test Data)"
  ) +
  theme_minimal() +
  coord_equal()

```

#combined data plot ROC
```{r}
library(dplyr)
library(ggplot2)

# Function to interpolate TPR at fixed FPR points (0 to 1)
interpolate_roc <- function(roc_df, fpr_points = seq(0, 1, length.out = 100)) {
  approx(x = roc_df$fpr, y = roc_df$tpr, xout = fpr_points, ties = mean)$y
}
# Assume combined_cv_roc or combined_test_roc is loaded as dataframe with columns:
# fpr, tpr, thresholds, split, model

combined_cv_roc <- readRDS("/scratch/cmp1g20/ML/combined_cv_roc.rds")
# OR
combined_test_roc <- readRDS("/scratch/cmp1g20/ML/combined_test_roc.rds")


# Example for combined CV ROC curves
roc_data <- combined_cv_roc  # or combined_test_roc

# Fixed FPR points for interpolation
fixed_fpr <- seq(0, 1, length.out = 100)

# Interpolate all ROC curves (each split/model) to fixed FPR grid
interp_list <- roc_data %>%
  group_by(split, model) %>%
  group_split() %>%
  lapply(function(df) {
    tpr_interp <- interpolate_roc(df, fixed_fpr)
    data.frame(fpr = fixed_fpr, tpr = tpr_interp, split = unique(df$split), model = unique(df$model))
  })

# Combine all interpolated curves into one dataframe
interp_df <- bind_rows(interp_list)
summary_df <- interp_df %>%
  group_by(model, fpr) %>%
  summarise(
    mean_tpr = mean(tpr, na.rm = TRUE),
    lower_ci = quantile(tpr, probs = 0.025, na.rm = TRUE),
    upper_ci = quantile(tpr, probs = 0.975, na.rm = TRUE),
    .groups = 'drop'
  )
summary_df <- interp_df %>%
  group_by(model, fpr) %>%
  summarise(
    mean_tpr = mean(tpr, na.rm = TRUE),
    lower_ci = quantile(tpr, probs = 0.025, na.rm = TRUE),
    upper_ci = quantile(tpr, probs = 0.975, na.rm = TRUE),
    .groups = 'drop'
  )
ggplot(summary_df, aes(x = fpr, y = mean_tpr, color = model, fill = model)) +
  geom_line(size = 1) +
  geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci), alpha = 0.2, color = NA) +
  labs(title = "Mean ROC Curve with 95% CI",
       x = "False Positive Rate",
       y = "True Positive Rate",
       color = "Model",
       fill = "Model") +
  theme_minimal() +
  coord_equal() +  # Make x and y scales equal
  theme(legend.position = "bottom")

```

