library(survival)
library(survminer)
## Loading required package: ggplot2
## Loading required package: ggpubr
## 
## Attaching package: 'survminer'
## The following object is masked from 'package:survival':
## 
##     myeloma
library(ggplot2)
data <- readRDS("~/YEPYENI/clinic_deconv_scc.rds")
data$vital_status <- as.character(data$vital_status)
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
data$overall_survival <- ifelse(data$vital_status == "Alive",
                                as.numeric(data$days_to_last_follow_up), 
                                as.numeric(data$days_to_death))
data$overall_survival <- as.numeric(data$overall_survival)

data$ajcc_pathologic_stage <- as.character(data$ajcc_pathologic_stage)

data$combined_stage <- NA

data$combined_stage <- ifelse(data$ajcc_pathologic_stage %in% c("Stage I", "Stage IA", "Stage IB"), "Stage I",
                              ifelse(data$ajcc_pathologic_stage %in% c("Stage IIA", "Stage IIB", "Stage IIIB"), "Stage II",
                                     ifelse(data$ajcc_pathologic_stage %in% c("Stage IIIA", "Stage IIIB"), "Stage III",
                                            ifelse(data$ajcc_pathologic_stage == "Stage IV", "Stage IV", NA))))

data$combined_stage <- factor(data$combined_stage, levels = c("Stage I", "Stage II", "Stage III", "Stage IV"))

table(data$combined_stage)
## 
##   Stage I  Stage II Stage III  Stage IV 
##       227       166        58         6
# Define columns for each domain
b_selected_columns <- c("Plasma cell", "Plasma cell dividing", "Mast cell", "B cell", "B cell dividing")
t_selected_columns <- c("T cell CD4", "T cell CD4 dividing", "T cell CD8 activated", "T cell CD8 dividing", 
                         "T cell CD8 effector memory", "T cell CD8 naive", "T cell CD8 terminally exhausted", 
                         "T cell NK-like", "T cell regulatory")
m_selected_columns <- c("Macrophage", "Myeloid dividing", "Macrophage alveolar", "cDC1", "cDC2", "pDC", 
                         "DC mature", "Monocyte classical", "Monocyte non-classical", "Myeloid dividing")
ct_selected_columns <- c("Alveolar cell type 1", "Alveolar cell type 2", "Cancer cells", "Ciliated", "Club",
                          "Endothelial cell arterial", "Endothelial cell capillary", "Endothelial cell lymphatic",
                          "Endothelial cell venous", "Fibroblast adventitial", "Fibroblast alveolar", 
                          "Fibroblast peribronchial", "Mesothelial", "Neutrophils", "NK cell", "NK cell dividing", 
                          "Pericyte", "ROS1+ healthy epithelial", "Smooth muscle cell", "Stromal dividing", 
                          "Transitional Club/AT2")

create_score_matrix <- function(data, selected_columns) {
  matrix <- data[, selected_columns]
  matrix <- as.matrix(matrix)
  row_sums <- rowSums(matrix, na.rm = TRUE)
  matrix_df <- as.data.frame(matrix)
  matrix_df$Sum <- row_sums
  return(matrix_df)
}


B_matrix_df <- create_score_matrix(data, b_selected_columns)
T_matrix_df <- create_score_matrix(data, t_selected_columns)
M_matrix_df <- create_score_matrix(data, m_selected_columns)
CT_matrix_df <- create_score_matrix(data, ct_selected_columns)

score_df_scc <- data.frame(
  B_Sum = B_matrix_df$Sum,
  T_Sum = T_matrix_df$Sum,
  M_Sum = M_matrix_df$Sum,
  Cold_Tumor_Sum = CT_matrix_df$Sum
)
rownames(score_df_scc) <- rownames(B_matrix_df)
saveRDS(score_df_scc, "~/YEPYENI/score_df_scc.rds")

data <- cbind(data, score_df_scc)
res.cut_combined <- surv_cutpoint(data, time = "overall_survival", event = "deceased", variables = c("B_Sum", "T_Sum", "M_Sum", "Cold_Tumor_Sum"))
global_cutoff <- res.cut_combined[["cutpoint"]][["cutpoint"]]

data$b_group <- ifelse(data$B_Sum > global_cutoff, "HIGH", "LOW")
## Warning in data$B_Sum > global_cutoff: longer object length is not a multiple
## of shorter object length
data$t_group <- ifelse(data$T_Sum > global_cutoff, "HIGH", "LOW")
## Warning in data$T_Sum > global_cutoff: longer object length is not a multiple
## of shorter object length
data$m_group <- ifelse(data$M_Sum > global_cutoff, "HIGH", "LOW")
## Warning in data$M_Sum > global_cutoff: longer object length is not a multiple
## of shorter object length
data$ct_group <- ifelse(data$Cold_Tumor_Sum > global_cutoff, "HIGH", "LOW")
## Warning in data$Cold_Tumor_Sum > global_cutoff: longer object length is not a
## multiple of shorter object length
# Perform Kaplan-Meier survival analysis for each group
fit_b <- survfit(Surv(overall_survival, deceased) ~ b_group, data = data)
fit_t <- survfit(Surv(overall_survival, deceased) ~ t_group, data = data)
fit_m <- survfit(Surv(overall_survival, deceased) ~ m_group, data = data)
fit_ct <- survfit(Surv(overall_survival, deceased) ~ ct_group, data = data)

b_df <- data.frame(Time = fit_b$time, Surv = fit_b$surv, Group = "B")
t_df <- data.frame(Time = fit_t$time, Surv = fit_t$surv, Group = "T")
m_df <- data.frame(Time = fit_m$time, Surv = fit_m$surv, Group = "M")
ct_df <- data.frame(Time = fit_ct$time, Surv = fit_ct$surv, Group = "CT")

plot_df <- rbind(b_df, t_df, m_df, ct_df)

p_luad <- ggplot(plot_df, aes(x = Time, y = Surv, color = Group)) +
  geom_line() +
  labs(title = "Kaplan-Meier Survival Curves with Global Cutoff for LUSC", x = "Time", y = "Survival Probability") +
  theme_minimal()

print(p_luad)

library(GGally)
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
ggpairs(score_df_scc[, c("B_Sum", "T_Sum", "M_Sum", "Cold_Tumor_Sum")]) +
  theme_minimal()

# Separate plots for each domain
p_b <- ggsurvplot(fit_b, data = data, pval = TRUE, risk.table = TRUE, title = "B Domain Survival")
p_t <- ggsurvplot(fit_t, data = data, pval = TRUE, risk.table = TRUE, title = "T Domain Survival")
p_m <- ggsurvplot(fit_m, data = data, pval = TRUE, risk.table = TRUE, title = "M Domain Survival")
p_ct <- ggsurvplot(fit_ct, data = data, pval = TRUE, risk.table = TRUE, title = "CT Domain Survival")

# Print plots
print(p_b)

print(p_t)

print(p_m)

print(p_ct)

library(ggplot2)
library(reshape2)

melted_data <- melt(score_df_scc[, c("B_Sum", "T_Sum", "M_Sum", "Cold_Tumor_Sum")])
## No id variables; using all as measure variables
ggplot(melted_data, aes(x = value, fill = variable)) +
  geom_density(alpha = 0.5) +
  labs(title = "Density Plot of All Domains", x = "Value", y = "Density") +
  scale_fill_manual(name = "Domain", values = c("B_Sum" = "blue", "T_Sum" = "red", "M_Sum" = "green", "Cold_Tumor_Sum" = "purple")) +
  theme_minimal()

fit_stage <- survfit(Surv(overall_survival, deceased) ~ combined_stage, data = data)
library(survminer)
ggsurvplot(fit_stage, 
           data = data, 
           pval = TRUE, 
           risk.table = TRUE, 
           legend.title = "Stage",
           title = "Kaplan-Meier Survival Curves by Stage LUSC")

library(ggplot2)

ggplot(score_df_scc, aes(x = T_Sum, y = B_Sum)) +
  geom_point(aes(size = M_Sum, color = M_Sum), shape = 16, alpha = 0.7) +
  scale_size_continuous(name = "M_Sum", range = c(2, 10)) +
  scale_color_gradient(low = "lightgreen", high = "darkgreen", name = "M_Sum") +
  labs(title = "Plot of B_Sum vs T_Sum with M_Sum as Points",
       x = "T_Sum",
       y = "B_Sum") +
  theme_minimal(base_size = 15) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "right",
    axis.title.x = element_text(face = "bold"),
    axis.title.y = element_text(face = "bold")
  )

library(ggplot2)

ggplot(score_df_scc, aes(x = M_Sum, y = B_Sum)) +
  geom_point(aes(size = T_Sum, color = T_Sum), shape = 16, alpha = 0.7) +
  scale_size_continuous(name = "T_Sum", range = c(2, 10)) +
  scale_color_gradient(low = "lightgreen", high = "darkgreen", name = "T_Sum") +
  labs(title = "Plot of M_Sum vs B_Sum with T_Sum as Points LUSC",
       x = "M_Sum",
       y = "B_Sum") +
  theme_minimal(base_size = 15) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "right",
    axis.title.x = element_text(face = "bold"),
    axis.title.y = element_text(face = "bold")
  )

library(ggplot2)

ggplot(score_df_scc, aes(x = M_Sum, y = Cold_Tumor_Sum)) +
  geom_point(aes(size = B_Sum, color = B_Sum), shape = 16, alpha = 0.7) +
  scale_size_continuous(name = "T_Sum", range = c(2, 10)) +
  scale_color_gradient(low = "lightgreen", high = "darkgreen", name = "B_Sum") +
  labs(title = "Plot of M_Sum vs Cold_Tumor_Sum with B_Sum as Points LUSC ",
       x = "M_Sum",
       y = "Cold_Tumor_Sum") +
  theme_minimal(base_size = 15) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "right",
    axis.title.x = element_text(face = "bold"),
    axis.title.y = element_text(face = "bold")
  )

library(ggplot2)

ggplot(score_df_scc, aes(x = B_Sum, y = T_Sum)) +
  geom_point(aes(size = Cold_Tumor_Sum, color = Cold_Tumor_Sum), shape = 16, alpha = 0.7) +
  scale_size_continuous(name = "T_Sum", range = c(2, 10)) +
  scale_color_gradient(low = "lightgreen", high = "darkgreen", name = "Cold_Tumor_Sum") +
  labs(title = "Plot of B_Sum vs T_Sum with Cold_Tumor_Sum as Points LUSC",
       x = "B_Sum",
       y = "T_Sum") +
  theme_minimal(base_size = 15) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "right",
    axis.title.x = element_text(face = "bold"),
    axis.title.y = element_text(face = "bold")
  )

library(ggplot2)

ggplot(score_df_scc, aes(x = B_Sum, y = M_Sum)) +
  geom_point(aes(size = Cold_Tumor_Sum, color = Cold_Tumor_Sum), shape = 16, alpha = 0.7) +
  scale_size_continuous(name = "T_Sum", range = c(2, 10)) +
  scale_color_gradient(low = "lightgreen", high = "darkgreen", name = "Cold_Tumor_Sum") +
  labs(title = "Plot of B_Sum vs M_Sum with Cold_Tumor_Sum as Points LUSC",
       x = "B_Sum",
       y = "M_Sum") +
  theme_minimal(base_size = 15) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "right",
    axis.title.x = element_text(face = "bold"),
    axis.title.y = element_text(face = "bold")
  )