bioinfoer»论坛 › 生信基础 › 可视化 › 线性回归画出双坐标轴的图

0 收听	12 听众	390 主题

线性回归画出双坐标轴的图

发表于 6 天前 | 查看: 34| 回复: 1

背景

线性回归，画出双坐标轴的图

出自文章PMID: 30027617

应用场景

同时计算和展示3组数据的相关性，A vs. B 和 A vs. C。

例如例文中同时展示B7-H3表达跟mir-29a的相关性，以及跟甲基化的相关性。由此看出：In IDH‐wild grade II gliomas, B7‐H3 expression seemed to be downregulated by both methylation on cg10586317 and microRNA‐29a, rationalizing the paradox in grade II gliomas with wild‐type IDH.

环境设置

使用国内镜像安装包

options("repos"= c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
install.packages("ggplot2")

加载包

library(ggplot2)
library(gtable) #方法二要用到
library(grid) #方法二要用到

Sys.setenv(LANGUAGE = "en") #显示英文报错信息
options(stringsAsFactors = FALSE) #禁止chr转成factor

参数设置

## 两个坐标轴的颜色
s1Colour="#6B9EC2"
s2Colour="#EE861A"

## xlab
xtitle="mRNA expression"
## ylab
ylefttitle="miRNA"
yrighttitle="cg10586317"

##Aspect ration of picture，图的宽高比例
AspectRatio = 5/2

输入文件

包含三列：对应X轴和两个Y轴

分别对应例文中的X轴：mRNA表达量，可以是FPKM等值。

左侧Y轴：miRNA的表达量

右侧Y轴：DNA甲基化水平

inputData <- read.table("easy_input.txt", header = T)
head(inputData)
#inputData$Y2 <- inputData$Y1 * 0.00001
#All secondary axes must be based on a one-to-one transformation of the primary axes.
#默认右侧Y轴跟左侧Y轴范围是一样的，然而示例数据的Y1和Y2范围差距很大，因此需要用Y1和Y2的数值计算两个Y轴的比例
(scaleFactor <- max(inputData$Y1) / max(inputData$Y2)) 
#在画右侧Y轴、X和Y2的散点和拟合线时都要用这个比例来换算
#除了用乘法，你还可以用加法

#前面参数设置了图的宽高比例，保存在AspectRatio里
#用Y1和Y2的数值换算坐标轴的宽高比，在后面的画图代码coord_fixed中要用到
(ratioValues <- ((max(inputData$X)-min(inputData$X))/(max(inputData$Y1)-min(inputData$Y1))) / AspectRatio)

计算相关系数

#写一个计算相关系数的函数，同时写成图中的格式
correlation <- function(data,X,Y){
  #计算相关系数
  R <- cor(X,Y,method = "pearson") #或"pearson", "kendall"
  
  #把相关系数写成图中展示的格式
  eq <- substitute(~~italic(R)~"="~r, list(r = format(round(R, 4),nsmall= 4)))
  Req <- as.character(as.expression(eq))
  return(Req)
}

#用上面的函数分别计算mRNA跟miRNA、DNA甲基化之间的相关系数
(R1cor <- correlation(inputData,inputData$X,inputData$Y1))
(R2cor <- correlation(inputData,inputData$X,inputData$Y2))

开始画图

此处提供两种方法画图：

方法一，优点：简单，易于理解；缺点：右侧Y轴的label的文字方向与原文相反。
方法二，优点：完全复现原文的图；缺点：用到gtable、grid，小伙伴可能不熟悉。

方法一

p <- ggplot(inputData, aes(x=inputData$X )) + 
  geom_point(aes(y=inputData$Y1), col=s1Colour) + #画X和Y1的散点
  geom_smooth(aes(x=inputData$X, y=inputData$Y1), 
              method="lm", se=F, col=s1Colour) + #添加线性拟合线
  
  geom_point(aes(y=(inputData$Y2) * scaleFactor), col=s2Colour) + #画X和Y2的散点
  geom_smooth(aes(x=inputData$X, y=inputData$Y2 * scaleFactor), 
              method="lm", se=F, colour=s2Colour) + #添加线性拟合线
  
  scale_y_continuous(name = ylefttitle, 
                     #breaks = c(5000, 10000, 15000, 20000, 25000), #左侧Y轴刻度
                     sec.axis=sec_axis(trans = ~./scaleFactor, name = yrighttitle)) + #画右侧Y轴
  coord_fixed(ratio = ratioValues) + #坐标轴的宽高比
  
  labs(x=xtitle) + theme_bw() + #去除背景色
  theme(panel.grid = element_blank(), #去除网格线
        axis.title.x = element_text(size=15), #X轴题目字号
        axis.text.x = element_text(size=12), #X轴刻度字号
        axis.title.y.right = element_text(size=15, color = s2Colour), #右侧Y轴刻度字体的颜色
        axis.text.y.right = element_text(size=12, color = s2Colour), #右侧Y轴title颜色
        axis.title.y.left = element_text(size=15, colour = s1Colour), #左侧Y轴刻度字体的颜色
        axis.text.y.left = element_text(size=12, color = s1Colour)) #左侧Y轴title颜色
p

把前面计算出的相关系数R1cor和R2cor标在图上

# 相关系数在图上的位置，需根据你自己的数据灵活调整
x <- max(inputData$X)*(0.95)
y <- max(inputData$Y1)*(0.97)
y2 <- max(inputData$Y1)*(0.9)

# 标相关系数
p1 <- p + annotate("text", x=x, y=y, label=R1cor, size=5, parse = TRUE, col=s1Colour ) +
  annotate("text", x=x, y=y2, label=R2cor, size=5, parse = TRUE, col=s2Colour ) 
p1

# 保存到文件
ggsave("TwoPlotOfLinearRegression.pdf", width = 6, height = 3)

右侧y轴的title头朝右，原文头朝左。我们用方法二来实现对原文完全复现。

方法二

先用ggplot2拿到两套坐标的plot object，修改之后重新画图。

用ggplot拿到两套坐标的plot object

#X和Y1的散点、拟合线、相关系数文字
p1 <- ggplot(inputData, aes(x=inputData$X, y=inputData$Y1)) +
  geom_point(col=s1Colour) +
  geom_smooth(aes(x=inputData$X, y=inputData$Y1), method="lm", se=F, col=s1Colour)+
  annotate("text", x=x, y=y , label=R1cor,size=5, parse = TRUE, col=s1Colour ) +
  annotate("text", x=x, y=y2, label=R2cor,size=5, parse = TRUE, col=s2Colour )+
  
  labs(x=xtitle, y=ylefttitle) + theme_bw()+
  theme(panel.grid = element_blank(), #去除网格线
        axis.title.x = element_text(size=15),
        axis.text.x = element_text(size=12),
        axis.title.y = element_text(size=15, colour = s1Colour),
        axis.text.y = element_text(size=12, color = s1Colour))+
  coord_fixed(ratio = ratioValues)

#X和Y2的散点、拟合线、背景透明
p2 <- ggplot(inputData, aes(x=inputData$X, y=inputData$Y2)) +
  geom_point(col=s2Colour) + 
  geom_smooth(aes(x=inputData$X, y=inputData$Y2), method="lm", se=F, colour=s2Colour)+

  labs(y=yrighttitle) + theme_bw() +
  theme(panel.background = element_rect(fill = NA), 
        panel.grid = element_blank(), #去除网格线
        axis.title.y = element_text(size=15, color = s2Colour),
        axis.text.y = element_text(size=12, color = s2Colour)) +
  coord_fixed(ratio = ratioValues)

#提取plot object，保存到gtable数据结构里
g1 <- ggplot_gtable(ggplot_build(p1))
g2 <- ggplot_gtable(ggplot_build(p2))

#查看一下提取出的信息长啥样
head(g1$layout)
head(g2$layout)

重新画图

**提示：**选择从grid.newpage()到任意一个grid.draw(g)之间的部分一起运行，就会看到这部分代码的画图效果。

#pdf("TwoPlotOfLinearRegression_plus.pdf", width = 6, height = 3)
grid.newpage()

#画X和Y1的散点、拟合线、相关系数文字
pp <- c(subset(g1$layout, name == "panel", se = t:r))
g <- gtable_add_grob(g1, g2$grobs[[which(g2$layout$name == "panel")]], pp$t, 
                     pp$l, pp$b, pp$l)
#head(g)
grid.draw(g)

#画X和Y2的散点、拟合线
ia <- which(g2$layout$name == "axis-l")
ga <- g2$grobs[[ia]]
ax <- ga$children[[2]]

ax$widths <- rev(ax$widths)
ax$grobs <- rev(ax$grobs)

g <- gtable_add_cols(g, g2$widths[g2$layout[ia, ]$l], length(g$widths) - 1)

grid.draw(g)

#画右侧Y轴刻度标签
g <- gtable_add_grob(g, ax, pp$t, length(g$widths) - 1, pp$b)
grid.draw(g)

#画右侧Y轴title
yLabTitle <- which(g2$layout$name == "ylab-l")
gaYLT <- g2$grobs[[yLabTitle]] 
axYLT <- gaYLT$children[[1]]

axYLT$widths <- rev(axYLT$widths)
axYLT$grobs <- rev(axYLT$grobs)

g <- gtable_add_cols(g, g2$widths[g2$layout[yLabTitle, ]$l], length(g$widths) - 1)
g <- gtable_add_grob(g, axYLT, pp$t, length(g$widths) - 1, pp$b)

grid.draw(g)

#dev.off()