TCGA数据库| 如何根据单个基因的表达对样本分组绘制KM生存曲线？

生信技能树

发布于 2025-06-30 08:53:45

1.1K0

文章被收录于专栏：生信技能树生信技能树

关于TCGA数据库，我们介绍了：

样本barcode的详细组成：TCGA样本barcode详细介绍
TCGA数据库| 如何将表达矩阵与样本临床数据进行合并？

今天来看看单个基因的生存分析。

此外这期内容也进行了了直播，我们开了一个新的微信交流群，准备做100期数据实战直播，带领大家加强代码训练，数据包括：GEO芯片+bulk转录组+单细胞转录组+空间转录组的分析。

如果你希望进群与大家一起交流，可以加我微信：Biotree123，发18.8的拉群费用拉你进群~

前面已经开展了5期啦，直播回放集合：https://www.bilibili.com/video/BV1eEMJznEog

本次讲解案例为高分文献肺癌耐药关键gene筛选经典思路：

01期数据集GSE231938-bulkRNA-Seq，包括上下游分析代码讲解（01期-03期代码文件：02-Human-8-NSCLC-Trans，链接: https://pan.baidu.com/s/19SezImI4PzgmlIfP21MX4Q?pwd=ukwb）
02期数据集GSE7670-Affymetrix芯片，包括超详细的芯片数据分析代码讲解
03期数据集GSE32863-Ilumina芯片数据处理
04期数据集TCGA数据库：数据下载，预处理，与临床信息整合，差异表达分析~（02-Human-8-NSCLC-Trans-TCGA.zip，链接: https://pan.baidu.com/s/10pOgfydqr_KmDqiJYDCnvw?pwd=y2v8）
05期数据集TCGA数据库：临床信息预处理，目标基因表达高低分组KM生存曲线绘制，目标基因在不同病理stage中表达比较~~（02-Human-8-NSCLC-Trans-TCGA-2.zip 链接: https://pan.baidu.com/s/1BCFchszjP3g2NwStHUpFWQ?pwd=swab）

读取表达矩阵

读取整理好的表达矩阵进来：tcga_mrna_fpkm_symbol.rds

这里使用R包 TCGAbiolinks 去TCGA官网下载数据，参考：TCGA数据库| 如何将表达矩阵与样本临床数据进行合并？

这个文件也可以在这里找到（https://pan.baidu.com/s/1BCFchszjP3g2NwStHUpFWQ?pwd=swab）

rm(list=ls())
# 加载包
library(survival)
library(survminer)

## 加载表达矩阵
exp <- readRDS("TCGA/tcga_mrna_fpkm_symbol.rds")
exp <- as.data.frame(exp)
rownames(exp) <- exp[,1]
exp <- exp[,-1]
exp[1:5,1:5]
exp["IGF2BP3",1:5]

读取临床信息

读取整理好的临床信息，并进行处理：

## 临床信息
clinical <- readRDS(file = "TCGA/TCGA-LUAD.clinical_patient.rds")
colnames(clinical)


## 提取对应的生存信息：生存结局，随访时间，死亡时间
clinical <- clinical[, c("bcr_patient_barcode", "vital_status", "days_to_death", "days_to_last_followup") ]
head(clinical)

## 去掉 没有生存结局的样本
clinical <- clinical[!is.na(clinical$vital_status), ]
clinical <- clinical[clinical$vital_status!="", ]
clinical$vital_status <- ifelse(clinical$vital_status=="Alive", 0, 1)
table(clinical$vital_status)

## 新增一列总生存时间OS
clinical$OS <- if_else( is.na(clinical$days_to_death), clinical$days_to_last_followup, clinical$days_to_death)
## 去掉没有生存时间的样本
clinical <- clinical[ !is.na(clinical$OS), ]
## 去掉小于30天的样本
clinical <- clinical[clinical$OS>=30, ]
# 单位变成月
clinical$OS <- clinical$OS / 30
rownames(clinical) <- clinical$bcr_patient_barcode
str(clinical)

这里生成了一个总生存时间OS，并去掉了小于30天的病人：

合并临床信息和基因表达

合并这部分的信息为绘图做准备：

## 整合表达矩阵与 生存信息
table(str_sub(colnames(exp),14,16))
exp_tumor <- as.data.frame(exp[, str_sub(colnames(exp),14,16)=="01A"])
head(colnames(exp_tumor))
head(clinical)
# 修改表达矩阵的列名为patient id
colnames(exp_tumor) <- str_sub(colnames(exp_tumor), 1,12)
head(colnames(exp_tumor))

# 获取 肿瘤样本与有临床信息样本的交集
comid <- intersect(colnames(exp_tumor), clinical$bcr_patient_barcode)
length(comid)

# 合并
data <- data.frame(patientid=comid,gene=as.numeric(t(exp_tumor["IGF2BP3", comid])), clinical[comid, ])
head(data)

生存曲线绘制

使用基因表达将样本分为高低两组，并比较这两组的总生存时间差异：

####### 生存曲线
# 将样本 按照基因表达分组
data$group <-  if_else(data$gene > median(data$gene), "high", "low")
table(data$group)

# 绘制基础生存曲线
fit <- survfit(Surv(data$OS, data$vital_status) ~ group, data = data)
fit
p <- ggsurvplot(fit, data = data, 
                pval.method = T,
                pval = TRUE,    # 添加log rank检验的p值
                conf.int= TRUE, # 添加置信区间
                palett = c("#fe7272","#4071a1"), # 曲线的颜色
                xlab = "Time (Months)",   # 设置X轴标题
                ylab = "Overall survival"
                )
p$plot

# 保存为矢量pdf文件
ggsave(filename = "Fig_i.pdf",width = 5,height = 4.5,plot = p$plot,bg="white")

结果如下：