背景
大学院の系統学の講義で系統のソーティング(lineage sorting)を扱いました. その時出された課題が, 種分化が生じてからの時間と種分化時の集団サイズが系統のソーティングにどのような影響を及ぼすのかをシミュレーションによって調べることでした. 系統樹の作成にはMCcoalを使っていて, ここではパラメーターだけ指定すれば自動的に系統樹のファイルを出力してくれます.
ここまでは良いのですが, 生成された系統樹の中で各種がどのような系統関係になっているかは目でチェックしなければならず非常に面倒. そこで, 系統樹のファイルから系統群を自動的に判別してくれるプログラムを考えました.
アルゴリズムを考える
MCcoalで生成されるファイルはNewickフォーマットと呼ばれる括弧とコンマを使って書かれたもので, このファイルをTree ViewやFigTreeなどのソフトウェアで開けば系統樹を描いてくれます. そこで, このNewickフォーマットをうまく処理して系統群の分類を行うことを考えていきます.今回は, 普段使い慣れているRを用いてプログラムを組んでいきたいと思います.
Newickフォーマットではある要素がいくつの括弧で囲まれているのかが枝分かれの数になります. 例えば, (((A,A),A),(((B,B),(B,C)),C))では一番左のAは3つの括弧に囲まれているため3回枝分かれをしており, 一番右のCは2つの括弧に囲まれているため2回枝分かれをしていることになります. また, 括弧がより内側に存在しているほど最近枝分かれしたことになります.
この辺をうまく利用して系統群の分類を行います. まず, それぞれの要素に対して最も近い共通祖先(MRCA)を探します. これは, ある要素が全て含まれる(例えば, 系統樹に含まれる全てのAが出現するような)最も内側の括弧を探すことで実現できます. 仮にAの場合, これをMRCA(A,A)とおきます. 次に, 異なる要素同士の全ての組み合わせに対してMRCAを探します. 例えば,
(((A,A),A),(((B,B),(B,C)),C))ではAとB, AとC, BとCの組み合わせになります. これは, それぞれの組み合わせの要素が全て含まれる最も内側の括弧を探すことで実現できます. 仮にAとBの組み合わせの場合, これをMRCA(A,B)とおきます. ここで,
MRCA(A,A)-MRCA(A,B)>0かつMRCA(A,A)-MRCA(A,C)>0になるならば, Aは単系統群と判断できます. そして単系統群にならない場合, その要素のMRCAが含む他の要素(Bに注目した場合のC)が全て単系統群である場合, 側系統群と判断できます. そして, 残ったものは多系統群と判断します.
実装
R
###系統群
#必要なパッケージ
library(magrittr)
###系統群を判定する関数. 引数はnewickフォーマットのデータファイル
judge.taxa = function(TreeData) {
original.trees = scan(TreeData, what = character(), sep = "\n", blank.lines.skip = F) #newickフォーマットのデータを読み込み
tree = rep(0, length(original.trees)) #それぞれの系統樹のデータを格納するための変数
result = list() #結果を格納するための変数
for (j in 1:length(tree)) {
#####括弧の対応関係と包含関係を調べる#####
tree[j] = gsub("[0-9\\.:; ]", "", original.trees[j])
f.table = gsub("[\\(\\)]", "", tree[j]) %>% strsplit(., ",") %>% table()
kakko = gsub("[^\\(\\)]", "", tree[j])
parenth = NULL
for (i in 1:nchar(kakko)) {
parenth[i] = substr(kakko, i, i)
}
bgkakko = grep("\\(", parenth)
edkakko = NULL
#はじめ括弧に対応する終わり括弧が何文字目にあるかを検索
c = 0
for (i in 1:length(bgkakko)) {
for (k in bgkakko[i]:length(parenth)) {
if (parenth[k] == "(") {
c = c - 1
} else {
c = c + 1
}
if (c == 0) {
edkakko[i] = k
break
}
}
}
n = nchar(kakko)
bgparenth = gregexpr("\\(", tree[j])[[1]][1:length(gregexpr("\\(", tree[j])[[1]])]
edparenth = gregexpr("\\)", tree[j])[[1]][1:length(gregexpr("\\)", tree[j])[[1]])]
corrkakko = data.frame(bg = rank(bgkakko), ed = rank(edkakko))
phasekakko = data.frame(bg = bgkakko, ed = edkakko, delta = edkakko - bgkakko)
corrparenth = data.frame(bg = bgparenth, ed = rep(0), delta = edkakko - bgkakko)
for (i in 1:(n/2)) {
corrparenth[i, 2] = edparenth[which(rank(edparenth[1:(n/2)]) == corrkakko[corrkakko$bg == i, ]$ed)]
}
order.parenth = corrparenth[order(phasekakko$delta), ]
#####全ての括弧の組み合わせに対して各要素の頻度と括弧の包含関係をまとめる#####
dgroup = data.frame()
for (i in 1:(nchar(kakko)/2)) {
test <- substr(tree[j], order.parenth[i, 1], order.parenth[i, 2]) %>%
gsub("[\\(\\)]", "", .) %>%
strsplit(., ",")
dgroup <- factor(test[[1]], levels = names(f.table)) %>%
table() %>%
as.data.frame() %>%
cbind(., phase = order.parenth$delta[i], id = i) %>%
rbind(dgroup, .)
}
names(dgroup)[1] = "group"
mat = matrix(0, length(names(f.table)), length(names(f.table))) #各要素どうしの距離を格納するための行列を用意
rownames(mat) = names(f.table)
colnames(mat) = names(f.table)
comb.homo = data.frame(r = 0, c = 0) #同じ要素同士の組み合わせ
for (i in 1:ncol(combn(names(f.table), 1))) {
comb.homo[i, 1] = combn(names(f.table), 1)[, i]
comb.homo[i, 2] = combn(names(f.table), 1)[, i]
}
comb.hetero = data.frame(r = 0, c = 0) #異なる二組の要素同士の組み合わせ
for (i in 1:ncol(combn(names(f.table), 2))) {
comb.hetero[i, 1] = combn(names(f.table), 2)[1, i]
comb.hetero[i, 2] = combn(names(f.table), 2)[2, i]
}
for (i in 1:nrow(comb.homo)) { #データに含まれるある要素のすべてが表れる括弧の最小値(単系統の判別用)
mat[comb.homo[i, 1], comb.homo[i, 2]] = min(dgroup[dgroup$group ==
comb.homo[i, 1], 3][dgroup[dgroup$group ==
comb.homo[i, 1], 2] - f.table[comb.homo[i, 1]] == 0])
}
for (i in 1:nrow(comb.hetero)) { #異なる二組の要素が枝分かれした最短の括弧(MRCA)を求める
for (k in 1:(nchar(kakko)/2)) {
if (dgroup[dgroup$id == k & dgroup$group == comb.hetero[i, 1], 2] * dgroup[dgroup$id == k &
dgroup$group == comb.hetero[i, 2], 2] == 0) {
} else {
mat[comb.hetero[i, 1], comb.hetero[i, 2]] = dgroup[dgroup$id == k, ]$phase[1]
break
}
}
}
result[[j]] = data.frame(group = names(f.table), taxa = -1, stringsAsFactors = FALSE) #それぞれの要素を-1で初期化
judge = NULL #判定用変数
pairset = data.frame() #単系統群にならなかった場合, その組み合わせを格納するための変数
#####単系統群になるかどうか判定#####
for (i in 1:nrow(comb.homo)) {
if (f.table[comb.homo[i, 1]] == 1) {
result[[j]][result[[j]]$group == comb.homo[i, 1], ]$taxa = 1 #単系統群の場合resultに1を代入
break
}
for (k in 1:(length(names(f.table)) - 1)) {
m.r = comb.hetero[comb.hetero$r == comb.homo[i, 1] | comb.hetero$c == comb.homo[i, 1], ]$r[k]
m.c = comb.hetero[comb.hetero$r == comb.homo[i, 1] | comb.hetero$c == comb.homo[i, 1], ]$c[k]
judge[k] = mat[m.r, m.c] - mat[comb.homo[i, 1], comb.homo[i, 2]]
if (judge[k] <= 0) { #単系統群でない場合
pairset = rbind(pairset, data.frame(comb.homo[i, 1], c(m.r, m.c)[c(m.r, m.c) != comb.homo[i, 1]]))
}
}
if (sum(judge > 0) == length(judge)) {
result[[j]][result[[j]]$group == comb.homo[i, 1], ]$taxa = 1
}
}
#####側系統群か多系統群かを判定#####
if (length(pairset)) {
colnames(pairset) = c("me", "you")
polyorpara = NULL
for (i in result[[j]][result[[j]]$taxa == -1, 1]) {
for (k in 1:nrow(pairset[pairset$me == i,])) {
polyorpara[k] = result[[j]][result[[j]]$group == pairset[pairset$me == i, 2][1], 2]
}
if (sum(polyorpara) == length(pairset$me == i)) {
result[[j]][result[[j]]$group == i, 2] = 0 #側系統群の場合0を代入
}
}
}
}
return(result)
}
コメントをお書きください