1.准备输入数据
load("TCGA-LUAD_sur_model.Rdata")
library(survival)
library(survminer)
library(ggplotify)
library(cowplot)
library(Hmisc)
options(scipen=200)
# http://dni-institute.in/blogs/cox-regression-interpret-result-and-predict/
library(pheatmap)
2.挑选感兴趣的基因构建coxph模型
e=t(exprSet[c('hsa-mir-31','hsa-mir-196b','hsa-mir-766','hsa-mir-519a-1','hsa-mir-375','hsa-mir-187','hsa-mir-331','hsa-mir-101-1'),])
colnames(e)=c('miR31','miR196b','miR766','miR519a1','miR375','miR187','miR331','miR101')
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" "miR31" "miR196b"
#> [13] "miR766" "miR519a1" "miR375" "miR187"
#> [17] "miR331" "miR101"
用survival::coxph()函数构建模型
library(survival)
library(survminer)
s=Surv(time, event) ~ miR31+miR196b+miR766+miR519a1+miR375+miR187+miR331+miR101
model <- coxph(s, data = dat )
summary(model,data=dat)
#> Call:
#> coxph(formula = s, data = dat)
#>
#> n= 513, number of events= 125
#>
#> coef exp(coef) se(coef) z Pr(>|z|)
#> miR31 0.07437 1.07720 0.04644 1.601 0.1093
#> miR196b 0.04605 1.04713 0.03752 1.227 0.2197
#> miR766 -0.03884 0.96190 0.12795 -0.304 0.7615
#> miR519a1 -0.03553 0.96509 0.05922 -0.600 0.5485
#> miR375 -0.11247 0.89363 0.05459 -2.060 0.0394 *
#> miR187 -0.03989 0.96089 0.04793 -0.832 0.4052
#> miR331 -0.08267 0.92066 0.12387 -0.667 0.5045
#> miR101 -0.42612 0.65304 0.13917 -3.062 0.0022 **
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> exp(coef) exp(-coef) lower .95 upper .95
#> miR31 1.0772 0.9283 0.9835 1.1799
#> miR196b 1.0471 0.9550 0.9729 1.1270
#> miR766 0.9619 1.0396 0.7485 1.2361
#> miR519a1 0.9651 1.0362 0.8593 1.0839
#> miR375 0.8936 1.1190 0.8029 0.9945
#> miR187 0.9609 1.0407 0.8747 1.0555
#> miR331 0.9207 1.0862 0.7222 1.1736
#> miR101 0.6530 1.5313 0.4971 0.8578
#>
#> Concordance= 0.662 (se = 0.029 )
#> Likelihood ratio test= 27.81 on 8 df, p=0.0005
#> Wald test = 28.49 on 8 df, p=0.0004
#> Score (logrank) test = 29.28 on 8 df, p=0.0003
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.66232618 0.02855476
library(Hmisc)
options(scipen=200)
with(dat,rcorr.cens(fp,Surv(time, event)))
#> C Index Dxy S.D. n missing
#> 0.33767382 -0.32465235 0.05710951 513.00000000 0.00000000
#> uncensored Relevant Pairs Concordant Uncertain
#> 125.00000000 40702.00000000 13744.00000000 221924.00000000
# 若要找到最佳模型，我们可以进行变量选择，可以采用逐步回归法进行分析
3.自己预测自己
查看?predict.coxph的帮助文档会发现：type参数有不同的取值,type不同则得出的预测值意义也不相同 https://mp.weixin.qq.com/s/hOgiWzvkWZLNKX1LDvrAaw
#fp <- predict(model,dat,type="risk");boxplot(fp)
#fp <- predict(model,dat,type="expected") ;boxplot(fp)
#names(fp) = rownames(meta)
fp <- predict(model,dat) ;boxplot(fp)
4.三图联动
这里的难点在于三个图的横坐标顺序是一致的。本强迫症还想让三个图严格对齐，然而无功而返，由于热图拼图本来就比较难弄，导致patchwork、cowplot都搞不定，多次尝试后只能通过对grid.arrange的调整对齐一下，应该是有别的办法的，如有后续，将在简书中更新
fp_dat=data.frame(patientid=1:length(fp),
fp=as.numeric(sort(fp)))
sur_dat=data.frame(patientid=1:length(fp),
time=meta[order(fp),'time'] ,
event=meta[order(fp),'event'])
#改变点的颜色（思路：先将0和1替换为alive和death，然后将因子的顺序改为death，alive）
sur_dat$event=ifelse(sur_dat$event==0,'alive','death')
sur_dat$event=factor(sur_dat$event,levels = c("death","alive"))
exp_dat = t(e[order(fp),])
###第一个图----
p1=ggplot(fp_dat,aes(x=patientid,y=fp))+geom_point()+theme_bw();p1
#第二个图
p2=ggplot(sur_dat,aes(x=patientid,y=time))+geom_point(aes(col=event))+theme_bw();p2
#第三个图
mycolors <- colorRampPalette(c("green","black", "red"), bias = 1.2)(100)
p3=pheatmap::pheatmap(exp_dat,
col= mycolors,
scale = "row",
breaks = seq(-1,1,length.out = 100),
show_colnames = F,
cluster_cols = F)
p3
#拼图实现三图联动
library(ggplotify)
plots = list(p1,p2,as.ggplot(as.grob(p3)))
library(gridExtra)
lay1 = rbind(c(rep(1,7),NA),
c(rep(2,8)),
c(rep(3,8))) #布局矩阵
grid.arrange(grobs = plots,
layout_matrix = lay1,
heigths = c(1, 2,3))
从上向下三个图分别表示:
1.每个病人的预测值，按照从小到大排序 2.每个病人的生存时间，颜色区分生死 3.热图，挑出的基因在每个样本中的表达量
*生信技能树课程笔记
