怎么获得Hazard Ratio？

生信喵

对于医学生来说，尤其是你的研究方向与肿瘤相关，那么避免不了的就是生存分析。生存分析的目的就是为了有个策略可以指示患者的生存期，以及寻找靶标继而作用于这些靶标是否可以延长患者的生存期。

学了统计的同学都知道RR（Relative Risk）是相对危险度，是指 2 个人群发病率的比值，通常为暴露人群的发病率和非暴露人群的发病率之比，计算公式就是[RR＝暴露组的发病或死亡率/非暴露组的发病或死亡率]。举个经典的例子就是吸烟与肺癌的关系，比如1991年随机调查了某社区10000名居民的吸烟情况，其中3000人吸烟，7000人不吸烟，且假定吸烟的人永远不会戒烟，不吸烟的人永远不会吸烟。随访30年后，2021年吸烟人群中有300人患了肺癌，不吸烟的人群中只有70人患肺癌。因此，与非吸烟人群相比，吸烟人群发生肺癌的相对危险度RR=（300/3000）/(70/7000)=10，即可以认为吸烟人群 30 年内发生肺癌的风险是非吸烟人群的 10 倍。

HR(Hazard Ratio)指的风险比，主要用于队列研究的生存分析，由 Cox 风险比例模型衍生出来。 HR 的解释与 RR 相似，即表示暴露组患病的概率为非暴露组的多少倍。区别在于 RR 只考虑结局是否发生，而 HR 还考虑了结局发生的时间，因此可以认为 HR 是考虑了时间因素的RR。

顺便再看一个大家很容易混淆的两个概念，即Kaplan-Meier 曲线(简称K-M曲线) 和Cox回归。

1.当我们需要描述生存分析数据时，我们常常使用生存曲线来展示，这个时候需要估计每个时间点位的生存率，用的就是K-M方法。因此准确来说，K-M方法是一种统计描述方法，就好比用饼状图来展示比例，用箱图来展示连续变量。K-M曲线还可以很直观地表现出两组或多组的生存率或死亡率，非常适合在文章中进行展示。但仅有描述是不够的，我们还需要进行统计推断，也就是比较。因为生存数据往往都是非正态分布，因此常使用非参数的检验方法，也就是log-rank检验，这就类似于对于非正态连续变量比较使用的方差分析。因此在实际写作中，K-M曲线都会搭配log-rank检验的P值，这样才算完整。

2.Cox回归是一种回归模型，它没有直接使用生存时间，而是使用了危险度（hazard ratio）作为因变量，该模型不用于估计生存率，而是用于因素分析，也就是找到某一个危险因素对结局事件发生的贡献度。

另外cox回归也可以考虑多因素回归，而K-M的log-rank检验只能进行单因素分析，这也是重要的一点区别。

本篇主要内容分为两部分，一是单因素cox回归，另一个是对结果可视化的森林图。

开始主线任务：

1. library(survival)
   
2. library(survminer) #载入R包

复制代码

1. load('1.Rda')
   
2. candidate <- colnames(rt)[3:10] #作为示例就选前8个基因吧，自定义更换自己喜欢的基因
   
3. rt <- rt[,c('futime','fustat',candidate)] #选中生存信息和基因列
   
4. rt[1:3,] #查看数据结构
   

复制代码

1. uniTab=data.frame() #建立空白数据框用于存结果
   
2. for(i in colnames(rt[,3:ncol(rt)])){
   
3.   cox <- coxph(Surv(futime, fustat) ~ rt[,i], data = rt)
   
4.   coxSummary = summary(cox)
   
5.   uniTab=rbind(uniTab,
   
6.                cbind(id=i,
   
7.                      HR=coxSummary$conf.int[,"exp(coef)"],
   
8.                      HR.95L=coxSummary$conf.int[,"lower .95"],
   
9.                      HR.95H=coxSummary$conf.int[,"upper .95"],
   
10.                      pvalue=coxSummary$coefficients[,"Pr(>|z|)"])
    
11.   )
    
12. } #循环计算每个基因的单因素cox
    
13. write.table(uniTab,file="uniCox.txt",sep="\t",row.names=F,quote=F) #保存结果
    
14. GetFactors_uni <- uniTab$id[which(uniTab$pvalue < 0.2)] %>% as.character() ## 选择单因素Cox结果中，P值＜0.2的变量可进一步纳入多因素Cox分析（此篇文章不介绍），这里8个基因就留下了5个，当然你也可以p值设为0.05。
    

复制代码

1. bioForest=function(coxFile=null,forestFile=null,forestCol=null){
   
2. 
   
3.   rt <- read.table(coxFile,header=T,sep="\t",row.names=1,check.names=F)
   
4.   gene <- rownames(rt)
   
5.   hr <- sprintf("%.3f",rt[        DISCUZ_CODE_3        ]quot;HR")
   
6.   hrLow  <- sprintf("%.3f",rt[        DISCUZ_CODE_3        ]quot;HR.95L")
   
7.   hrHigh <- sprintf("%.3f",rt[        DISCUZ_CODE_3        ]quot;HR.95H")
   
8.   Hazard.ratio <- paste0(hr,"(",hrLow,"-",hrHigh,")")
   
9.   pVal <- ifelse(rt$pvalue<0.001, "<0.001", sprintf("%.3f", rt$pvalue))
   
10. 
    
11.   pdf(file=forestFile, width = 9,height = 6.45)
    
12.   n <- nrow(rt)
    
13.   nRow <- n+1
    
14.   ylim <- c(1,nRow)
    
15.   layout(matrix(c(1,2),nc=2),width=c(3,2.5))
    
16. 
    
17.   xlim = c(0,3)
    
18.   par(mar=c(4,2,2,1))
    
19.   plot(1,xlim=xlim,ylim=ylim,type="n",axes=F,xlab="",ylab="")
    
20.   text.cex=0.8
    
21.   text(0,n:1,gene,adj=0,cex=text.cex)
    
22.   text(1.5-0.5*0.2,n:1,pVal,adj=1,cex=text.cex);text(1.5-0.5*0.2,n+1,'pvalue',cex=text.cex,font=2,adj=1)
    
23.   text(3,n:1,Hazard.ratio,adj=1,cex=text.cex);text(3,n+1,'Hazard ratio',cex=text.cex,font=2,adj=1,)
    
24. 
    
25.   par(mar=c(4,1,2,1),mgp=c(2,0.5,0))
    
26.   xlim = c(0,max(as.numeric(hrLow),as.numeric(hrHigh)))
    
27.   plot(1,xlim=xlim,ylim=ylim,type="n",axes=F,ylab="",xaxs="i",xlab="Hazard ratio")
    
28.   arrows(as.numeric(hrLow),n:1,as.numeric(hrHigh),n:1,angle=90,code=3,length=0.05,col="darkblue",lwd=3.5)
    
29.   abline(v=1,col="black",lty=2,lwd=2)
    
30.   boxcolor = ifelse(as.numeric(hr) > 1, forestCol, forestCol)
    
31.   points(as.numeric(hr), n:1, pch = 15, col = boxcolor, cex=2)
    
32.   axis(1)
    
33.   dev.off()
    
34. }#画图的函数
    
35. 
    
36. bioForest(coxFile="uniCox.txt",forestFile="uniForest.pdf",forestCol="red")#一步画图

复制代码

下面就是展示最终结果啦，喜欢的不要忘记分享此篇推文出去哦！
pdf结果我就截图放出来大家看看吧。

生信喵

代码中需要用到的输入数据：生存信息和表达（表达的数据是log2fpkm）。获取示例数据请在公众号"生信喵实验柴"后台回复“20211112。