生存分析

Surv：用于创建生存数据对象
survfit：创建KM生存曲线或是Cox调整生存曲线
survdiff：用于不同组的统计检验

作图survfit

使用survfit+ggsurvplot：只需要三个参数：分组、time、event

ggsurvplot(survfit(Surv(time, event)~group, data=phe), conf.int=F, pval=TRUE)

统计数据

使用survdiff：同样只需要三个参数：分组、time、event

data.survdiff<- survdiff(Surv(time, event)~group, data=phe)
p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
####### 充分理解p.val计算方法
data.survdiff$chisq
data.survdiff$n
length(data.survdiff$n)  #### 在进行多组计算时，这个值会发生变化，容易报错，或者计算错误
p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
p.val=1-pchisq(4,1)

length(data.survdiff$n) #### 在进行多组计算时，这个值会发生变化，容易报错，或者计算错误

按照上调75%，或者下调25%进行log-rank生存曲线计算的循环代码

gene_name<- 'COL4A2'


if(gene_name%in%row.names(exprSet)){
  phe$group=ifelse(exprSet[gene_name,]>quantile(exprSet[gene_name,],0.75),'high',ifelse(exprSet[gene_name,]<quantile(exprSet[gene_name,],0.25),'low',NA))
  phe_loop<- phe[!is.na(phe$group),]
  table(phe_loop$group)
  data.survdiff<- survdiff(Surv(time, event)~group, data=phe_loop)
  p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)  #### 同时这里计算了P值
  p.val
  ggsurvplot(survfit(Surv(time, event)~group, data=phe_loop), conf.int=F, pval=TRUE,  surv.median.line = "hv",legend.title = gene_name, font.main = c(16, "bold", "darkblue"),
             font.x = c(14, "bold.italic", "red"),
             font.y = c(14, "bold.italic", "darkred"),
             font.tickslab = c(12, "plain", "darkgreen"))
}

想怎么调整就怎么调整

合并作图

基于以下几句代码

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_ggsurvplots(splots, print = TRUE,  ncol = 2, nrow = 1, risk.table.height = 0.4)

最终呈现

############################################################  合并作图
gene_name<- 'PLAU'
if(gene_name%in%row.names(exprSet)){
  phe$group=ifelse(exprSet[gene_name,]>quantile(exprSet[gene_name,],0.6),'high',ifelse(exprSet[gene_name,]<quantile(exprSet[gene_name,],0.4),'low',NA))
  phe_loop<- phe[!is.na(phe$group),]
  table(phe_loop$group)
  data.survdiff<- survdiff(Surv(time, event)~group, data=phe_loop)
  p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
  p.val
  ggsurvplot(survfit(Surv(time, event)~group, data=phe_loop), conf.int=F, pval=TRUE)
  sfit1<- survfit(Surv(time, event)~group, data=phe_loop)
}

# dev.off()

splots <- list()
splots[[1]] <- ggsurvplot(sfit1,pval =TRUE, data = phe, risk.table = FALSE,legend.title = gene_name, 
                          font.main = c(14, "plain", "black"),
                          font.x = c(14, "plain", "black"),
                          font.y = c(14, "plain", "black"),
                          font.tickslab = c(14, "plain", "black"))
############
gene_name<- 'SPP1'
if(gene_name%in%row.names(exprSet)){
  phe$group=ifelse(exprSet[gene_name,]>quantile(exprSet[gene_name,],0.6),'high',ifelse(exprSet[gene_name,]<quantile(exprSet[gene_name,],0.4),'low',NA))
  phe_loop<- phe[!is.na(phe$group),]
  table(phe_loop$group)
  data.survdiff<- survdiff(Surv(time, event)~group, data=phe_loop)
  p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
  p.val
  ggsurvplot(survfit(Surv(time, event)~group, data=phe_loop), conf.int=F, pval=TRUE)
  sfit2<- survfit(Surv(time, event)~group, data=phe_loop)
}

# dev.off()
splots[[2]] <- ggsurvplot(sfit2,pval =TRUE, data = phe,legend.title = gene_name, risk.table = FALSE, 
                          font.main = c(14, "plain", "black"),
                          font.x = c(14, "plain", "black"),
                          font.y = c(14, "plain", "black"),
                          font.tickslab = c(14, "plain", "black"))
############
gene_name<- 'BGN'
if(gene_name%in%row.names(exprSet)){
  phe$group=ifelse(exprSet[gene_name,]>quantile(exprSet[gene_name,],0.6),'high',ifelse(exprSet[gene_name,]<quantile(exprSet[gene_name,],0.4),'low',NA))
  phe_loop<- phe[!is.na(phe$group),]
  table(phe_loop$group)
  data.survdiff<- survdiff(Surv(time, event)~group, data=phe_loop)
  p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
  p.val
  ggsurvplot(survfit(Surv(time, event)~group, data=phe_loop), conf.int=F, pval=TRUE)
  sfit3<- survfit(Surv(time, event)~group, data=phe_loop)
}

# dev.off()

splots[[3]] <- ggsurvplot(sfit3,pval =TRUE, data = phe,legend.title = gene_name, risk.table = FALSE, font.main = c(14, "plain", "black"),
                          font.x = c(14, "plain", "black"),
                          font.y = c(14, "plain", "black"),
                          font.tickslab = c(14, "plain", "black"))
# Arrange multiple ggsurvplots and print the output
res<- arrange_ggsurvplots(splots, print = TRUE,  ncol = 3, nrow = 1,surv.plot.height=0.75, risk.table.height = FALSE)
ggsave('./figure/survival.pdf',res, width = 45, height = 15, units = "cm")

最终效果