library(survival)
library(reactable)

select <- c("ENSG00000135365", "ENSG00000162747", "ENSG00000180871", "ENSG00000186529", "ENSG00000105835",
            "ENSG00000081237", "ENSG00000090104", "ENSG00000211772", "ENSG00000117091","ENSG00000145649",
            "ENSG00000106066", "ENSG00000140968", "ENSG00000100079", "ENSG00000127951", "ENSG00000128815",
            "ENSG00000112149", "ENSG00000167604", "ENSG00000166501", "ENSG00000184588", "ENSG00000173193",
            "ENSG00000158050", "ENSG00000111796", "ENSG00000188042", "ENSG00000101109", "ENSG00000277734")



Adeno_bulk <- readRDS("~/new_data/Adeno_bulk.rds")
bulk <- Adeno_bulk
bulk <- bulk[select,]
bulk <- t(bulk)
bulk <- as.data.frame(bulk)
LUAD_BULK_CLINIC <- readRDS("~/new_data/LUAD_BULK_CLINIC.rds")
select_col <- c("overall_survival", "vital_status")
clinic <- LUAD_BULK_CLINIC[,select_col]
all_clinic <- clinic
data <- cbind(all_clinic,bulk)
data$vital_status <- ifelse(data$vital_status == "Alive", 1, 0)
data_clean <- na.omit(data)
univ_cox <- function(variable) {
  cox_model <- coxph(Surv(overall_survival, vital_status) ~ variable, data = data_clean)
  cox_summary <- summary(cox_model)
  
  c(
    coef = cox_summary$coefficients[1],  
    HR = exp(cox_summary$coefficients[1]), 
    p_value = cox_summary$coefficients[5], 
    CI_lower = exp(cox_summary$conf.int[ ,3]), 
    CI_upper = exp(cox_summary$conf.int[ ,4])
  )
}
univ_results <- lapply(data_clean[, 3:ncol(data_clean)], univ_cox) 
univ_results_df <- do.call(rbind, univ_results)
univ_results_df <- as.data.frame(univ_results_df)
rownames(univ_results_df) <- colnames(data_clean[, 3:ncol(data_clean)])
reactable(univ_results_df)
univ_results_df$cell_type[rownames(univ_results_df) == "ENSG00000135365"] <- "neutrophils"

univ_results_df$cell_type[rownames(univ_results_df) == "ENSG00000162747"] <- "neutrophils"

univ_results_df$cell_type[rownames(univ_results_df) == "ENSG00000180871"] <- "neutrophils"

univ_results_df$cell_type[rownames(univ_results_df) == "ENSG00000186529"] <- "neutrophils"

univ_results_df$cell_type[rownames(univ_results_df) == "ENSG00000105835"] <- "neutrophils"

#"ENSG00000105835" "ENSG00000135365" "ENSG00000162747" "ENSG00000180871" "ENSG00000186529"

univ_results_df$cell_type[rownames(univ_results_df) == "ENSG00000081237"] <- "t_cell_CD8_effector_memory"

univ_results_df$cell_type[rownames(univ_results_df) == "ENSG00000117091"] <- "t_cell_CD8_effector_memory"

univ_results_df$cell_type[rownames(univ_results_df) == "ENSG00000090104"] <- "t_cell_CD8_effector_memory"

univ_results_df$cell_type[rownames(univ_results_df) == "ENSG00000211772"] <- "t_cell_CD8_effector_memory"

univ_results_df$cell_type[rownames(univ_results_df) == "ENSG00000145649"] <- "t_cell_CD8_effector_memory"

#"ENSG00000081237"  "ENSG00000090104" "ENSG00000211772"  "ENSG00000145649", ENSG00000117091

univ_results_df$cell_type[rownames(univ_results_df) == "ENSG00000158050"] <- "t_cell_CD8_activated"

univ_results_df$cell_type[rownames(univ_results_df) == "ENSG00000111796"] <- "t_cell_CD8_activated"

univ_results_df$cell_type[rownames(univ_results_df) == "ENSG00000188042"] <- "t_cell_CD8_activated"

univ_results_df$cell_type[rownames(univ_results_df) == "ENSG00000101109"] <- "t_cell_CD8_activated"

univ_results_df$cell_type[rownames(univ_results_df) == "ENSG00000277734"] <- "t_cell_CD8_activated"

#"ENSG00000158050","ENSG00000111796", "ENSG00000188042", ENSG00000101109, ENSG00000277734
univ_results_df$cell_type[rownames(univ_results_df) == "ENSG00000112149"] <- "B_cell"

univ_results_df$cell_type[rownames(univ_results_df) == "ENSG00000167604"] <- "B_cell"

univ_results_df$cell_type[rownames(univ_results_df) == "ENSG00000166501"] <- "B_cell"

univ_results_df$cell_type[rownames(univ_results_df) == "ENSG00000184588"] <- "B_cell"

univ_results_df$cell_type[rownames(univ_results_df) == "ENSG00000173193"] <- "B_cell"


# "ENSG00000112149", "ENSG00000081237", "ENSG00000166501", "ENSG00000184588", "ENSG00000173193"
univ_results_df$cell_type[rownames(univ_results_df) == "ENSG00000106066"] <- "cDC1"

univ_results_df$cell_type[rownames(univ_results_df) == "ENSG00000140968"] <- "cDC1"

univ_results_df$cell_type[rownames(univ_results_df) == "ENSG00000100079"] <- "cDC1"

univ_results_df$cell_type[rownames(univ_results_df) == "ENSG00000127951"] <- "cDC1"

univ_results_df$cell_type[rownames(univ_results_df) == "ENSG00000128815"] <- "cDC1"

# "ENSG00000106066", "ENSG00000140968", "ENSG00000100079", "ENSG00000127951", "ENSG00000128815"

univ_results_df$ln_B <- abs(log(univ_results_df$coef))

#"ENSG00000105835" "ENSG00000135365" "ENSG00000162747" "ENSG00000180871" "ENSG00000186529":neu
univ_results_df$si[rownames(univ_results_df) == "ENSG00000105835"] <- 1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000162747"] <- 1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000180871"] <- 1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000186529"] <- 1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000135365"] <- 1
#"ENSG00000081237" "ENSG00000117091" "ENSG00000090104" "ENSG00000211772"  "ENSG00000145649": eff

univ_results_df$si[rownames(univ_results_df) == "ENSG00000081237"] <- -1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000090104"] <- -1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000145649"] <- -1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000211772"] <- -1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000117091"] <- -1

#"ENSG00000158050", "ENSG00000111796", "ENSG00000188042", "ENSG00000101109", "ENSG00000277734": act
univ_results_df$si[rownames(univ_results_df) == "ENSG00000158050"] <- 1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000111796"] <- 1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000188042"] <- 1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000101109"] <- 1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000277734"] <- 1

# "ENSG00000112149", "ENSG00000167604", "ENSG00000166501", "ENSG00000184588", "ENSG00000173193": B
univ_results_df$si[rownames(univ_results_df) == "ENSG00000112149"] <- -1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000167604"] <- -1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000166501"] <- -1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000184588"] <- -1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000173193"] <- -1

# "ENSG00000106066", "ENSG00000140968", "ENSG00000100079", "ENSG00000127951", "ENSG00000128815": cDC1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000106066"] <- -1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000140968"] <- -1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000100079"] <- -1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000127951"] <- -1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000128815"] <- -1

univ_results_df$beta <- univ_results_df$si * univ_results_df$ln_B
univ_results_df$ens <- rownames(univ_results_df)
library(biomaRt)
mart <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
ensembl_ids <- rownames(univ_results_df)
gene_conversion <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol"), 
                         filters = "ensembl_gene_id", 
                         values = ensembl_ids, 
                         mart = mart)
univ_results_df$Gene_Name <- gene_conversion$hgnc_symbol[match(rownames(univ_results_df), gene_conversion$ensembl_gene_id)]
univ_results_df <- as.data.frame(univ_results_df)
univ_results_df <- univ_results_df[!is.na(univ_results_df$Gene_Name), ] 
rownames(univ_results_df) <- univ_results_df$Gene_Name  
univ_results_df$Gene_Name <- NULL
reactable(univ_results_df)
final <- data_clean[, c("ENSG00000135365", "ENSG00000162747", "ENSG00000180871", "ENSG00000186529", "ENSG00000105835",
            "ENSG00000081237", "ENSG00000090104", "ENSG00000211772", "ENSG00000117091","ENSG00000145649",
            "ENSG00000106066", "ENSG00000140968", "ENSG00000100079", "ENSG00000127951", "ENSG00000128815",
            "ENSG00000112149", "ENSG00000167604", "ENSG00000166501", "ENSG00000184588", "ENSG00000173193",
            "ENSG00000158050", "ENSG00000111796", "ENSG00000188042", "ENSG00000101109", "ENSG00000277734")]

si <- univ_results_df$si
ln_B <- univ_results_df$ln_B
final$risk_scores <- apply(final, 1, function(x) sum(x * si * ln_B))
head(final$risk_scores)
## [1] -221725.164  -56358.127  -92280.171  187195.936  -57500.873    8218.054
library(survival)
library(survminer)
## Loading required package: ggplot2
## Loading required package: ggpubr
## 
## Attaching package: 'survminer'
## The following object is masked from 'package:survival':
## 
##     myeloma
final$overall_survival <- data_clean$overall_survival
final$vital_status <- data_clean$vital_status
final$deceased <- ifelse(final$vital_status == "1", FALSE, TRUE)
final$overall_survival <- as.numeric(final$overall_survival)
plot_list <- list()
res.cut <- surv_cutpoint(final, time = "overall_survival", event = "deceased", 
                         variables = c("risk_scores"))
cutpoint <- res.cut[["cutpoint"]][["cutpoint"]]
final$risk_scores <- ifelse(final$risk_scores > cutpoint, "HIGH", "LOW")
fit_risk_scores <- survfit(Surv(overall_survival, deceased) ~ risk_scores, data = final)
ggsurvplot(fit_risk_scores,
           final,
           pval = TRUE, 
           )