八、cox-forest
2022-01-22 本文已影响0人
白米饭睡不醒
1.准备输入数据
load("TCGA-KIRC_sur_model.Rdata")
2.挑选感兴趣的基因构建coxph模型
出自文章Integrated genomic analysis identifies subclasses and prognosis signatures of kidney cancer中，五个miRNA是miR-21,miR-143,miR-10b,miR-192,miR-183
将他们从表达矩阵中取出来，就得到了5个基因在522个肿瘤样本中的表达量，可作为列添加在meta表格的后面,组成的数据框赋值给dat。
#输入数据为几个目的基因的表达矩阵转置一下，再和临床信息粘贴在一起
e=t(exprSet[c('hsa-mir-21','hsa-mir-143','hsa-mir-10b','hsa-mir-192','hsa-mir-183'),])
colnames(e)=c('miR21','miR143','miR10b','miR192','miR183')#列名中不能有特殊字符，所以做了一步简化去掉减号
dat=cbind(meta,e)
dat$gender=factor(dat$gender)
dat$stage=factor(dat$stage)
colnames(dat)
#> [1] "ID" "event" "death" "last_followup"
#> [5] "race" "age" "gender" "stage"
#> [9] "time" "age_group" "miR21" "miR143"
#> [13] "miR10b" "miR192" "miR183"
用survival::coxph()函数构建模型
library(survival)
library(survminer)
s=Surv(time, event) ~ gender + stage + age + miR21+miR143+miR10b+miR192+miR183
#s=Surv(time, event) ~ miR21+miR143+miR10b+miR192+miR183
model <- coxph(s, data = dat )
3.模型可视化-森林图
options(scipen=1)
ggforest(model, data =dat,
main = "Hazard ratio",
cpositions = c(0.10, 0.22, 0.4),
fontsize = 1.0,
refLabel = "1", noDigits = 4)
4.模型预测
fp <- predict(model)
summary(model)$concordance
#> C se(C)
#> 0.79000355 0.01796942
library(Hmisc)
options(scipen=200)
with(dat,rcorr.cens(fp,Surv(time, event)))
#> C Index Dxy S.D. n missing
#> 0.20999645 -0.58000710 0.03593883 522.00000000 0.00000000
#> uncensored Relevant Pairs Concordant Uncertain
#> 158.00000000 90192.00000000 18940.00000000 181754.00000000
# 若要找到最佳模型，我们可以进行变量选择，可以采用逐步回归法进行分析
这里只是举个栗子，自己预测自己的C-index是1-with(dat,rcorr.cens(fp,Surv(time, event)))[[1]],实战应该是拿另一个数据集来预测，或者将一个数据集分两半，一半构建模型，一半验证,可以使用机器学习的R包切割数据。
C-index用于计算生存分析中的COX模型预测值与真实之间的区分度（discrimination），也称为Harrell’s concordanceindex。C-index在0.5-1之间。0.5为完全不一致,说明该模型没有预测作用,1为完全一致,说明该模型预测结果与实际完全一致。
5.切割数据构建模型并预测
5.1 切割数据
用R包caret切割数据，生成的结果是一组代表列数的数字，用这些数字来给表达矩阵和meta取子集即可。
library(caret)
set.seed(12345679)
sam<- createDataPartition(meta$event, p = .5,list = FALSE)
head(sam)
#> Resample1
#> [1,] 5
#> [2,] 9
#> [3,] 13
#> [4,] 17
#> [5,] 19
#> [6,] 22
可查看两组一些临床参数切割比例
train <- exprSet[,sam]
test <- exprSet[,-sam]
train_meta <- meta[sam,]
test_meta <- meta[-sam,]
prop.table(table(train_meta$stage))
#>
#> i ii iii iv
#> 0.4636015 0.1072797 0.2796935 0.1494253
prop.table(table(test_meta$stage))
#>
#> i ii iii iv
#> 0.5249042 0.1111111 0.1954023 0.1685824
prop.table(table(test_meta$race))
#>
#> asian black or african american white
#> 0.01171875 0.08593750 0.90234375
prop.table(table(train_meta$race))
#>
#> asian black or african american white
#> 0.01937984 0.13953488 0.84108527
5.2 切割后的train数据集建模
和上面的建模方法一样。
e=t(train[c('hsa-mir-21','hsa-mir-143','hsa-mir-10b','hsa-mir-192','hsa-mir-183'),])
colnames(e)=c('miR21','miR143','miR10b','miR192','miR183')
dat=cbind(train_meta,e)
dat$gender=factor(dat$gender)
dat$stage=factor(dat$stage)
colnames(dat)
#> [1] "ID" "event" "death" "last_followup"
#> [5] "race" "age" "gender" "stage"
#> [9] "time" "age_group" "miR21" "miR143"
#> [13] "miR10b" "miR192" "miR183"
s=Surv(time, event) ~ gender + stage + age + miR21+miR143+miR10b+miR192+miR183
#s=Surv(time, event) ~ miR21+miR143+miR10b+miR192+miR183
model <- coxph(s, data = dat )
5.3 模型可视化
options(scipen=1)
ggforest(model, data =dat,
main = "Hazard ratio",
cpositions = c(0.10, 0.22, 0.4),
fontsize = 1.0,
refLabel = "1", noDigits = 4)
5.4 用切割后的数据test数据集验证模型
e=t(test[c('hsa-mir-21','hsa-mir-143','hsa-mir-10b','hsa-mir-192','hsa-mir-183'),])
colnames(e)=c('miR21','miR143','miR10b','miR192','miR183')
test_dat=cbind(test_meta,e)
fp <- predict(model)
summary(model)$concordance
#> C se(C)
#> 0.75981243 0.02639671
library(Hmisc)
options(scipen=200)
with(test_dat,rcorr.cens(fp,Surv(time, event)))
#> C Index Dxy S.D. n missing
#> 0.45294331 -0.09411337 0.07194385 261.00000000 0.00000000
#> uncensored Relevant Pairs Concordant Uncertain
#> 74.00000000 22016.00000000 9972.00000000 45842.00000000
# C-index为1减0.45294331
C-index不到0.6， 模型全靠猜了。
*生信技能树课程笔记
