这是关于KM生存分析的notebook，记录了关于KM的基本方法，使用案例，基本结果，误差调整方案，统计检验方法等基本内容。希望彻底把这一块搞清楚。
还包括了一些基本的绘图，生存曲线，统计量，特别是样本分组表这种描述性的结果。
1 2 3 4 5 library(tidyverse) library(survival) library(survivalROC) library(survminer) library(tidyfst)
1 2 3 expset <- read.table(file = "./data/TCGA-COAD-fpkm.tsv",check.names = F,sep = "\t",header = T) meta <- read.csv(file = "./data/TCGA-COAD-meta.csv",check.names = F,header = T) surv <- read.table(file = "./data/TCGA-COAD.survival.tsv",check.names = F,sep = "\t",header = T)
1 2 3 row.names(expset) <- expset$Ensembl_ID expset <- expset[,-1] names(meta)[1] <- "sample"
1 2 3 4 #查看数据形状 dim(expset) dim(meta) dim(surv)
1 2 3 4 boxplot(expset[,1:5]) #查看缺失情况 Amelia::missmap(meta) Amelia::missmap(surv)
表达矩阵有6W+个基因，512个病人，meta中有120项meta，571个病人，而生存数据则有539个病人，4个特征，需要取最小交集后续分析。
此外，表达矩阵是带版本号的ensembl基因ID，理论上也可以最后再进行注释，但真正分析中肯定会在一开始进行注释，为了后续其他的分析，在此使用genecode.v22.probemap文件进行注释。
数据完整性上，表达矩阵和生存数据没有缺失，但meta信息中存在大量的缺失数据，需要进行缺失处理。
因此，数据处理的流程应当是：
1.处理meta中的缺失;
2.获得最小sample的交集作为全部的样本集;
3.表达矩阵的注释;
4.数据检查，固定外键为sample;
5.干净数据的存储持久化;
1 2 3 4 5 6 7 8 9 #处理meta中的缺失，阈值选择为5%(571*5%=30) meta[meta==""] <- NA coad.meta <- meta %>% tidyfst::delete_na_cols(n = 30) %>% tidyfst::delete_na_rows(n = 1) #获取最小sampleid的集合 commonsample = intersect(surv$sample,intersect(colnames(expset),coad.meta$sample))
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 #借助probemap进行基因注释，同时使用同版本的gtf文件注释编码基因 probemap <- read.table(file = "./data/gencodev22probeMap",sep = "\t",header = T,check.names = F) expset$id <- row.names(expset) coad.exp <- probemap %>% select(.,id,gene) %>% left_join(.,expset,by = "id") %>% subset(.,select = -id) gtfv22 <- rtracklayer::import("./data/gencode.v22.annotation.gtf") %>% as.data.frame(.) %>% select(.,c("gene_id","gene_type", "gene_name")) %>% subset(., gene_type == "protein_coding") %>% distinct(.,gene_name,.keep_all = T) coad.exp <- coad.exp %>% filter(.,gene %in% all_of(gtfv22$gene_name)) %>% distinct(.,gene,.keep_all = T) row.names(coad.exp) <- coad.exp$gene
1 2 3 4 5 #固定外键为sample，完成三者meta的全匹配 coad.exp <- coad.exp %>% select(.,all_of(commonsample)) coad.meta <- coad.meta %>% filter(.,sample %in% all_of(commonsample)) coad.surv <- surv %>% filter(.,sample %in% all_of(commonsample)) coad.surv <- coad.surv[,-3]
1 2 3 4 #干净数据的持久化 write.csv(coad.exp,file = "./data/coad.exp.csv") write.csv(coad.meta,file = "./data/coad.meta.csv") write.csv(coad.surv,file = "./data/coad.surv.csv")
1 2 3 #清理变量，释放内存 rm(expset,meta,surv,gtfv22,probemap) gc()
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 coad.exp.tmp <- as.data.frame(t(coad.exp)) coad.surv.tmp <- coad.surv %>% mutate(.,CXCL5 = coad.exp.tmp$CXCL5) survgrouptmp <- ifelse(coad.surv.tmp$CXCL5>median(coad.surv.tmp$CXCL5),'CXCL5High','CXCL5low') survivaltmp <-survfit(Surv(coad.surv.tmp$OS.time, coad.surv.tmp$OS)~survgrouptmp, data=coad.surv.tmp,) ggsurvplot(survivaltmp, conf.int=T, pval=T, pval.method = T, risk.table = F, ggtheme = theme_bw(), surv.median.line = "hv",) ggsave(filename = "survtmp.pdf",height = 8,width = 8)
1 2 3 rm(coad.exp.tmp,coad.surv.tmp) rm(survgrouptmp) gc()
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 #组合成完整的数据表 coad.exp <- as.data.frame(t(coad.exp)) #以FPKM>1作为表达阈值，需要在至少70%样本内有表达视为有效基因 coad.exp.use<- coad.exp[colSums(coad.exp>=1)>=314] #如果还需要对样本至少表达10000基因进行筛选，形式类似 #coad.exp.tmp <- coad.exp[rowSums(coad.exp>=1)>=10000] surv.allgene <- cbind(coad.exp.use,coad.surv) surv.allgene[1:5,1:5] #方法1:手写循环 mySurv=with(surv.allgene,Surv(OS.time, OS)) log_rank_p <- apply(surv.allgene, 2 , function(gene){ surv.allgene$group=ifelse(gene>median(gene),'high','low') data.survdiff=survdiff(mySurv~group,data=surv.allgene) p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1) return(p.val) }) log_rank_p=sort(log_rank_p) log_rank_p_sig <- as.data.frame(log_rank_p[log_rank_p < 0.05]) kmsiggene <- row.names(log_rank_p_sig) kmsiggene <- kmsiggene[c(-1,-2)]
可以看到，最后我们获得了740个log_rank_p显著的基因，可以使用这种方法从大量的基因中缩小基因群体，整个运算速度也很快，后续可以接批量的绘图就能得到全部显著的单基因km了。当然，也可以根据功能、DEG、相关性等进行进一步的基因筛选，最后得到最佳的功能基因集。此外，此处只考虑了编码基因，如果要进行调控网络分析，还可以纳入非编码的基因。
1 2 3 4 5 6 7 8 9 #方法2，tinyarray::surv_km，既然有手写，那肯定有大佬已经帮忙写好了，这里就推荐一个最佳选择，tinyarray中批量生存的功能。surv_km只需要有表达矩阵和匹配的meta即可，返回log-rank-p值 require(tinyarray) #tinyarray数据格式需要将生存的名字改为event和time coad.exp.use <- as.data.frame(t(coad.exp.use)) names(coad.surv) <- c("sample","event","time") log_rank_p <- surv_KM(coad.exp.use,coad.surv,cut.point = FALSE) log_rank_p_sig <- as.data.frame(log_rank_p[log_rank_p < 0.05]) kmsiggene <- row.names(log_rank_p_sig)
结果一致，也是740个基因，基因也一致，表明这两种方法的确是等价的，但是看surv_km的源码，它用的是for循环，但是手写部分用的apply，就R来说应该选择前者，虽然用tinyarray的方法看上去更简洁。
变体问题1，统计检验方法问题，log-rank-p的适用性和检验power
log-rank-p的检验对等比例风险假定没有特别要求，但是根据搜索，有文章提到如果不满足PH假设，检验效能会受到影响而削弱，而在满足PH假定的前提下，log-rank检验的效能是最强的。此外，如果生存曲线有交叉，则存在混杂因素，理论上应该选择cox多因素分析之后去除混杂因素。log rank检验对于远期变化敏感，各时间点的存活死亡权重一样，具体到结果里就是，看上去近期曲线分得不是很开，但log-rank-p仍然显著。
然而这并不是说为了让p值显著，就应该“近处显著选近，远处显著选远，都不显著选中”，一切基于具体问题，比如研究短期疗效，就应该选择Breslow-Wilcoxon检验近期效能，要研究长期效应，就应该选择log-rank
变体问题2，log-rank、Breslow-Wilcoxon与Tarone-Ware
在survdiff中，有一个参数是rho，根据文档：
就是rho=0就是默认的log-rank检验，rho=1则是广义Wilcoxon（Gehan-Wilcoxon和Breslow-Wilcoxon就是一个东西），一个远期一个近期，具体到代码里，只需要在survdiff的那一行加上rho=0/1即可。
变体问题3，KM调整，最佳cutoff与模型调整
很多时候，以median、mean或者其他常见统计量作为截断点得到的生存并不显著，这既和检验方法有关，也和具体数据的分布有关，检验方法的问题上面已经说了，长期远期要选对；其二就是数据分布的问题，mean、median等统计量其实只代表数据的总体特征，但不代表这种特征和其他变量存在关联，因此并不一定合适，这时就可以借用survminer的surv_cutpoint函数，就可以返回所谓的最佳截断点，这个数值代表了所选指标在统计学上的最优位置，简单来说就是截断点数值的benchmark然后返回p值或者某种评价指标最佳的数值，如果最佳截断点分群也无法让分组差异变得明显，那就说明在当前数据下，该数值模型本身就是不显著，要么增删新的变量，要么调整建模方法，要么接受结果。例如：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 #选几个不怎么显著的基因 surv.tmp <- subset(surv.allgene,select = c("OS.time","OS","HSPA1A","XPC")) #使用median,不怎么显著 surv.tmp$group <- ifelse(surv.tmp$XPC>median(surv.tmp$XPC),'high','low') survdiff(Surv(OS.time,OS)~group,data=surv.tmp) #使用mean，更不显著了 surv.tmp$group <- ifelse(surv.tmp$XPC>mean(surv.tmp$XPC),'high','low') survdiff(Surv(OS.time,OS)~group,data=surv.tmp) #使用最佳截断,虽然还是不显著，但比最初好太多了.....(其实意思就是没得救了) res.cut <- surv_cutpoint(surv.tmp, time = "OS.time", event = "OS",variables = "XPC") surv.tmp$group <- ifelse(surv.tmp$XPC>mean(res.cut$cutpoint$cutpoint),'high','low') survdiff(Surv(OS.time,OS)~group,data=surv.tmp) #看一下数据分布 plot(res.cut, "XPC", palette = "npg")
大致能看到cutoff的原理，以及选择的cutoff处于数据中的位置，这里的cutoff明显是偏得比较多了，这种情况下即便能够让模型变得显著也最好别用。
因此能够得到对KM优化的基本原则：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 surv.tmp$group <- ifelse(surv.tmp$HSPA1A>median(surv.tmp$HSPA1A),'high','low') surv.fit <- survfit(Surv(OS.time,OS)~group,data=surv.tmp) #简单曲线,一行即可 #plot(surv.fit)开玩笑的，这种线算个屁的可视化 ggsurvplot(surv.fit,pval.method = T,pval = T) #标准曲线,根据医学文献，标准的km曲线应当包含统计量，统计方法和风险表。 ggsurvplot(surv.fit, pval.method = T,#是否添加统计检验的方法 pval = T,#是否添加p值 conf.int = T,#是否展示置信区间 risk.table = T,#是否展示风险进度表格 surv.median.line = "hv",#是否绘制50%生存线 tables.theme = theme_classic2(),#风险进度表的主题 ggtheme = theme_classic2(),#生存曲线图的主题 ) #进阶修改 ggsurvplot(surv.fit, palette = "uchicago",#改配色，个人觉得这个配色最好 pval.method = T,#添加统计检验的方法 pval = T,#添加p值 conf.int = T,#展示置信区间 conf.int.style = "step",#置信区间风格,可选ribbon和step #conf.int.fill设置置信区间填充的颜色 #conf.int.alpha指定置信区间填充颜色的透明度； risk.table = T,#展示风险进度表格 surv.median.line = "hv",#是否绘制50%生存线 tables.theme = theme_classic2(),#风险进度表的主题 ggtheme = theme_classic2(),#生存曲线图的主题 add.all = TRUE,#是否添加全部患者的生存情况，一般这条线应该在中间 legend="right", #指定图例的位置 legend.title = "Group By Gene", #指定图例组的名字 legend.labs = c("All","HighGroup","LowGroup"), #指定图例标签的字符向量 axes.offset = FALSE#指定坐标轴从0开始 #fun = "cumhaz" 使用该fun则会绘制累计风险曲线 )
ggsurvplot的内置参数可谓相当地多，基本所有图上面的东西，从字体到位置全部能直接调整，唯一的问题是，ggsurvplot返回的不是ggplot2对象，反而不能适配一些ggplot2调整的额外包，如果想要像ggplot2一样用的话，使用print(p)即可进行组合或者pdf保存。
可视化的部分，最后的进阶可视化其实没啥必要，也不是最重要的部分，重要的是前面统计的内容。
基本上KM的内容就到此为止，更加进阶的KM，比如分层分时其实没什么必要，那种数据用COX或者其他机器学习的模型，加上多因素的筛选之类的综合方法来做更合适，不要强求KM。
另外，说一个最重要的点，KM曲线除了p值之外，另外一个重要的东西就是曲线是否交叉，理论上KM曲线应该是不会交叉的，如果交叉说明当前的分组存在问题，多变量之间有混杂效应，而如果只有一个自变量分组依然存在交叉，说明PH假设有问题，这种情况可以考虑进行数据截断。
最后再次提醒，KM是一个相对简单的计算方法，不要指望能够完成什么惊世骇俗的结果，真正的结果也就一张图和一个统计结果而已。
最后的最后，推荐一下《总之就是非常可爱》动漫～