仅仅是改变了统计学显著性呢？还是说改变了其本性

生信技能树

发布于 2021-07-29 11:25:46

7450

发布于 2021-07-29 11:25:46

文章被收录于专栏：生信技能树生信技能树

前面我们分享过：学徒数据挖掘之谁说生存分析一定要按照表达量中位值或者平均值分组呢？，可以很容易对一个基因，根据表达量分组，然后进行生存分析，判断它是风险因子或者是保护因子，当然了前提是具有统计学显著性啦。

然后很多粉丝留言说，如果并不是按照表达量中位值或者平均值分组，而是取巧使用了surv_cutpoint这样的函数，得到的结果并不好解释，认为这样的的数据处理方式简直是黑白颠倒！

我看到了这样的留言，确实思考了一下，使用了surv_cutpoint这样的函数是仅仅是改变了统计学显著性指标呢？还是说，确实有可能造成黑白颠倒呢？风险因子会变成保护因子？？？

布置它为学徒作业吧

前面自己把表达量矩阵和临床信息准备好，得到 exprSet这个表达量矩阵，以及 meta 这个临床信息，然后就可以使用了surv_cutpoint这样的函数对指定基因做生存分析啦，代码如下所示：

library(survival)
library(survminer)

gs=c('BCL2','BCL2L1','MCL1')

splots <- lapply(gs, function(g){
  meta$gene = t(exprSet[g,])
  sur.cut <- surv_cutpoint(meta,   time= 'time',
                           event = 'event' ,
                           variables = 'gene' )
  sur.cat <- surv_categorize(sur.cut)
  sfit1=survfit(Surv(time, event)~gene, data=sur.cat)
  ggsurvplot(sfit1,pval =TRUE, data = meta)
}) 
arrange_ggsurvplots(splots, print = TRUE,  
                    ncol = 2, nrow = 2, risk.table.height = 0.4)

如下所示；

这个 surv_cutpoint 函数确实效果给力，基本上没有它搞不定的统计学显著性。

表达量矩阵和临床信息的准备过程

需要自己去tcga所在的ucsc的xena网站 ，去下载 TCGA.LGG.sampleMap_HiSeqV2.gz 文件以及 TCGA.LGG.sampleMap_LGG_clinicalMatrix 文件，然后下面的代码就可以复制粘贴直接运行啦！

1.准备数据

rm(list = ls())
library(stringr)
exp = read.table('../TCGA.LGG.sampleMap_HiSeqV2.gz',header = T,row.names = 1,check.names = F)
exp = as.matrix(exp)
clinical = data.table::fread('../TCGA.LGG.sampleMap_LGG_clinicalMatrix')
dim(exp)
#> [1] 20530   530
clinical[1:4,1:4]
#>           sampleID    _INTEGRATION     _PATIENT                       _cohort
#> 1: TCGA-CS-4938-01 TCGA-CS-4938-01 TCGA-CS-4938 TCGA Lower Grade Glioma (LGG)
#> 2: TCGA-CS-4941-01 TCGA-CS-4941-01 TCGA-CS-4941 TCGA Lower Grade Glioma (LGG)
#> 3: TCGA-CS-4942-01 TCGA-CS-4942-01 TCGA-CS-4942 TCGA Lower Grade Glioma (LGG)
#> 4: TCGA-CS-4943-01 TCGA-CS-4943-01 TCGA-CS-4943 TCGA Lower Grade Glioma (LGG)

2.数据过滤

#过滤一下
dim(exp)
#> [1] 20530   530
table(str_sub(colnames(exp),14,15))
#> 
#>  01  02 
#> 516  14
# 01肿瘤，02复发，11正常
group_list = ifelse(as.numeric(str_sub(colnames(exp),14,15)) < 10,'tumor','normal')
group_list = factor(group_list,levels = c("normal","tumor"))
table(group_list)
#> group_list
#> normal  tumor 
#>      0    530
#这里是根据较小样本（normal）取值，当然可以有别的取值
exp = exp[apply(exp, 1, function(x) sum(x > 1) > 0), ]
dim(exp)
#> [1] 19713   530
exp[1:4,1:4]
#>           TCGA-E1-5319-01 TCGA-HT-7693-01 TCGA-CS-6665-01 TCGA-S9-A7J2-01
#> ARHGEF10L         10.1208         11.8084         10.1872         10.9576
#> HIF3A             10.4506         10.5173          3.6466          4.8099
#> RNF17              0.0000          0.0000          0.0000          0.4657
#> RNF10             11.2171         11.8974         11.5369         11.2675
#save(exp,group_list,clinical,file = '../data/01TCGA-LGG.Rdata')

3.生存分析前准备

这里可以选取自己需要的列，

### 1.简化临床信息，选出要用的列
tmp = data.frame(colnames(clinical))
#colnames(clinical)
#write.csv(clinical,file = 'clinical.csv')
clinical = clinical[,c(
  'bcr_patient_barcode',
  'vital_status',
  "days_to_last_followup"  ,
  'days_to_death',
  'days_to_birth'
)]

dim(clinical)
#> [1] 530   5
#rownames(clinical) <- NULL
clinical = clinical[!duplicated(clinical$bcr_patient_barcode),]
rownames(clinical) <- clinical$bcr_patient_barcode
clinical[1:4,1:4]
#>    bcr_patient_barcode vital_status days_to_last_followup days_to_death
#> 1:        TCGA-CS-4938       LIVING                  3574            NA
#> 2:        TCGA-CS-4941     DECEASED                    NA           234
#> 3:        TCGA-CS-4942     DECEASED                    NA          1335
#> 4:        TCGA-CS-4943     DECEASED                   481          1106

meta = clinical
exprSet=exp[,group_list=='tumor']
table(group_list)
#> group_list
#> normal  tumor 
#>      0    530
exprSet = log2(edgeR::cpm(exprSet)+1)

#简化meta的列名
colnames(meta)=c('ID','event','last_followup','death','age')

#空着的值改为NA
meta[meta==""]=NA

### 2.实现表达矩阵与临床信息的匹配 
# 以病人为中心，表达矩阵按病人ID去重复
k = !duplicated(str_sub(colnames(exprSet),1,12));table(k)
#> k
#> FALSE  TRUE 
#>    14   516
exprSet = exprSet[,k]
exprSet = as.data.frame(exprSet)

#调整meta的ID顺序与exprSet列名一致
exprSet = exprSet[,str_sub(colnames(exprSet),1,12) %in% meta$ID]
meta=meta[match(str_sub(colnames(exprSet),1,12),meta$ID),]
identical(meta$ID,str_sub(colnames(exprSet),1,12))
#> [1] TRUE

4.整理生存分析输入数据

#1.1计算生存时间(月)
meta$time = ifelse(meta$event=="LIVING",
                   meta$last_followup,
                   meta$death)
meta$time = as.numeric(meta$time)/30
#1.2 根据生死定义event，活着是0，死的是1
meta$event=ifelse(meta$event=='LIVING',
                  0,
                  1)
table(meta$event)
#> 
#>   0   1 
#> 389 126

#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 
#>     246     268

#save(meta,exprSet,file = "../data/TCGA-LGG_sur_model.Rdata")

得到 exprSet这个表达量矩阵，以及 meta 这个临床信息，后续的全部生存分析相关，都是依赖于它哦！

作业

对全部的基因，首先使用表达量中位值进行分组后，批量进行 cox 分析，拿到HR值和P值，输出成为表格。

然后对基因根据surv_cutpoint函数进行分组后，再批量cox分析，拿到HR值和P值，输出成为表格。

比较两个表格，看看是否有基因的HR值的方向冲突了，还是说，仅仅是统计学指标的改变。

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2021-07-15，如有侵权请联系 cloudcommunity@tencent.com 删除

html

本文分享自生信技能树微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

html

登录后参与评论

0 条评论

热度