自学R语言（十六）-相关性分析和具体基因也组别之间的相关性热图的绘制

迷途中小小书童 2019-04-24

展开全文

这幅图的缘起至少一个月以前了，在easychart的群里反反复复被提及过很多次，一直想尝试解决这个图，一直各种拖延。上周三左右开始动手，周天出图，持续了四天，如此之低效率主要不是在画这幅图上，而是我用了最傻的box()函数去测试corrplot包的相关系数图的绘图区域，当发现没有框住全部绘图区域的时候，自以为是的认为原作者定义绘图区域的时候用了什么黑科技，导致我不能准确定位图中的每个点的坐标。这时，我没有仔细的去阅读源代码，而是选择了自己造个类似的“假轮子”，然后再用这个轮子添加图中的其它元素。“假轮子”造了两天，基本重现了corrplot包中的相关系数热图，问题也解决了，一个偶然的机会想起问下原包作者魏大佬怎么定义坐标的，他让我用text(1:n, 1:n, 1:n)测试下，这才恍然大悟我前面的方法多么愚蠢。

原图结构

画比较复杂的图之前，要尝试图分解成比较简单的块，然后分别绘制。从原图我们很容易发现，主要有三部分：右上角是类似于corrplot包中的上三角相关系数图；下三角是一组点之间的连接线（作者用了弧线，直线也能达到同样的效果）；剩余部分主要是图例等其它辅助绘图元素。接下来将分别讨论如何绘制这些元素。

相关系数图

原图的上三角部分主要是调用corrplot包中的corrplot()函数进行绘制，原包是基于基础绘图系统开发，能满足多样需求，详细了解该包用法可以查看vignette("corrplot-intro")。画这幅图还有个棘手的问题是没有原始数据，我自己也不懂相关模型，只能根据图的样子用R自带数据集mtcars进行测试。

第一步非常简单，method = 'square'表示使用正方形符号或者"circle"表示使用圆形的符号，type = 'upper'表示只画上三角区域， "lower"表示显示右下角部分，"full"表示显示全部的相关性热点图。除此之外的其它参数几乎不用管

library(corrplot)par(omi = c(0.3, 0.3, 0.3, 0.3),    cex = 1.2,    family = 'Arial') # windows系统可能需要安装其他字体包M <- cor(mtcars) #以mtcars这个已有的数据为例，计算相关系数矩阵

corrplot(M, method = "pie", type = 'upper') #使用corrplot包绘制相关系数热点图，符号使用圆形的，只显示右上部分

mpg cyl disp hp drat wt

mpg 1.0000000 -0.8521620 -0.8475514 -0.7761684 0.68117191 -0.8676594

cyl -0.8521620 1.0000000 0.9020329 0.8324475 -0.69993811 0.7824958

disp -0.8475514 0.9020329 1.0000000 0.7909486 -0.71021393 0.8879799

hp -0.7761684 0.8324475 0.7909486 1.0000000 -0.44875912 0.6587479

drat 0.6811719 -0.6999381 -0.7102139 -0.4487591 1.00000000 -0.7124406

wt -0.8676594 0.7824958 0.8879799 0.6587479 -0.71244065 1.0000000

qsec 0.4186840 -0.5912421 -0.4336979 -0.7082234 0.09120476 -0.1747159

vs 0.6640389 -0.8108118 -0.7104159 -0.7230967 0.44027846 -0.5549157

am 0.5998324 -0.5226070 -0.5912270 -0.2432043 0.71271113 -0.6924953

gear 0.4802848 -0.4926866 -0.5555692 -0.1257043 0.69961013 -0.5832870

carb -0.5509251 0.5269883 0.3949769 0.7498125 -0.09078980 0.4276059

qsec vs am gear carb

mpg 0.41868403 0.6640389 0.59983243 0.4802848 -0.55092507

cyl -0.59124207 -0.8108118 -0.52260705 -0.4926866 0.52698829

disp -0.43369788 -0.7104159 -0.59122704 -0.5555692 0.39497686

hp -0.70822339 -0.7230967 -0.24320426 -0.1257043 0.74981247

drat 0.09120476 0.4402785 0.71271113 0.6996101 -0.09078980

wt -0.17471588 -0.5549157 -0.69249526 -0.5832870 0.42760594

qsec 1.00000000 0.7445354 -0.22986086 -0.2126822 -0.65624923

vs 0.74453544 1.0000000 0.16834512 0.2060233 -0.56960714

am -0.22986086 0.1683451 1.00000000 0.7940588 0.05753435

gear -0.21268223 0.2060233 0.79405876 1.0000000 0.27407284

carb -0.65624923 -0.5696071 0.05753435 0.2740728 1.00000000

连接线

没有原始数据，我随便模拟生成了三组，分别是“Group01”、“Group02”和“Group03”，因为每组都要和每个变量连线，所以线条的数量是组数和相关系数矩阵行数的乘机（这里是3*11 = 33个）。

R语言本质是向量化的，基础绘图函数也基本支持使用向量作为位置、点的形状、大小、颜色、线条类型、宽度颜色等的参数值。所以我们数据处理部分将相关参数的值统一整理到数据框中，方便后面调用。

library(dplyr)library(corrplot)# 准备数据，这一部分的数据在实际应用中应该是自己计算各个组别与具体某几个基因的相关性包括p值和r值去表征相关性
set.seed(20190420)n <- ncol(mtcars)grp <- c('Group01', 'Group02', 'Group03') # 分组名称sp <- c(rep(0.0008, 6), rep(0.007, 2), rep(0.03, 3), rep(0.13, 22)) # P值gx <- c(-4.5, -2.5, 1) # 分组的X坐标gy <- c(n-1, n-5, 2.5) # 分组的Y坐标df <- data.frame(  grp = rep(grp, each = n), # 分组名称，每个重复n次  gx = rep(gx, each = n), # 组X坐标，每个重复n次  gy = rep(gy, each = n), # 组Y坐标，每个重复n次  x = rep(0:(n - 1) - 0.5, 3), # 变量连接点X坐标  y = rep(n:1, 3), # 变量连接点Y坐标  p = sample(sp), # 对人工生成p值进行随机抽样  r = sample(c(rep(0.8, 4), rep(0.31, 7), rep(0.12, 22)))   # 对人工生成r值进行随机抽样) #以上数据中的P值和r值本质是根据具体的实际实验当中的数据进行计算得到的具有意义的数值，所以在以后的具体实战中，只需要修改p值和r即可，gx，gy，x，y代表的是连接线的线段的初始点和终点不用调整

# 这一部分代码是按照原图图例说明处理线条宽度和颜色映射df <- df %>% #此处使用的是dplyr包中的管道符将df传递给下一个分析函数作为输入值。
  mutate( #此处使用的plyr包的mutate函数在原有数据框的基础上，对已有数据进行汇总，并且添加至新的列，以便于下游绘制不同p值和不同r值得线条颜色和线条粗细使用。    lcol = ifelse(p <= 0.001, '#1B9E77', NA), #使用ifelse函数根据p值和r值对线条的粗细和颜色进行因子化处理，并且添加至新的列便于下游分析使用。    # p值小于0.001时，颜色为绿色，下面依次类推    lcol = ifelse(p > 0.001 & p <= 0.01, '#88419D', lcol),    lcol = ifelse(p > 0.01 & p <= 0.05, '#A6D854', lcol),    lcol = ifelse(p > 0.05, '#B3B3B3', lcol),    lwd = ifelse(r >= 0.5, 14, NA),    # r >= 0.5 时，线性宽度为14，下面依次类推    lwd = ifelse(r >= 0.25 & r < 0.5, 7, lwd),    lwd = ifelse(r < 0.25, 1, lwd)  )

经过汇总，和处理之后就会得到如下的数据框，方便我们接下来绘制连接线

X grp gx gy x y p r lcol lwd

1 Group01 -4.5 10.0 -1 11 0.0070 0.12 #88419D 2

2 Group01 -4.5 10.0 0 10 0.0008 0.12 #1B9E77 2

3 Group01 -4.5 10.0 1 9 0.0008 0.12 #1B9E77 2

4 Group01 -4.5 10.0 2 8 0.0008 0.12 #1B9E77 2

5 Group01 -4.5 10.0 3 7 0.1300 0.12 #B3B3B3 2

6 Group01 -4.5 10.0 4 6 0.0240 0.80 #A6D854 11

7 Group01 -4.5 10.0 5 5 0.1300 0.12 #B3B3B3 2

8 Group01 -4.5 10.0 6 4 0.0008 0.12 #1B9E77 2

9 Group01 -4.5 10.0 7 3 0.1300 0.31 #B3B3B3 6

10 Group01 -4.5 10.0 8 2 0.0260 0.12 #A6D854 2

11 Group01 -4.5 10.0 9 1 0.3400 0.31 #B3B3B3 6

12 Group02 -2.5 6.0 -1 11 0.0360 0.12 #A6D854 2

13 Group02 -2.5 6.0 0 10 0.0580 0.12 #B3B3B3 2

14 Group02 -2.5 6.0 1 9 0.1300 0.12 #B3B3B3 2

15 Group02 -2.5 6.0 2 8 0.0070 0.12 #88419D 2

16 Group02 -2.5 6.0 3 7 0.1300 0.12 #B3B3B3 2

17 Group02 -2.5 6.0 4 6 0.1300 0.80 #B3B3B3 11

18 Group02 -2.5 6.0 5 5 0.1300 0.80 #B3B3B3 11

19 Group02 -2.5 6.0 6 4 0.0008 0.80 #1B9E77 11

20 Group02 -2.5 6.0 7 3 0.1300 0.12 #B3B3B3 2

21 Group02 -2.5 6.0 8 2 0.0240 0.12 #A6D854 2

22 Group02 -2.5 6.0 9 1 0.0150 0.31 #A6D854 6

23 Group03 1.0 2.5 -1 11 0.1300 0.31 #B3B3B3 6

24 Group03 1.0 2.5 0 10 0.1300 0.31 #B3B3B3 6

25 Group03 1.0 2.5 1 9 0.0300 0.12 #A6D854 2

26 Group03 1.0 2.5 2 8 0.0300 0.12 #A6D854 2

27 Group03 1.0 2.5 3 7 0.0008 0.31 #1B9E77 6

28 Group03 1.0 2.5 4 6 0.1300 0.12 #B3B3B3 2

29 Group03 1.0 2.5 5 5 0.1300 0.12 #B3B3B3 2

30 Group03 1.0 2.5 6 4 0.0025 0.31 #88419D 6

31 Group03 1.0 2.5 7 3 0.0260 0.12 #A6D854 2

32 Group03 1.0 2.5 8 2 0.1300 0.12 #B3B3B3 2

33 Group03 1.0 2.5 9 1 0.0140 0.12 #A6D854 2

可以发现，把每个图形元素及其属性参数整理成一个数据框，画图过程简单很多。很多时候我们觉得基础绘图系统很复杂，一个简单的图可能需要很长的代码才能解决，其实也和我们没有很好的利用R向量化运算的特点，没有去寻找最简洁的方案有关系。

segments(df$gx, df$gy, df$x, df$y, lty = 'solid', lwd = df$lwd,          col = df$lcol, xpd = TRUE) # 使用segments函数绘制连接线，df$gx, df$gy表示线段的起点，df$x, df$y表示线段的终点。points(df$gx, df$gy, pch = 24, col = 'blue', bg = 'blue', cex = 3, xpd = TRUE) # 组标记点，绘制每个组的标记点。df$gx, df$
gy表示组标记点的坐标，pch表示标记点的符号类型，24表示三角形，text(df$gx - 0.5, df$gy, labels = df$grp, adj = c(1, 0.5), cex = 1.5, xpd = TRUE)# 表示每个组的名称，首先限定名称的显示位置，然后就是显示内容。

添加连接线之后的图片

添加标记点之后的图片

这一部分主要在前面基础图的基础上确定每个元素标记位置，出图之后根据细节进行微调，没有太多复杂的内容。

labels01 <- c('<= 0.001','0.001 < x <= 0.01','0.01 < x <= 0.05','> 0.05')labels02 <- c('>= 0.5', '0.25 - 0.5', '< 0.25')labels_x <- rep(-6, 4)labels_y <- seq(4.6, 2.6, length.out = 4)text(-6.5, 5.2, 'P-value', adj = c(0, 0.5), cex = 1.2, font = 2, xpd = TRUE)text(labels_x, labels_y, labels01, adj = c(0, 0.5), cex = 1.2, xpd = TRUE)points(labels_x - 0.5, labels_y, pch = 20, col = c('#1B9E77', '#88419D','#A6D854', '#B3B3B3'),       cex = 3, xpd = TRUE)lines_x <- c(-6.5, -3, 0.5)lines_y <- rep(1.2, 3)text(-6.5, 1.9, "Mantel's r", adj = c(0, 0.5), cex = 1.2, font = 2, xpd = TRUE)text(lines_x + 1.5, lines_y, labels02, adj = c(0, 0.5), cex = 1.2, xpd = TRUE)segments(lines_x, lines_y, lines_x + 1, lines_y, lwd = c(14, 7, 2.5), lty = 'solid',          col = '#B3B3B3', xpd = TRUE)
## 图例框框，这一部分就是绘制图注信息外面框，本质就是一条线段一条线段的拼接起来的，严格按照坐标信息标记即可。segments(-6.9, 5.6, -2.8, 5.6, lty = 'solid', lwd = 1.2,          col = 'grey50', xpd = TRUE)segments(-2.8, 5.6, -2.8, 1.8, lty = 'solid', lwd = 1.2,          col = 'grey50', xpd = TRUE)segments(-2.8, 1.8, 3.6, 1.8, lty = 'solid', lwd = 1.2,          col = 'grey50', xpd = TRUE)segments(3.6, 1.8, 3.6, 0.7, lty = 'solid', lwd = 1.2,          col = 'grey50', xpd = TRUE)segments(3.6, 0.7, -6.9, 0.7, lty = 'solid', lwd = 1.2,          col = 'grey50', xpd = TRUE)segments(-6.9, 0.7, -6.9, 5.6, lty = 'solid', lwd = 1.2,          col = 'grey50', xpd = TRUE)