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( )」で出来る!
  1. library(ggtree)
  2. library(treeio)
  3. # nwkデータの取り込み
  4. tree <- read.newick("ML_stimpsonii.nwk")
  5. # ラベル (OTU名) の確認
  6. tree$tip.label
長い、、、

  1. # 上記のデータを対象に、"_"で区切り、その2番目を「tree$tip.label」に入れる
  2. library(stringr)
  3. # "_"ごとに区切って2番目を取り出す
  4. a<-str_split(tree$tip.label, pattern = "_", simplify = TRUE)
  5. tree$tip.label <- a[,2]
  6. # ラベル (OTU名) の確認
  7. tree$tip.label
短くなった。

あとは、グループごとに色付けして出力。
  1. # OTU名に基づきグループを作成する
  2. cls <- list(Kumejima=c("St-55", "St-56"), Okinawa=c("St-57", "St-58"・・・))
  3. # nwkデータにグループ名を紐つけ
  4. tree <- groupOTU(tree, cls)
  5. # 描画
  6. library("colorspace") #色に必要?
  7. ggtree(tree, aes(color=group), layout='circular') +
  8. scale_color_manual(values=c("black", "azure4", "red", "magenta", "royal blue", "#00552e")) +
  9. geom_tiplab()