library(survival)
library(survminer)
## Loading required package: ggplot2
## Loading required package: ggpubr
##
## Attaching package: 'survminer'
## The following object is masked from 'package:survival':
##
## myeloma
signature_sq_matrix <- readRDS("~/final_results/14_02_2025_signature_sq_matrix.rds")
sq <- readRDS("~/final_results/14_02_2025_deconvolution_sq.rds")
lusc_clinic_for_bulk <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
deconvolution_sq_scaled <- sq / rowSums(sq) * 100
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
res.cut <- surv_cutpoint(data, time = "overall_survival", event = "deceased",
variables = c("NK_cell_dividing"))
cutpoint <- res.cut[["cutpoint"]][["cutpoint"]]
data$NK_cell_dividing <- ifelse(data$NK_cell_dividing > cutpoint, "HIGH", "LOW")
fit_NK_cell_div <- survfit(Surv(overall_survival, deceased) ~ NK_cell_dividing, data = data)
ggsurvplot(fit_NK_cell_div,
data,
pval = TRUE,
risk.table = TRUE)
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
colnames(data) <- gsub("-", "_", colnames(data))
coxa <- coxph(Surv(overall_survival, deceased) ~ NK_cell_dividing, data = data)
coxa
## Call:
## coxph(formula = Surv(overall_survival, deceased) ~ NK_cell_dividing,
## data = data)
##
## coef exp(coef) se(coef) z p
## NK_cell_dividing -0.05737 0.94425 0.06027 -0.952 0.341
##
## Likelihood ratio test=0.93 on 1 df, p=0.3353
## n= 144, number of events= 143
## (253 observations deleted due to missingness)
library(survival)
library(reactable)
signature_sq_matrix <- as.data.frame(signature_sq_matrix)
Alveolar_exp <- signature_sq_matrix$`NK cell dividing`
other_cells_exp <- rowMeans(signature_sq_matrix[, colnames(signature_sq_matrix) != "NK_cell_dividing"])
Alveolar_lusc_filtered_genes <- signature_sq_matrix[(Alveolar_exp - other_cells_exp) / other_cells_exp > 0.90, ]
sq_bulk <- readRDS("~/final_results/lusc_bulk.rds")
Alveolar <- intersect(rownames(Alveolar_lusc_filtered_genes), rownames(sq_bulk))
Alveolar <- sq_bulk[Alveolar,]
Alveolar <- t(Alveolar)
Alveolar <- as.data.frame(Alveolar)
lusc_BULK_CLINIC <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
lusc_BULK_CLINIC <- as.data.frame(lusc_BULK_CLINIC)
select_col <- c("overall_survival", "vital_status")
clinic <- lusc_BULK_CLINIC[,select_col]
data <- cbind(clinic, Alveolar)
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
library(survival)
exclude_cols <- c("overall_survival", "vital_status", "deceased")
gene_cols <- setdiff(colnames(data), exclude_cols)
cox_results <- list()
for (gene in gene_cols) {
formula <- as.formula(paste("Surv(overall_survival, deceased) ~", gene))
cox_model <- coxph(formula, data = data)
cox_results[[gene]] <- summary(cox_model)
}
negative_genes <- c()
negative_coeffs <- c()
for (gene in names(cox_results)) {
coef_value <- cox_results[[gene]]$coefficients[1, "coef"]
if (!is.na(coef_value) && coef_value < 0) {
negative_genes <- c(negative_genes, gene)
negative_coeffs <- c(negative_coeffs, coef_value)
}
}
library(biomaRt)
ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl", version = 105)
gene_annotation <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol"),
filters = "ensembl_gene_id",
values = negative_genes,
mart = ensembl)
ort <- intersect(negative_genes, rownames(signature_sq_matrix))
sig <- signature_sq_matrix[ort, ]
NK_cell_div_negative_gene_df <- data.frame(
ENSG_ID = negative_genes,
Gene_Name = gene_annotation$hgnc_symbol,
Coefficient = negative_coeffs,
Cell_Type = "NK_cell_div"
)
NK_cell_div_negative_gene_df$exp.level <- sig$`NK cell dividing`
reactable(NK_cell_div_negative_gene_df)
library(survival)
library(survminer)
signature_sq_matrix <- readRDS("~/final_results/14_02_2025_signature_sq_matrix.rds")
sq <- readRDS("~/final_results/14_02_2025_deconvolution_sq.rds")
lusc_clinic_for_bulk <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
deconvolution_sq_scaled <- sq / rowSums(sq) * 100
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
res.cut <- surv_cutpoint(data, time = "overall_survival", event = "deceased",
variables = c("T_cell_CD8_dividing"))
cutpoint <- res.cut[["cutpoint"]][["cutpoint"]]
data$T_cell_CD8_dividing <- ifelse(data$T_cell_CD8_dividing > cutpoint, "HIGH", "LOW")
fit_T_cell_CD8_dividing <- survfit(Surv(overall_survival, deceased) ~ T_cell_CD8_dividing, data = data)
ggsurvplot(fit_T_cell_CD8_dividing,
data,
pval = TRUE,
risk.table = TRUE)
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
colnames(data) <- gsub("-", "_", colnames(data))
coxa <- coxph(Surv(overall_survival, deceased) ~ T_cell_CD8_dividing, data = data)
coxa
## Call:
## coxph(formula = Surv(overall_survival, deceased) ~ T_cell_CD8_dividing,
## data = data)
##
## coef exp(coef) se(coef) z p
## T_cell_CD8_dividing -0.0844 0.9191 0.1137 -0.742 0.458
##
## Likelihood ratio test=0.58 on 1 df, p=0.4466
## n= 144, number of events= 143
## (253 observations deleted due to missingness)
library(survival)
library(reactable)
signature_sq_matrix <- as.data.frame(signature_sq_matrix)
Alveolar_exp <- signature_sq_matrix$`T cell CD8 dividing`
other_cells_exp <- rowMeans(signature_sq_matrix[, colnames(signature_sq_matrix) != "T_cell_CD8_dividing"])
Alveolar_lusc_filtered_genes <- signature_sq_matrix[(Alveolar_exp - other_cells_exp) / other_cells_exp > 0.90, ]
sq_bulk <- readRDS("~/final_results/lusc_bulk.rds")
Alveolar <- intersect(rownames(Alveolar_lusc_filtered_genes), rownames(sq_bulk))
Alveolar <- sq_bulk[Alveolar,]
Alveolar <- t(Alveolar)
Alveolar <- as.data.frame(Alveolar)
lusc_BULK_CLINIC <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
lusc_BULK_CLINIC <- as.data.frame(lusc_BULK_CLINIC)
select_col <- c("overall_survival", "vital_status")
clinic <- lusc_BULK_CLINIC[,select_col]
data <- cbind(clinic, Alveolar)
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
library(survival)
exclude_cols <- c("overall_survival", "vital_status", "deceased")
gene_cols <- setdiff(colnames(data), exclude_cols)
cox_results <- list()
for (gene in gene_cols) {
formula <- as.formula(paste("Surv(overall_survival, deceased) ~", gene))
cox_model <- coxph(formula, data = data)
cox_results[[gene]] <- summary(cox_model)
}
negative_genes <- c()
negative_coeffs <- c()
for (gene in names(cox_results)) {
coef_value <- cox_results[[gene]]$coefficients[1, "coef"]
if (!is.na(coef_value) && coef_value < 0) {
negative_genes <- c(negative_genes, gene)
negative_coeffs <- c(negative_coeffs, coef_value)
}
}
library(biomaRt)
ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl", version = 105)
gene_annotation <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol"),
filters = "ensembl_gene_id",
values = negative_genes,
mart = ensembl)
ort <- intersect(negative_genes, rownames(signature_sq_matrix))
sig <- signature_sq_matrix[ort, ]
T_cell_CD8_dividing_negative_gene_df <- data.frame(
ENSG_ID = negative_genes,
Gene_Name = gene_annotation$hgnc_symbol,
Coefficient = negative_coeffs,
Cell_Type = "T_cell_CD8_dividing"
)
T_cell_CD8_dividing_negative_gene_df$exp.level <- sig$`T cell CD8 dividing`
reactable(T_cell_CD8_dividing_negative_gene_df)
library(survival)
library(survminer)
signature_sq_matrix <- readRDS("~/final_results/14_02_2025_signature_sq_matrix.rds")
sq <- readRDS("~/final_results/14_02_2025_deconvolution_sq.rds")
lusc_clinic_for_bulk <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
deconvolution_sq_scaled <- sq / rowSums(sq) * 100
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
res.cut <- surv_cutpoint(data, time = "overall_survival", event = "deceased",
variables = c("B_cell"))
cutpoint <- res.cut[["cutpoint"]][["cutpoint"]]
data$B_cell <- ifelse(data$B_cell > cutpoint, "HIGH", "LOW")
fit_B_cell<- survfit(Surv(overall_survival, deceased) ~ B_cell, data = data)
ggsurvplot(fit_B_cell,
data,
pval = TRUE,
risk.table = TRUE)
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
colnames(data) <- gsub("-", "_", colnames(data))
coxa <- coxph(Surv(overall_survival, deceased) ~ B_cell, data = data)
coxa
## Call:
## coxph(formula = Surv(overall_survival, deceased) ~ B_cell, data = data)
##
## coef exp(coef) se(coef) z p
## B_cell 0.4112 1.5086 0.2191 1.877 0.0606
##
## Likelihood ratio test=2.98 on 1 df, p=0.08412
## n= 144, number of events= 143
## (253 observations deleted due to missingness)
library(survival)
library(reactable)
signature_sq_matrix <- as.data.frame(signature_sq_matrix)
Alveolar_exp <- signature_sq_matrix$`B cell`
other_cells_exp <- rowMeans(signature_sq_matrix[, colnames(signature_sq_matrix) != "B_cell"])
Alveolar_lusc_filtered_genes <- signature_sq_matrix[(Alveolar_exp - other_cells_exp) / other_cells_exp > 0.90, ]
sq_bulk <- readRDS("~/final_results/lusc_bulk.rds")
Alveolar <- intersect(rownames(Alveolar_lusc_filtered_genes), rownames(sq_bulk))
Alveolar <- sq_bulk[Alveolar,]
Alveolar <- t(Alveolar)
Alveolar <- as.data.frame(Alveolar)
lusc_BULK_CLINIC <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
lusc_BULK_CLINIC <- as.data.frame(lusc_BULK_CLINIC)
select_col <- c("overall_survival", "vital_status")
clinic <- lusc_BULK_CLINIC[,select_col]
data <- cbind(clinic, Alveolar)
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
library(survival)
exclude_cols <- c("overall_survival", "vital_status", "deceased")
gene_cols <- setdiff(colnames(data), exclude_cols)
cox_results <- list()
for (gene in gene_cols) {
formula <- as.formula(paste("Surv(overall_survival, deceased) ~", gene))
cox_model <- coxph(formula, data = data)
cox_results[[gene]] <- summary(cox_model)
}
negative_genes <- c()
negative_coeffs <- c()
for (gene in names(cox_results)) {
coef_value <- cox_results[[gene]]$coefficients[1, "coef"]
if (!is.na(coef_value) && coef_value > 0) {
negative_genes <- c(negative_genes, gene)
negative_coeffs <- c(negative_coeffs, coef_value)
}
}
library(biomaRt)
ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl", version = 105)
gene_annotation <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol"),
filters = "ensembl_gene_id",
values = negative_genes,
mart = ensembl)
ort <- intersect(negative_genes, rownames(signature_sq_matrix))
sig <- signature_sq_matrix[ort, ]
B_cell_negative_gene_df <- data.frame(
ENSG_ID = negative_genes,
Gene_Name = gene_annotation$hgnc_symbol,
Coefficient = negative_coeffs,
Cell_Type = "B_cell"
)
B_cell_negative_gene_df$exp.level <- sig$`B cell`
reactable(B_cell_negative_gene_df)
library(survival)
library(survminer)
signature_sq_matrix <- readRDS("~/final_results/14_02_2025_signature_sq_matrix.rds")
sq <- readRDS("~/final_results/14_02_2025_deconvolution_sq.rds")
lusc_clinic_for_bulk <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
deconvolution_sq_scaled <- sq / rowSums(sq) * 100
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
res.cut <- surv_cutpoint(data, time = "overall_survival", event = "deceased",
variables = c("DC_mature"))
cutpoint <- res.cut[["cutpoint"]][["cutpoint"]]
data$DC_mature <- ifelse(data$DC_mature > cutpoint, "HIGH", "LOW")
fit_DC_mature <- survfit(Surv(overall_survival, deceased) ~ DC_mature, data = data)
ggsurvplot(fit_DC_mature,
data,
pval = TRUE,
risk.table = TRUE)
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
colnames(data) <- gsub("-", "_", colnames(data))
coxa <- coxph(Surv(overall_survival, deceased) ~ DC_mature, data = data)
coxa
## Call:
## coxph(formula = Surv(overall_survival, deceased) ~ DC_mature,
## data = data)
##
## coef exp(coef) se(coef) z p
## DC_mature 0.159 1.172 0.101 1.575 0.115
##
## Likelihood ratio test=2.31 on 1 df, p=0.1285
## n= 144, number of events= 143
## (253 observations deleted due to missingness)
library(survival)
library(reactable)
signature_sq_matrix <- as.data.frame(signature_sq_matrix)
Alveolar_exp <- signature_sq_matrix$`DC mature`
other_cells_exp <- rowMeans(signature_sq_matrix[, colnames(signature_sq_matrix) != "DC_mature"])
Alveolar_lusc_filtered_genes <- signature_sq_matrix[(Alveolar_exp - other_cells_exp) / other_cells_exp > 0.90, ]
sq_bulk <- readRDS("~/final_results/lusc_bulk.rds")
Alveolar <- intersect(rownames(Alveolar_lusc_filtered_genes), rownames(sq_bulk))
Alveolar <- sq_bulk[Alveolar,]
Alveolar <- t(Alveolar)
Alveolar <- as.data.frame(Alveolar)
lusc_BULK_CLINIC <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
lusc_BULK_CLINIC <- as.data.frame(lusc_BULK_CLINIC)
select_col <- c("overall_survival", "vital_status")
clinic <- lusc_BULK_CLINIC[,select_col]
data <- cbind(clinic, Alveolar)
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
library(survival)
exclude_cols <- c("overall_survival", "vital_status", "deceased")
gene_cols <- setdiff(colnames(data), exclude_cols)
cox_results <- list()
for (gene in gene_cols) {
formula <- as.formula(paste("Surv(overall_survival, deceased) ~", gene))
cox_model <- coxph(formula, data = data)
cox_results[[gene]] <- summary(cox_model)
}
negative_genes <- c()
negative_coeffs <- c()
for (gene in names(cox_results)) {
coef_value <- cox_results[[gene]]$coefficients[1, "coef"]
if (!is.na(coef_value) && coef_value > 0) {
negative_genes <- c(negative_genes, gene)
negative_coeffs <- c(negative_coeffs, coef_value)
}
}
library(biomaRt)
ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl", version = 105)
gene_annotation <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol"),
filters = "ensembl_gene_id",
values = negative_genes,
mart = ensembl)
ort <- intersect(negative_genes, rownames(signature_sq_matrix))
sig <- signature_sq_matrix[ort, ]
DC_mature_negative_gene_df <- data.frame(
ENSG_ID = negative_genes,
Gene_Name = gene_annotation$hgnc_symbol,
Coefficient = negative_coeffs,
Cell_Type = "DC_mature"
)
DC_mature_negative_gene_df$exp.level <- sig$`DC mature`
reactable(DC_mature_negative_gene_df)
library(survival)
library(survminer)
signature_sq_matrix <- readRDS("~/final_results/14_02_2025_signature_sq_matrix.rds")
sq <- readRDS("~/final_results/14_02_2025_deconvolution_sq.rds")
lusc_clinic_for_bulk <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
deconvolution_sq_scaled <- sq / rowSums(sq) * 100
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
res.cut <- surv_cutpoint(data, time = "overall_survival", event = "deceased",
variables = c("Fibroblast_alveolar"))
cutpoint <- res.cut[["cutpoint"]][["cutpoint"]]
data$Fibroblast_alveolar <- ifelse(data$Fibroblast_alveolar > cutpoint, "HIGH", "LOW")
fit_Fibroblast_alveolar <- survfit(Surv(overall_survival, deceased) ~ Fibroblast_alveolar, data = data)
ggsurvplot(fit_Fibroblast_alveolar,
data,
pval = TRUE,
risk.table = TRUE)
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
colnames(data) <- gsub("-", "_", colnames(data))
coxa <- coxph(Surv(overall_survival, deceased) ~ Fibroblast_alveolar, data = data)
coxa
## Call:
## coxph(formula = Surv(overall_survival, deceased) ~ Fibroblast_alveolar,
## data = data)
##
## coef exp(coef) se(coef) z p
## Fibroblast_alveolar 0.08755 1.09150 0.05749 1.523 0.128
##
## Likelihood ratio test=2.13 on 1 df, p=0.1448
## n= 144, number of events= 143
## (253 observations deleted due to missingness)
library(survival)
library(reactable)
signature_sq_matrix <- as.data.frame(signature_sq_matrix)
Alveolar_exp <- signature_sq_matrix$`Fibroblast alveolar`
other_cells_exp <- rowMeans(signature_sq_matrix[, colnames(signature_sq_matrix) != "Fibroblast_alveolar"])
Alveolar_lusc_filtered_genes <- signature_sq_matrix[(Alveolar_exp - other_cells_exp) / other_cells_exp > 0.90, ]
sq_bulk <- readRDS("~/final_results/lusc_bulk.rds")
Alveolar <- intersect(rownames(Alveolar_lusc_filtered_genes), rownames(sq_bulk))
Alveolar <- sq_bulk[Alveolar,]
Alveolar <- t(Alveolar)
Alveolar <- as.data.frame(Alveolar)
lusc_BULK_CLINIC <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
lusc_BULK_CLINIC <- as.data.frame(lusc_BULK_CLINIC)
select_col <- c("overall_survival", "vital_status")
clinic <- lusc_BULK_CLINIC[,select_col]
data <- cbind(clinic, Alveolar)
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
library(survival)
exclude_cols <- c("overall_survival", "vital_status", "deceased")
gene_cols <- setdiff(colnames(data), exclude_cols)
cox_results <- list()
for (gene in gene_cols) {
formula <- as.formula(paste("Surv(overall_survival, deceased) ~", gene))
cox_model <- coxph(formula, data = data)
cox_results[[gene]] <- summary(cox_model)
}
negative_genes <- c()
negative_coeffs <- c()
for (gene in names(cox_results)) {
coef_value <- cox_results[[gene]]$coefficients[1, "coef"]
if (!is.na(coef_value) && coef_value > 0) {
negative_genes <- c(negative_genes, gene)
negative_coeffs <- c(negative_coeffs, coef_value)
}
}
library(biomaRt)
ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl", version = 105)
gene_annotation <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol"),
filters = "ensembl_gene_id",
values = negative_genes,
mart = ensembl)
ort <- intersect(negative_genes, rownames(signature_sq_matrix))
sig <- signature_sq_matrix[ort, ]
Fibroblast_alveolar_negative_gene_df <- data.frame(
ENSG_ID = negative_genes,
Gene_Name = gene_annotation$hgnc_symbol,
Coefficient = negative_coeffs,
Cell_Type = "NK_cell_div"
)
Fibroblast_alveolar_negative_gene_df$exp.level <- sig$`Fibroblast alveolar`
reactable(Fibroblast_alveolar_negative_gene_df)
library(survival)
library(survminer)
signature_sq_matrix <- readRDS("~/final_results/14_02_2025_signature_sq_matrix.rds")
sq <- readRDS("~/final_results/14_02_2025_deconvolution_sq.rds")
lusc_clinic_for_bulk <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
deconvolution_sq_scaled <- sq / rowSums(sq) * 100
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
res.cut <- surv_cutpoint(data, time = "overall_survival", event = "deceased",
variables = c("Fibroblast_peribronchial"))
cutpoint <- res.cut[["cutpoint"]][["cutpoint"]]
data$Fibroblast_peribronchial <- ifelse(data$Fibroblast_peribronchial > cutpoint, "HIGH", "LOW")
fit_Fibroblast_peribronchial <- survfit(Surv(overall_survival, deceased) ~ Fibroblast_peribronchial, data = data)
ggsurvplot(fit_Fibroblast_peribronchial,
data,
pval = TRUE,
risk.table = TRUE)
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
colnames(data) <- gsub("-", "_", colnames(data))
coxa <- coxph(Surv(overall_survival, deceased) ~ Fibroblast_peribronchial, data = data)
coxa
## Call:
## coxph(formula = Surv(overall_survival, deceased) ~ Fibroblast_peribronchial,
## data = data)
##
## coef exp(coef) se(coef) z p
## Fibroblast_peribronchial 0.008238 1.008272 0.010770 0.765 0.444
##
## Likelihood ratio test=0.56 on 1 df, p=0.4534
## n= 144, number of events= 143
## (253 observations deleted due to missingness)
library(survival)
library(reactable)
signature_sq_matrix <- as.data.frame(signature_sq_matrix)
Alveolar_exp <- signature_sq_matrix$`Fibroblast peribronchial`
other_cells_exp <- rowMeans(signature_sq_matrix[, colnames(signature_sq_matrix) != "Fibroblast_peribronchial"])
Alveolar_lusc_filtered_genes <- signature_sq_matrix[(Alveolar_exp - other_cells_exp) / other_cells_exp > 0.90, ]
sq_bulk <- readRDS("~/final_results/lusc_bulk.rds")
Alveolar <- intersect(rownames(Alveolar_lusc_filtered_genes), rownames(sq_bulk))
Alveolar <- sq_bulk[Alveolar,]
Alveolar <- t(Alveolar)
Alveolar <- as.data.frame(Alveolar)
lusc_BULK_CLINIC <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
lusc_BULK_CLINIC <- as.data.frame(lusc_BULK_CLINIC)
select_col <- c("overall_survival", "vital_status")
clinic <- lusc_BULK_CLINIC[,select_col]
data <- cbind(clinic, Alveolar)
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
library(survival)
exclude_cols <- c("overall_survival", "vital_status", "deceased")
gene_cols <- setdiff(colnames(data), exclude_cols)
cox_results <- list()
for (gene in gene_cols) {
formula <- as.formula(paste("Surv(overall_survival, deceased) ~", gene))
cox_model <- coxph(formula, data = data)
cox_results[[gene]] <- summary(cox_model)
}
negative_genes <- c()
negative_coeffs <- c()
for (gene in names(cox_results)) {
coef_value <- cox_results[[gene]]$coefficients[1, "coef"]
if (!is.na(coef_value) && coef_value > 0) {
negative_genes <- c(negative_genes, gene)
negative_coeffs <- c(negative_coeffs, coef_value)
}
}
library(biomaRt)
ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl", version = 105)
gene_annotation <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol"),
filters = "ensembl_gene_id",
values = negative_genes,
mart = ensembl)
ort <- intersect(negative_genes, rownames(signature_sq_matrix))
sig <- signature_sq_matrix[ort, ]
Fibroblast_peribronchial_negative_gene_df <- data.frame(
ENSG_ID = negative_genes,
Gene_Name = gene_annotation$hgnc_symbol,
Coefficient = negative_coeffs,
Cell_Type = "Fibroblast_peribronchial"
)
Fibroblast_peribronchial_negative_gene_df$exp.level <- sig$`Fibroblast peribronchial`
reactable(Fibroblast_peribronchial_negative_gene_df)
library(survival)
library(survminer)
signature_sq_matrix <- readRDS("~/final_results/14_02_2025_signature_sq_matrix.rds")
sq <- readRDS("~/final_results/14_02_2025_deconvolution_sq.rds")
lusc_clinic_for_bulk <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
deconvolution_sq_scaled <- sq / rowSums(sq) * 100
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
res.cut <- surv_cutpoint(data, time = "overall_survival", event = "deceased",
variables = c("Macrophage"))
cutpoint <- res.cut[["cutpoint"]][["cutpoint"]]
data$Macrophage <- ifelse(data$Macrophage > cutpoint, "HIGH", "LOW")
fit_Macrophage <- survfit(Surv(overall_survival, deceased) ~ Macrophage, data = data)
ggsurvplot(fit_Macrophage,
data,
pval = TRUE,
risk.table = TRUE)
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
colnames(data) <- gsub("-", "_", colnames(data))
coxa <- coxph(Surv(overall_survival, deceased) ~ Macrophage, data = data)
coxa
## Call:
## coxph(formula = Surv(overall_survival, deceased) ~ Macrophage,
## data = data)
##
## coef exp(coef) se(coef) z p
## Macrophage 0.03033 1.03079 0.02377 1.276 0.202
##
## Likelihood ratio test=1.53 on 1 df, p=0.2154
## n= 144, number of events= 143
## (253 observations deleted due to missingness)
library(survival)
library(reactable)
signature_sq_matrix <- as.data.frame(signature_sq_matrix)
Alveolar_exp <- signature_sq_matrix$Macrophage
other_cells_exp <- rowMeans(signature_sq_matrix[, colnames(signature_sq_matrix) != "Macrophage"])
Alveolar_lusc_filtered_genes <- signature_sq_matrix[(Alveolar_exp - other_cells_exp) / other_cells_exp > 0.90, ]
sq_bulk <- readRDS("~/final_results/lusc_bulk.rds")
Alveolar <- intersect(rownames(Alveolar_lusc_filtered_genes), rownames(sq_bulk))
Alveolar <- sq_bulk[Alveolar,]
Alveolar <- t(Alveolar)
Alveolar <- as.data.frame(Alveolar)
lusc_BULK_CLINIC <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
lusc_BULK_CLINIC <- as.data.frame(lusc_BULK_CLINIC)
select_col <- c("overall_survival", "vital_status")
clinic <- lusc_BULK_CLINIC[,select_col]
data <- cbind(clinic, Alveolar)
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
library(survival)
exclude_cols <- c("overall_survival", "vital_status", "deceased")
gene_cols <- setdiff(colnames(data), exclude_cols)
cox_results <- list()
for (gene in gene_cols) {
formula <- as.formula(paste("Surv(overall_survival, deceased) ~", gene))
cox_model <- coxph(formula, data = data)
cox_results[[gene]] <- summary(cox_model)
}
negative_genes <- c()
negative_coeffs <- c()
for (gene in names(cox_results)) {
coef_value <- cox_results[[gene]]$coefficients[1, "coef"]
if (!is.na(coef_value) && coef_value > 0) {
negative_genes <- c(negative_genes, gene)
negative_coeffs <- c(negative_coeffs, coef_value)
}
}
library(biomaRt)
ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl", version = 105)
gene_annotation <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol"),
filters = "ensembl_gene_id",
values = negative_genes,
mart = ensembl)
ort <- intersect(negative_genes, rownames(signature_sq_matrix))
sig <- signature_sq_matrix[ort, ]
Macrophage_negative_gene_df <- data.frame(
ENSG_ID = negative_genes,
Gene_Name = gene_annotation$hgnc_symbol,
Coefficient = negative_coeffs,
Cell_Type = "Macrophage"
)
Macrophage_negative_gene_df$exp.level <- sig$Macrophage
reactable(Macrophage_negative_gene_df)
library(survival)
library(survminer)
signature_sq_matrix <- readRDS("~/final_results/14_02_2025_signature_sq_matrix.rds")
sq <- readRDS("~/final_results/14_02_2025_deconvolution_sq.rds")
lusc_clinic_for_bulk <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
deconvolution_sq_scaled <- sq / rowSums(sq) * 100
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
res.cut <- surv_cutpoint(data, time = "overall_survival", event = "deceased",
variables = c("Macrophage_alveolar"))
cutpoint <- res.cut[["cutpoint"]][["cutpoint"]]
data$Macrophage_alveolar <- ifelse(data$Macrophage_alveolar > cutpoint, "HIGH", "LOW")
fit_Macrophage_alveolar <- survfit(Surv(overall_survival, deceased) ~ Macrophage_alveolar, data = data)
ggsurvplot(fit_Macrophage_alveolar,
data,
pval = TRUE,
risk.table = TRUE)
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
colnames(data) <- gsub("-", "_", colnames(data))
coxa <- coxph(Surv(overall_survival, deceased) ~ Macrophage_alveolar, data = data)
coxa
## Call:
## coxph(formula = Surv(overall_survival, deceased) ~ Macrophage_alveolar,
## data = data)
##
## coef exp(coef) se(coef) z p
## Macrophage_alveolar 0.06940 1.07187 0.03033 2.288 0.0221
##
## Likelihood ratio test=4.5 on 1 df, p=0.03398
## n= 144, number of events= 143
## (253 observations deleted due to missingness)
library(survival)
library(reactable)
signature_sq_matrix <- as.data.frame(signature_sq_matrix)
Alveolar_exp <- signature_sq_matrix$`Macrophage alveolar`
other_cells_exp <- rowMeans(signature_sq_matrix[, colnames(signature_sq_matrix) != "Macrophage_alveolar"])
Alveolar_lusc_filtered_genes <- signature_sq_matrix[(Alveolar_exp - other_cells_exp) / other_cells_exp > 0.90, ]
sq_bulk <- readRDS("~/final_results/lusc_bulk.rds")
Alveolar <- intersect(rownames(Alveolar_lusc_filtered_genes), rownames(sq_bulk))
Alveolar <- sq_bulk[Alveolar,]
Alveolar <- t(Alveolar)
Alveolar <- as.data.frame(Alveolar)
lusc_BULK_CLINIC <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
lusc_BULK_CLINIC <- as.data.frame(lusc_BULK_CLINIC)
select_col <- c("overall_survival", "vital_status")
clinic <- lusc_BULK_CLINIC[,select_col]
data <- cbind(clinic, Alveolar)
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
library(survival)
exclude_cols <- c("overall_survival", "vital_status", "deceased")
gene_cols <- setdiff(colnames(data), exclude_cols)
cox_results <- list()
for (gene in gene_cols) {
formula <- as.formula(paste("Surv(overall_survival, deceased) ~", gene))
cox_model <- coxph(formula, data = data)
cox_results[[gene]] <- summary(cox_model)
}
negative_genes <- c()
negative_coeffs <- c()
for (gene in names(cox_results)) {
coef_value <- cox_results[[gene]]$coefficients[1, "coef"]
if (!is.na(coef_value) && coef_value > 0) {
negative_genes <- c(negative_genes, gene)
negative_coeffs <- c(negative_coeffs, coef_value)
}
}
library(biomaRt)
ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl", version = 105)
gene_annotation <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol"),
filters = "ensembl_gene_id",
values = negative_genes,
mart = ensembl)
ort <- intersect(negative_genes, rownames(signature_sq_matrix))
sig <- signature_sq_matrix[ort, ]
Macrophage_alveolar_negative_gene_df <- data.frame(
ENSG_ID = negative_genes,
Gene_Name = gene_annotation$hgnc_symbol,
Coefficient = negative_coeffs,
Cell_Type = "Macrophage_alveolar"
)
Macrophage_alveolar_negative_gene_df$exp.level <- sig$`Macrophage alveolar`
reactable(Macrophage_alveolar_negative_gene_df)
library(survival)
library(survminer)
signature_sq_matrix <- readRDS("~/final_results/14_02_2025_signature_sq_matrix.rds")
sq <- readRDS("~/final_results/14_02_2025_deconvolution_sq.rds")
lusc_clinic_for_bulk <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
deconvolution_sq_scaled <- sq / rowSums(sq) * 100
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
res.cut <- surv_cutpoint(data, time = "overall_survival", event = "deceased",
variables = c("Mast_cell"))
cutpoint <- res.cut[["cutpoint"]][["cutpoint"]]
data$Mast_cell <- ifelse(data$Mast_cell > cutpoint, "HIGH", "LOW")
fit_Mast_cell <- survfit(Surv(overall_survival, deceased) ~ Mast_cell, data = data)
ggsurvplot(fit_Mast_cell,
data,
pval = TRUE,
risk.table = TRUE)
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
colnames(data) <- gsub("-", "_", colnames(data))
coxa <- coxph(Surv(overall_survival, deceased) ~ Mast_cell, data = data)
coxa
## Call:
## coxph(formula = Surv(overall_survival, deceased) ~ Mast_cell,
## data = data)
##
## coef exp(coef) se(coef) z p
## Mast_cell 0.06052 1.06239 0.45609 0.133 0.894
##
## Likelihood ratio test=0.02 on 1 df, p=0.895
## n= 144, number of events= 143
## (253 observations deleted due to missingness)
library(survival)
library(reactable)
signature_sq_matrix <- as.data.frame(signature_sq_matrix)
Alveolar_exp <- signature_sq_matrix$`NK cell dividing`
other_cells_exp <- rowMeans(signature_sq_matrix[, colnames(signature_sq_matrix) != "Mast_cell"])
Alveolar_lusc_filtered_genes <- signature_sq_matrix[(Alveolar_exp - other_cells_exp) / other_cells_exp > 0.90, ]
sq_bulk <- readRDS("~/final_results/lusc_bulk.rds")
Alveolar <- intersect(rownames(Alveolar_lusc_filtered_genes), rownames(sq_bulk))
Alveolar <- sq_bulk[Alveolar,]
Alveolar <- t(Alveolar)
Alveolar <- as.data.frame(Alveolar)
lusc_BULK_CLINIC <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
lusc_BULK_CLINIC <- as.data.frame(lusc_BULK_CLINIC)
select_col <- c("overall_survival", "vital_status")
clinic <- lusc_BULK_CLINIC[,select_col]
data <- cbind(clinic, Alveolar)
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
library(survival)
exclude_cols <- c("overall_survival", "vital_status", "deceased")
gene_cols <- setdiff(colnames(data), exclude_cols)
cox_results <- list()
for (gene in gene_cols) {
formula <- as.formula(paste("Surv(overall_survival, deceased) ~", gene))
cox_model <- coxph(formula, data = data)
cox_results[[gene]] <- summary(cox_model)
}
negative_genes <- c()
negative_coeffs <- c()
for (gene in names(cox_results)) {
coef_value <- cox_results[[gene]]$coefficients[1, "coef"]
if (!is.na(coef_value) && coef_value > 0) {
negative_genes <- c(negative_genes, gene)
negative_coeffs <- c(negative_coeffs, coef_value)
}
}
library(biomaRt)
ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl", version = 105)
gene_annotation <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol"),
filters = "ensembl_gene_id",
values = negative_genes,
mart = ensembl)
ort <- intersect(negative_genes, rownames(signature_sq_matrix))
sig <- signature_sq_matrix[ort, ]
Mast_cell_negative_gene_df <- data.frame(
ENSG_ID = negative_genes,
Gene_Name = gene_annotation$hgnc_symbol,
Coefficient = negative_coeffs,
Cell_Type = "Mast_cell"
)
Mast_cell_negative_gene_df$exp.level <- sig$`Mast cell`
reactable(Mast_cell_negative_gene_df)
library(survival)
library(survminer)
signature_sq_matrix <- readRDS("~/final_results/14_02_2025_signature_sq_matrix.rds")
sq <- readRDS("~/final_results/14_02_2025_deconvolution_sq.rds")
lusc_clinic_for_bulk <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
deconvolution_sq_scaled <- sq / rowSums(sq) * 100
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
res.cut <- surv_cutpoint(data, time = "overall_survival", event = "deceased",
variables = c("Monocyte_classical"))
cutpoint <- res.cut[["cutpoint"]][["cutpoint"]]
data$Monocyte_classical <- ifelse(data$Monocyte_classical > cutpoint, "HIGH", "LOW")
fit_Monocyte_classical <- survfit(Surv(overall_survival, deceased) ~ Monocyte_classical, data = data)
ggsurvplot(fit_Monocyte_classical,
data,
pval = TRUE,
risk.table = TRUE)
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
colnames(data) <- gsub("-", "_", colnames(data))
coxa <- coxph(Surv(overall_survival, deceased) ~ Monocyte_classical, data = data)
coxa
## Call:
## coxph(formula = Surv(overall_survival, deceased) ~ Monocyte_classical,
## data = data)
##
## coef exp(coef) se(coef) z p
## Monocyte_classical 0.2020 1.2238 0.1278 1.581 0.114
##
## Likelihood ratio test=2.27 on 1 df, p=0.1321
## n= 144, number of events= 143
## (253 observations deleted due to missingness)
library(survival)
library(reactable)
signature_sq_matrix <- as.data.frame(signature_sq_matrix)
Alveolar_exp <- signature_sq_matrix$`Monocyte classical`
other_cells_exp <- rowMeans(signature_sq_matrix[, colnames(signature_sq_matrix) != "Monocyte_classical"])
Alveolar_lusc_filtered_genes <- signature_sq_matrix[(Alveolar_exp - other_cells_exp) / other_cells_exp > 0.90, ]
sq_bulk <- readRDS("~/final_results/lusc_bulk.rds")
Alveolar <- intersect(rownames(Alveolar_lusc_filtered_genes), rownames(sq_bulk))
Alveolar <- sq_bulk[Alveolar,]
Alveolar <- t(Alveolar)
Alveolar <- as.data.frame(Alveolar)
lusc_BULK_CLINIC <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
lusc_BULK_CLINIC <- as.data.frame(lusc_BULK_CLINIC)
select_col <- c("overall_survival", "vital_status")
clinic <- lusc_BULK_CLINIC[,select_col]
data <- cbind(clinic, Alveolar)
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
library(survival)
exclude_cols <- c("overall_survival", "vital_status", "deceased")
gene_cols <- setdiff(colnames(data), exclude_cols)
cox_results <- list()
for (gene in gene_cols) {
formula <- as.formula(paste("Surv(overall_survival, deceased) ~", gene))
cox_model <- coxph(formula, data = data)
cox_results[[gene]] <- summary(cox_model)
}
negative_genes <- c()
negative_coeffs <- c()
for (gene in names(cox_results)) {
coef_value <- cox_results[[gene]]$coefficients[1, "coef"]
if (!is.na(coef_value) && coef_value > 0) {
negative_genes <- c(negative_genes, gene)
negative_coeffs <- c(negative_coeffs, coef_value)
}
}
library(biomaRt)
ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl", version = 105)
gene_annotation <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol"),
filters = "ensembl_gene_id",
values = negative_genes,
mart = ensembl)
ort <- intersect(negative_genes, rownames(signature_sq_matrix))
sig <- signature_sq_matrix[ort, ]
Monocyte_classical_negative_gene_df <- data.frame(
ENSG_ID = negative_genes,
Gene_Name = gene_annotation$hgnc_symbol,
Coefficient = negative_coeffs,
Cell_Type = "Monocyte_classical"
)
Monocyte_classical_negative_gene_df$exp.level <- sig$`Monocyte classical`
reactable(Monocyte_classical_negative_gene_df)
library(survival)
library(survminer)
signature_sq_matrix <- readRDS("~/final_results/14_02_2025_signature_sq_matrix.rds")
sq <- readRDS("~/final_results/14_02_2025_deconvolution_sq.rds")
lusc_clinic_for_bulk <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
deconvolution_sq_scaled <- sq / rowSums(sq) * 100
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
res.cut <- surv_cutpoint(data, time = "overall_survival", event = "deceased",
variables = c("Myeloid_dividing"))
cutpoint <- res.cut[["cutpoint"]][["cutpoint"]]
data$Myeloid_dividing <- ifelse(data$Myeloid_dividing > cutpoint, "HIGH", "LOW")
fit_Myeloid_dividing <- survfit(Surv(overall_survival, deceased) ~ Myeloid_dividing, data = data)
ggsurvplot(fit_Myeloid_dividing,
data,
pval = TRUE,
risk.table = TRUE)
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
colnames(data) <- gsub("-", "_", colnames(data))
coxa <- coxph(Surv(overall_survival, deceased) ~ Myeloid_dividing, data = data)
coxa
## Call:
## coxph(formula = Surv(overall_survival, deceased) ~ Myeloid_dividing,
## data = data)
##
## coef exp(coef) se(coef) z p
## Myeloid_dividing 0.02782 1.02821 0.01958 1.421 0.155
##
## Likelihood ratio test=1.94 on 1 df, p=0.1637
## n= 144, number of events= 143
## (253 observations deleted due to missingness)
library(survival)
library(reactable)
signature_sq_matrix <- as.data.frame(signature_sq_matrix)
Alveolar_exp <- signature_sq_matrix$`Myeloid dividing`
other_cells_exp <- rowMeans(signature_sq_matrix[, colnames(signature_sq_matrix) != "Myeloid_dividing"])
Alveolar_lusc_filtered_genes <- signature_sq_matrix[(Alveolar_exp - other_cells_exp) / other_cells_exp > 0.90, ]
sq_bulk <- readRDS("~/final_results/lusc_bulk.rds")
Alveolar <- intersect(rownames(Alveolar_lusc_filtered_genes), rownames(sq_bulk))
Alveolar <- sq_bulk[Alveolar,]
Alveolar <- t(Alveolar)
Alveolar <- as.data.frame(Alveolar)
lusc_BULK_CLINIC <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
lusc_BULK_CLINIC <- as.data.frame(lusc_BULK_CLINIC)
select_col <- c("overall_survival", "vital_status")
clinic <- lusc_BULK_CLINIC[,select_col]
data <- cbind(clinic, Alveolar)
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
library(survival)
exclude_cols <- c("overall_survival", "vital_status", "deceased")
gene_cols <- setdiff(colnames(data), exclude_cols)
cox_results <- list()
for (gene in gene_cols) {
formula <- as.formula(paste("Surv(overall_survival, deceased) ~", gene))
cox_model <- coxph(formula, data = data)
cox_results[[gene]] <- summary(cox_model)
}
negative_genes <- c()
negative_coeffs <- c()
for (gene in names(cox_results)) {
coef_value <- cox_results[[gene]]$coefficients[1, "coef"]
if (!is.na(coef_value) && coef_value > 0) {
negative_genes <- c(negative_genes, gene)
negative_coeffs <- c(negative_coeffs, coef_value)
}
}
library(biomaRt)
ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl", version = 105)
gene_annotation <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol"),
filters = "ensembl_gene_id",
values = negative_genes,
mart = ensembl)
ort <- intersect(negative_genes, rownames(signature_sq_matrix))
sig <- signature_sq_matrix[ort, ]
Myeloid_dividing_negative_gene_df <- data.frame(
ENSG_ID = negative_genes,
Gene_Name = gene_annotation$hgnc_symbol,
Coefficient = negative_coeffs,
Cell_Type = "Myeloid_dividing"
)
Myeloid_dividing_negative_gene_df$exp.level <- sig$`Myeloid dividing`
reactable(Myeloid_dividing_negative_gene_df)
library(survival)
library(survminer)
signature_sq_matrix <- readRDS("~/final_results/14_02_2025_signature_sq_matrix.rds")
sq <- readRDS("~/final_results/14_02_2025_deconvolution_sq.rds")
lusc_clinic_for_bulk <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
deconvolution_sq_scaled <- sq / rowSums(sq) * 100
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
res.cut <- surv_cutpoint(data, time = "overall_survival", event = "deceased",
variables = c("pDC"))
cutpoint <- res.cut[["cutpoint"]][["cutpoint"]]
data$pDC <- ifelse(data$pDC > cutpoint, "HIGH", "LOW")
fit_pDC <- survfit(Surv(overall_survival, deceased) ~ pDC, data = data)
ggsurvplot(fit_pDC,
data,
pval = TRUE,
risk.table = TRUE)
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
colnames(data) <- gsub("-", "_", colnames(data))
coxa <- coxph(Surv(overall_survival, deceased) ~ pDC, data = data)
coxa
## Call:
## coxph(formula = Surv(overall_survival, deceased) ~ pDC, data = data)
##
## coef exp(coef) se(coef) z p
## pDC 0.3391 1.4036 0.7356 0.461 0.645
##
## Likelihood ratio test=0.2 on 1 df, p=0.6575
## n= 144, number of events= 143
## (253 observations deleted due to missingness)
library(survival)
library(reactable)
signature_sq_matrix <- as.data.frame(signature_sq_matrix)
Alveolar_exp <- signature_sq_matrix$pDC
other_cells_exp <- rowMeans(signature_sq_matrix[, colnames(signature_sq_matrix) != "pDC"])
Alveolar_lusc_filtered_genes <- signature_sq_matrix[(Alveolar_exp - other_cells_exp) / other_cells_exp > 0.90, ]
sq_bulk <- readRDS("~/final_results/lusc_bulk.rds")
Alveolar <- intersect(rownames(Alveolar_lusc_filtered_genes), rownames(sq_bulk))
Alveolar <- sq_bulk[Alveolar,]
Alveolar <- t(Alveolar)
Alveolar <- as.data.frame(Alveolar)
lusc_BULK_CLINIC <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
lusc_BULK_CLINIC <- as.data.frame(lusc_BULK_CLINIC)
select_col <- c("overall_survival", "vital_status")
clinic <- lusc_BULK_CLINIC[,select_col]
data <- cbind(clinic, Alveolar)
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
library(survival)
exclude_cols <- c("overall_survival", "vital_status", "deceased")
gene_cols <- setdiff(colnames(data), exclude_cols)
cox_results <- list()
for (gene in gene_cols) {
formula <- as.formula(paste("Surv(overall_survival, deceased) ~", gene))
cox_model <- coxph(formula, data = data)
cox_results[[gene]] <- summary(cox_model)
}
negative_genes <- c()
negative_coeffs <- c()
for (gene in names(cox_results)) {
coef_value <- cox_results[[gene]]$coefficients[1, "coef"]
if (!is.na(coef_value) && coef_value > 0) {
negative_genes <- c(negative_genes, gene)
negative_coeffs <- c(negative_coeffs, coef_value)
}
}
library(biomaRt)
ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl", version = 105)
gene_annotation <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol"),
filters = "ensembl_gene_id",
values = negative_genes,
mart = ensembl)
ort <- intersect(negative_genes, rownames(signature_sq_matrix))
sig <- signature_sq_matrix[ort, ]
pDC_negative_gene_df <- data.frame(
ENSG_ID = negative_genes,
Gene_Name = gene_annotation$hgnc_symbol,
Coefficient = negative_coeffs,
Cell_Type = "pDC"
)
pDC_negative_gene_df$exp.level <- sig$pDC
reactable(pDC_negative_gene_df)
library(survival)
library(survminer)
signature_sq_matrix <- readRDS("~/final_results/14_02_2025_signature_sq_matrix.rds")
sq <- readRDS("~/final_results/14_02_2025_deconvolution_sq.rds")
lusc_clinic_for_bulk <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
deconvolution_sq_scaled <- sq / rowSums(sq) * 100
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
res.cut <- surv_cutpoint(data, time = "overall_survival", event = "deceased",
variables = c("Plasma_cell"))
cutpoint <- res.cut[["cutpoint"]][["cutpoint"]]
data$Plasma_cell <- ifelse(data$Plasma_cell > cutpoint, "HIGH", "LOW")
fit_Plasma_cell <- survfit(Surv(overall_survival, deceased) ~ Plasma_cell, data = data)
ggsurvplot(fit_Plasma_cell,
data,
pval = TRUE,
risk.table = TRUE)
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
colnames(data) <- gsub("-", "_", colnames(data))
coxa <- coxph(Surv(overall_survival, deceased) ~ Plasma_cell, data = data)
coxa
## Call:
## coxph(formula = Surv(overall_survival, deceased) ~ Plasma_cell,
## data = data)
##
## coef exp(coef) se(coef) z p
## Plasma_cell 0.005021 1.005033 0.005369 0.935 0.35
##
## Likelihood ratio test=0.85 on 1 df, p=0.3569
## n= 144, number of events= 143
## (253 observations deleted due to missingness)
library(survival)
library(reactable)
signature_sq_matrix <- as.data.frame(signature_sq_matrix)
Alveolar_exp <- signature_sq_matrix$`Plasma cell`
other_cells_exp <- rowMeans(signature_sq_matrix[, colnames(signature_sq_matrix) != "Plasma_cell"])
Alveolar_lusc_filtered_genes <- signature_sq_matrix[(Alveolar_exp - other_cells_exp) / other_cells_exp > 0.90, ]
sq_bulk <- readRDS("~/final_results/lusc_bulk.rds")
Alveolar <- intersect(rownames(Alveolar_lusc_filtered_genes), rownames(sq_bulk))
Alveolar <- sq_bulk[Alveolar,]
Alveolar <- t(Alveolar)
Alveolar <- as.data.frame(Alveolar)
lusc_BULK_CLINIC <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
lusc_BULK_CLINIC <- as.data.frame(lusc_BULK_CLINIC)
select_col <- c("overall_survival", "vital_status")
clinic <- lusc_BULK_CLINIC[,select_col]
data <- cbind(clinic, Alveolar)
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
library(survival)
exclude_cols <- c("overall_survival", "vital_status", "deceased")
gene_cols <- setdiff(colnames(data), exclude_cols)
cox_results <- list()
for (gene in gene_cols) {
formula <- as.formula(paste("Surv(overall_survival, deceased) ~", gene))
cox_model <- coxph(formula, data = data)
cox_results[[gene]] <- summary(cox_model)
}
negative_genes <- c()
negative_coeffs <- c()
for (gene in names(cox_results)) {
coef_value <- cox_results[[gene]]$coefficients[1, "coef"]
if (!is.na(coef_value) && coef_value > 0) {
negative_genes <- c(negative_genes, gene)
negative_coeffs <- c(negative_coeffs, coef_value)
}
}
library(biomaRt)
ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl", version = 105)
gene_annotation <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol"),
filters = "ensembl_gene_id",
values = negative_genes,
mart = ensembl)
ort <- intersect(negative_genes, rownames(signature_sq_matrix))
sig <- signature_sq_matrix[ort, ]
Plasma_cell_negative_gene_df <- data.frame(
ENSG_ID = negative_genes,
Gene_Name = gene_annotation$hgnc_symbol,
Coefficient = negative_coeffs,
Cell_Type = "Plasma_cell"
)
Plasma_cell_negative_gene_df$exp.level <- sig$`Plasma cell`
reactable(Plasma_cell_negative_gene_df)
library(survival)
library(survminer)
signature_sq_matrix <- readRDS("~/final_results/14_02_2025_signature_sq_matrix.rds")
sq <- readRDS("~/final_results/14_02_2025_deconvolution_sq.rds")
lusc_clinic_for_bulk <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
deconvolution_sq_scaled <- sq / rowSums(sq) * 100
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
res.cut <- surv_cutpoint(data, time = "overall_survival", event = "deceased",
variables = c("T_cell_CD8_naive"))
cutpoint <- res.cut[["cutpoint"]][["cutpoint"]]
data$T_cell_CD8_naive <- ifelse(data$T_cell_CD8_naive > cutpoint, "HIGH", "LOW")
fit_T_cell_CD8_naive <- survfit(Surv(overall_survival, deceased) ~ T_cell_CD8_naive, data = data)
ggsurvplot(fit_T_cell_CD8_naive,
data,
pval = TRUE,
risk.table = TRUE)
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
colnames(data) <- gsub("-", "_", colnames(data))
coxa <- coxph(Surv(overall_survival, deceased) ~ T_cell_CD8_naive, data = data)
coxa
## Call:
## coxph(formula = Surv(overall_survival, deceased) ~ T_cell_CD8_naive,
## data = data)
##
## coef exp(coef) se(coef) z p
## T_cell_CD8_naive 0.12939 1.13813 0.06223 2.079 0.0376
##
## Likelihood ratio test=4.13 on 1 df, p=0.04213
## n= 144, number of events= 143
## (253 observations deleted due to missingness)
library(survival)
library(reactable)
signature_sq_matrix <- as.data.frame(signature_sq_matrix)
Alveolar_exp <- signature_sq_matrix$`T cell CD8 naive`
other_cells_exp <- rowMeans(signature_sq_matrix[, colnames(signature_sq_matrix) != "T_cell_CD8_naive"])
Alveolar_lusc_filtered_genes <- signature_sq_matrix[(Alveolar_exp - other_cells_exp) / other_cells_exp > 0.90, ]
sq_bulk <- readRDS("~/final_results/lusc_bulk.rds")
Alveolar <- intersect(rownames(Alveolar_lusc_filtered_genes), rownames(sq_bulk))
Alveolar <- sq_bulk[Alveolar,]
Alveolar <- t(Alveolar)
Alveolar <- as.data.frame(Alveolar)
lusc_BULK_CLINIC <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
lusc_BULK_CLINIC <- as.data.frame(lusc_BULK_CLINIC)
select_col <- c("overall_survival", "vital_status")
clinic <- lusc_BULK_CLINIC[,select_col]
data <- cbind(clinic, Alveolar)
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
library(survival)
exclude_cols <- c("overall_survival", "vital_status", "deceased")
gene_cols <- setdiff(colnames(data), exclude_cols)
cox_results <- list()
for (gene in gene_cols) {
formula <- as.formula(paste("Surv(overall_survival, deceased) ~", gene))
cox_model <- coxph(formula, data = data)
cox_results[[gene]] <- summary(cox_model)
}
negative_genes <- c()
negative_coeffs <- c()
for (gene in names(cox_results)) {
coef_value <- cox_results[[gene]]$coefficients[1, "coef"]
if (!is.na(coef_value) && coef_value > 0) {
negative_genes <- c(negative_genes, gene)
negative_coeffs <- c(negative_coeffs, coef_value)
}
}
library(biomaRt)
ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl", version = 105)
gene_annotation <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol"),
filters = "ensembl_gene_id",
values = negative_genes,
mart = ensembl)
ort <- intersect(negative_genes, rownames(signature_sq_matrix))
sig <- signature_sq_matrix[ort, ]
T_cell_CD8_naive_negative_gene_df <- data.frame(
ENSG_ID = negative_genes,
Gene_Name = gene_annotation$hgnc_symbol,
Coefficient = negative_coeffs,
Cell_Type = "T_cell_CD8_naive"
)
T_cell_CD8_naive_negative_gene_df$exp.level <- sig$`T cell CD8 naive`
reactable(T_cell_CD8_naive_negative_gene_df)
library(survival)
library(survminer)
signature_sq_matrix <- readRDS("~/final_results/14_02_2025_signature_sq_matrix.rds")
sq <- readRDS("~/final_results/14_02_2025_deconvolution_sq.rds")
lusc_clinic_for_bulk <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
deconvolution_sq_scaled <- sq / rowSums(sq) * 100
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
res.cut <- surv_cutpoint(data, time = "overall_survival", event = "deceased",
variables = c("T_cell_CD8_terminally_exhausted"))
cutpoint <- res.cut[["cutpoint"]][["cutpoint"]]
data$T_cell_CD8_terminally_exhausted <- ifelse(data$T_cell_CD8_terminally_exhausted > cutpoint, "HIGH", "LOW")
fit_T_cell_CD8_terminally_exhausted <- survfit(Surv(overall_survival, deceased) ~ T_cell_CD8_terminally_exhausted, data = data)
ggsurvplot(fit_T_cell_CD8_terminally_exhausted,
data,
pval = TRUE,
risk.table = TRUE)
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
colnames(data) <- gsub("-", "_", colnames(data))
coxa <- coxph(Surv(overall_survival, deceased) ~ T_cell_CD8_terminally_exhausted, data = data)
coxa
## Call:
## coxph(formula = Surv(overall_survival, deceased) ~ T_cell_CD8_terminally_exhausted,
## data = data)
##
## coef exp(coef) se(coef) z p
## T_cell_CD8_terminally_exhausted 0.05342 1.05488 0.04410 1.211 0.226
##
## Likelihood ratio test=1.4 on 1 df, p=0.237
## n= 144, number of events= 143
## (253 observations deleted due to missingness)
library(survival)
library(reactable)
signature_sq_matrix <- as.data.frame(signature_sq_matrix)
Alveolar_exp <- signature_sq_matrix$`T cell CD8 terminally exhausted`
other_cells_exp <- rowMeans(signature_sq_matrix[, colnames(signature_sq_matrix) != "T_cell_CD8_terminally_exhausted"])
Alveolar_lusc_filtered_genes <- signature_sq_matrix[(Alveolar_exp - other_cells_exp) / other_cells_exp > 0.90, ]
sq_bulk <- readRDS("~/final_results/lusc_bulk.rds")
Alveolar <- intersect(rownames(Alveolar_lusc_filtered_genes), rownames(sq_bulk))
Alveolar <- sq_bulk[Alveolar,]
Alveolar <- t(Alveolar)
Alveolar <- as.data.frame(Alveolar)
lusc_BULK_CLINIC <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
lusc_BULK_CLINIC <- as.data.frame(lusc_BULK_CLINIC)
select_col <- c("overall_survival", "vital_status")
clinic <- lusc_BULK_CLINIC[,select_col]
data <- cbind(clinic, Alveolar)
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
library(survival)
exclude_cols <- c("overall_survival", "vital_status", "deceased")
gene_cols <- setdiff(colnames(data), exclude_cols)
cox_results <- list()
for (gene in gene_cols) {
formula <- as.formula(paste("Surv(overall_survival, deceased) ~", gene))
cox_model <- coxph(formula, data = data)
cox_results[[gene]] <- summary(cox_model)
}
negative_genes <- c()
negative_coeffs <- c()
for (gene in names(cox_results)) {
coef_value <- cox_results[[gene]]$coefficients[1, "coef"]
if (!is.na(coef_value) && coef_value > 0) {
negative_genes <- c(negative_genes, gene)
negative_coeffs <- c(negative_coeffs, coef_value)
}
}
library(biomaRt)
ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl", version = 105)
gene_annotation <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol"),
filters = "ensembl_gene_id",
values = negative_genes,
mart = ensembl)
ort <- intersect(negative_genes, rownames(signature_sq_matrix))
sig <- signature_sq_matrix[ort, ]
T_cell_CD8_terminally_exhausted_negative_gene_df <- data.frame(
ENSG_ID = negative_genes,
Gene_Name = gene_annotation$hgnc_symbol,
Coefficient = negative_coeffs,
Cell_Type = "T_cell_CD8_terminally_exhausted"
)
T_cell_CD8_terminally_exhausted_negative_gene_df$exp.level <- sig$`T cell CD8 terminally exhausted`
reactable(T_cell_CD8_terminally_exhausted_negative_gene_df)
library(survival)
library(survminer)
signature_sq_matrix <- readRDS("~/final_results/14_02_2025_signature_sq_matrix.rds")
sq <- readRDS("~/final_results/14_02_2025_deconvolution_sq.rds")
lusc_clinic_for_bulk <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
deconvolution_sq_scaled <- sq / rowSums(sq) * 100
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
res.cut <- surv_cutpoint(data, time = "overall_survival", event = "deceased",
variables = c("T_cell_regulatory"))
cutpoint <- res.cut[["cutpoint"]][["cutpoint"]]
data$T_cell_regulatory <- ifelse(data$T_cell_regulatory > cutpoint, "HIGH", "LOW")
fit_T_cell_regulatory <- survfit(Surv(overall_survival, deceased) ~ T_cell_regulatory, data = data)
ggsurvplot(fit_T_cell_regulatory,
data,
pval = TRUE,
risk.table = TRUE)
data <- cbind(lusc_clinic_for_bulk, deconvolution_sq_scaled)
colnames(data) <- gsub(" ", "_", colnames(data))
colnames(data) <- gsub("-", "_", colnames(data))
coxa <- coxph(Surv(overall_survival, deceased) ~ T_cell_regulatory, data = data)
coxa
## Call:
## coxph(formula = Surv(overall_survival, deceased) ~ T_cell_regulatory,
## data = data)
##
## coef exp(coef) se(coef) z p
## T_cell_regulatory 0.16993 1.18522 0.09721 1.748 0.0804
##
## Likelihood ratio test=2.84 on 1 df, p=0.0919
## n= 144, number of events= 143
## (253 observations deleted due to missingness)
library(survival)
library(reactable)
signature_sq_matrix <- as.data.frame(signature_sq_matrix)
Alveolar_exp <- signature_sq_matrix$`T cell regulatory`
other_cells_exp <- rowMeans(signature_sq_matrix[, colnames(signature_sq_matrix) != "T_cell_regulatory"])
Alveolar_lusc_filtered_genes <- signature_sq_matrix[(Alveolar_exp - other_cells_exp) / other_cells_exp > 0.90, ]
sq_bulk <- readRDS("~/final_results/lusc_bulk.rds")
Alveolar <- intersect(rownames(Alveolar_lusc_filtered_genes), rownames(sq_bulk))
Alveolar <- sq_bulk[Alveolar,]
Alveolar <- t(Alveolar)
Alveolar <- as.data.frame(Alveolar)
lusc_BULK_CLINIC <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
lusc_BULK_CLINIC <- as.data.frame(lusc_BULK_CLINIC)
select_col <- c("overall_survival", "vital_status")
clinic <- lusc_BULK_CLINIC[,select_col]
data <- cbind(clinic, Alveolar)
data$deceased <- ifelse(data$vital_status == "Alive", FALSE, TRUE)
library(survival)
exclude_cols <- c("overall_survival", "vital_status", "deceased")
gene_cols <- setdiff(colnames(data), exclude_cols)
cox_results <- list()
for (gene in gene_cols) {
formula <- as.formula(paste("Surv(overall_survival, deceased) ~", gene))
cox_model <- coxph(formula, data = data)
cox_results[[gene]] <- summary(cox_model)
}
negative_genes <- c()
negative_coeffs <- c()
for (gene in names(cox_results)) {
coef_value <- cox_results[[gene]]$coefficients[1, "coef"]
if (!is.na(coef_value) && coef_value > 0) {
negative_genes <- c(negative_genes, gene)
negative_coeffs <- c(negative_coeffs, coef_value)
}
}
library(biomaRt)
ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl", version = 105)
gene_annotation <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol"),
filters = "ensembl_gene_id",
values = negative_genes,
mart = ensembl)
ort <- intersect(negative_genes, rownames(signature_sq_matrix))
sig <- signature_sq_matrix[ort, ]
T_cell_regulatory_negative_gene_df <- data.frame(
ENSG_ID = negative_genes,
Gene_Name = gene_annotation$hgnc_symbol,
Coefficient = negative_coeffs,
Cell_Type = "T_cell_regulatory"
)
T_cell_regulatory_negative_gene_df$exp.level <- sig$`T cell regulatory`
reactable(T_cell_regulatory_negative_gene_df)
library(survival)
library(survminer)
library(ggplot2)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:biomaRt':
##
## select
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
df_list <- list(
NK_cell_div_negative_gene_df,
T_cell_CD8_dividing_negative_gene_df,
B_cell_negative_gene_df,
DC_mature_negative_gene_df,
Fibroblast_alveolar_negative_gene_df,
Fibroblast_peribronchial_negative_gene_df,
Macrophage_alveolar_negative_gene_df,
Macrophage_negative_gene_df,
Mast_cell_negative_gene_df,
Monocyte_classical_negative_gene_df,
Myeloid_dividing_negative_gene_df,
pDC_negative_gene_df,
Plasma_cell_negative_gene_df,
T_cell_CD8_naive_negative_gene_df,
T_cell_CD8_terminally_exhausted_negative_gene_df,
T_cell_regulatory_negative_gene_df)
df_names <- c(
"NK_cell_div_negative_gene_df",
"T_cell_CD8_dividing_negative_gene_df",
"B_cell_negative_gene_df",
"DC_mature_negative_gene_df",
"Fibroblast_alveolar_negative_gene_df",
"Fibroblast_peribronchial_negative_gene_df",
"Macrophage_alveolar_negative_gene_df",
"Macrophage_negative_gene_df",
"Mast_cell_negative_gene_df",
"Monocyte_classical_negative_gene_df",
"Myeloid_dividing_negative_gene_df",
"pDC_negative_gene_df",
"Plasma_cell_negative_gene_df",
"T_cell_CD8_naive_negative_gene_df",
"T_cell_CD8_terminally_exhausted_negative_gene_df",
"T_cell_regulatory_negative_gene_df"
)
library(dplyr)
select <- unique(unlist(lapply(df_list, function(df) {
if (!is.null(df)) {
df %>%
arrange(desc(exp.level)) %>%
head(5) %>%
pull(ENSG_ID)
}
})))
sq <- readRDS("~/final_results/lusc_bulk.rds") |> as.data.frame() |> t() |> as.data.frame()
lusc_clinic_for_bulk <- readRDS("~/final_results/lusc_clinic_for_bulk.rds")
data <- cbind(lusc_clinic_for_bulk, sq)
data$deceased <- ifelse(data$vital_status == "Alive", 0, 1)
ensembl_ids <- select
mart <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl", version = 105)
gene_map <- getBM(
attributes = c("ensembl_gene_id", "external_gene_name"),
filters = "ensembl_gene_id",
values = ensembl_ids,
mart = mart
)
colnames(data)[colnames(data) %in% gene_map$ensembl_gene_id] <- gene_map$external_gene_name
data <- cbind(lusc_clinic_for_bulk, data[, gene_map$external_gene_name])
data$overall_survival <- as.numeric(data$overall_survival)
univ_cox <- function(variable) {
formula <- reformulate(variable, response = "Surv(overall_survival, deceased)")
cox_model <- coxph(formula, data = data)
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])
)
}
colnames(data) <- gsub("-", "_", colnames(data))
univ_results <- lapply(colnames(data[, 8:ncol(data)]), univ_cox)
univ_results_df <- do.call(rbind, univ_results) |> as.data.frame()
rownames(univ_results_df) <- colnames(data[, 8:ncol(data)])
gen_expression <- data[, 8:ncol(data)]
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)
data$risk_scores <- risk_scores$risk_scores
res.cut <- surv_cutpoint(data, time = "overall_survival", event = "deceased",
variables = c("risk_scores"))
cutpoint <- res.cut[["cutpoint"]][["cutpoint"]]
data$risk_scores <- ifelse(data$risk_scores > cutpoint, "HIGH", "LOW")
fit_risk_scores <- survfit(Surv(overall_survival, deceased) ~ risk_scores, data = data)
ggsurvplot(fit_risk_scores, data, pval = TRUE, risk.table = TRUE)
coxa <- coxph(Surv(overall_survival, deceased) ~ risk_scores, data = data)
coxa
## Call:
## coxph(formula = Surv(overall_survival, deceased) ~ risk_scores,
## data = data)
##
## coef exp(coef) se(coef) z p
## risk_scoresLOW -0.8689 0.4194 0.1834 -4.739 2.15e-06
##
## Likelihood ratio test=21.72 on 1 df, p=3.152e-06
## n= 144, number of events= 143
## (253 observations deleted due to missingness)