2020年3月26日木曜日

OTU名の変更

午前
・修論手伝(標本整理)
・標本整理

午後
・データ解析

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()