十二、三图联动

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.热图，挑出的基因在每个样本中的表达量
*生信技能树课程笔记