三阴性乳腺癌预后预测：基于TCGA数据库的生存分析

library(cBioPortalData)
library(survival)
library(survminer)
library(TCGAbiolinks)
library(dplyr)
library(ggplot2)


cbio <- cBioPortal(
  hostname = "www.cbioportal.org",
  protocol = "https",
  api. = "/api"  # 更新 API 路径
)

clinical <- cBioPortalData::clinicalData(cbio, studyId = "brca_tcga")

# 筛选 TNBC 样本
tnbc_samples <- clinical[
  clinical$ER_STATUS_BY_IHC == "Negative" &
    clinical$PR_STATUS_BY_IHC == "Negative" &
    clinical$IHC_HER2 == "Negative",
]

# 查看 TNBC 样本数量
# 去除 patientId 为 NA 的行
clin_data <- tnbc_samples[!is.na(tnbc_samples$patientId), ]

save(clin_data,file = "TNBC_clin.Rdata")

# 创建查询对象
query <- GDCquery(
  project = "TCGA-BRCA",
  data.category = "Transcriptome Profiling",
  data.type = "Gene Expression Quantification",
  workflow.type = "STAR - Counts")

# 启用多线程下载
GDCdownload(query,  
            files.per.chunk = 10,
            directory = "GDCdata")

data <- GDCprepare(query)
expre_matrix <- assay(data)

save(expre_matrix,file = "BRCA_expre_matrix.Rdata")

# 过滤出 TNBC 样本的表达数据
# 提取简化的 patientId
colnames(expre_matrix) <- substr(colnames(expre_matrix), 1, 12)

# 查看转换后的列名
colnames(expre_matrix)[1:4]

# 筛选出在 clin_data 中存在的样本
matched_samples <- colnames(expre_matrix) %in% clin_data$patientId
exp_matrix <- expre_matrix[, matched_samples]

# 查看筛选后的样本数量
print(paste("匹配的样本数量:", ncol(exp_matrix)))

exp_mat <- exp_matrix[,!duplicated(colnames(exp_matrix))] ##去除重复样本

# 加载 biomaRt 包
library(biomaRt)

# 连接到 Ensembl 数据库
ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")

# 获取基因 ID 列表
rownames(exp_mat) <- sub("\\.[0-9]+$", "", rownames(exp_mat))
gene_ids <- rownames(exp_mat)

# 定义需要获取的属性
attributes <- c("ensembl_gene_id", "hgnc_symbol")

# 使用 getBM() 进行转换
gene_info <- getBM(
  attributes = attributes,
  filters = "ensembl_gene_id",
  values = gene_ids,
  mart = ensembl
)

# 将 gene_info 与 exp_mat 合并
exp_mat <- as.data.frame(exp_mat)
exp_mat$ensembl_gene_id <- rownames(exp_mat)

# 使用 left_join() 合并数据
exp_mat <- left_join(exp_mat, gene_info, by = "ensembl_gene_id")

exp_mat=exp_mat[!duplicated(exp_mat$hgnc_symbol),]#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果s

# 去除 hgnc_symbol 列中的 NA 值
exp_mat <- exp_mat[complete.cases(exp_mat$hgnc_symbol), ]


rownames(exp_mat) <- exp_mat$hgnc_symbol

# 去除多余的列
exp_mat <- exp_mat[, !colnames(exp_mat) %in% c("ensembl_gene_id", "hgnc_symbol")]


# 查看转换后的表达矩阵
exp_mat[1:4, 1:4]
# 对表达量矩阵取 log2
log_exp_mat <- log2(exp_mat + 1)  # 加 1 是为了避免对 0 取对数

# 查看转换后的数据
range(log_exp_mat)

save(log_exp_mat,file = "gene_exp_mat.Rdata")

# 假设你的基因列表是一个字符向量
gene_list <- c("ANXA8","IL1B","HCLS1","STAT1","S100A16","C1QBP","TLR2","MRM3","CLDN4","DSP","EPPK1","JUP","G6PC3","SAMSN1","MTSS1","IL13RA1","IRGM","PDPN","SERPINB5","AMB3","LAMC2","SMN1")  # 替换为你的基因列表

expression_data <- log_exp_mat[rownames(log_exp_mat) %in% gene_list, ]

# 将表达数据转置，使样本在行中
expression_data <- t(expression_data)
expression_data <- as.data.frame(expression_data)
expression_data$patientId <- rownames(expression_data)

# 合并临床数据和表达数据
merged_data <- merge(clin_data[,c("patientId","OS_STATUS","OS_MONTHS")], expression_data, by = "patientId")

merged_data <- merge(clin_data[,c("patientId","DFS_STATUS","DFS_MONTHS")], expression_data, by = "patientId")

# 去除 NA 值
merged_data <- merged_data[complete.cases(merged_data$DFS_STATUS), ]

# # 从临床数据中提取生存时间和生存状态,创建生存对象
merged_data$DFS_STATUS <- ifelse(merged_data$DFS_STATUS == "1:Recurred/Progressed", 1, 0)
merged_data$DFS_MONTHS <- as.numeric(merged_data$DFS_MONTHS)
surv_obj <- Surv(time = merged_data$DFS_MONTHS, event = merged_data$DFS_STATUS)

![image-20250312111556963](TNBC%E7%94%9F%E5%AD%98%E5%88%86%E6%9E%90%EF%BC%9A%E6%9D%A5%E8%87%AATCGA%E7%9A%84%E9%A2%84%E5%90%8E%E5%88%86%E6%9E%90.assets/image-20250312111556963.png)# 加载必要的包
library(survival)
library(survminer)
library(ggplot2)

# 定义输出文件夹
output_dir <- "TNBC_output"
if (!dir.exists(output_dir)) {
  dir.create(output_dir)
}

# 初始化 significant_genes
significant_genes <- list()

for (gene in gene_list) {
# 检查基因是否存在于 merged_data 中
if (!gene %in% colnames(merged_data)) {
    message(paste("跳过基因", gene, "，因为它不存在于 merged_data 中。"))
    next
  }

# 检查基因表达数据是否为数值型
if (!is.numeric(merged_data[[gene]])) {
    message(paste("跳过基因", gene, "，因为它的表达数据不是数值型。"))
    next
  }

# 检查基因表达数据中是否有缺失值
if (any(is.na(merged_data[[gene]]))) {
    message(paste("基因", gene, "的表达数据中存在缺失值，已忽略缺失值。"))
  }

# 计算中位数并分组
  median_value <- median(merged_data[[gene]], na.rm = TRUE)
  merged_data$group <- ifelse(merged_data[[gene]] > median_value, "High", "Low")

# 检查分组是否成功
if (length(unique(merged_data$group)) < 2) {
    message(paste("基因", gene, "的分组失败，可能因为表达量全部相同。"))
    next
  }

# 拟合生存曲线
  fit <- survfit(surv_obj ~ group, data = merged_data)

# 计算 p 值（使用 log-rank 检验）
  logrank_test <- survdiff(surv_obj ~ group, data = merged_data)
  p_value <- 1 - pchisq(logrank_test$chisq, df = 1)  # 提取 p 值

# 将基因和 p 值保存到列表中
  significant_genes[[gene]] <- p_value

# 绘制生存曲线
  surv_plot <- ggsurvplot(
    fit, 
    data = merged_data, 
    pval = TRUE, 
    risk.table = TRUE, 
    title = paste("Survival Curve for", gene),
    legend.labs = c("High Expression", "Low Expression")
  )

# 定义保存路径
  output_file <- file.path(output_dir, paste0("survival_curve_", gene, ".png"))

# 保存图像
  ggsave(output_file, surv_plot$plot, width = 8, height = 6, dpi = 300)

# 打印保存路径
  message(paste("生存曲线已保存到：", output_file))
}