用散点可视化一个矩阵
当然是ggplot2画,才好和ggtree兼容,你首先问自己,抛开树,你能右边吗?本身你有个矩阵,row的顺序要跟着树来,你可以先不管,先随意,这是y轴。那x轴呢?你矩阵里的数字,打点你总该会的,而点连线,同样column的数字打出来的点,就连在一起,你先按这个思路画出来,再来问后续的。
这是我在「知识星球」中的回答,然而小伙伴表示还是不懂。那么我来详细讲解一下:
首先我们读入一个FASTA数据:
require(treeio)
fasta <- system.file("examples/FluA_H3_AA.fas", package="ggtree")
aa <- read.fasta(fasta)
数据用treeio
包里的read.fasta
去读这个氨基酸序列,它有个好处是你直接as.character(aa[[i]])
就是个字母的向量,方便比较。
n = length(aa)
d = matrix(NA, ncol=n, nrow=n)
nm = labels(aa)
rownames(d) = colnames(d) = nm
for (i in seq_len(n)) {
for (j in seq_len(i)) {
d[i, j] = mean(as.character(aa[[i]]) != as.character(aa[[j]]))
d[j, i] = d[i, j]
}
}
这段代码生成了pairwise distance矩阵,你可以搞其它的统计量,都OK的。大家可以想一想,一般这种矩阵怎么可视化?用热图!热图用数字来填充颜色,而现在我们把数字当成x轴上的变量,然后打点,仅此而已,我们有没有办法画出来,当然可以。同样row的数值会在同一y轴上,如果是热图,同一column的数值会在同一x轴上,方便横向和纵向比较,而现在同一column的数值不在同一直线一了,因为我们打的点的位置使用了矩阵中的数值,同一column的点被拉扯到不同的位置上去了,所以为了能够比较或者是看清楚同一column上的点的走势,我们把同一column上的点用线条连接在一起。这就是小伙伴提问的图,就是这么画出来的,非常简单。我在「知识星球」里其实已经解答得很清楚了。
画图实战
require(tibble)
dd <- as_data_frame(d)
dd$seq1 <- nm
require(tidyr)
td <- gather(dd,seq2, dist, -seq1)
这里宽变长,这种tidy格式适合于我们用ggplot2
可视化。
> td
# A tibble: 5,776 x 3
seq1 seq2 dist
<fct> <chr> <dbl>
1 A/Swine/Binh_Duong/03_08/2010 A/Swine/Binh_Duong/03_08/2010 0.
2 A/Swine/Binh_Duong/03_10/2010 A/Swine/Binh_Duong/03_08/2010 0.00178
3 A/Swine/HK/2503/2011 A/Swine/Binh_Duong/03_08/2010 0.0267
4 A/Swine/HK_NS2439/2011 A/Swine/Binh_Duong/03_08/2010 0.0267
5 A/Swine/GX_NS2783/2010 A/Swine/Binh_Duong/03_08/2010 0.0249
6 A/Swine/GD/3767/2011 A/Swine/Binh_Duong/03_08/2010 0.0320
7 A/Swine/GX/2242/2011 A/Swine/Binh_Duong/03_08/2010 0.0338
8 A/Swine/GX/2803/2011 A/Swine/Binh_Duong/03_08/2010 0.0249
9 A/Swine/GX_NS3106/2011 A/Swine/Binh_Duong/03_08/2010 0.0285
10 A/Swine/GX_NS3108/2011 A/Swine/Binh_Duong/03_08/2010 0.0249
# ... with 5,766 more rows
看了上面的数据,看懂了我的解释,画图似乎是显而易见的。这里legend太多,所以用theme
去掉不显示,实际操作中,比如提问者的图中,颜色没这么多,一般分clade来上色,或者如果你单纯想画这个图,你可以聚类,然后cut tree,通过子树来分类上色,或者使用kmeans聚类,直接上k种颜色。
require(ggplot2)
ggplot(td, aes(dist, seq1, group=seq2, color=seq2)) +
geom_point() + geom_line() +
theme(legend.position="none")
那么如何与树关联,首先我们要可视化相应的树,然后找出taxa在这颗树上的排列(可视化的角度,从上到下的顺序),我在ggtree
里提供了get_taxa_name
这个函数,方便用户得到这个排列,这样方便用户按照这个排列去做别的图,然后可以去拼图,这样可以保证图中的顺序与树是对应的。你可以用别的软件做图,你可以在AI里拼图,只要你能够按照顺序对应就行。
require(ggtree)
beast_file <- system.file("examples/MCC_FluA_H3.tree", package="ggtree")
beast_tree <- read.beast(beast_file)
p <- ggtree(beast_tree) + geom_tiplab() + xlim(NA, 30)
lv <- get_taxa_name(p)
td$seq1 = factor(td$seq1, levels=rev(lv))
我们排序后,重新画,会发现规律就出来了,其实你如果是单纯要画这个图,用聚类算法来排个序,也能出来挺好的效果,如果减少了颜色的数目,如我前面所说,规律应该更加清楚。当然我这个演示的,没有用户提问的图那么清楚漂亮,你知道原因吗?因为我这颗树的树型是梯子状的,而用户的那颗树是平衡的。可以说不看树,从这个矩阵点图,我们也可以大致猜出来树的形状!
p2 <- ggplot(td, aes(dist, as.numeric(seq1), group=seq2, color=seq2)) +
geom_point() + geom_line() + theme(legend.position="none")
好了,最后是拼图,用cowplot
,要使用align="h"
参数,好对齐,ggtree
画的图,y轴是数字型的,这也是为什么我在生成p2
的时候,要使用as.numeric(seq1)
把y轴变为数字型,大家都是数字型的,好对齐。
p2 <- p2 + theme(axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
axis.line.y=element_blank()) +
xlab(NULL) + ylab(NULL)
require(cowplot)
plot_grid(p, p2, align="h", rel_widths=c(1, .8))