library(tidyverse)
library(tidyr)
library(Rtsne)
library(readxl)
library(ggrepel)


# Data --------------------------------------------------------------------

#methylation profiles
#2/3 split of methylation data (training set)
meth <-  readRDS(file="~/")
#1/3 split of methylation data (validation set)
meth2 <- readRDS(file="~/")
# Methylation data of cCUP samples
CUP <- readRDS(file="~/")


#cg sites
cgs <- readRDS(file="~/")

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


# Figure 1 ----------------------------------------------------------------

train <- meth[cgs,]
test <- meth2[cgs,]
testtrain <- cbind(test,train)
cup <- CUP[cgs,]
tsne.data <-cbind(testtrain,cup)
tsne.data <- t(tsne.data)

set.seed(11)
rtsne.result <- Rtsne(tsne.data, perplexity = 50)
tsneData <- rtsne.result$Y

colnames(tsneData) <- c("t-SNE dim1", "t-SNE dim2")
rownames(tsneData) <- rownames(tsne.data)
data <- as_tibble(tsneData) %>% mutate(,file_name = rownames(tsneData))
data1<-left_join(data,ann)
data1 <- data1[complete.cases(data1[,6]),]
data1 <- data1[!(data1$project=="Valencia"),]
data_plot1a <- data1[!(data1$Primary_Site=="CUP"),]
data_plot1a <- data_plot1a[!(data_plot1a$Primary_Site=="Hypopharynx"),]

as.factor(data_plot1a$Primary_Site) %>% levels()
data_plot1a$Primary_Site<-factor(data_plot1a$Primary_Site,levels=c("Larynx","Oral_Cavity","Oropharynx","Esophagus","Lung","Lymphoma","Skin")  )

Plot1a <- ggplot(data_plot1a, aes(`t-SNE dim1`,`t-SNE dim2`))  + 
  geom_point(size=1.4,aes(colour = `Primary_Site`),shape=1) +
  labs(colour = "Primary tumor site") +
  scale_color_manual(values = c("#004C99","#3399FF","#99CCFF","#ff8000","#00CC00","#FF3333","#FF00FF"))+
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
print(Plot1a)

# Figure 2 ----------------------------------------------------------------

data_plot2 <- data1 %>%
  mutate(Label = if_else(Primary_Site == "CUP", str_sub(SampleID, 1, 3), "")) %>% 
  arrange(Label)
plot <- data_plot2 %>%
  mutate(CUP_subtype = if_else(`t-SNE dim2` <= 10 & Primary_Site == "CUP", "skin-like", "")) %>% 
  arrange(CUP_subtype)
plot$CUP_subtype[plot$`t-SNE dim1`<= -10 & plot$Primary_Site == "CUP"] <- "oropharynx-like"
plot$CUP_subtype[plot$`t-SNE dim1`>= -10 & plot$`t-SNE dim2`>= 0 & plot$Primary_Site == "CUP"] <- "atypical"
plot$CUP_subtype[plot$CUP_subtype == ""] <- NA


data_plot2 <- data_plot2[!(data_plot2$Primary_Site=="Hypopharynx"),]

Plot2 <- ggplot(plot, aes(`t-SNE dim1`,`t-SNE dim2`))  + 
  geom_point(aes(colour=`CUP_subtype`),size=1.4,shape=1) +
  geom_text_repel(min.segment.length = 0,aes(label = Label),size=2,nudge_x = 2,max.overlaps = Inf) +
  scale_color_manual(values = c("red","blue","#66FF00"),limits= c("oropharynx-like","skin-like","atypical")) +
  theme_bw() +
  labs(colour = "CUP subtype") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
print(Plot2)


# Hierarchical clustering CUPs -------------------------------------------------

library(dendextend)
library(pheatmap)
library(BiocManager)
library(ComplexHeatmap)

plot2 <- plot[c(1104:1131), ]
CUP1<- t(CUP)
CUP1 <- CUP1[,cgs]
CUP2 <- as_tibble(CUP1) %>% mutate(file_name = rownames(CUP1))
CUP2 <- left_join(CUP2,plot2, by = 'file_name')
rownames(CUP1) <- CUP2$SampleID

ann_col <- data.frame(CUP2$CUP_subtype)
rownames(ann_col) <- CUP2$SampleID
colnames(ann_col) <- "CUP subtype"

CUP1 <- t(CUP1)
rownames(ann_col) <- colnames(CUP1)

pheatmap(CUP1,annotation_col = ann_col,show_rownames = F, treeheight_row =  0, 
         main = '')