Kaplan |
您所在的位置：网站首页 › 生存分析中删失是什么意思呀 › Kaplan |
|
1. Kaplan-Meier生存分析
Kaplan Meier是一种单因素生存分析，它可用于研究1个因素对于生存时间的影响。 参考：生存分析之Kaplan-Meier曲线都告诉我们什么 思路一般来说，我们做生存分析，会有（P0.05）两种结果。KM plot在生物医学中很常见，主要用来做预后分析，比如可以根据表达量把病人分成两组，然后比较哪组病人预后好，进而可以得出基因表达量高低与病人预后好坏相关性的结论。(但凡是能把病人分成两个组的，都可以拿来做生存分析，根据基因表达量，人种，性别，stage等等来分，都是可以的。) 画KM plot时，有时候会比较纠结怎样对病人进行分组，如何来设置分组的cutoff。一般来说常见的几种设置cutoff值得思路如下： 1：大多数情况下，根据表达量从低到高对样本进行排序，取前50%为低表达，后50%为高表达，然后画KM plot。 2：还有一些文章也会将样本表达量均分为三组或者四组。 3：一些文章也会选一些其它的cutoff，比如前1/3和后2/3，前25%和后25%（中间50%的数据去掉）。 参考：https://zhuanlan.zhihu.com/p/86307194#:~:text=KM%20plot 2. 结果解读左图：KM-plot的横坐标是生存时间（可以是OS，也可以是PFS等），纵坐标是生存率。起点是随访开始的时间，曲线下降代表患者死亡，曲线上的+号表示删失（仍然存活的病人的最后随访时间）。p值是log-rank test的p值，衡量的绘制的曲线之间生存率有没有显著差别。 右上：左图的美化版，曲线的阴影表示置信区间。 右中：每个时间节点有多少人在随访中 右下：删失表格，右上图中每出现一个+，也就是删失，右下对应的位置就会有一个条图。纵轴是在该时间点删失的数目。 3. 数据整理⚠️生存分析需要的输入数据有两个：表达矩阵数据和临床信息数据 生存分析只需要tumor数据，不要normal，将其去掉，新表达矩阵数据命名为exprSet； clinical信息需要进一步整理，成为生存分析需要的格式，新临床信息数据命名为meta。 由于不同癌症的临床信息表格组织形式不同，这里的代码需要根据实际情况修改。 rm(list=ls()) options(stringsAsFactors = F) load("TCGA-CHOL_gdc.Rdata") load("TCGA-CHOL_symbol_exp.Rdata") library(stringr) 3.1 整理表达矩阵与临床信息 整理表达矩阵 exprSet=exp[,Group=='tumor'] #只留下肿瘤样本 #前面的过滤标准比较宽松，这里可以再次过滤 k = apply(exprSet,1, function(x){sum(x==0) TCGA-W5-AA2Q TCGA-W5-AA2Q Alive #> TCGA-W5-AA2O TCGA-W5-AA2O Dead 640 #> TCGA-W5-AA39 TCGA-W5-AA39 Dead 170 #> TCGA-W6-AA0S TCGA-W6-AA0S Alive #> days_to_last_followup #> TCGA-W5-AA2Q 50 #> TCGA-W5-AA2O #> TCGA-W5-AA39 #> TCGA-W6-AA0S 352 #简化meta的列名 colnames(meta)=c('ID','event','death','last_followup','race','age','gender','stage') #空着的值改为NA meta[meta==""]=NA 3.2 实现表达矩阵与临床信息的匹配有的病人会有两个或两个以上的肿瘤样本，就有重复。两种可行的办法： （1）以病人为中心，对表达矩阵的列按照病人ID去重复，每个病人只保留一个样本。（本文档） （2）以样本为中心，把meta里的病人ID替换成样本ID，这样同一个病人的两个样本就会有两行完全一致的临床信息。（zz.R） # 以病人为中心，表达矩阵按病人ID去重复 k = !duplicated(str_sub(colnames(exprSet),1,12));table(k) # k # TRUE # 36 exprSet = exprSet[,k] #调整meta的ID顺序与exprSet列名一致 meta=meta[match(str_sub(colnames(exprSet),1,12),meta$ID),] identical(meta$ID,str_sub(colnames(exprSet),1,12)) # [1] TRUE 3.3 整理生存分析的输入数据 #1.1由随访时间和死亡时间计算生存时间(月) table(meta$event) # # Alive Dead # 20 16 meta$time = ifelse(meta$event=="Alive", meta$last_followup, meta$death) meta$time = as.numeric(meta$time)/30 #1.2 根据生死定义event，活着是0，死的是1 meta$event=ifelse(meta$event=='Alive', 0, 1) table(meta$event) # # 0 1 # 20 16 #1.3 年龄和年龄分组 meta$age=ceiling(abs(as.numeric(meta$age))/365) meta$age_group=ifelse(meta$age>median(meta$age,na.rm = T),'older','younger') table(meta$age_group) # # older younger # 18 18 #1.4 stage library(stringr) head(meta$stage) # [1] "6thStage IVT3N0M1" "7thStage IIIT3N0M0" "7thStage IIT2bNXMX" # [4] "7thStage IT1N0M0" "7thStage IIT2N0M0" "7thStage IT1N0M0" a = str_extract_all(meta$stage,"I|V");head(a) #提取分期信息 b = sapply(a,paste,collapse = "");head(b) #> [1] "IV" "III" "II" "I" "II" "I" meta$stage = b # 去掉生存信息不全或者生存时间小于0.1(月)的病人，样本纳排标准不唯一，且差别很大 k1 = meta$time>=0.1;table(k1) # k1 # FALSE TRUE # 1 35 k2 = !(is.na(meta$time)|is.na(meta$event));table(k2) # k2 # TRUE # 36 exprSet = exprSet[,k1&k2] meta = meta[k1&k2,] save(meta,exprSet,proj,file = paste0(proj,"_sur_model.Rdata")) 整理好的表达矩阵 整理好的临床信息矩阵 4. 绘制生存曲线 # 格式： sfit 25430 91 table(log_rank_p #> FALSE TRUE #> 24859 662 # 挑一个p值小的基因来做KM_plot g = names(log_rank_p)[1] meta$gene = ifelse(exprSet[g,]> median(exprSet[g,]),'high','low') sfit1=survfit(Surv(time, event)~gene, data=meta) ggsurvplot(sfit1,pval =TRUE, data = meta, risk.table = TRUE) 5.2 cox批量生存分析Cox 回归本质上是一种回归模型 ，它没有直接使用生存时间，而是使用了风险比（ hazard ratio ）作为因变量，该模型不用于估计生存率，而是用于因素分析也就是 找到某一个危险因素对结局事件发生的贡献度。 Cox 回归的重要统计指标： 风险比（ hazard ratio） • 当 HR>1 时，说明研究对象是一个危险因素。 • 当 HR |
| CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |