R语言批量计算OR（优势比）并可视化—科研工具箱

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

比值比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\')

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)