luad_lusc_risk_score

library(survival)
library(reactable)

select <- c("ENSG00000135365", "ENSG00000162747", "ENSG00000180871", "ENSG00000186529", "ENSG00000105835",
            "ENSG00000106066", "ENSG00000140968", "ENSG00000100079", "ENSG00000127951", "ENSG00000128815")
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)])
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) == "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

# "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",
            "ENSG00000106066", "ENSG00000140968", "ENSG00000100079", "ENSG00000127951", "ENSG00000128815")]


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))
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, 
           )

library(survival)
library(reactable)

select <- c("ENSG00000169429", "ENSG00000123689", "ENSG00000136689" ,"ENSG00000119535", "ENSG00000276070",
            "ENSG00000231389", "ENSG00000179344", "ENSG00000140968", "ENSG00000128815", "ENSG00000100079")

sq_bulk <- readRDS("~/new_data/Squamous_bulk.rds")
bulk <- sq_bulk
bulk <- bulk[select,]
bulk <- t(bulk)
bulk <- as.data.frame(bulk)
LUSC_BULK_CLINIC <- readRDS("~/new_data/LUSC_BULK_CLINIC.rds")
select_col <- c("overall_survival", "vital_status")
clinic <- LUSC_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)])
univ_results_df$cell_type[rownames(univ_results_df) == "ENSG00000169429"] <- "neutrophils"

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

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

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

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

#"ENSG00000169429", "ENSG00000123689", "ENSG00000136689" ,"ENSG00000119535", "ENSG00000276070"

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

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

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

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

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

# "ENSG00000231389", "ENSG00000179344", "ENSG00000140968", "ENSG00000128815", "ENSG00000100079"

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

#""ENSG00000169429", "ENSG00000123689", "ENSG00000136689" ,"ENSG00000119535", "ENSG00000276070":neu
univ_results_df$si[rownames(univ_results_df) == "ENSG00000169429"] <- 1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000123689"] <- 1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000136689"] <- 1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000119535"] <- 1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000276070"] <- 1

# "ENSG00000231389", "ENSG00000179344", "ENSG00000140968", "ENSG00000128815", "ENSG00000100079": cDC1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000231389"] <- -1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000179344"] <- -1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000140968"] <- -1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000128815"] <- -1
univ_results_df$si[rownames(univ_results_df) == "ENSG00000100079"] <- -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("ENSG00000169429", "ENSG00000123689", "ENSG00000136689" ,"ENSG00000119535", "ENSG00000276070",
            "ENSG00000231389", "ENSG00000179344", "ENSG00000140968", "ENSG00000128815", "ENSG00000100079")]
ln_B <- univ_results_df$ln_B
final$risk_scores <- apply(final, 1, function(x) sum(x * si * ln_B))
library(survival)
library(survminer)
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, 
           )