经历过整个过程，才会对最终结果弥足珍惜。————The journey is the reward.
就像之前说的，在每天发的推送中，我准备把我自己开发的 R 包所经历的过程来放进来。
我希望现在把 0.0.0.1 版本先上传上来，让大家可以看看它是如何一步一步成长起来的，让大家体验一下这个 R 包不断变强的过程。
既然我已经上传到了 GitHub，那么它就是开源的，大家可以尽情去看里面的源代码。在开源协议上，我选择了一个比较保守的协议（a copyleft license，GPL3），因为我希望自己的工作能够得到欣赏。如果有任何引用最好注明出处。
所以现在我一天最多会推送三篇文章：
-
第一篇是公众号简介部分讲的内容
-
第二篇是“每天简化地学习一个基因/家族”
-
第三篇是 R 包开发的动态
关注我的 github（https://github.com/yudalang3/PathwayI，或者gitee: https://gitee.com/yudalang3/projects），我新开了一个 github 账号，后面专门用来开发 R 包。下面是本文的目录：
-
安装该包
-
可视化转录因子，共表达基因和 TF 靶基因
-
快速使用
-
一些说明
-
具有生物学意义的例子
-
难道是一个玩具？
-
-
两种基本布局的进化树绘制
-
生成树
-
矩形布局
-
环形布局
-
定制的可视化，调整连接方式
-
随便挑一个 ggtree 里面的图进行重复
-
-
一个最简单的底物、产物图绘制
-
让我们来绘制热图吧
-
什么？这竟然也能绘制热图
-
挑一个 ComplexHeatmap 的热图的局部我们来进行重复
-
-
说明
安装该包
# install.packages("devtools")
devtools::install_github("yudalang3/PathwayI")
# 通过调用help查看一下是否已经安装成功。
help(package='PathwayIlluminator')
安装 github 上的包竟然如此简单。什么？你访问不了 github? 我也注册了 gitee，你可以直接访问 gitee 上的项目。
# install.packages("remotes")
# install.packages("git2r")
# 上面两个R包需要被安装
remotes::install_git("https://gitee.com/yudalang3/PathwayI.git")
下面介绍本文主要内容，功能：
可视化转录因子，共表达基因和 TF 靶基因
快速使用
library(PathwayIlluminator)
display_TF_targets_coexpress()
display_TF_targets_coexpress(buttomStyleBezier = F)
运行之后给出的两个图如下， buttomStyleBezier 参数影响的是黄色区域底部的展示方式。
这个图毋庸置疑是浅显易懂的，带箭头的线表示直接的 TF 和靶基因的关系，虚线表示共表达关系。下面的贝塞尔连接线表示的是对应关系。
看到没有，其实“富集的思想”已经蕴含在里面了。后面我会专门讲富集分析。这里就对应着一个生物学问题：“这些和 TF 共表达的基因是否显著富集在了 TF 的靶基因中”？或者也可以相反，TF 及其靶基因是否有共表达关系。
小技巧：在Rstudio中，有很多曲线的图最好输出中矢量图，例如`pdf`，否则默认的绘图设备是png，会有很强的锯齿感。
一些说明：
你可能有疑问：咦，数据呢？我什么都没有输入为什么会出来一个图？此时我们打开这个函数的说明：
奥奥，原来是这里给出了默认的数据啊。嗯，根据描述只要给出一个 list 即可，有三个 key 是固定的：TF, targets 和 coexpression。
看 Details 部分，我们要注意一下这里面输入的格式与相关的生物学意义。
-
coexprerssion 对应的字符串向量是和 TF 有共表达关系的基因，当然这些基因可能属于 TF 的靶基因，也可能不是。
-
targets 后面给出的字符串向量就是 TF 的所有靶基因。
-
要注意 coexpression 的内容中要讲属于 targets 的字符串放在前面部分，而不是乱放。（当然这个功能可以加，欢迎各位用户来提需求）
具有生物学意义的例子
尽管上面已经放了一个可视化的例子，但是缺乏实际的使用场景。因此我挑选了一个我们组研究的信号通路作为一个具体的例子。
WNT 信号通路的下游转录因子为 TCF/LEF，有四个同源基因可以编码这个 TF。它的靶基因有很多，我们挑选基因作为示例。因此我们可以组成如下的一个具体的例子：
set_global_pars('fontsize', 7)
data = list(
TF = "TCF/LEF",
targets = c('Axin2','TCF7','ST7','BZRAP1', 'PDK2', 'PIK3C2A', 'CLK1' , 'SERPINB1','SPEN','NUCB2','JMJD6','DGKD','MLLT10'),
coexpression = c('SPEN','NUCB2','JMJD6','DGKD','MLLT10','AXIN2','PIK3C2A','SERPINB1','TP53','CXCL1','INS')
)
display_TF_targets_coexpress(data = data,buttomStyleBezier = T)
我想这个 fontsize 参数是不言而喻的。
难道是一个玩具？
不过这好像是一个很简单的玩具啊，我们实际的组学分析里面，有成千上万个基因。共表达基因动辄几百几千个，靶基因也有这么多。我怎么用呢？
不要心急，你会看到这个包慢慢强大起来的。当然首先你自己要认识清楚，这么多基因肯定是不可能继续用这种方式来可视化的。所以，请期待后面的博客。
那么现在这个功能有什么用呢？我觉得还是有一个功能的可以用的。但是前提是你要自己对这些基因进行归类。
假如，你对共表达基因进行归类，你分为了 g1, g2, g3, g4 四类和 k1, k2, k3 ,k4 四类，共八类。这些是你根据你对自己问题的理解分的。接着你又对几千个靶基因进行了分类，分为了 g1~g10 十大类。那么你就可以用这个功能了。
set_global_pars('fontsize', 12)
data = list(
TF = "TF1",
targets = paste0("g", 1:10),
coexpression = c(paste0("g", 1:4), paste0("k", 1:4))
)
display_TF_targets_coexpress(data = data)
两种基本布局的进化树绘制
进化树？我没有听错吧，这个东西很进化树有什么关系。
你注意到了吗？这里面有很多的节点，这些节点很进化树其实可以是同一个数据结构呢！所以呢，我顺手就把进化树的可视化就做好了。我这里都可以叫做 GraphicTree。实际上我这里的继承关系是这样的：
ListTreeNode <-- GraphicNode <-- BioGraphicNode
生成树
那么在下一步进行之前我们如何生成树呢？最简单的的，我们直接 new 出来不就好了吗？（new 指的是面向对象编程里面的一个术语，就是创建一个对象的时候，我们会用 new 关键字。）然后根据这个拓扑结构的关系，将其组成一棵树。
手动生成一棵树
# 我们手动生成一棵树
node1 <- GraphicNode$new(1)
node2 <- GraphicNode$new(2)
node1.2 <- GraphicNode$new()
node1.2$addChild(node1)
node1.2$addChild(node2)
node3 <- GraphicNode$new(3)
node4 <- GraphicNode$new(4)
node3.4 <- GraphicNode$new()
node3.4$addChild(node3)
node3.4$addChild(node4)
root <- GraphicNode$new()
root$addChild(node1.2);
root$addChild(node3.4)
displayTheTree(root)
上面的代码就是人工手动产生一棵进化树结构的代码。它输出的结果是：
Current id is 1.NA Children: 2 id is 1.NA id is 1.NA
Current id is 1.NA Children: 2 ID 1 ID 2
Current ID 1 is leaf.
Current ID 2 is leaf.
Current id is 1.NA Children: 2 ID 3 ID 4
Current ID 3 is leaf.
Current ID 4 is leaf.
嗯？这个 1.NA 是个什么东西。我们还需要执行：
assignTheCGBID(root)
displayTheTree(root)
此时出来的结果就是一个以表格形式展示出来的树结构。
Current id is 1.3 Children: 2 id is 1.2 id is 3.4
Current id is 1.2 Children: 2 ID 1 ID 2
Current ID 1 is leaf.
Current ID 2 is leaf.
Current id is 3.4 Children: 2 ID 3 ID 4
Current ID 3 is leaf.
Current ID 4 is leaf.
以第一行为例，这个信息就是说，当前节点的名称叫做 1.3， 它有两个孩子， 第一个叫做 1.2 第二个叫做 3.4
你再仔细瞧瞧，这个方式表示的结构是不是很像 R 中的 “hclust" 对象和 “ape”对象啊。准确地说和"ape"包中的phylo对象很像，只不过它第一行是扩展成了两行，因为节点1.3有两个孩子，所以可以写成两行。
后面我们要专门出一期来讲讲如何在计算机编程中如何表示树对象。
好，你可能注意到了还有一个疑问，这个 assignTheCGBID方法是个什么东西？
他的作用就是生成这个 _1.3_， _1.2_ 的字符串。这个是我们在已发表的论文中提出来的一种标记内节点的方法。
详见：https://doi.org/10.1093/bib/bbab583 内容在附录里面。
直接根据层次聚类的结果生成树
你可能又会问：嗯？难道你们做的又是一个玩具吗？你刚才是一棵手动生成的树，我要一棵我们实际数据分析生成的树。
好，那么我们就用 R 中的层次聚类方法（hclust）生成的树。对了，值得一说，原来进化中的 UPGMA 方法就是其中层次聚类的一种，因为层次聚类有两个重要参数。在给定一个距离矩阵之后，我们需要有叶子（OUT）之间的距离计算方法和组别之间的计算方法。调好参数之后，UPGMA 就是其中的一种。我们从这个实际结果将其转换为一棵树。
# Get Tree From the hclust object
hc <- hclust(dist(USArrests), "ave")
dendrogram <- as.dendrogram(hc)
intuitiveTree <- process_dendrogram2intuitiveTree(dendrogram)
displayTheTree(intuitiveTree$rootNode)
其实我们这个树结构是比较直观的树结构，和 dendrogram 差不多，它是一个 list in list 结构，属性放到了每个节点的 attribute 里面。
调用process_dendrogram2intuitiveTree函数可以将 dendrogram 对象转化为我们的树对象。
你可能要问：为什么不直接用 dendrogram 这个结构，或者直接产生的 hclust 结构。不直接用 hclust 结构是因为这种像表格一样的记录方式不利于实现一些算法，你想想看，为什么信息科学里面二叉树，多叉树，红黑树，平衡树，这些数据结构不用一个多维数组来表示？
像 C 语言会用结构体，C++/JAVA 里面是直接用类来表示节点。因为这样表示之后可以很方便的实现一些算法。至于这个 dendrogram 结构么，主要是因为 R 中的 S3/S4 面向对象编程系统，它没有引用语义（reference sematics） 这真的是太要命了。所以每次赋值之后，都会触发 R 的 copy-on-modify机制。没办法，我只能借助于基于 environment 的 R6，OOP。
除了上面两种生成树的方式之外，我还写了一个模拟生成满多叉树的一个函数：
ret <- process_createNodes(2,numOfChildren = 2)
displayTheTree(ret$rootNode)
得到的结果 ret 变量是一个 GraphicTree 对象的实例，下面我们来进行可视化。
矩形布局
这是最简单的方式，我们用前面我们手动创建的进化树进行可视化，
# Create tree manually ----------------------------------------------------
node1 <- GraphicNode$new(1)
node2 <- GraphicNode$new(2)
node1.2 <- GraphicNode$new()
node1.2$addChild(node1)
node1.2$addChild(node2)
node3 <- GraphicNode$new(3)
node4 <- GraphicNode$new(4)
node3.4 <- GraphicNode$new()
node3.4$addChild(node3)
node3.4$addChild(node4)
root <- GraphicNode$new()
root$addChild(node1.2);
root$addChild(node3.4)
displayTheTree(root)
assignTheCGBID(root)
displayTheTree(root)
tree <- GraphicTree$new(root)
tree$assignAttributes()
# set the root node length to be 0
root$branchLength <- 0
rectLayoutDesigner<- RectangularLayoutDesigner$new()
rectLayoutDesigner$calculate(tree)
rectLayoutDesigner$plotTree(tree)
circLayoutDesigner <- CircularLayoutDesigner$new();
# adjust the start degree to make it more beautiful
tree$startDegree <- 45
circLayoutDesigner$calculate(tree)
circLayoutDesigner$plotTree(tree)
于是得到下面两张图片：
下面我们用 hclust 函数生成的分析之后形成的进化树来验证我们的程序：
# Get Tree From the hclust object ----------------------------------------------------
hc <- hclust(dist(USArrests), "ave")
dendrogram <- as.dendrogram(hc)
intuitiveTree <- process_dendrogram2intuitiveTree(dendrogram)
displayTheTree(intuitiveTree$rootNode)
intuitiveTree$innerRadius <- 1
circLayoutDesigner <- CircularLayoutDesigner$new();
circLayoutDesigner$calculate(intuitiveTree)
circLayoutDesigner$plotTree(intuitiveTree)
rectLayoutDesigner<- RectangularLayoutDesigner$new()
rectLayoutDesigner$calculate(intuitiveTree)
rectLayoutDesigner$plotTree(intuitiveTree)
上面是我们的绘图函数绘制的效果。
如果直接 plot(dendrogram, horiz = T)，给出的结果是：
只是转了一个方向，结果是一样的。如果大家有这个需求，我也可以绘制纵向的可视化。但是你要注意原始的 hclust 对象你绘制一下好像不长这个样子，as.dendrogram()函数默认会进行截断操作，通过调整 hang = 0.2 可以得到和hclust结构一样的树。
为什么这个 as.dendxx 函数会这样呢？我估计主要是人们关注的是聚类的类群所属关系，对这个枝长不感冒，所以倾向于将叶子对齐。但是在进化里面这个枝长是有意义的。
注意，这里有一个有趣的彩蛋哦。假如，outterRadius比innerRadius值小怎么办呢？
这个问题很好，要注意默认的outterRadius值为-1，当布局计算程序（circLayoutDesigner$calculate(intuitiveTree)） 碰到这个值小于0时，会重新根据当前画布来计算合适的半径。那么我们来试一试吧！
intuitiveTree$innerRadius <- 2
intuitiveTree$outterRadius <- 0.7
circLayoutDesigner <- CircularLayoutDesigner$new();
circLayoutDesigner$calculate(intuitiveTree)
circLayoutDesigner$plotTree(intuitiveTree)
嗯？这个根结点怎么会事？它怎么延伸到圆心去了。哈，想要正确处理它，就继续看下去。看到了吗？大小互换之后变成了一个外翻的圆，你可以想想这是为什么？有空我再出一期讲一讲原理。
环形布局
其实上面的例子已经有了，我们再调整参数使其更加美观一点：
# adjust the start degree to make it more beautiful
tree$startDegree <- 0
tree$extendDegree <- 90
tree$innerRadius <- 1
circLayoutDesigner$calculate(tree)
circLayoutDesigner$plotTree(tree)
于是得到了下面这个图：
猜一猜，这调整的几个参数，都是什么意思。
下面绘制一个饱满一点的树，用环形可视化来看看。
# Create tree by simulation ----------------------------------------------------
ret <- process_createNodes(3,numOfChildren = 6)
displayTheTree(ret$rootNode)
ret$innerRadius <- 1
circLayoutDesigner <- CircularLayoutDesigner$new();
circLayoutDesigner$calculate(ret)
circLayoutDesigner$plotTree(ret)
rectLayoutDesigner<- RectangularLayoutDesigner$new()
rectLayoutDesigner$calculate(ret)
rectLayoutDesigner$plotTree(ret)
确实饱满的更好看。
定制的可视化，调整连接方式
这里是很重要的部分，我们可视化往往是需要进行一大堆的定制的需求，那么如何实现多种定制的绘制方法。我看很多人还因为这个定制化写了几个 R 包，最后发到了比较有影响力的期刊上。
因此我推出了 treeCustomizedDrawer()函数，这个函数接受四个参数：
当传入 tree 对象之后，可以填写三个 drawer，这样你就可以绘制你要的任何一种图形了。后面根据这个需求，我们来了解如何使用这个函数。
从需求出发来绘图
不需要绘制节点，直接绘制最最简单的线条就行：
## the simplest tree drawer -------------------------
# First, generate the tree
hc <- hclust(dist(USArrests), "ave")
dendrogram <- as.dendrogram(hc)
intuitiveTree <- process_dendrogram2intuitiveTree(dendrogram)
# calculate the information for the rectangular layout
rectLayoutDesigner<- RectangularLayoutDesigner$new()
rectLayoutDesigner$calculate(intuitiveTree)
root_node_drawer <- NULL
node_drawer <- function(parent, node) {
x_parent <- parent$xAxis_or_radius
y_parent <- parent$yAxis_or_angle
x_node <- node$xAxis_or_radius
y_node <- node$yAxis_or_angle
grid.segments(
x0 = x_parent,
y0 = y_parent,
x1 = x_node,
y1 = y_node,
default.units = 'in'
)
}
treeCustomizedDrawer(tree = intuitiveTree,
root_node_drawer = root_node_drawer,
inner_node_drawer = node_drawer,
leaf_drawer = node_drawer)
# calculate the information for the circular layout
intuitiveTree$innerRadius <- 1
circLayoutDesigner <- CircularLayoutDesigner$new();
circLayoutDesigner$calculate(intuitiveTree)
xCenter <- circLayoutDesigner$xCenter;
yCenter <- circLayoutDesigner$yCenter;
node_drawer <- function(parent, node) {
radius_parent <- parent$xAxis_or_radius
angle_parent <- parent$yAxis_or_angle
radius_node <- node$xAxis_or_radius
angle_node <- node$yAxis_or_angle
parent_coor <-
polar2cartesianCoor(radius_parent, angle_parent, xCenter, yCenter)
node_coor <-
polar2cartesianCoor(radius_node, angle_node, xCenter, yCenter)
grid.segments(
x0 = parent_coor[1],
y0 = parent_coor[2],
x1 = node_coor[1],
y1 = node_coor[2],
default.units = 'in'
)
}
treeCustomizedDrawer(tree = intuitiveTree,
root_node_drawer = root_node_drawer,
inner_node_drawer = node_drawer,
leaf_drawer = node_drawer)
效果是这样的：
环形布局直接这样绘制效果不好，我们还是把 arc 补上吧。包中给出了一个 draw_arc 函数。
# pretty circular layout
node_drawer <- function(parent, node) {
radius_parent <- parent$xAxis_or_radius
angle_parent <- parent$yAxis_or_angle
radius_node <- node$xAxis_or_radius
angle_node <- node$yAxis_or_angle
parent_coor <-
polar2cartesianCoor(radius_parent, angle_parent, xCenter, yCenter)
node_coor <-
polar2cartesianCoor(radius_node, angle_node, xCenter, yCenter)
radius_x1 <- radius_parent
angle_y1 <- angle_node
point_coor <-
polar2cartesianCoor(radius_x1, angle_y1, xCenter, yCenter)
grid.segments(
x0 = point_coor[1],
y0 = point_coor[2],
x1 = node_coor[1],
y1 = node_coor[2],
default.units = 'in'
)
draw_arc(angle_node, angle_parent, radius_parent, xCenter, yCenter,default_unit = 'in')
}
treeCustomizedDrawer(tree = intuitiveTree,
root_node_drawer = root_node_drawer,
inner_node_drawer = node_drawer,
leaf_drawer = node_drawer)
回归正常。
需求：我记得也有人用贝塞尔曲线来绘制边
当然可以了，我们直接调整这个 node_drawer 就好了。这里给出矩形布局的绘制方式，环形布局也是一样的道理，期待你自己去实现，如果需求多，我就给出脚本。
## use bezier curve the draw the rectangular tree
# First, generate the tree
hc <- hclust(dist(USArrests), "ave")
dendrogram <- as.dendrogram(hc)
intuitiveTree <- process_dendrogram2intuitiveTree(dendrogram)
# calculate the information for the rectangular layout
rectLayoutDesigner<- RectangularLayoutDesigner$new()
rectLayoutDesigner$calculate(intuitiveTree)
root_node_drawer <- NULL
node_drawer <- function(parent, node) {
x_parent <- parent$xAxis_or_radius
y_parent <- parent$yAxis_or_angle
x_node <- node$xAxis_or_radius
y_node <- node$yAxis_or_angle
draw_bezier_twoPoints(
x0 = x_parent,
y0 = y_parent,
x1 = x_node,
y1 = y_node,
default_unit = 'in',
vertical = F
)
}
treeCustomizedDrawer(tree = intuitiveTree,
root_node_drawer = root_node_drawer,
inner_node_drawer = node_drawer,
leaf_drawer = node_drawer)
结果是：
我们提供了一个贝塞尔曲线快速绘制的函数：draw_bezier_twoPoints()
这种曲线一般比较适合，cladogram，什么是cladogram呢？就是每个节点的枝长都相同，我们模拟产生的那个树就是一个cladogram。
# First, generate the tree
intuitiveTree <- process_createNodes(2,numOfChildren = 6)
# calculate the information for the rectangular layout
rectLayoutDesigner<- RectangularLayoutDesigner$new()
rectLayoutDesigner$calculate(intuitiveTree)
root_node_drawer <- NULL
node_drawer <- function(parent, node) {
x_parent <- parent$xAxis_or_radius
y_parent <- parent$yAxis_or_angle
x_node <- node$xAxis_or_radius
y_node <- node$yAxis_or_angle
draw_bezier_twoPoints(
x0 = x_parent,
y0 = y_parent,
x1 = x_node,
y1 = y_node,
default_unit = 'in',
vertical = F
)
}
treeCustomizedDrawer(tree = intuitiveTree,
root_node_drawer = root_node_drawer,
inner_node_drawer = node_drawer,
leaf_drawer = node_drawer)
需求：我要显示每个叶子的名字
代码奉上：
## Dispaly the leaf name --------------------------------
# First, generate the tree
global_pars$blank_area_ratio$r <- 0.2
hc <- hclust(dist(USArrests), "ave")
dendrogram <- as.dendrogram(hc)
intuitiveTree <- process_dendrogram2intuitiveTree(dendrogram)
# calculate the information for the rectangular layout
rectLayoutDesigner<- RectangularLayoutDesigner$new()
rectLayoutDesigner$calculate(intuitiveTree)
root_node_drawer <- NULL
inner_node_drawer <- function(parent, node) {
x_parent <- parent$xAxis_or_radius
y_parent <- parent$yAxis_or_angle
x_node <- node$xAxis_or_radius
y_node <- node$yAxis_or_angle
grid.lines(x = c(x_parent, x_parent, x_node), y = c( y_parent, y_node, y_node),default.units = 'in')
}
leaf_drawer <- function(parent, node) {
inner_node_drawer(parent, node)
# draw the leaf label
x_node <- node$xAxis_or_radius
y_node <- node$yAxis_or_angle
grid.text(
label = node$label,
just = 'left',
x= x_node + 0.1,
y = y_node,
default.units = 'in',
gp = gpar(fontsize = 8)
)
}
treeCustomizedDrawer(tree = intuitiveTree,
root_node_drawer = root_node_drawer,
inner_node_drawer = inner_node_drawer,
leaf_drawer = leaf_drawer)
效果：
这里实现这个功能的主要函数是编写leaf_drawer()函数，但是这个函数其实相对于inner_node_drawer()对了一个绘制字符串的过程，所以直接调用了上面一个函数，再格外加了绘制字符串的代码。
什么？你还看不懂这个代码？没关系，学呗，就是 grid 包的一些基础用法。
如果它提示找不到函数，那么就将library(grid)代码加上去。
那么环状的结果如何呢？也类似不过字体需要旋转一下。如果读者多，我就写一下。
随便挑一个 ggtree 里面的图进行重复
这个包貌似宣传做得特别好，我们根据上面讲的定制可视化的方式来绘制一个 ggtree 的图。
P.S. 我们这些实现方式有点类似于低级函数，ggtree 能绘制的我们都能绘制，它不能绘制的，我们也能。
这里我将两幅比较有趣的图进行组合，把他们的特性都放进去来进行绘制。
挑选了这两幅图，第一幅的特性就是对特定的节点绘制特定的颜色，并高亮一些叶子而已。这里能看出来大意，叶子节点就是根据一定的物种分类标记颜色，内节点就是根据所包含叶子颜色最多的一种来确定颜色。内节点的染色实现有多种方案，因为本包的目的不是写一个多功能的进化树可视化，这个进化树可视化只是一个附带的产品，所以先不费时间去写了。比如我提供一个思路，在 GraphicNode 类里面建立一个 list 属性，记录下当前节点中包含叶子阶段所归属类群的信息。然后用后序遍历先计算好每个内节点，哪种类群最多，记录一个颜色属性。接着绘制颜色即可。我们来对叶子进行染色。
第二幅图就是把叶子线条补齐，label 加了几个颜色，就这么简单。我们用矩形布局做示范吧，当然我们不知道它的 nwk 数据是什么，所以我们还是用 hclust，用其他的数据集去模仿效果好了。
### Any ggtree picture can be drawed --------------------------------
# First, generate the tree
global_pars$blank_area_ratio$r <- 0.2
hc <- hclust(dist(USArrests), "ave")
dendrogram <- as.dendrogram(hc, hang = 0.1)
intuitiveTree <- process_dendrogram2intuitiveTree(dendrogram)
# calculate the information for the rectangular layout
rectLayoutDesigner<- RectangularLayoutDesigner$new()
rectLayoutDesigner$calculate(intuitiveTree)
root_node_drawer <- NULL
inner_node_drawer <- function(parent, node) {
x_parent <- parent$xAxis_or_radius
y_parent <- parent$yAxis_or_angle
x_node <- node$xAxis_or_radius
y_node <- node$yAxis_or_angle
grid.lines(x = c(x_parent, x_parent, x_node), y = c( y_parent, y_node, y_node),default.units = 'in')
}
colorMapper <- function(x){
if (startsWith(x = x, prefix = 'M')) {
return('wheat2')
}else if(startsWith(x = x, prefix = 'N')){
return('yellowgreen')
}else {
return('purple')
}
}
leaf_drawer <- function(parent, node) {
x_parent <- parent$xAxis_or_radius
y_parent <- parent$yAxis_or_angle
# draw the leaf label
x_node <- node$xAxis_or_radius
y_node <- node$yAxis_or_angle
## draw the leaf line color
color <- colorMapper(node$label)
grid.text(
label = node$label,
just = 'left',
x= unit(1 - global_pars$blank_area_ratio$r, 'npc') + unit(0.1,'in'),
y = y_node,
default.units = 'in',
gp = gpar(fontsize = 8, col = color)
)
# Fill in the missing lines
grid.segments(
x0 = x_node,
y0 = y_node,
x1 = unit(1 - global_pars$blank_area_ratio$r, 'npc'),
y1 = y_node,
default.units = 'in',
gp = gpar(lty = 'dashed',col = color)
)
grid.lines(x = c(x_parent, x_parent, x_node), y = c( y_parent, y_node, y_node),default.units = 'in',gp = gpar(col = color))
}
treeCustomizedDrawer(tree = intuitiveTree,
root_node_drawer = root_node_drawer,
inner_node_drawer = inner_node_drawer,
leaf_drawer = leaf_drawer)
完成。可以看到我在例子中用一个 colorMapper 模拟了用户的操作。我在这里用那 0.2 的右侧留白来定位叶子标签绘制的空间。
奥，对了，树和热图的组合图，会在后面的热图部分进行示范。
一个最简单的底物、产物图绘制
其实两个节点的树是可以用来很好的表示这个过程的，我来示例一下：
额，后来一想，我发现我有点陷进去了。为什么呢？因为绘制一个简单的三角形不就好了吗？何必用一棵树呢？
后面我会对这个代码的功能进行封装，叫做 simple_reaction_illustraction()?
让我们来绘制热图吧
热图原理：给定一个配色方案，指定颜色的冷暖程度，我们可以称其为 colorMapper。然后对于给定数据，利用 colorMapper 得到对应的颜色。最后画图即可。
什么这竟然也能绘制热图
难道你不觉得热图的原理实在是太简单了吗？现在流行的热图，无非就是环形的和矩形的。当然其它各种形状也行。
对于矩形的热图，这也实在太简单了，可能定义好这个 colorMapper 之后，两层 for 循环也就解决了。包中定义了两个 ColorMapper，一个是针对连续型随机变量的continuousMapper，和另一个针对离散型随机变量的discreteMapper。具体的应用可以看后面的代码。
P.S.: 这个 mapper 我是用 R 中的 S4 对象写的，和树结构的 R6 不同，还是要尝试一下 R 中原生的面向对象。
那么环形的热图呢？我们接上刚才的进化树绘制代码，也是很快解决。
下面我们来绘制一下。
矩形热图
我们直接上代码：
ret <- process_createNodes(2,numOfChildren = 3)
set_global_pars('blank_area_ratio', list(l = 0.05, r = 0.4, t = 0.1, b = 0.1))
rectLayoutDesigner<- RectangularLayoutDesigner$new()
rectLayoutDesigner$calculate(ret)
rectLayoutDesigner$plotTree(ret)
displayTheTree(ret$rootNode)
data <- matrix(rnorm(9 * 3), ncol = 3)
rownames(data) <- as.character(1:9)
mapper <- continuousMapper(colors = colorRampPalette(c("blue", "red"))( 50 ), range = range(data))
mapped_colors <- generate_colors(data, mapper)
rownames(mapped_colors) <- rownames(data)
root_node_drawer <- NULL
inner_node_drawer <- function(parent, node) {
x_parent <- parent$xAxis_or_radius
y_parent <- parent$yAxis_or_angle
x_node <- node$xAxis_or_radius
y_node <- node$yAxis_or_angle
grid.lines(x = c(x_parent, x_parent, x_node), y = c( y_parent, y_node, y_node),default.units = 'in')
}
leafHeight <- global_pars$leafIntervalLength
leaf_drawer <- function(parent, node) {
inner_node_drawer(parent, node)
# draw the leaf label
x_node <- node$xAxis_or_radius
y_node <- node$yAxis_or_angle
textGrob <- grid.text(
label = paste0('ID ',node$id),
just = 'left',
x= x_node + 0.1,
y = y_node,
default.units = 'in',
gp = gpar(fontsize = 12)
)
colo <- mapped_colors[paste0(node$id),]
xx <- unit(x_node + 0.35 * 1:length(colo), 'in') + grobWidth(textGrob)
grid.rect(
x = xx,
y = y_node,
width = 0.35,
height = leafHeight,
gp = gpar(fill = colo, col = NA),
hjust = 0,
vjust = 0.5,
default.units = 'in'
)
}
treeCustomizedDrawer(tree = ret,
root_node_drawer = root_node_drawer,
inner_node_drawer = inner_node_drawer,
leaf_drawer = leaf_drawer)
效果：
实现原理呢？其实很简单，重写 leaf_drawer 就行，我们在绘制叶子的时候再让它额外绘制一下矩形方块就行了。在示例中我设置了一个模拟数据，行名对应着叶子的 id。因为我没有赋予 Leaf label，当然你也可以行名是 leaf label。
中间用到了一个连续型数据到颜色的 mapper，就是这么简单的调用过去了。
环形热图
直接上代码：
# Circular heatmap --------------------------------------------------------
# intuitiveTree <- process_createNodes(2,numOfChildren = 3)
hc <- hclust(dist(USArrests), "ave")
dendrogram <- as.dendrogram(hc)
intuitiveTree <- process_dendrogram2intuitiveTree(dendrogram)
## convert the colors
color = colorRampPalette(rev(RColorBrewer::brewer.pal(n = 7, name ="RdYlBu")))(100)
mapper <- continuousMapper(colors = color, range = range(USArrests))
mapped_colors <- generate_colors(as.matrix(USArrests), mapper)
rownames(mapped_colors) <- rownames(USArrests)
## define the designer
circLayoutDesigner <- CircularLayoutDesigner$new();
intuitiveTree$innerRadius <- 0.5;
intuitiveTree$outterRadius <- 2
# get the right degree not easy
intuitiveTree$startDegree <- 30
intuitiveTree$extendDegree <- 270
circLayoutDesigner$calculate(intuitiveTree)
# circLayoutDesigner$plotTree(intuitiveTree)
half_leaf_degree <- get_global_pars()$leafIntervalLength * 0.5
xCenter <- circLayoutDesigner$xCenter;
yCenter <- circLayoutDesigner$yCenter;
# pretty circular layout
node_drawer <- function(parent, node) {
radius_parent <- parent$xAxis_or_radius
angle_parent <- parent$yAxis_or_angle
radius_node <- node$xAxis_or_radius
angle_node <- node$yAxis_or_angle
parent_coor <-
polar2cartesianCoor(radius_parent, angle_parent, xCenter, yCenter)
node_coor <-
polar2cartesianCoor(radius_node, angle_node, xCenter, yCenter)
radius_x1 <- radius_parent
angle_y1 <- angle_node
point_coor <-
polar2cartesianCoor(radius_x1, angle_y1, xCenter, yCenter)
grid.segments(
x0 = point_coor[1],
y0 = point_coor[2],
x1 = node_coor[1],
y1 = node_coor[2],
default.units = 'in'
)
draw_arc(angle_node, angle_parent, radius_parent, xCenter, yCenter,default_unit = 'in')
}
leaf_drawer <- function(parent, node) {
node_drawer(parent,node)
radius_parent <- parent$xAxis_or_radius
angle_parent <- parent$yAxis_or_angle
radius_node <- node$xAxis_or_radius
angle_node <- node$yAxis_or_angle
# parent_coor <-
# polar2cartesianCoor(radius_parent, angle_parent, xCenter, yCenter)
node_coor <-
polar2cartesianCoor(radius_node, angle_node, xCenter, yCenter)
colo <- mapped_colors[node$label, ]
for (i in 1:length(colo)) {
inner_circle_radius <- 0.1 *i
draw_sector(
xCenter = xCenter,
yCenter = yCenter,
startDegree = angle_node - half_leaf_degree,
endDegree = angle_node + half_leaf_degree,
circle_inner_radius = radius_node + inner_circle_radius,
circle_out_radius = radius_node + inner_circle_radius + 0.1,
fill_colo = colo[i]
)
}
}
treeCustomizedDrawer(tree = intuitiveTree,
root_node_drawer = NULL,
inner_node_drawer = node_drawer,
leaf_drawer = leaf_drawer)
是不是很方便啊，灵活度很高。那么我们看看结果：
你可能会问：
*. 咦？我在 Rstudio 中画出来的弧线没有那么圆润？
这是因为，默认 Rstuido 的绘图设备是 png，所以是个位图，我上图的截图是用矢量图绘制的。直接pdf() + dev.off()。
*. 咦？为什么那个根结点没有线段呢？应该有线段啊
那是因为我没有顶一个 root_node_drawer，我们设置一个好了：
root_node_drawer <-
function(rootNode) {
radius_node <- rootNode$xAxis_or_radius
angle_node <- rootNode$yAxis_or_angle
node_coor <-
polar2cartesianCoor(radius_node, angle_node, xCenter, yCenter)
grid.circle(
x = node_coor[1],
y = node_coor[2],
r = unit(3,'pt'),
gp = gpar(fill = 'black', col = NA),
default.units = 'in'
)
grid.segments(
x0 = xCenter,
y0 = yCenter,
x1 = node_coor[1],
y1 = node_coor[2],
default.units = 'in'
)
}
treeCustomizedDrawer(tree = intuitiveTree,
root_node_drawer = root_node_drawer,
inner_node_drawer = node_drawer,
leaf_drawer = leaf_drawer)
我标注了一个偌大的黑点。
挑一个 ComplexHeatmap 的热图的局部我们来进行重复
挑选这张图来绘制吧。其实这就是一个矩形布局。它在热图的顶部，右侧和底部分别绘制很多的其它图形。人们看着这个图好像很酷，其实就是很多东西的叠加而已，那我现在简单实现一下，毕竟我们这个包的这个功能只是一个附属功能。
所以这里先实现一些主要的特性，这里其实就是**分区域绘制多个热图**而已，在我看来，理论上其实是不复杂的。因为我们还没有将绘制的代码封装到一个 viewport （grid 包一个很重要的功能），所以现在实现其局部的特性，全局的图其实就是直接拼一拼。
我们主要挑选其中间的一个热图吧。
就是这个，它主要的特性其实就是叶子节点那么多了很多的矩形方块。当然这对应基因的很多不同注释。
# ComplexHeatmap partial mimic -----------------------------------------------------
ret <- process_createNodes(2,numOfChildren = 3)
set_global_pars('blank_area_ratio', list(l = 0.05, r = 0.6, t = 0.1, b = 0.1))
rectLayoutDesigner<- RectangularLayoutDesigner$new()
rectLayoutDesigner$calculate(ret)
data <- matrix(rnorm(9 * 3), ncol = 3)
rownames(data) <- as.character(1:9)
mapper <- continuousMapper(colors = colorRampPalette(c("navy", "white", "firebrick3"))( 50 ), range = range(data))
mapper2 <- continuousMapper(colors = c("blue", "red"), range = range(data))
mapped_colors <- generate_colors(data, mapper)
rownames(mapped_colors) <- rownames(data)
mapped_colors2 <- generate_colors(data, mapper2)
rownames(mapped_colors2) <- rownames(data)
root_node_drawer <- NULL
inner_node_drawer <- function(parent, node) {
x_parent <- parent$xAxis_or_radius
y_parent <- parent$yAxis_or_angle
x_node <- node$xAxis_or_radius
y_node <- node$yAxis_or_angle
grid.lines(x = c(x_parent, x_parent, x_node), y = c( y_parent, y_node, y_node),default.units = 'in')
}
leafHeight <- global_pars$leafIntervalLength
leaf_drawer <- function(parent, node) {
inner_node_drawer(parent, node)
# draw the leaf label
x_node <- node$xAxis_or_radius
y_node <- node$yAxis_or_angle
textGrob <- grid.text(
label = paste0('ID ',node$id),
just = 'left',
x= x_node + 0.1,
y = y_node,
default.units = 'in',
gp = gpar(fontsize = 12)
)
colo <- mapped_colors[paste0(node$id),]
xx <- unit(x_node + 0.5 * 1:length(colo), 'in') + grobWidth(textGrob)
grid.rect(
x = xx,
y = y_node,
width = 0.35,
height = leafHeight,
gp = gpar(fill = colo, col = NA),
hjust = 0,
vjust = 0.5,
default.units = 'in'
)
colo <- mapped_colors2[paste0(node$id),]
xx <- unit(x_node + 0.5 * length(colo) + 0.2 + 0.35 * 1:length(colo), 'in') + grobWidth(textGrob)
grid.rect(
x = xx,
y = y_node,
width = 0.35,
height = leafHeight,
gp = gpar(fill = colo, col = NA),
hjust = 0,
vjust = 0.5,
default.units = 'in'
)
}
treeCustomizedDrawer(tree = ret,
root_node_drawer = root_node_drawer,
inner_node_drawer = inner_node_drawer,
leaf_drawer = leaf_drawer)
# 最后再添加个图例
grid.text('Type', x = 0.95, y = 0.85)
grid.rect(x = 0.95 , y = 0.8, width = 0.05,height = 0.05, gp = gpar(col =NA, fill = 'blue'))
grid.rect(x = 0.95 , y = 0.75, width = 0.05,height = 0.05, gp = gpar(col =NA, fill = 'red'))
其实绘制的主要原理就是这样，像ggtree和ComplexHeatmap，主要的就在于节省人们的时间和精力吧，毕竟有直接能用的工具为什么不直接用呢？当然，我们可以有理论上的创新，我借用进化大牛 Joseph Felsenstein 在《Inferring Phylogenies》一书中的描述来看待这个写可视化 R 包与实现可视化的关系。
当然他说得比较书面化，直白地说人话就是：你只是做一个好用的工具，而且绘制出来还符合一定的美学，其实不算科研。如果你能够做到把其中的信息用一种新型的方式更好地展示出来，那是科研。你还能提出一派理论，更好的把专业知识和可视化元素结合起来，最后实现你的理论，那也是科研。
最后，我们当然是要进行理论上的创新，例如你开发了一个新的展示方式（布局），而不是把别人本来写 200 行代码的图，你用 5 行绘制出来了。
说明
此为每日更新的大型写实电视连续剧《小小计算生物学》第一季，From-2023-08-20。
欢迎大家阅读，点赞，转发，关注。你们的欣赏是我不断坚持的动力。
欢迎关注：
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
