Posted on

比值比(Odds ratio, OR):又名机会比,优势比,交叉乘积比( Cross-product Ratio),相对比值( Relative Odds ),两个比值的比。暴露因素与疾病的关系 在进行暴露因素和疾病关系的研究时,暴露和疾病的关系可以总结为下列四格表:

暴露无暴露
病例ab
对照cd

比值比OR=ad/bc。

①在病例对照研究中,比值比指病例组中暴露与非暴露人数的比值(a/b)和对照组中暴露与非暴露人数的比值 c/d) 的比,得出OR=ad/bc,所以又叫交叉乘积比,该值用作相对危险度的估计值。

②在队列研究中,指的是暴露组中患病与非患病者的比值(a/c)和非暴露组中患病与非患病者的比值(b/d)的比。也用作相对危险度的估计值。在队列研究中,可以计算相对危险度,所以一般不计算比值比,但是有的时候根据需要也应用OR作为联系的强度的指标,例如在应用logistic回归模型对队列研究的资料进行多因素分析时,即应用OR值。

本文将介绍OR的批量计算以及结果可视化,以下代码经供参考,并不是最好的实现方式,欢迎大家留言讨论,一起学习。

aa<- read.csv('I:/COR-NGS-ALL.csv')
C:\Users\JINWAN~1\AppData\Local\Temp\1646299808(1).png
glm1<- glm(formula = ASXL1 ~ TET2,
           data=aa,
           family = binomial)
glm2<- summary(glm1);glm2
#计算OR并保留两位数
OR<-round(exp(coef(glm1)),2)

将上述计算包装成函数进行批量分析:

Uni_glm_model<- 
  function(x,y,aa){
    #拟合结局和变量
    FML<-as.formula(paste0(x,"~",y))
    #glm()逻辑回归
    glm1<-glm(FML,data=aa,family = binomial)
    #提取所有回归结果放入glm2中
    glm2<-summary(glm1)
    #1-计算OR值保留两位小数
    OR<-round(exp(coef(glm1)),3)
    #2-提取SE
    SE<-glm2$coefficients[,3]
    #3-计算CI保留两位小数并合并
    CI5<-round(exp(coef(glm1)-1.96*SE),3)
    CI95<-round(exp(coef(glm1)+1.96*SE),3)
    CI<-paste0(CI5,'-',CI95)
    #4-提取P值
    P<-round(glm2$coefficients[,4],3)
    #5-将变量名、OR、CI、P合并为一个表,删去第一行
    Uni_glm_model <- data.frame('Characteristics'=x,
                                'OR' = OR,
                                'CI' = CI,
                                'P' = P)[-1,]
    #返回循环函数继续上述操作                     
    return(Uni_glm_model)
  }  

for循环批量计算基因两两OR值:

#这里我选择了本数据的第2-14列的变量
#把它们放入variable.names中
variable.names<- colnames(aa)[c(2:14)];variable.names 
#变量带入循环函数
results_or<-data.frame(matrix(NA,13,13))
colnames(results_or)<-variable.names
rownames(results_or)<-variable.names
results_p<-results_or
for (i in (1:(length(variable.names)-1))) {
  for (j in (i+1):(length(variable.names))) {
    results_or[i,j]<-Uni_glm_model(variable.names[i],variable.names[j],aa)$OR
    results_p[i,j]<-Uni_glm_model(variable.names[i],variable.names[j],aa)$P
  }
}
results_or<-results_or[-13,-1]
results_p<-results_p[-13,-1]

results_or$x<-rownames(results_or)
results_p$x<-rownames(results_p)

ggplot2绘制热图,其它函数都可以:

library(dplyr)
library(ggplot2)
library(grDevices)
results_or %>% 
  reshape2::melt(id.vars="x",variable.name="y") %>% 
  na.omit() -> dftmp
results_p %>% 
  reshape2::melt(id.vars="x",variable.name="y") %>% 
  na.omit() -> dftmp_p
dftmp$p<-dftmp_p$value
dftmp$plab<-ifelse( dftmp_p$value<0.001,"***",ifelse(dftmp_p$value<0.01,"**",ifelse(dftmp_p$value<0.05,"*","")))

dftmp$x<-factor(dftmp$x,
                levels = variable.names)
dftmp$y<-factor(dftmp$y,
                levels = variable.names)

cols<-c(colorRampPalette(c("#6BAED6","#FFFFFF"))(10),
        colorRampPalette(c("#FFFFFF","#FB6A4A"),bias=2)(90))

ggplot(data=dftmp,aes(x,y,fill=value))+
  geom_tile(colour="grey80", 
            #linewidth=2, 
            width=.9, 
            height=.9)+
  theme_minimal()+
  theme(panel.grid = element_blank(),axis.text.x = element_text(angle=90, hjust=-0.02))+
  scale_x_discrete(position = "top")+
  scale_y_discrete(position = "left")+
  labs(x=NULL,y=NULL)+scale_fill_gradientn(colours= cols, limits=c(0, 10),
                                           breaks=c(0,1,10), 
                                          na.value=rgb(246, 246, 246, max=255),
                                          labels=c("0", "1", "10"),
                                          guide=guide_colourbar(ticks=T,barwidth=10))+
  theme(legend.position = c(0.7, 0.2),legend.direction = "horizontal")+
  geom_text(aes(label = plab),col='black',cex=3)

发表评论

邮箱地址不会被公开。 必填项已用*标注