研究メモブログの目次

ultrametric treeを作ってみよう@RStudio

はじめに注意ですが、何かの解析のインプットとして用いる場合、同じUltrametricTreeでも、作り方によってその後の解析の結果の精度が変わるようです(例:GMYC解析)。したがって、先行研究を見て、どんな方法でultrametric treeを作っているのかを確認した方がいいです。

ちなみに、GMYC解析の場合はBEASTBEAST Software - Bayesian Evolutionary Analysis Sampling Trees | BEAST Documentation)で作ったultrametric treeのほうが精度がいいみたいなので、BEASTのほうでも計算が収束する程度のデータサイズならばBEASTでつくるのが無難です。自分は使うデータの系統間距離が大きすぎてBEASTでの計算が収束しなくてしぶしぶRで作る方法を模索した、という感じなので、はじめはBEASTでトライしてみるのをお勧めします。BEASTで作る方法のチュートリアルへは以下のURLから飛べます。

https://botany.natur.cuni.cz/hodnocenidat/Lesson_06_tutorial.pdf

 

ここからはultarametric tree@Rstudio の説明に入っていきます。

 

@Rstudio
参考ページはこちら:Update: 'new' functions for generating starting trees for BEAST or *BEAST in R(https://justinbagley.rbind.io/2015/08/30/update-new-functions-for-generating-starting-trees-for-beast-or-starbeast-in-r/

最初に作業するフォルダを読み込む
Session>Set Working Directory>作業するフォルダの選択

ape package のインストール
Tools>install packages...>
install.packages(ape)
※R が管理者でないと入れないフォルダに保存されていてうまくPackage がインストールできない時は、Rstudio に入るときにアイコン右クリックで「管理者として実行」して入る
参考:https://qiita.com/fujino-fpu/items/543a05044f58cde37eb3

ape package をインストールできたら
library(ape)

fishtree <- read.tree(file='自分のnewick.tree') #拡張子が.tree でない場合変えておく

ultrametric tree を作る方法は以下二通り(これ以外の方法も追記したので、これでダメならこのメモの一番下のほうを見てください)
1.timetree = chronos(fishtree) と打ち込む
2.以下の2行を打ち込む
real_calib1 = makeChronosCalib(fishtree, node = 'root', age.min = 1, age.max = 1, interactive = FALSE, soft.bounds = FALSE)
timetree2 = chronos(fishtree, lambda = 0, model = 'relaxed', calibration = real_calib1)

最後に木の書き出し
write.tree(定義したtree, file="ultrametric.tree")
例:write.tree(timetree2, file="ultrametric2.tree")

choronoPL→こっちのほうがGMYC(種境界解析)の精度が高いらしい
https://www.biostars.org/p/4576/

 


これ以下は11月に自分でやってみて上記の方法が上手くいかなかった時の記録です。
結果から言うと、インプットのMLツリーをRAXMLで作ったら成功しました。

ちなみに、今まではMEGA7(MEGA: Molecular Evolutionary Genetics Analysis software (tmu.ac.jp))で作っていました。
解決方法だけ知りたい人はそういうことなので、以下の記録は読まなくても大丈夫です。
エラーが出た時の試行錯誤方法の一例として一応載せておきます。


timetree を作った時に以下の警告が表示された
Setting initial dates...
Fitting in progress... get a first set of estimates
(Penalised) log-lik = -1.195483e+12
Optimising rates... dates... -1.195483e+12

log-Lik = -1e+100
PHIIC = 2e+100
警告メッセージ:
false convergence (8)

また、timetree2 を作った時には以下の警告が表示された
Setting initial dates...
Fitting in progress... get a first set of estimates
(Penalised) log-lik = -5727357948
Optimising rates... dates... -5727357948
Optimising rates... dates... -5727357945

log-Lik = -5727357945
PHIIC = 11454717055
警告メッセージ:
false convergence (8)

ただ、いずれもtree の書き出しはできた。謎エラー。警告だから大丈夫なのか。
Research Gate でこのエラーが何なのかについてアルバータ大学のCristián Gutiérrez-Ibáñez さんが「この問題は、対数尤度推定値が無限大または0になる傾向があるときに急増すると思います。 これが起こっているかどうかを確認する方法は、ここで説明されています: http://www.mpcm-evolution.org/OPM/Chapter5_OPM/OPM_chap5.pdf
本質的には、ラムダが 0 から 1 のときに PGLS を実行し、対数尤度とラムダをグラフ化します。いくつかのケースでは、ラムダ 0 または 1 に近づくにつれて対数尤度が無限大になることがわかります。」と説明していました。尤度がバカデカくなるか0になるときに出る警告らしいです。

ultrametric tree 作成に悩む人々のページ
https://hcliedtke.github.io/R-scrapheap/be_ultrametric.html
https://www.biostars.org/p/64179/

ultrametrictree 作成代替案
第一候補
1.Mesquite というソフトウェア→ダメだった
使い方
File>Open>Tree(Nexus形式)>新しいウインドウが開く>左手のImported Treesを右クリック>View Trees
ImportedTree タブのTreeタブ>Alter/TransformBranchLength>arbitrarily ultramize

第二候補 別の方法@Rstudioで作るultrametric tree

##や#は上記のお悩みな人々のサイトに書かれていたコメントをそのまま書いています。

>の後ろがRで実際に打ち込む部分になります。
>library(splits)
>library(adephylo)
>library(phytools)
>tree <- read.tree('nwk.tree')
>is.ultrametric(tree)
[1] FALSE
ここではまだultrametric ではないようです
>tip.heights<-distRoot(tree)
50 件以上の警告がありました (最初の 50 個の警告を見るには warnings() を使って下さい)
警告は全部これ→In getEdge(x, e) : Not all nodes are descendants in this tree

>heights.summary<-table(tip.heights)
>options(digits=最大だったものの半分を指定?)

# set to maximum allowed digits
>real.tree.height<-as.numeric(names(which.max(heights.summary)))
>over.under<-tip.heights-real.tree.height

## we can now paint the branches that were problematic using phytools()
>painted.tree<-paintBranches(tree,which(round(over.under,5)!=0),"2")

## extract all terminal edges for tips that do not have the final height we want:
>tip.ids <- tree$edge[, 2] <= Ntip(tree)
>terminal.edges <- tree$edge.length[tip.ids]

## add/subtract the extra length from the terminal branches
>corrected.terminal.edges<-terminal.edges-over.under

## change the termnial edges in the phylo object
>tree$edge.length[tip.ids]<-corrected.terminal.edges
>plot(tree, show.tip.label = F, main="its ultrametric!!")

最後に出来上がったツリーがウルトラメトリックなのか確認
>is.ultrametric(tree)
[1] TRUE

木の書き出し
>write.tree(定義したtree, file="ultrametric.tree")

そして早速GMYC
write.tree(tree, file="ultrametric.tree")
> library(splits)
> tr <- read.tree('ultrametric.tree')
> gmyc.result <- gmyc(tr)
des[j, ] <- as.integer(tr$edge[, 2])[as.integer(tr$edge[, 1]) == でエラー:
置き換えるべき項目数が、置き換える数の倍数ではありませんでした
ええ…困りました
とりあえず、ここまできたらインプットファイルが悪い気がしてきた
MEGA から最初のMLツリーのファイルを開く
Blanch LengthとNode 両方保存(今まで枝長しか保存していなかった)
is.ultrametric(tree)
.is.ultrametric_ape(phy, tol, option, length(phy$tip.label)) でエラー:
the tree has no branch lengths
枝長がないって言われました
なんで~~~~~!!!!!

 

RAXMLでインプット作成やってみる→成功した

 

おしまい。