玩转 TCGA 数据库 - 生存分析（三）

生信菜鸟团

发布于 2025-05-21 15:22:34

1.8K0

生存分析：事件的结果和出现这一结果所经历的事件结合起来分析的一种方法。通常情况下，我们设定一个入组时间区间，在这个区间内搜寻患者的第一次患病时间称为起始时间，当过了这个入组时间区间，我们就不再收集患者了，当患者发生终点事件（比如死亡）时，我们记录此事件为终点时间。同时实验会设置一个实验截止时间，实验终止后患者仍未发生终点事件我们将实验截止时间记录为这个患者的终点时间，但终点事件记录为删失。

起始时间：患者的入组时间。

终点时间：患者发生结局事件的时间 / 研究终止时间。

生存事件：起始时间和终点时间的差值。

生存结局：是否发生结局事件。删失/死亡

删失 (censoring)：删失也细分为多种类型。

右删失：研究截止时，感兴趣终点事件未出现。
左删失：患者在入组时间起始之前就被确认为患病。
区间删失：无法确定准确的患病时间，但患病时间是在入组时间区间内。

生存率，也称累积生存概率，是在给定的事件点，所研究的个体继续存活的概率。

生存分析方法：

描述分析：根据样本生存资料估计总体生存率或者其他有关指标。（K-M 生存曲线）。首先要算出生存率及其标准误，在此基础上得到生存曲线，进而得到中位生存时间（线性内差法/生存曲线法）。
比较分析：对不同组生存率进行比较分析（log-Rank 检验和 Breslow 检验）。log-rank/mantel-cox 检验对所有时间点权重相同，对远期差异敏感，如果生存曲线的分离发生在研究的后期，可以选择这种方法。breslow/wilcoxon 检验是根据每个时间点进行加权，按时间点死亡例数加权，对近期差异敏感。
影响因素分析：通过生存分析模型来探讨影响生存时间的因素（COX 比例风险模型，挖掘影响生存结局的潜在风险，并得到矫正后的风险比）。模型中以生存时间和生存结局为应变量，感兴趣的因素为自变量，分为保护因素和有害因素。

TCGA 数据库有多种类型的癌症数据，我们可以直接下载临床数据和基因表达谱数据，做生存分析。

加载包

library(tidyverse)
library(TCGAbiolinks)
library(SummarizedExperiment)

下载 TCGA 数据

TCGAbiolinks 下载 project 数据

project <- "TCGA-PAAD" #⭐️ 关注的 project
TCGAbiolinks:::getProjectSummary(project)
query <- GDCquery(project = project,
                  data.category = "Transcriptome Profiling",
                  data.type = "Gene Expression Quantification",
                  workflow.type = "STAR - Counts"
                 )
GDCdownload(query)
data <- GDCprepare(query = query)

整理临床信息和基因表达数据

clinical <- colData(data)

data@assays@data@listData %>% names()
counts <- assay(data, "unstranded")
fpkm <- assay(data, "tpm_unstrand")
tpm <- assay(data, "tpm_unstrand")

使用 TCGAbiolinks 包下载 & 准备的数据是属于 SummarizedExperimen 结构的，在这个下载对象中包含了多个表达矩阵。

> data@assays@data@listData %>% names()
[1] "unstranded"       "stranded_first"   "stranded_second"  "tpm_unstrand"     "fpkm_unstrand"   
[6] "fpkm_uq_unstrand"

unstranded：表示 RNA-Seq 数据在处理时没有考虑链的方向性。
stranded_first：表示 RNA-Seq 数据是链特异性的，并且第一条链（通常是正义链）被考虑。
stranded_second：表示 RNA-Seq 数据是链特异性的，并且第二条链（通常是反义链）被考虑。
tpm_unstrand：表示未区分链的 TPM（Transcripts Per Million）值。这是另一种标准化基因表达的方法，考虑了转录本长度和测序深度。
fpkm_unstrand：表示未区分链的 FPKM（Fragments Per Kilobase of transcript per Million mapped reads）值，用于标准化基因表达。
fpkm_uq_unstrand：表示未区分链的 Upper Quartile-normalized FPKM 值，这是一种通过上四分位数标准化的 FPKM，用于减少样本间的技术变异。

表达矩阵行名 ID 转换

将 ENSEMBL ID 转换为 SYMBOL ID

library(tinyarray)
exp <- tpm # 这里我选择的是tpm作为表达谱方式
exp <- trans_exp_new(exp)
exp[1:4,1:4]

基因过滤

nrow(exp)
# 仅保留在一半以上样本里表达的基因
exp <- exp[apply(exp, 1, function(x) sum(x > 0) > 0.5*ncol(exp)), ]
nrow(exp)

# 仅去除在所有样本里表达量都为零的基因
# exp1 <- exp[rowSums(exp)>0,]
# nrow(exp1)

生存分析数据整理

分组信息获取

Group <- make_tcga_group(exp)
table(Group)

去除 norm 样本

这里可以根据生存分析的定义思考一下为什么要去除非癌症样本呀。

table(Group)
exprSet <- exp[, Group=='tumor']
ncol(exp)
ncol(exprSet)

基因过滤

因为我们去掉了一部分样本，这里就再次过滤一下吧。

nrow(exprSet)
exprSet <- exprSet[apply(exprSet,1, function(x){sum(x>0)>0.5*ncol(exprSet)}),]
nrow(exprSet)

使用logCPM或logTPM数据

因为不同样本间进行基因表达量比较使用统一的归一化手段才有可比性，我使用的是TPM。

# logTPM
exprSet <- log2(exprSet+1)
exprSet[1:4,1:4]

整理生存分析和临床信息

根据生存分析的定义可知，我们需要生存时间和生存结局这两个临床信息，如果需要根据基因表达量对样本分组，我们还需要基因表达谱数据。

TCGA 临床信息整理

meta <- as.data.frame(clinical)
nrow(meta)
length(unique(meta$sample))
meta <- distinct(meta,sample,.keep_all = T)

样本处理

这一步在不同的数据库中整理的临床信息字段是不一样的，所以自己要拿出临床信息字段去看一下什么意思再去定义。可以去官网探索一下临床信息意义。https://gdc.cancer.gov/about-data/gdc-data-processing/clinical-data-standardization

meta <- meta %>%
  mutate(overall_survival = ifelse(vital_status == "Dead", days_to_death, days_to_last_follow_up)) %>%
  mutate(event = vital_status)

去掉生存信息不全的样本或者自己不想要的样本

# 去掉生存信息不全或者生存时间小于30天的样本
k1 <- meta$overall_survival >= 30;table(k1)
k2 <- !(is.na(meta$overall_survival)|is.na(meta$overall_survival));table(k2)
meta <- meta[k1 & k2,]

简化、规范化变量

结局事件：生存分析的输入数据里，要求结局事件必须用0和1表示，1表示阳性结局。

table(meta$event)
meta$event <- ifelse(meta$event == "Dead", 1, 0)

生存时间认清生存时间的单位（通常是月，也可以用天和年）。生存分析中一般最开始是使用天进行统计的，要转化为月或者年时，精细一点来说，一个月是 30.4375 天，一年是 365.25 天。但也可粗糙计算，影响不会很大。

range(meta$overall_survival)
meta$overall_survival <- meta$overall_survival/30.4375
range(meta$overall_survival)

表达矩阵和临床信息匹配

head(rownames(meta))
head(colnames(exprSet))

s <- intersect(rownames(meta),colnames(exprSet));length(s)
exprSet <- exprSet[,s]
meta <- meta[s,]

dim(exprSet)
dim(meta)
identical(rownames(meta),colnames(exprSet))

生存分析

我们常说的生存分析其实是生存分析比较分析

library(survival)
library(survminer)

分类变量生存分析

# 根据性别分组
sfit <- survfit(Surv(overall_survival, event)~gender, data = meta)
ggsurvplot(sfit,pval=TRUE)
ggsurvplot(sfit,
           palette = "jco",
           risk.table =TRUE,
           pval =TRUE,
           conf.int =TRUE)

连续变量生存分析

g <- "TP53"
meta$gene <- ifelse(exprSet[g,]> median(exprSet[g,]),'high','low')
table(meta$gene)
sfit <- survfit(Surv(overall_survival, event)~gene, data = meta)
ggsurvplot(sfit, pval =TRUE, data = meta, risk.table = TRUE)

批量生存分析

如果想知道在我们数据中哪些基因高低表达是显著影响生存率的话，可以批量进行生存分析，然后看有没有关注的。

geneKM <- function(gene){
  meta$group <- ifelse(gene>median(gene),'high','low')  
  data.survdiff <- survdiff(Surv(overall_survival, event)~group,data=meta)
  p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
  return(p.val)
}

log_rank_p <- apply(exprSet, 1, geneKM) %>% sort()
head(log_rank_p)
table(log_rank_p<0.01) 
table(log_rank_p<0.05)
log_rank_p %>% as.data.frame() %>% rownames_to_column() %>% view()

lr <- names(log_rank_p)[log_rank_p<0.01];length(lr)

批量单因素 COX

也可以做一下生存分析影响因素分析，看哪些基因属于风险因素或保护因素。

genecox <- function(gene){
  meta$gene <- gene
  #可直接使用连续型变量
  m <- coxph(Surv(overall_survival, event) ~ gene, data =  meta)
  #也可使用二分类变量
  #meta$group=ifelse(gene>median(gene),'high','low') 
  #meta$group = factor(meta$group,levels = c("low","high"))
  #m=coxph(Surv(overall_survival, event) ~ group, data =  meta)

  beta <- coef(m) # 提取Cox模型的系数。
  se <- sqrt(diag(vcov(m))) # 计算系数的标准误差
  HR <- exp(beta) # 计算风险比（Hazard Ratio），通过对系数取指数得到。
  HRse <- HR * se # 计算HR的标准误差

  #summary(m)
  tmp <- round(cbind(coef = beta, 
                     se = se, # 系数的标准误差
                     z = beta/se, # z统计量，系数除以标准误差
                     p = 1 - pchisq((beta/se)^2, 1), # p值，基于z统计量的卡方分布。
                     HR = HR, # 风险比
                     HRse = HRse, # 风险比的标准误差
                     HRz = (HR - 1) / HRse, # HR的z值
                     HRp = 1 - pchisq(((HR - 1)/HRse)^2, 1), # HR的p值
                     HRCILL = exp(beta - qnorm(.975, 0, 1) * se), # HR的95%置信区间下限
                     HRCIUL = exp(beta + qnorm(.975, 0, 1) * se)), 3) #HR的95%置信区间上限

  return(tmp['gene',]) 
  #return(tmp['grouphigh',])#二分类变量
}


cox_results <-apply(exprSet, 1, genecox)
cox_results <- as.data.frame(t(cox_results))

table(cox_results$p<0.01)
table(cox_results$p<0.05)

cox <- rownames(cox_results)[cox_results$p<0.01];length(cox)

length(intersect(lr,cox))

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2025-05-19，如有侵权请联系 cloudcommunity@tencent.com 删除

数据库