原文链接:http://tecdat.cn/?p=23891
能够应用环状图形展现基因数据比拟。能够增加多种图展信息,如热图、散点图等。
本文指标:
- 可视化基因组数据
制作环形热图
环形热图很漂亮。能够通过 R 来实现环形热图。
首先,让咱们生成一个随机矩阵,并将其随机分成五组。
mat1 = rbind(cbind(matrix(rnorm(50*5, mean = 1), nr = 50),
matrix(rnorm(50*5, mean = -1), nr = 50)),
cbind(matrix(rnorm(50*5, mean = -1), nr = 50),
matrix(rnorm(50*5, mean = 1), nr = 50))
)
上面的图是热图的失常布局。
Heatmap(mat1, row_split = split)
在接下来的章节中,我将演示如何将其可视化。
输出数据
heatmap()的输出应该是一个矩阵(或者一个将被转换为单列矩阵的向量)。如果矩阵被宰割成组,必须用 split 参数指定一个分类变量。留神 spilt 的值应该是一个字符向量或一个因子。如果它是一个数字向量,它将被转换为字符。
色彩是矩阵中数值的重要美学映射。用户必须用用户定义的色彩模式指定 col 参数。如果矩阵是间断数字,如果矩阵是字符,col 的值应该是一个命名的色彩向量。
上面的图是之前热图的圆形版本。请留神,矩阵的行沿圆形方向散布,矩阵的列沿径向方向散布。在上面的图中,圆形被分成五个局部,每个局部对应一个行组。
heatmap(mat1col_fun1)
有一件事十分重要,那就是在创立圆形热图之后,你必须齐全删除布局。
如果没有指定 split,就只有一个大的扇区蕴含残缺的热图。
环形布局
与生成的其余圆形图相似,环形布局能够在制作图之前由 par()管制。
热图轨道的参数能够在 circos()函数中管制,如 track.height(轨道的高度)和 bg.border(轨道的边界)。
在上面的例子中,通过设置 show.sector.labels 参数,减少了扇区的标签。扇区的程序是 c(“a”, “b”, “c”, “d”, “e”),按时钟方向排列。你能够在上面的图中看到,a 扇区从 \(\theta = 90^{\circ}\)开始。
heatmap(bg.border)
如果 split 参数的值是一个因子,那么因子程度的顺序控制热图的程序。如果 split 是一个简略的向量,热图的程序是 unique(split)。
# 留神,因为在前一个图中调用了 circos.clear()。# 当初布局从 theta = 0 开始(第一个扇区是 'e')。heatmap(levels = c("e", "d", "c", "b", "a))
树状图和行名
默认状况下,数字矩阵是按行聚类的,因而,有聚类产生的树状图。side 参数管制树状图绝对于热图轨道的地位。留神,树枝图是在一个拆散的轨道上。
heatmap(dend.side = "inside")
树状图的高度是由 dend.track.height 参数管制的。
矩阵的行名能够通过设置 rownames.side 参数来绘制。行名也会被绘制在一个拆散的轨道中。
heatmap(rownames.side = "inside")
矩阵的行名和树状图能够同时绘制。当然,它们不能在热图轨道的同一侧。
dend.side = "inside",
rownames.side = "outside"
行名的图形参数能够设置为标量或向量,长度与矩阵中的行数雷同。
heatmap(col = col_fun1, rownames.side = "outside")
树状图的图形参数能够通过回调函数间接渲染树状图来设置,这一点将在前面演示。
聚类
默认状况下,数字矩阵是按行聚类的。cluster 参数能够设置为 FALSE 来敞开聚类。
当然,当 cluster 被设置为 FALSE 时,即便 dend.side 被设置,也不会绘制树状图。
聚类办法和间隔办法由 clustering.method 和 distance.method 参数管制。
请留神 heatmap()不间接反对对矩阵列的聚类。你应该在应用 heatmap()之前利用列的从新排序,例如。
hclust(dist(t(mat1)))$order
对树状图的回调
聚类产生树状图。回调函数能够在每个树状图生成后利用于相应的类。回调函数能够编辑树状图,例如:1. 重新排列树状图,或者 2. 给树状图着色。
在 circos.heatmap()中,一个用户定义的函数应该被设置为 callback 参数。该用户定义的函数应该有三个参数。
dend
: 以后扇区的树状图。m
: 与以后扇区绝对应的子矩阵。si
: 以后扇区的扇区索引(或扇区名称)。
默认的回调函数定义如下,它通过对矩阵行均值加权来重新排列树状图。
reorder(dend, rowMeans(m))
上面的例子通过 dendsort()对每个扇区的树状图从新排序。
heatmap( col = col_fun1, dend.side = "inside",
dendsort(dend)
}
咱们能够应用 color()来渲染树状图的边缘。例如,为五个区的树枝图调配不同的色彩。这里,树枝图轨道的高度由 height 参数减少。
den = function(dend, m, si) {
# 当 k = 1 时,它为整个树状图渲染一种雷同的色彩
color\_branches(dend, k = 1, col = dend\_col\[si\])
或者如果矩阵没有被宰割,咱们能够给子树状图调配不同的色彩。
color_branches(dend, k = 4, col = 2:5)
多个热图轨迹
如果你制作的环状图只蕴含一个热图轨迹,应用 heatmap()是非常简单的。如果你制作一个蕴含多个轨道的更简单的环状图,你应该理解对于 heatmap()的更多细节。
heatmap()的第一次调用实际上是初始化布局,即利用聚类和拆分矩阵。树状图和宰割变量是外部存储的。这就是为什么你应该明确地调用 clear()来删除所有的外部变量,这样能够确保当你制作一个新的圆形热图时,heatmap()的第一次调用是在一个新的环境中。
heatmap()的第一次调用决定了所有轨道的行程序(循环方向的程序),因而,接下来的轨道中的矩阵共享与第一个轨道中雷同的行程序。另外,前面轨道中的矩阵也会依据第一个 heatmap 轨道中的宰割状况进行宰割。
如果在第一个热图轨道中没有利用聚类,则应用行的天然排序(即 c(1,2,…,n))。
mat1\[sample(100, 100), \] # 按行随机排列 mat1
heatmap(mat1, split, col_fun1, dend.side = "outside")
如果我切换两个轨道,你能够看到当初的聚类是由第一个热图轨道控制的,也就是蓝 - 红热图轨道。
heatmap(mat1, col = col_fun2)
你可能想问,如果我不心愿聚类是由第一个轨道决定的,而第二个或第三个轨道呢?解决办法很简略。实际上,初始化能够通过明确调用 initialize()函数来手动实现。
在 initialize()中,你指定你想利用聚类的任何矩阵以及宰割变量,而后,上面的 heatmap()调用都共享这个布局。
在上面的例子中,全局布局是由 mat1 决定的,它在第二个轨道中被可视化。我在第一个轨道中设置了 side = “outside”,实际上你能够发现树状图实际上是依据第二个轨道中的矩阵生成的。
circos.heatmap.initialize(mat1, split = split)
在下一个例子中,热图布局是由 mat1 生成的,而两个热图轨道别离只蕴含五列。
initialize(mat1, split = split)
与其余轨道 整合
其余非热图轨道整合。在环形布局中,x 轴和 y 轴上的值只是数字索引。假如在一个扇形区域内有 nr 行和 nc 列的热图,热图行的绘制距离为(0,1),c(1,2),...,c(nr-1,nr),热图列也相似。同时,原始矩阵也被从新排序。如果减少更多的轨道,须要思考所有这些影响,以确保与热图轨道有正确的对应关系。
热图布局实现后,轨道 / 扇区 / 单元的额定信息能够通过非凡变量 CELL_META 来检索。单元 / 扇区的附加元数据列举如下,它们对于正确对应热图轨道十分重要。
CELL_META$row_dend 或简称 CELL_META$dend:以后扇区的树状图。如果没有进行聚类,则该值为 NULL。
CELL_META$row_order 或简称 CELL_META$order:聚类后以后扇区中子矩阵的行排序。如果没有进行聚类,其值为 c(1, 2, ...,)。
CELL_META$subset。原始残缺矩阵中指数的子集。这些值的排序是递增的。
以下是 CELL\_META$row\_dend、CELL\_META$row\_order 和 CELL_META$subset 在循环热图例子中的第一个扇区的输入。
row_order
subset
在上面的例子中,我增加了一个轨迹,它将 mat1 中前五列的行平均值可视化。我增加了 cell.padding = c(0.02, 0, 0.02, 0),这样最大和最小的点就不会与单元格的高低边界重叠了。
track(ylim = range(row_mean{lines(CELL_META$cell.xlim, c(0, 0),
points(seq_along(y) - 0.5, y, )
同样地,如果把点数轨道作为第一条轨道,则应当时对布局进行初始化。
initialize(mat1, split = split)
# 这与后面的例子雷同
track(ylim = range(row_mean), panel.fun = function(x, y) {circose(y > 0, "red", "blue"))
}, cell.padding = c(0.02, 0, 0.02, 0))
# 不须要在这里指定 "宰割"。
Boxplots 被用来对应矩阵的行。
circos.heatmap(mat1, split = split, col = col_fun1)
circos.track(ylim = range(mat1), panel.fun = function(x, y) {m = mat1\[CELL_META$subset, 1:5, drop = FALSE\]
m = m\[CELL\_META$row\_order, , drop = FALSE\]
n = nrow(m)
# circos.boxplot is applied on matrix columns, so here we transpose it.
circos.boxplot(t(m), pos = 1:n - 0.5, pch = 16, cex = 0.3)
circos.lines(CELL_META$cell.xlim, c(0, 0), lty = 2, col = "grey")
}, cell.padding = c(0.02, 0, 0.02, 0))
增加正文
能够通过设置 abels = TRUE 来增加区块的标签,然而,这并不提供对标签的任何定制。用户能够通过定义 fun 函数来定制本人的标签,演示如下。
heatmap(col = col_fun1)
circos.track
adj = c(0.5, 0), niceFacing = TRUE)
heatmap()不间接反对矩阵的列名,但也能够通过自行定义 fun 函数轻松增加。在上面的例子中,我通过 par()中的 after 参数在最初一个扇区(第五扇区)后设置了较大的空间(10 度,用户通常须要尝试几个值来获得最佳空间),之后我在 fun 中绘制了最初一个扇区中的列名。
par(gap.after = c(2, 2, 2, 2, 10))
heatmap(mat1, split , col = col_fun1)
track(track.index = 1}, bg.border = NA)
下一个例子增加了矩形和标签来显示矩阵中的两组列。fun 外面的代码很简略。convert_x()将 x 方向上的单位转换为环形坐标系中测量的适当数值。
par(gap.after = c(2, 2, 2, 2, 10))
heatmap(mat1, split , col = col_fun1)
track(track.index = 1
circos.text(cell.xlim\[2\] + convert_x(3, "mm"), 7.5,
}, bg.border = NA)
circlize 不生成图例,然而图例能够由 Legend()函数手动生成并增加到圆形图中。上面是一个增加图例的简略例子。在下一节中,你能够找到一个增加许多图例的更简单的例子。
heatmap(mat1, split = split)
clear()
grid.draw(lgd)
一个简单的圆形热图的例子
在本节中,我将演示如何制作简单的圆形热图。下图是失常布局的热图,当初我将用圆形布局扭转它们。
热图直观地显示了 DNA 甲基化、基因表白和其余基因组程度信息之间的相关性。
原始热图是用随机数据集生成的。
与原始热图相似,通过对甲基化矩阵(mat_meth)的行进行 k -means 聚类,将所有热图的行分成 5 组。
km = kmeans(mat_meth, centers = 5)$cluster
当初有以下矩阵 / 向量须要被可视化为热图。
mat:一个矩阵,其中各行对应不同的甲基化区域(DMRs)。矩阵中的值是每个样本中 DMR 的均匀甲基化程度。
expr:一个矩阵,其中的行对应于与 DMR 相干的基因(即与 DMR 最近的基因)。矩阵中的值是每个样本中每个基因的表白程度。每个基因在不同样本中的表白量是有比例的。
direction:甲基化变动的方向(hyper 示意肿瘤样本中甲基化程度较高,hypo 示意肿瘤样本中甲基化水平较低)。
pvalue:甲基化和相干基因表白之间的相干测试的 P 值。数值是 -log10 转换的。
type:基因的类型(如蛋白质编码基因或林肯 RNAs)。
gene:对基因模型的正文(即基因间、基因内或转录起始地位(TSS))。
dist:从 DMRs 到被鉴定基因的 TSS 的间隔。
enhancer:每个 DMR 中与增强器重叠的局部。
在这些变量中,mat\_meth、mat\_expr、cor\_pvalue、dist 和 anno\_enhancer 是数字变量,我为它们设置了色彩映射函数。对于其余变量,我设置了命名的色彩向量。
在上面的代码中,我在 heatmap()的第一次调用中指定了决裂,这是甲基化热图。轨道的高度是手动调整的。
col_meth = colorRamp2(c(0, 0.5, 1), c("blue", "white", "red"))
heatmap(mat, split col_meth, track.height = 0.12)
循环热图看起来很美! 因为矩阵中的行是基因组区域(差别甲基化区域),如果咱们能在一些区域之间建立联系,例如三维染色体构造中的物理相互作用,那么这个图就会更丑陋、更有用。
在上面的代码中,我在 DMRs 之间生成一些随机的相互作用。df_link 中的每一行意味着有一个从第 i 个 DMR 到第 j 个 DMR 的互动。
data.frame(
from_index
to_index
在圆形热图上找到这些 DMR 的地位是很辣手的。请看上面代码中的正文。留神这里的子集和行序元数据是通过 get.data()函数明确指定扇形索引来检索的。
for(i in seq_len(nrow)) {# 让咱们把索引为 df\_link$from\_index\[i\]的 DMR 称为 DMR1。# 另一个索引为 df\_link$to\_index\[i\]的为 DMR2。# DMR1 所在的扇区。group1 = km\[from_index\[i\] \]。# DMR2 所在的扇区。group2 = km\[to_index\[i\] \]。# 区块 \`group1\` 中的 DMRs 子集(来自 mat_meth 的行指数)。data("subset", sector.index = group1)
# 扇区 \`group1\` 中的行排序。data("row_order", sector.index = group1)
# 这是 DMR1 在 \`group1\` 热图中的地位。which(subset1\[row\_order1\] == from\_index\[i\])
# 区块 \`group2\` 中的 DMRs 子集(来自 mat_meth 的行指数)。data("subset", sector.index = group2)
# 扇区 \`group2\` 中的行排序。ret.data("r sector.indexoup2)
# 这是 DMR2 在 \`group2\` 热图中的地位。index\[i\])
# 咱们取两头的点,在 D 和 DMR2 之间画一个链接
link(group1, x1 - 0.5, grup2,col = rcoor(1))
我实现了一个函数。当初绘制矩阵行之间的链接更简略。
for(i in seq\_len(nrow(df\_link))) {heatmap.link(from\[i\],
to_index\[i\],
color(1))
}
增加链接后,绘图看起来更丑陋了!
图例对于了解热图十分重要。
绘制圆形图的函数只是后面代码的一个封装,没有任何批改。
图例对于了解热图十分重要。依照该链接的阐明,咱们须要一个绘制圆形图的函数和一个 Legends 对象。
当初咱们应用 gridBase 来联合根底图形和网格图形。
h = dev.size()\[2\]
draw(lgd_list, x = circle, just = "left")
最受欢迎的见解
1.R 语言动态图可视化:如何、创立具备精美动画的图
2.R 语言生存剖析可视化剖析
3.Python 数据可视化 -seaborn Iris 鸢尾花数据
4.r 语言对布丰投针(蒲丰投针)试验进行模仿和动静
5.R 语言生存剖析数据分析可视化案例
6.r 语言数据可视化剖析案例:摸索 brfss 数据数据分析
7.R 语言动静可视化:制作历史寰球平均温度的累积动静折线图动画 gif 视频图
8.R 语言高维数据的主成分 pca、t-SNE 算法降维与可视化剖析案例报告
9.python 主题 LDA 建模和 t -SNE 可视化