R里面实现生存分析非常简单

代码来自老大github：https://github.com/jmzeng1314/tcga_example

用my.surv <- surv(OS_MONTHS,OS_STATUS=='DECEASED')构建生存曲线。

用kmfit2 <- survfit(my.surv~TUMOR_STAGE_2009)来做某一个因子的KM生存曲线。

用 survdiff(my.surv~type, data=dat)来看看这个因子的不同水平是否有显著差异，其中默认用是的logrank test 方法。

用coxph(Surv(time, status) ~ ph.ecog + tt(age), data=lung) 来检测自己感兴趣的因子是否受其它因子(age,gender等等)的影响。

rm(list=ls())
options(stringsAsFactors = F)

Rdata_dir='../Rdata/'
Figure_dir='../figures/'
# 加载上一步从RTCGA.miRNASeq包里面提取miRNA表达矩阵和对应的样本临床信息。
load( file = 
        file.path(Rdata_dir,'TCGA-KIRC-miRNA-example.Rdata')
)
dim(expr)
dim(meta)
# 可以看到是 537个病人，但是有593个样本，每个样本有 552个miRNA信息。
# 当然，这个数据集可以下载原始测序数据进行重新比对，可以拿到更多的miRNA信息

# 这里需要解析TCGA数据库的ID规律，来判断样本归类问题。
group_list=ifelse(as.numeric(substr(colnames(expr),14,15)) < 10,'tumor','normal')

table(group_list)
exprSet=na.omit(expr)
dim(expr)

library(survival)
library(survminer)

# 这里做生存分析，已经不需要正常样本的表达矩阵了，所以需要过滤。
# 而且临床信息，有需要进行整理。
### survival analysis only for patients with tumor.
if(F){
  exprSet=na.omit(expr)
  exprSet=exprSet[,group_list=='tumor']
  
  head(meta)
  colnames(meta)
  meta[,3][is.na(meta[,3])]=0
  meta[,4][is.na(meta[,4])]=0
  meta$days=as.numeric(meta[,3])+as.numeric(meta[,4])
  meta=meta[,c(1:2,5:9)]
  colnames(meta)
  colnames(meta)=c('ID','event','race','age','gender','stage',"days")

  library(survival)
  library(survminer)
  meta$event=ifelse(meta$event=='alive',0,1)
  meta$age=as.numeric(meta$age)
  library(stringr) 
  meta$stage=str_split(meta$stage,' ',simplify = T)[,2]
  table(meta$stage)
  boxplot(meta$age)
  meta$age_group=ifelse(meta$age>median(meta$age),'older','younger')
  table(meta$race)
  meta$time=meta$days/30
  phe=meta
  
  head(phe)
  phe$ID=toupper(phe$ID) 
  phe=phe[match(substr(colnames(exprSet),1,12),phe$ID),]
  head(phe)
  exprSet[1:4,1:4]
  
  save(exprSet,phe,
       file = 
         file.path(Rdata_dir,'TCGA-KIRC-miRNA-survival_input.Rdata')
      )
}
# 上面被设置为if(F)的代码，就是在整理临床信息和生存分析的表达矩阵。
# 整理好的数据，直接加载即可
load(file = 
         file.path(Rdata_dir,'TCGA-KIRC-miRNA-survival_input.Rdata')
)

整理出的exprSet

image-20200103000945032

整理出的phe

image-20200103000905201

利用ggsurvplot快速绘制漂亮的生存曲线图

# 利用ggsurvplot快速绘制漂亮的生存曲线图
sfit <- survfit(Surv(time, event)~stage, data=phe)
sfit
summary(sfit)
ggsurvplot(sfit, conf.int=F, pval=TRUE)
## more complicate figures.
ggsurvplot(sfit,palette = c("#E7B800", "#2E9FDF",'#FFC1C1','#6A5ACD'),
           risk.table =TRUE,pval =TRUE,
           conf.int =TRUE,xlab ="Time in months", 
           ggtheme =theme_light(), 
           ncensor.plot = TRUE)

image-20200103001445792

## more complicate figures.
ggsurvplot(sfit,palette = c("#E7B800", "#2E9FDF",'#FFC1C1','#6A5ACD'),
           risk.table =TRUE,pval =TRUE,
           conf.int =TRUE,xlab ="Time in months", 
           ggtheme =theme_light(), 
           ncensor.plot = TRUE)

image-20200102235244702

多个 ggsurvplots作图生存曲线代码合并

sfit1=survfit(Surv(time, event)~gender, data=phe)
sfit2=survfit(Surv(time, event)~age_group, data=phe)
splots <- list()
splots[[1]] <- ggsurvplot(sfit1,pval =TRUE, data = phe, risk.table = TRUE)
splots[[2]] <- ggsurvplot(sfit2,pval =TRUE, data = phe, risk.table = TRUE)
# Arrange multiple ggsurvplots and print the output
arrange_ggsurvplots(splots, print = TRUE,  ncol = 2, nrow = 1, risk.table.height = 0.4)
dev.off()

用的arrange_ggsurvplots函数

image-20200103001906332

可以很明显看到，肿瘤病人的生存受着诊断癌症的年龄的影响，却与性别无关。

在相对年长的时候诊断的癌症患者通常会死的快一点。

挑选感兴趣的基因做生存分析

• 来自于文章：2015-TCGA-ccRCC-5-miRNAs-signatures

  Integrated genomic analysis identifies subclasses and prognosis signatures of kidney cancer miR-21,miR-143,miR-10b,miR-192,miR-183

tmp=as.data.frame(rownames(exprSet))
g1='hsa-mir-21' # p value = 0.0059
g2='hsa-mir-143' # p value = 0.0093
g3='hsa-mir-192' # p value = 0.00073
g4='hsa-mir-183' # p value = 0.00092
g5='hsa-mir-10b' # p value < 0.0001
gs=c('hsa-mir-21','hsa-mir-143','hsa-mir-192',
     'hsa-mir-183','hsa-mir-10b') 
splots <- lapply(gs, function(g){
  phe$gene=ifelse(exprSet[g,]>median(exprSet[g,]),'high','low')
  table(phe$gene)
  sfit1=survfit(Surv(time, event)~gene, data=phe)
  ggsurvplot(sfit1,pval =TRUE, data = phe, risk.table = TRUE)
}) 
arrange_ggsurvplots(splots, print = TRUE,  
                    ncol = 2, nrow = 3, risk.table.height = 0.4)
dev.off()

这里还是用的arrange_ggsurvplots函数，可以是结果出现在一张绘图结果中

image-20200103003156875

批量生存分析 使用 logrank test 方法

mySurv=with(phe,Surv(time, event))
log_rank_p <- apply(exprSet , 1 , function(gene){
  # gene=exprSet[1,]
  phe$group=ifelse(gene>median(gene),'high','low')  
  data.survdiff=survdiff(mySurv~group,data=phe)
  p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
  return(p.val)
})
log_rank_p=sort(log_rank_p)
head(log_rank_p)
tail(log_rank_p)
boxplot(log_rank_p)  
table(log_rank_p<0.01)
log_rank_p[log_rank_p<0.01]

上面还是好好理解下，就是把高表达和低表达的基因分开，来做生存分析，用p值来看对生存的影响，但是应该不是都是高表达的才对生存有影响，某些高表达对生存就没有影响；也有可能某些低表达对生存有影响。

文章里面挑选出来的生存分析相关的miRNA基因，在我们的分析里面都是显著的。

> c('hsa-mir-21','hsa-mir-143','hsa-mir-192',
+   'hsa-mir-183','hsa-mir-10b')  %in% names(log_rank_p[log_rank_p<0.01])
[1] TRUE TRUE TRUE TRUE TRUE