risk_score

library(survival)

library(reactable)
T_cell_CD8_effector_memory <- readRDS("~/new_data/T_cell_CD8_effector_memory.rds")
reactable(T_cell_CD8_effector_memory)

T_cell_CD8_effector_memory$cell_type <- "CD8 effector memory"

#ENSG00000180644: PRF1

library(reactable)
b_cell_dividing <- readRDS("~/new_data/b_cell_dividing.rds")
reactable(b_cell_dividing)

b_cell_dividing$cell_type <- "B cell dividing"

#ENSG00000148773: MKI67

library(reactable)
cdc1 <- readRDS("~/new_data/cdc1.rds")
reactable(cdc1)

cdc1$cell_type <- "cDC1"

#ENSG00000140968: IRF8

library(reactable)
t_cell_cd8_activated <- readRDS("~/new_data/t_cell_cd8_activated.rds")
t_cell_cd8_activateda <- t_cell_cd8_activated[order(t_cell_cd8_activated$coef), ]
reactable(t_cell_cd8_activateda)

t_cell_cd8_activated$cell_type <- "CD8 activated"

#ENSG00000110848: CD69

library(reactable)
neutrophils <- readRDS("~/new_data/neutrophil.rds")
reactable(neutrophils)

neutrophils$cell_type <- "neutrophils"

#ENSG00000162747: FCGR3B

select <- c("ENSG00000180644", "ENSG00000148773", "ENSG00000140968", "ENSG00000110848", "ENSG00000162747")

Adeno_bulk <- readRDS("~/new_data/Adeno_bulk.rds")
Squamous_bulk <- readRDS("~/new_data/Squamous_bulk.rds")
bulk <- cbind(Adeno_bulk,Squamous_bulk)

bulk <- bulk[select,]
bulk <- t(bulk)
bulk <- as.data.frame(bulk)
LUSC_BULK_CLINIC <- readRDS("~/new_data/LUSC_BULK_CLINIC.rds")
LUAD_BULK_CLINIC <- readRDS("~/new_data/LUAD_BULK_CLINIC.rds")
select_col <- c("overall_survival", "vital_status")

clinic <- LUAD_BULK_CLINIC[,select_col]
clinic1 <- LUSC_BULK_CLINIC[,select_col]
all_clinic <- rbind(clinic,clinic1)
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)])


combine <- combine[select,]
univ_results_df$gene <- combine$Gene_Name
univ_results_df$cell_type <- combine$cell_type

library(ggplot2)

ggplot(univ_results_df, aes(x = gene, y = coef, fill = factor(coef > 0))) +
  geom_col(width = 0.6) +
  scale_fill_manual(values = c("TRUE" = "lightblue", "FALSE" = "salmon"),
                    name = "Coef Value",
                    labels = c("Negative", "Positive")) +
  coord_cartesian(ylim = c(min(univ_results_df$coef),max(univ_results_df$coef))) + 
  labs(title = "Coef Values by Gene Names",
       x = "Gene Names",
       y = "Coef") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top") +
  geom_text(aes(label = cell_type), vjust = -0.5, size = 3, color = "black")

final <- data_clean[, c("ENSG00000180644", "ENSG00000148773", "ENSG00000140968", "ENSG00000110848", "ENSG00000162747")]
gen_expression <- final
coef_values <- univ_results_df$coef
risk_scores <- apply(gen_expression, 1, function(x) sum(x * coef_values))
risk_scores <- as.data.frame(risk_scores)
final$risk_scores <- risk_scores$risk_scores
saveRDS(final, "~/new_data/final.rds")
saveRDS(data_clean, "~/new_data/data_clean.rds")

final <- readRDS("~/new_data/final.rds")
data_clean <- readRDS("~/new_data/data_clean.rds")

library(survival)
library(survminer)

## Loading required package: ggpubr

## 
## Attaching package: 'survminer'

## The following object is masked from 'package:survival':
## 
##     myeloma

final <- na.omit(final)

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)