生存分析有必要把连续值依据中位值进行高低分组变成分类变量吗

前面的教程：estimate或者CIBERSORT结果真的是很好的临床预后指标吗，我们针对 estimate 的StromalSignature 和  ImmuneSignature 这样的打分值进行了生存分析。

但是呢，我们其实是根据每个癌症内部自己的 estimate 的StromalSignature 和  ImmuneSignature的打分的中位值，首先分成为了高低两个组，然后进行生存分析看是否有统计学显著。estimate 的打分本身是超级简单， 如果你还不懂就去看前面的教程：不同癌症内部按照estimate的两个打分值高低分组看蛋白编码基因表达量差异 ：

全部的癌症批量就可以跑完生存分析，然后我们查看了BRCA这个癌症的结果， estimate 算法得到的StromalSignature 和  ImmuneSignature都是可以区分生存，因为p值都是0.05附近，结合生存分析的图表，可以看到：

• 其中StromalSignature 高死得快，是风险因子。
• 而ImmuneSignature高死的慢，是保护因子。

然后有小伙伴就留言了，为什么要把连续值依据中位值进行高低分组变成分类变量，然后使用survdiff来做两个组的统计检验呢，既然是连续值，可以直接cox方法啊！

ox方法的操作代码如下所示：

rm(list = ls()) 
load('../expression/estimate_results/ssgseaScore-for-ssGSEA_for_seurat.Rdata')
head(ssgseaScore)
ssgseaScore$pid=substring(rownames(ssgseaScore),1,12)
head(ssgseaScore)

library(survminer) 
library(survival)
load(file = 'phe_stage.Rdata')
phe[1:4,1:4] 
phe$OS.time =  phe$OS.time/365
tp=unique(phe$type)
tp
options(digits = 2)

cox_results <- do.call(rbind,
                      lapply(tp, function(x){
                        # x=tp[1]
                        this_phe=phe[phe$type==x,]
                        # 这里先看 os 
                        survival_dat=as.data.frame(this_phe[,c('bcr_patient_barcode','OS','OS.time')])
                        head(survival_dat) 
                        colnames(survival_dat)=c('pid','event','time')
                        survival_dat=merge(survival_dat,ssgseaScore,by='pid')  
                        colnames(survival_dat)
                        attach(survival_dat)
                        s=Surv(time, event) ~    StromalSignature+ImmuneSignature 
                        model <- coxph(s, data = survival_dat )
                        sm=summary(model,data=survival_dat) 
                        return( c(sm$coefficients[1,c(5,2)],
                                  sm$coefficients[2,c(5,2)]))
                      }))
rownames(cox_results)=tp
round(cox_results,2)
save(cox_results,file = 'cox_results_for_estimate.Rdata')

load(file = 'km_results_for_estimate.Rdata')
tmp=cbind(cox_results,km_results)
colnames(tmp)=NULL
round(tmp,2)


可以看到两次生存分析结果还是有比较好的一致性，为了节省空间，下面的表格结合了cox和km的两种生存分析结果，都是  stromal_p.val,stromal_HR,  immune_p.val,immune_HR的顺序。

前面的4列是cox结果，后面的4列是km的结果。可以看到cox的生存分析把打分当做是连续变量，计算得到的HR值非常的大，但是km方法把打分根据中位值进行了高低分组，得到的HR整体低很多！

> round(tmp,2)
     [,1]     [,2] [,3]    [,4] [,5] [,6] [,7] [,8]
ACC  0.68     5.45 0.28    0.03 0.70 0.86 0.03 0.43
BLCA 0.00    20.16 0.00    0.07 0.03 1.38 0.85 0.97
BRCA 0.00    18.18 0.00    0.09 0.06 1.31 0.04 0.75
CESC 0.13     5.88 0.01    0.07 0.47 0.84 0.06 0.64
CHOL 0.91     1.48 0.45    0.11 0.15 0.54 0.04 0.39
COAD 0.14     5.51 0.28    0.22 0.41 1.17 0.81 1.05
DLBC 0.70     0.07 0.94    0.62 0.87 0.90 0.63 1.35
ESCA 0.95     1.08 0.84    0.72 0.66 1.11 0.68 1.10
GBM  0.26     9.45 0.74    0.62 0.20 1.24 0.21 1.24
HNSC 0.59     1.41 0.01    0.17 0.23 0.86 0.04 0.77
KICH 0.16   217.22 0.49    0.06 0.76 1.19 0.64 0.76
KIRC 0.06     0.15 0.01    5.62 0.78 0.96 0.13 1.24
KIRP 0.00   135.75 0.12    0.14 0.06 1.70 0.83 0.94
LAML 0.01     0.00 0.00 7617.42 0.80 1.05 0.01 1.70
LGG  0.03    81.55 0.75    1.60 0.00 1.77 0.01 1.58
LIHC 0.71     0.64 0.91    1.12 0.65 1.07 0.75 1.05
LUAD 0.71     1.38 0.17    0.33 0.16 0.82 0.07 0.78
LUSC 0.26     2.30 0.84    1.15 0.03 1.32 0.18 1.19
MESO 0.01    46.32 0.08    0.10 0.01 1.83 0.45 0.84
OV   0.04     4.64 0.12    0.35 0.63 1.06 0.36 0.89
PAAD 0.40     3.12 0.92    0.87 0.54 1.13 0.80 0.95
PCPG 0.98     0.87 0.46    0.00 0.17 0.34 0.64 0.72
PRAD 0.36     0.01 0.91    0.50 0.22 0.45 0.59 0.71
READ 0.44     7.13 0.74    0.35 0.16 1.65 0.80 0.91
SARC 0.85     0.79 0.24    0.34 0.06 0.68 0.01 0.59
SKCM 0.06     4.66 0.00    0.07 0.04 0.76 0.00 0.57
STAD 0.00    10.25 0.16    0.34 0.10 1.30 0.33 1.17
TGCT 0.04 46133.57 0.12  137.01 0.04  Inf 0.19 3.95
THCA 0.01 11888.26 0.01    0.00 1.00 1.00 0.75 0.86
THYM 0.84     1.96 0.51   12.86 0.92 1.07 0.73 0.79
UCEC 0.84     1.28 0.02    0.08 0.48 0.87 0.01 0.60
UCS  0.79     0.60 0.79    1.61 0.68 1.15 0.73 0.89
UVM  0.82     3.06 0.01  164.61 0.01 2.86 0.00 4.38


对上面的表格进行简单的统计：

> table(tmp[,1] > 0.05 ,tmp[,5] > 0.05)
# cox和km对 StromalSignature 指标判断为显著，重合的癌症就3个，各自有自己的特殊性       
        FALSE TRUE
  FALSE     4    6
  TRUE      3   20
> table(tmp[,3] > 0.05 ,tmp[,7] > 0.05)
       
        FALSE TRUE
  FALSE     6    4
  TRUE      4   19
> table(tmp[,2] > 1 ,tmp[,6] > 1)
# 可以看到   StromalSignature 在绝大部分癌症都是风险因子，因为HR值大于1     
        FALSE TRUE
  FALSE     5    3
  TRUE      7   18
> table(tmp[,4] > 1 ,tmp[,8] > 1)
# 可以看到   ImmuneSignature 在绝大部分癌症都是风险因子，因为HR值小于1          
# 而且两个方法异一致性还行
        FALSE TRUE
  FALSE    19    5
  TRUE      2    7


可以看到，无论是从p值0.05这个阈值的角度来看，cox和km的生存分析是否有统计学意义的一致性都还行吧！

另外，从HR值角度看 cox和km对该因素的风险因子和保护因子的判断也是勉强可以的！

主要是可以看看卡方检验的p值，如下所示：

> chisq.test(table(tmp[,1] > 0.05 ,tmp[,5] > 0.05)) 
X-squared = 2, df = 1, p-value = 0.2 
> chisq.test(table(tmp[,3] > 0.05 ,tmp[,7] > 0.05))  
X-squared = 4, df = 1, p-value = 0.04 
> chisq.test(table(tmp[,2] > 1 ,tmp[,6] > 1)) 
X-squared = 2, df = 1, p-value = 0.2 
> chisq.test(table(tmp[,4] > 1 ,tmp[,8] > 1)) 
X-squared = 7, df = 1, p-value = 0.009


这就有点麻烦，在0.05附近蹦跶，很难下结论。

最后你会发现其实上面的数值不直观

所以，需要一个很好的可视化！我随便写了一点：

library(ggpubr) 
library(ggplot2)
library(patchwork)

check <- function(df){
  colnames(df)=c('p.val','HR')
  df$v=log(df$HR)
  df$sig = ifelse(df$p.val < 0.05,'sig','no')
  table(df$sig)
  df$cancer=rownames(df)
  ggdotchart(df, y = "v", x = "cancer",
             color = "sig",       
             sorting = 'none',
             add = "segments", 
             rotate = TRUE,
             xlab=""
  ) 
}
df=as.data.frame(tmp[,1:2])
p1=check(df)
df=as.data.frame(tmp[,5:6])
p2=check(df)

p1+p2

出图如下：

 

可以看到ACC这个癌症就是cox和km算出来的的HR值反过来了，对stromal来说。