开篇先看个风险森林图吧~~

image.png

1.准备输入数据

load("TCGA-KIRCsur_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。

#首先把这5个基因的表达量挑选出来，一行是一个基因，一列是一个样本，根据表达矩阵按行取子集[,]
e=t(exprSet[c('hsa-mir-21','hsa-mir-143','hsa-mir-10b','hsa-mir-192','hsa-mir-183'),])
e=log2(e)
#补充一个小小知识点，就是把上面的中划线替换成下划线
#library(stringr)
#str_replace("hsa-mir-21","-","_")   
#上面是替换一个，如果是要替换所有，可以加上all
#str_replace_all("hsa-mir-21","-","_")   
colnames(e)=c('miR21','miR143','miR10b','miR192','miR183')  #改名
dat=cbind(meta,e)
dat$gender=factor(dat$gender)
dat$stage=factor(dat$stage)
colnames(dat)

用survival::coxph()函数构建模型

library(survival)
library(survminer)
s=Surv(time, event) ~ gender + stage + age + miR21+miR143+miR10b+miR192+miR183   #越往前就是越基础，gender,stage,age然后才是miR，不要乱改代码
#s=Surv(time, event) ~ miR21+miR143+miR10b+miR192+miR183
model <- coxph(s, data = dat )

3.模型可视化-森林图

ggforest 风险森林图

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
library(Hmisc)
options(scipen=200)
with(dat,rcorr.cens(fp,Surv(time, event)))
# 若要找到最佳模型，我们可以进行变量选择，可以采用逐步回归法进行分析

这里只是举个栗子，自己预测自己的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)

可查看两组一些临床参数切割比例

train <- exprSet[,sam]
test <- exprSet[,-sam]
train_meta <- meta[sam,]
test_meta <- meta[-sam,]

prop.table(table(train_meta$stage))
prop.table(table(test_meta$stage)) 
prop.table(table(test_meta$race)) 
prop.table(table(train_meta$race)) 

5.2 切割后的train数据集建模

和上面的建模方法一样。

e=t(train[c('hsa-mir-21','hsa-mir-143','hsa-mir-10b','hsa-mir-192','hsa-mir-183'),])
e=log2(e)
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)
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'),])
e=log2(e)
colnames(e)=c('miR21','miR143','miR10b','miR192','miR183')
test_dat=cbind(test_meta,e)

fp <- predict(model)
summary(model)$concordance
library(Hmisc)
options(scipen=200)
with(test_dat,rcorr.cens(fp,Surv(time, event)))

C-index不到0.6， 模型全靠猜了。

内部评价：训练集
外部评价：测试集
在没有拆分数据的时候我们只有训练集，就自己预测自己

莎莎写于2020.4.30 今天是噶宝生日~🎂

1.准备输入数据

输入数据是肿瘤样本表达矩阵exprSet、临床信息meta

load("sur_model_KIRC.Rdata")

library(randomForest)
library(ROCR)
library(genefilter)
library(Hmisc)

2.构建随机森林模型

输入数据是表达矩阵(仅含tumor样本)和对应的生死。

x=t(log2(exprSet+1))
y=meta$event
#构建模型，一个叫randomForest的函数,运行时间很长,存Rdata跳过
tmp_rf="TCGA_KIRC_miRNA_rf_output.Rdata"
if(!file.exists(tmp_rf)){
  rf_output=randomForest(x=x, y=y,importance = TRUE, ntree = 10001, proximity=TRUE )   #这一步运行要20多分钟，所以设置了跳过，如果再运行自己的数据的时候再修改
  save(rf_output,file = tmp_rf)
}
load(file = tmp_rf)
#top30的基因
varImpPlot(rf_output, type=2, n.var=30, scale=FALSE, 
           main="Variable Importance (Gini) for top 30 predictors",cex =.7)
rf_importances=importance(rf_output, scale=FALSE)
head(rf_importances)
#解释量top30的基因,和图上是一样的，从下往上看。
choose_gene=rownames(tail(rf_importances[order(rf_importances[,2]),],30))   #用tail函数取了后30个

image.png
image.png
image.png
这样就实现了倒序

3.模型预测和评估

3.1自己预测自己

rf.prob <- predict(rf_output, x)
re=cbind(y ,rf.prob)
head(re)

3.2 箱线图

对预测结果进行可视化。以实际的生死作为分组，画箱线图整体上查看预测结果。

re=as.data.frame(re)
colnames(re)=c('event','prob')
re$event=as.factor(re$event)
library(ggpubr) 
p1 = ggboxplot(re, x = "event", y = "prob",
               color = "event", palette = "jco",
               add = "jitter")+ stat_compare_means()
p1

对比lasso回归，这个似乎更好一些？

3.3 ROC曲线

library(ROCR)
#library(caret)

pred <- prediction(re[,2], re[,1])
auc = performance(pred,"auc")@y.values[[1]]
auc

自己预测自己，auc值竟然是1。

4.切割数据构建模型并预测

4.1 切割数据

用R包caret切割数据，生成的结果是一组代表列数的数字，用这些数字来给表达矩阵和meta取子集即可。

library(caret)
set.seed(12345679)
sam<- createDataPartition(meta$event, p = .5,list = FALSE)

train <- exprSet[,sam]
test <- exprSet[,-sam]
train_meta <- meta[sam,]
test_meta <- meta[-sam,]

4.2 切割后的train数据集建模

和上面的建模方法一样。

#建模
x = t(log2(train+1))
y = train_meta$event
tmp_rf="TCGA_KIRC_miRNA_t_rf_output.Rdata"
if(!file.exists(tmp_rf)){
  rf_output=randomForest(x=x, y=y,importance = TRUE, ntree = 10001, proximity=TRUE )
  save(rf_output,file = tmp_rf)
}
load(file = tmp_rf)

choose_gene=rownames(tail(rf_importances[order(rf_importances[,2]),],30))
length(choose_gene)

5.模型预测

用训练集构建模型，预测测试集的生死。

x = t(log2(test+1))
y = test_meta$event
rf.prob <- predict(rf_output, x)
re=cbind(y ,rf.prob)
re=as.data.frame(re)
colnames(re)=c('event','prob')
re$event=as.factor(re$event)
library(ggpubr) 
p1 = ggboxplot(re, x = "event", y = "prob",
               color = "event", palette = "jco",
               add = "jitter")+ stat_compare_means()
p1

再看AUC值。

library(ROCR)
# 训练集模型预测测试集
pred <- prediction(re[,2], re[,1])
auc= performance(pred,"auc")@y.values[[1]]
auc

还行。

心理感受：怎么说呢，这部分内容对我来说太抽象了，我不是很懂，具体哪里不懂也不知道，就是觉得很难懂，很难理解，后面也不知道怎么用，还有一个支持向量机，还有timeROC和风险因子评估的也没有搞懂，代码就不贴了，慢慢学习吧~~~😔

小洁老师画的图，真好看