15、使用ggtree实现进化树的可视化和注释（转载）

本文作者：余光创，目前就读于香港大学公共卫生系，开发过多个R/Bioconductor包，包括ChIPseeker, clusterProfiler, DOSE,ggtree，GOSemSim和ReactomePA。

 

进化树看起来和层次聚类很像。有必要解释一下两者的一些区别。

 

层次聚类的侧重点在于分类，把距离近的聚在一起。而进化树的构建虽然也可以说是一个聚类过程，但侧重点在于推测进化关系和进化距离(evolutionary distance)。

 

层次聚类的输入是距离，比如euclidean或manhattan距离。把距离近的聚在一起。而进化树推断是从生物序列（DNA或氨基酸）的比对开始。最简单的方法是计算一下序列中不匹配的数目，称之为hamming distance（通常用序列长度做归一化），使用距离当然也可以应用层次聚类的方法。进化树的构建最简单的方法是非加权配对平均法（Unweighted Pair Group Method with Arithmetic Mean, UPGMA），这其实是使用average linkage的层次聚类。这种方法在进化树推断上现在基本没人用。更为常用的是邻接法（neighbor joining），两个节点距离其它节点都比较远，而这两个节点又比较近，它们就是neighbor，可以看出neighbor不一定是距离最近的两个节点。真正做进化的人，基本不用这些基于距离的方法。现在主流的方法是最大似然法(Maximum likelihood, ML)，通过进化模型（evolutionary model)估计拓朴结构和分支长度，估计的结果具有最高的概率能够产生观测数据（多序列比对）。另外还有最大简约法和贝叶斯推断等方法用于构建进化树。

Newick是最常用的存储进化树的文件格式，如上面这个树，拓朴结构用newick格式可以表示为：

(B,(A,C,E),D);

括号最外层是根节点，它有三个子节点，B, (A,C,E)和D，而节点(A,C,E)也有三个子节点A，C和E。

 

加上分支长度，使用 : 来分隔：

(B:6.0,(A:5.0,C:3.0,E:4.0):5.0,D:11.0);

比如A:5.0代表的是A与其父节点的距离是5.0。

 

内部节点也可以有label，写在相应的括号外面，如下所示：

(B:6.0,(A:5.0,C:3.0,E:4.0)Ancestor1:5.0,D:11.0);

这是最为广泛支持的文件格式，很多进化树可视软件只支持newick格式。

 

ggtree的开发源自于我需要在树上做注释，发现并没有软件可以很容易地实现，通常情况下我们把统计信息加到节点的label上来展示，比如CodeML的dN/dS分析，输出文件里就给用户准备了newick树文本，把dN/dS (ω) 加于节点label之上:

 

codeml_file <- system.file("extdata/PAML_Codeml/mlc", package="ggtree")tree_text <- readLines(codeml_file)[375:376]tree_text

## [1] "w ratios as labels for TreeView:"                            ## [2] "(K #0.0224 , N #0.0095 , (D #0.0385 , (L #0.0001 , (J #0.0457 , (G #0.1621 , ((C #0.0461 , (E #0.0641 , O #0.0538 ) #0.0001 ) #0.0395 , (H #0.1028 , (I #0.0001 , (B #0.0001 , (A #0.0646 , (F #0.2980 , M #0.0738 ) #0.0453 ) #0.0863 ) #1.5591 ) #0.0001 ) #0.0001 ) #0.0549 ) #0.0419 ) #0.0001 ) #0.0964 ) #0.0129 );"

这种做法只能展示一元信息，而且修改节点lab