Newick形式で系統樹の樹形が同じか確かめる

複数のNewick形式の系統樹から樹形が同じものを抽出する必要に駆られたので備忘録を兼ねたまとめです。

当たり前ですが、

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

の2つのNewick形式は同じ樹形を表しており、

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

の2つのNewick形式は異なる樹形を表しています。

問題は、文字列としては異なる2つのNewick形式で樹形が同一かを判定しなければならない点にあります。枝長を無視して判定したいような場合には、さらに処理がややこしくなります。この辺りで頭が痛くなってきました。

という訳で既存のコードを探します。

調べた結果、Rのapeパッケージにシンプルに樹形比較をしてくれる関数があるようなのでこれを利用します。

インストール方法はごく普通です。

install.packages("ape")

使い方もシンプルです。同じくapeのread.tree関数でNewick形式を読み込んで、all.equal.phylo関数で樹形が同じかを判定します。今回は文字列をそのまま読み込ませていますが、ファイルのパスを渡すこともできます。

> t1 <- read.tree(text='((A,B),(C,D));')
> t2 <- read.tree(text='((C,D),(A,B));')
> all.equal.phylo(t1, t2)
[1] TRUE

> t1 <- read.tree(text='((A,B),(C,D));')
> t3 <- read.tree(text='((A,C),(B,D));')
> all.equal.phylo(t1, t3)
[1] FALSE

枝長を無視して樹形比較を行いたい場合は、use.edge.length引数をFALSEにします。デフォルトではTRUEになっています。

以上をまとめて今回作った、Newick形式の系統樹ファイルを2つ受け取って、樹形が同じかを表示するプログラムを下に載せました。

#!/usr/bin/Rscript

library('ape')
args <- commandArgs(trailingOnly = TRUE)
nwk1 <- args[1]
nwk2 <- args[2]
tree1 <- read.tree(nwk1)
tree2 <- read.tree(nwk2)
is.same.topology <- all.equal.phylo(tree1, tree2, use.edge.length=FALSE)

if (is.same.topology == TRUE) {
    cat("identical\n")
} else {
    cat("different\n")
}

本当はPythonで作りたかったんですけどね。Rも勉強しないといけないってことですかね。