library(tidyverse)
library(caret)


# Data --------------------------------------------------------------------
#2/3 split of methylation data (training set)
beta <- readRDS(file="~/")

#annotation sheet
ann <- readRDS(file="~/")


# top 3000 CpG sites ---------------------------------------------------------------------
colnames(beta)[736:length(colnames(beta))] <- str_sub(colnames(beta)[736:length(colnames(beta))],end=10L)
sd <- apply(beta,1,sd)
sds <- sd > 0.25
o <- order(sd, decreasing = TRUE)
top_n <- 3000
top <- beta[o[1:top_n],]
sd[o[1:top_n]] %>% hist()
cgs <- rownames(top)

# Preparing test and training sets---------------------------------------------------------------------
df <- beta[cgs,]
meth3k<-t(df)
data <- as_tibble(meth3k) %>% mutate(,file_name = rownames(meth3k))
meth_ann <- left_join(data,ann)
table(meth_ann$Primary_Site)
# Esophagus   Larynx     Lung         Lymphoma   Oral_Cavity  Oropharynx  Skin 
# 63          78         270          32         213          54          25 

test <- readRDS(file="~/")
colnames(test)[369:length(colnames(test))] <- str_sub(colnames(test)[369:length(colnames(test))],end=10L)
test <- test[cgs,]
test2 <- as_tibble(t(test)) %>% mutate(,file_name = colnames(test))
test_ann <- left_join(test2,ann)
table(test_ann$Primary_Site)
# Esophagus   Larynx     Lung         Lymphoma   Oral_Cavity  Oropharynx  Skin 
# 32          39         135          16         107          27          12 

methTest <- subset(test_ann,select=-c(SampleID,Source,Entity,project,file_name)) 
methTrain <- subset(meth_ann,select=-c(SampleID,Source,Entity,project,file_name))

# support vector machine --------------------------------------------------

trctrl <- trainControl(method = "repeatedcv", number = 10, repeats = 3)
set.seed(12)

svm_Linear <- train(Primary_Site ~., data = methTrain, method = "svmLinear",
                    trControl=trctrl,
                    preProcess = c("center", "scale"),
                    tuneLength = 10)

svm <- svm_Linear
pred <- predict(svm, methTest)
pred
methTest$Primary_Site <- as.factor(methTest$Primary_Site)
confusionMatrix(data = pred, reference = methTest$Primary_Site)

# Random forest  ----------------------------------------------------------------------
library(randomForest)
library(renv)

tuneGrid <- expand.grid(.mtry = c(1 : 10))
fitControl <- trainControl(method = "cv",
                           number = 5,
                           search = 'grid',
                           classProbs = TRUE,
                           savePredictions = "final")
ntree <- 500 
nodesize <- 1

rf_model <- train(Primary_Site~.,
                  data = methTrain,
                  method = "rf",
                  importance=TRUE,
                  metric= "Accuracy",
                  tuneGrid = tuneGrid,
                  trControl = fitControl,
                  ntree = ntree,
                  nodesize = nodesize)
forest <- rf_model
pred_rf <- predict(forest, methTest)
pred_rf
methTest$Primary_Site <- as.factor(methTest$Primary_Site)
confusionMatrix(data = pred_rf, reference = methTest$Primary_Site)

# for loop CG# ------------------------------------------------------------

library(ggplot2)
library(dplyr, warn.conflicts = FALSE)
library(officer)
library(rvg)

results <- list()
results2 <- list()

#specify trainControl for machine learning
trctrl <- trainControl(method = "repeatedcv", number = 10, repeats = 3)
#enter file name for powerpoint with plots
fileout <- "tSNE-plots-13Sep22.pptx"
#specify slide layout in powerpoint
read_pptx() %>%
  add_slide(layout = "Title Only", master = "Office Theme") %>%
  ph_with(fig_vg, location = ph_location(left = .5, top = 1.3, height = 5, width = 8)) %>%
  print(fileout)
#fill in cg sizes
cg_sizes <- c(30,100,300,1000,3000,5000,7000,10000,15000,10,200)
# fill in repeats
repeats <-1

df1 <- t(beta)
#for (sample in c(1:repeats)){
#df3 <- samples[[sample]]
for (n in c(10:11)) {
  o <- order(sd, decreasing = TRUE)
  top_n <- cg_sizes[n]
  top <- df1[,o[1:top_n]]
  sd[o[1:top_n]] %>% hist()
  print(dim(top))
  cgs <- colnames(top)
  methTrain1 <- methTrain[,c(cgs,"Primary_Site")]
  methTest1  <- methTest[,c(cgs,"Primary_Site")]
  svm <- train(Primary_Site ~., data = methTrain1, method = "svmLinear",
               trControl=trctrl,
               preProcess = c("center", "scale"),
               tuneLength = 10)
  pred <- predict(svm, methTest1)
  methTest$Primary_Site <- as.factor(methTest1$Primary_Site)
  res <- confusionMatrix(data = pred, reference = methTest1$Primary_Site)
  results[[n]] <- res[3]
}
results2[[length(results2)+1]] <- results
saveRDS(results,file="~/Documents/Data_Paper/RDS/results_SVMaccuracyVScgs#.rds")  
#}
res3 <- lapply(results2, function(x) lapply(x, function(y) y$overall["Accuracy"])) # unlist(x))
tab <- bind_rows(res3) %>% unlist() %>% matrix(ncol = 11, byrow = TRUE)
cg_sizes
tab <- as_tibble(t(tab))%>% mutate("Number of CGs" = cg_sizes[1:11])

tab2 <- pivot_longer(tab,cols=1:1,values_to = "Accuracy of SVM")
tab2$`Number of CGs` <- as.factor(tab2$`Number of CGs`)
ggplot(tab2,aes(`Number of CGs`,`Accuracy of SVM`))+
  geom_point(size=1.4) + expand_limits(y = c(0,1)) +
  theme_bw()

# test classifier with CUP samples ----------------------------------------
CUP3k  <- readRDS(file="~/")
CUP3k <- CUP4k[cgs,]

is <- intersect(colnames(CUP3k),colnames(meth3k))
sub <- t(CUP3k)

pred2 <- predict(svm, sub)

cup <- colnames(CUP4k) %>% as_tibble()
cup2 <- mutate(cup,Prediction = pred2)
cup2 <- mutate(cup2, sampleID = ann$SampleID)