・修論手伝(標本整理)
・標本整理
午後
・データ解析
ggtreeを使った系統樹描画。OTU名の変更も簡単だった。昨日、うまく行かなかったのは、単純なデータの勘違いだった。
流れは、下記方法で「nwk」データを「tree」に読み込んだ場合、OTU名は「tree$tip.label」に入っているので、「tree$tip.label <- ***」で好きな名前に変更すれば良い。
なぜ、そんなことをする必要があったのかというと、KAKUSAN4を使って同じ配列データを一つのハプロタイプにまとめたので、サンプル数が多いハプロタイプは名前が非常に長かったため。例えば、こんなのとか。
"Amami-38_St-6_109_Amami-38_St-6_110_Uke-1_St-6_119_Uke-1_St-6_120_Uke-1_St-6_121_Miyazaki-1_St-6_132_Amami-6_St-6_160_Amami-6_St-6_163"
これは、各シーケンスデータに「locality_haplotptype id_sample id」が書かれていて、それが9サンプルが合わさったため。
でも、結果オーライなのだが、"_"で区切った2番目にはハプロタイプ番号(St-6)が記されているので、それをOTUに使えばOK!
テキストの一部を切り出すのは「library(stringr)」の「str_split( )」で出来る!
library(ggtree) library(treeio) # nwkデータの取り込み tree <- read.newick("ML_stimpsonii.nwk") # ラベル (OTU名) の確認 tree$tip.label長い、、、
# 上記のデータを対象に、"_"で区切り、その2番目を「tree$tip.label」に入れる library(stringr) # "_"ごとに区切って2番目を取り出す a<-str_split(tree$tip.label, pattern = "_", simplify = TRUE) tree$tip.label <- a[,2] # ラベル (OTU名) の確認 tree$tip.label短くなった。
あとは、グループごとに色付けして出力。
# OTU名に基づきグループを作成する cls <- list(Kumejima=c("St-55", "St-56"), Okinawa=c("St-57", "St-58"・・・)) # nwkデータにグループ名を紐つけ tree <- groupOTU(tree, cls) # 描画 library("colorspace") #色に必要? ggtree(tree, aes(color=group), layout='circular') + scale_color_manual(values=c("black", "azure4", "red", "magenta", "royal blue", "#00552e")) + geom_tiplab()