【中药网络药理学】筛选细胞衰老和预后相关基因（附分类代码和画图代码）

1、衰老相关基因

从HAGR和msigdb数据获取细胞衰老相关基因，将两者取交集后构建基因蛋白互作网络

HAGR数据库

该库本身提供了下载链接，我在下载后对其进行了清洗

msigdb数据库

以"aging"作为关键词，Search Filters中collection设置为"all collections"，source species设置为"Homo sapiens",contributors设置为"all contributions"

在msigdb数据库中共得到56个衰老相关基因集

2、小细胞肝癌预后相关基因_TCGA数据库

library("survival")
library("survminer")

使用TCGA数据库通过TCGA-LIHC数据集(2023年1月,n=377)，下载clinical（临床信息），点击Metadata（样本信息），cart（基因文件）[/data/TCGA数据库]获取转录组测序数据与临床数据，并排除了临床数据不完整的患者队列

我使用了Harmony软件包进行批次效应调整，

接下来，基于拟合多因素Cox回归和lasso回归，构建预后效果评价分类模型

以"stranded_first"作为特定基因的表达水平评价量，临床信息中选择vital_status：生存状态（例如，存活、死亡）和days_to_death：诊断后至死亡的天数，作为评价的临床指标，最后整合得到如下表格：

病例ID	基因名	生存时间	结局
ID1	GeneA	时间1	结果1
ID2	GeneB	时间2	结果2
ID3	GeneC	时间3	结果3

描述性分析

首先使用生存分析中的Cox比例风险模型建立一个风险评分系统，然后利用风险评分中位数对其进行分组，将其分为“高风险组”和”低风险组“，并绘制生存分析图（针对筛选出的基因）

# 生存分析
load("files/exp_clinical")install.packages("survival") 
install.packages("survminer")
install.packages("glmnet")
install.packages("caret")library(survival)
library(survminer)
library(glmnet)
library(dplyr)
library(caret)# 绘制Kaplan-Meier曲线
exp_clinical <- as.data.frame(t(exp_clinical))
exp_clinical$time <- as.numeric(exp_clinical$time)
exp_clinical <- exp_clinical[!is.na(exp_clinical$time), ]# 观察小范围内数据框
View(exp_clinical[(1:50), (1:50)])
View(gene_data[(1:50), (1:50)])
View(normalized_data[(1:50), (1:50)])gene_data <- exp_clinical[(1:(nrow(exp_clinical)-5)), ]
clinical_data <- exp_clinical[((nrow(exp_clinical)-4) : nrow(exp_clinical)), ]# 选择至少在一个样本中表达量高于该基因90%分位数的基因
sapply(gene_data, class)gene_data[is.na(gene_data)] <- 0
gene_data[] <- lapply(gene_data, function(x) as.numeric(as.character(x)))gene_data_scale <- as.data.frame(scale(gene_data))normalize_minmax <- function(x) {return((x - min(x, na.rm = TRUE)) / (max(x, na.rm = TRUE) - min(x, na.rm = TRUE)))
}gene_data <- as.data.frame(lapply(gene_data_scale, normalize_minmax))names(gene_data) <- names(clinical_data)
data <- rbind(clinical_data, gene_data)
data <- as.data.frame(t(data))rownames(data)[6:nrow(data)] <- rownames(exp_clinical)[1:(nrow(exp_clinical)-5)]
save(data, file = "files/gene_clinical_data")# 做相关性分析
load("files/gene_clinical_data")View(data[(1:50), (1:50)])

rm(list = ls())load("files/exp_clinical")
load("files/agging_HCC_tcm_gene")library(ggplot2)
library(ggpubr)
library(survival)
library(survminer)
library(glmnet)
library(mice)
library(dplyr)
library(car)
library(randomForestSRC)
library(caret)## 预处理
rownames(exp_clinical) <- data_rownames
exp_clinical <- replace(exp_clinical, is.na(exp_clinical), 0)
temp <- as.data.frame(t(exp_clinical))temp <- data.frame(lapply(temp, as.numeric))
temp <- replace(temp, is.na(temp), 0)colnames(temp) <- data_rownames
colnames(temp)gene_data <- as.data.frame(temp[, unlist(aging_HCC_tcm_gene)])
temp_gene_data <- as.data.frame(lapply(gene_data, as.numeric))str(temp_gene_data)temp_gene_data <- scale(temp_gene_data)
gene_data <- as.data.frame(temp_gene_data)
# temp <- scale(gene_data)rownames(gene_data) <- t(colnames(exp_clinical))exp_clinical <- as.data.frame(t(exp_clinical))
clinical_data <- as.data.frame(cbind(exp_clinical$time, exp_clinical$flag))colnames(clinical_data) <- c("time", "flag")clinical_data <- as.data.frame(lapply(clinical_data, as.numeric))
clinical_data[is.na(clinical_data)] <- 0sum(clinical_data$flag)rownames(clinical_data) <- rownames(gene_data)rows_to_remove <- rownames(clinical_data[clinical_data$time == 0, ])gene_data <- gene_data[!rownames(gene_data) %in% rows_to_remove, ]
clinical_data <- clinical_data[!rownames(clinical_data) %in% rows_to_remove, ]
# clinical_data <- clinical_data[!rownames(clinical_data) %in% rows_to_remove, ]# rownames(gene_data) <- t(colnames(exp_clinical))# clinical_data <- as.matrix(temp[, c("time", "flag")])
# clinical_data <- as.data.frame(clinical_data)
# 
# clinical_data <- as.data.frame(lapply(clinical_data, as.numeric))
# sum(clinical_data$flag)
# 
data <- cbind(clinical_data, gene_data)
# # temp_data <- cbind(clinical_data, temp)
# 
# sum(is.na(clinical_data))
# sum(is.na(gene_data))# gene_data <- preProcess(gene_data, method = "range")save(data, gene_data, clinical_data, file = "files/surv_pre_data")## 风险回归
rm(list = ls())
load("files/surv_pre_data")str(clinical_data)
str(gene_data)custom_control <- coxph.control(eps = 1e-09, iter.max = 10000, toler.chol = 1e-10)surv_model <- coxph(Surv(time, flag) ~ ., data = data, control = custom_control)# temp_surv_model <- coxph(
#  Surv(time, flag) ~ ., 
#  data = temp_data, 
#  control = coxph.control(iter.max = 10000)
#)surv_model$coefficients# surv_model <- temp_surv_modelcoefficients <- as.data.frame(surv_model$coefficients)
coefficients[is.na(coefficients)] <- 0
coefficientsgene_data <- data.frame(gene_data)
gene_coef_data <- rbind(gene_data, t(coefficients))
rownames(gene_coef_data)[nrow(gene_coef_data)] <- c("coef")surv_results <- sweep(gene_data, 2, t(coefficients), `*`)
surv_row_sums <- rowSums(surv_results)surv_row_sums[is.na(surv_row_sums)] <- 0
surv_row_sumssurv_median_value <- median(surv_row_sums)
surv_median_value# 分类并添加分组
risk_group <- ifelse(surv_row_sums > surv_median_value,"High Risk","Low Risk")
surv_results$risk_group <- risk_group
final_surv_results <- surv_results[, !apply(surv_results, 2,function(x) all(x == 0))]save(data, risk_group,gene_data, clinical_data, coefficients, final_surv_results, file = "files/surv_final_data")# 绘制生存曲线
rm(list = ls())load("files/surv_final_data")colnames(final_surv_results)temp_data <- data
temp_data$time <- ceiling(temp_data$time/30)
# temp_data <- lapply(temp_data, as.numeric)# 剔除大于5年的
rows_to_remove <- rownames(temp_data[temp_data$time > 60, ])
temp_data <- temp_data[!rownames(temp_data) %in% rows_to_remove, ]
clinical_data <- clinical_data[!rownames(clinical_data) %in% rows_to_remove, ]
risk_group <- as.data.frame(risk_group)
risk_group <- risk_group[!rownames(risk_group) %in% rows_to_remove, ]temp_surv_obj <- Surv(temp_data$time, temp_data$flag)
temp_data <- cbind(temp_data, risk_group)
temp_surv_fit <- survfit(temp_surv_obj ~ risk_group, data = temp_data)
temp_max <- max(temp_data$time)
ggsurvplot(temp_surv_fit, data = temp_data,risk.table = TRUE, # 显示风险表pval = TRUE, # 显示P值conf.int = TRUE, # 显示置信区间xlim = c(0, temp_max), # 可选：设置X轴的范围xlab = "Times(months)", # 设置X轴标签ylab = "Survival probability", # 设置Y轴标签
)coefficients$gene_name <- rownames(coefficients)
colnames(coefficients)[1] <- c("coef")
writexl::write_xlsx(coefficients, "data/washed_data/coef.xlsx")surv_obj <- Surv(data$time, data$flag)
surv_fit <- survfit(surv_obj ~ risk_group, data = data)
max <- max(data$time)
ggsurvplot(surv_fit, data = final_surv_results,risk.table = TRUE, # 显示风险表pval = TRUE, # 显示P值conf.int = TRUE, # 显示置信区间xlim = c(0, max), # 可选：设置X轴的范围xlab = "Times(days)", # 设置X轴标签ylab = "Survival probability", # 设置Y轴标签title = "Kaplan-Meier 生存曲线" # 设置图标题
)