R语言绘制圈图、环形热图可视化基因组实战:展示基因数据比较

简介: R语言绘制圈图、环形热图可视化基因组实战:展示基因数据比较

可以使用环状图形展示基因数据比较。可以添加多种图展信息,如热图、散点图等。

本文目标:

  • 可视化基因组数据

制作环形热图

环形热图很漂亮。可以通过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}θ=90开始。

heatmap(
    bg.border )

如果split参数的值是一个因子,那么因子水平的顺序控制热图的顺序。如果split是一个简单的向量,热图的顺序是unique(split)。



点击标题查阅往期内容


R语言k-means聚类、层次聚类、主成分(PCA)降维及可视化分析鸢尾花iris数据集


左右滑动查看更多

01

02

03

04




# 注意,因为在前一个图中调用了 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")

相关文章
|
3月前
|
数据采集 机器学习/深度学习 数据可视化
R语言从数据到决策:R语言在商业分析中的实践
【9月更文挑战第1天】R语言在商业分析中的应用广泛而深入,从数据收集、预处理、分析到预测模型构建和决策支持,R语言都提供了强大的工具和功能。通过学习和掌握R语言在商业分析中的实践应用,我们可以更好地利用数据驱动企业决策,提升企业的竞争力和盈利能力。未来,随着大数据和人工智能技术的不断发展,R语言在商业分析领域的应用将更加广泛和深入,为企业带来更多的机遇和挑战。
|
4月前
|
数据可视化 数据挖掘 数据处理
R语言高级可视化技巧:使用Plotly与Shiny制作互动图表
【8月更文挑战第30天】通过使用`plotly`和`shiny`,我们可以轻松地创建高度互动的数据可视化图表。这不仅增强了图表的表现力,还提高了用户与数据的交互性,使得数据探索变得更加直观和高效。本文仅介绍了基本的使用方法,`plotly`和`shiny`还提供了更多高级功能和自定义选项,等待你去探索和发现。希望这篇文章能帮助你掌握使用`plotly`和`shiny`制作互动图表的技巧,并在你的数据分析和可视化工作中发挥更大的作用。
|
3月前
R语言基于表格文件的数据绘制具有多个系列的柱状图与直方图
【9月更文挑战第9天】在R语言中,利用`ggplot2`包可绘制多系列柱状图与直方图。首先读取数据文件`data.csv`,加载`ggplot2`包后,使用`ggplot`函数指定轴与填充颜色,并通过`geom_bar`或`geom_histogram`绘图。参数如`stat`, `position`, `alpha`等可根据需要调整,实现不同系列的图表展示。
|
3月前
|
数据采集 数据可视化 数据挖掘
R语言在金融数据分析中的深度应用:探索数据背后的市场智慧
【9月更文挑战第1天】R语言在金融数据分析中展现出了强大的功能和广泛的应用前景。通过丰富的数据处理函数、强大的统计分析功能和优秀的可视化效果,R语言能够帮助金融机构深入挖掘数据价值,洞察市场动态。未来,随着金融数据的不断积累和技术的不断进步,R语言在金融数据分析中的应用将更加广泛和深入。
|
4月前
|
数据可视化
R语言可视化设计原则:打造吸引力十足的数据可视化
【8月更文挑战第30天】R语言可视化设计是一个综合性的过程,需要综合运用多个设计原则来创作出吸引力十足的作品。通过明确目标、选择合适的图表类型、合理运用色彩与视觉层次、明确标注与引导视线以及引入互动性与动态效果等原则的应用,你可以显著提升你的数据可视化作品的吸引力和实用性。希望本文能为你提供一些有益的启示和帮助。
|
4月前
|
数据采集 存储 数据可视化
R语言时间序列分析:处理与建模时间序列数据的深度探索
【8月更文挑战第31天】R语言作为一款功能强大的数据分析工具,为处理时间序列数据提供了丰富的函数和包。从数据读取、预处理、建模到可视化,R语言都提供了灵活且强大的解决方案。然而,时间序列数据的处理和分析是一个复杂的过程,需要结合具体的应用场景和需求来选择合适的方法和模型。希望本文能为读者在R语言中进行时间序列分析提供一些有益的参考和启示。
|
2月前
|
数据挖掘 C语言 C++
R语言是一种强大的统计分析工具,提供了丰富的函数和包用于时间序列分析。
【10月更文挑战第21天】时间序列分析是一种重要的数据分析方法,广泛应用于经济学、金融学、气象学、生态学等领域。R语言是一种强大的统计分析工具,提供了丰富的函数和包用于时间序列分析。本文将介绍使用R语言进行时间序列分析的基本概念、方法和实例,帮助读者掌握R语言在时间序列分析中的应用。
50 3
|
7月前
|
数据可视化 数据挖掘 API
【R语言实战】聚类分析及可视化
【R语言实战】聚类分析及可视化
|
7月前
|
机器学习/深度学习 数据可视化
R语言逻辑回归logistic模型ROC曲线可视化分析2例:麻醉剂用量影响、汽车购买行为2
R语言逻辑回归logistic模型ROC曲线可视化分析2例:麻醉剂用量影响、汽车购买行为
|
4月前
|
机器学习/深度学习 数据采集 数据可视化
R语言在数据科学中的应用实例:探索与预测分析
【8月更文挑战第31天】通过上述实例,我们展示了R语言在数据科学中的强大应用。从数据准备、探索、预处理到建模与预测,R语言提供了完整的解决方案和丰富的工具集。当然,数据科学远不止于此,随着技术的不断发展和业务需求的不断变化,我们需要不断学习和探索新的方法和工具,以更好地应对挑战,挖掘数据的潜在价值。 未来,随着大数据和人工智能技术的普及,R语言在数据科学领域的应用将更加广泛和深入。我们期待看到更多创新的应用实例,为各行各业的发展注入新的动力。