ラベル R の投稿を表示しています。 すべての投稿を表示
ラベル R の投稿を表示しています。 すべての投稿を表示

2024年10月28日月曜日

geom_strip( )でイタリック

午前
・データ解析
・Rゼミ

午後
・野外調査
・データ整理

系統樹の描画には、library(ggtree)を使っている。

Cladeを示すOTUのannotationは、geom_strip( )で行っているが、学名をイタリックにすることに苦戦したのでメモ。

geom_text( )をベースとしている(?)ようなので、同じ方法でやれば良いみたい。

微妙なコマンドミスがあったようでエラーが出まくって諦めかけたけど、どうにか正解に辿りつけた。

本来はこんな感じで書く(必要な部分のみ)。
  geom_strip(label = "B. ocellatus") +
  geom_strip(label = "B. cf. purpura"))
すると学名が立体になる。
イタリックにするにはitalic( )を使うが、paste( )とparse=Tを合わせて使用する。また、テキスト部分はc(' ')で括る必要がある。

途中に立字を入れる場合、面倒だけど「italic( ), 立字, italic ( )」とする。
  geom_strip(label = c('paste(italic("B. ocellatus"))'), parse=T) +
  geom_strip(label = c('paste(italic("B. "), "cf. ", italic("purpura"))'), parse=T)
参考にしたサイト
・How to make R and p italic in geom_text()
・13 Gallery of Reproducible Examples

2024年10月24日木曜日

read.xlsx( )

午前
・卒論手伝い
・データ解析
・個人ゼミ

午後
・データ解析
・個人ゼミ
・研究相談

Rでxlsxファイルを読み込んだ時に生じた不具合メモ。

これまではcsvをread.csv( )やread_csv( )で読み込んでいた。

ただ、データがどうしても一つのファイルにまとめられない場合、csvだと複数のファイルとして保存せざるを得ない。2〜3つなら良いけど、それよりも多い場合は、xlsxのシートの方が整理しやすい。

ということで、最近はxlsxから直接読み込むことも多くなった。xlsxファイルを読み込みは、library(openxlsx)のread.xlsx( )を使っている。

使い方はread.csv( )と同じ。
library(openxlsx)
read.xlsx(data)
ただ、この方法で読み込んだデータをdplyr::filter( )で処理をすると「! Can't transform a data frame with duplicate names.」とエラーが出ることがある。

同じ名前が複数あったらダメなの?とか考えたけど、そいうことでもないもみたい。

このサイトによると、data.frame( )を介せば解消されるそうだ。

私も同様に改善された。ということで、最近は
library(openxlsx)
data.frame(read.xlsx(data))
として読み込んでいる。

参考にしたサイト
・can't transform a data frame with duplicate names mutate() function in dplyr

2024年9月24日火曜日

NMDSに環境データをバイプロット

午前
・卒論手伝い(野外調査)
・データ解析

午後
・データ解析
・論文書き

NMDSに環境データをバイプロットする際にハマったのでメモ。

metaMDS( )でNMDSを実施し、envfit( )で関連のある環境データを抽出し、その座標データを「vector$arrows」で抽出してggpolt( )で作図すると、ordiplot( )とplot( )で作図した場合と環境データの表示が異なることが分かった。

下図は、Rによる群集組成の解析のデータを使わせて頂き、ordiplot( )とplot( )で作成した図。
環境データのベクトルよりも上部に点がある。そして、「vector$arrows」で抽出した座標データを用いてggplot( )で作成したのが下図。
環境データのベクトルが点よりも上部まで伸びている。多分、長さは全て同じなのだろう。

朧げな理解しかできていないが、「vector$arrows」は関係性の(座標上での)方向性のみの情報を有しており、関係性の強さが考慮されてないみたい。下記サイトによると、この座標データは「direction cosines」と記されている。

関係性の強さを加えて補正するには、r2のルートを書ければ良いそうだ。

もしくは、scores( )を使って座標を抽出すれば良い。下図はscores( )で抽出したデータを用いて作成した図。ordiplot( )とplot( )と同じ図になった!
参考にしたサイト
Plotting envfit vectors (vegan package) in ggplot2
Rによる群集組成の解析

2024年9月19日木曜日

四分位数

午前
・ハンミョウ調査

午後
・資料作成
・論文書き

論文用に箱ヒゲ図をチェックしている最中にビックリしたこと。

例えば、「4,5,11,23,27,44,75,87,90,100」で箱ヒゲ図を描き、プロットを載せると下記のような図になる。
ちなみにコマンドは以下の通り。
a <- c(4,5,11,23,27,44,75,87,90,100)
b <- rep("B", 10)
ab <- data.frame(a,b)

ggplot(data=ab, aes(x=b, y=a)) +
  geom_boxplot() +
  geom_point(size=5, col="red")
データ数が10個なので、中央値は下(or上)から5番目と6番目の値の平均値となる。

上記のデータでは「27, 44」の平均値で「35.5」となる。図から正確な数値は読み取れないが、中央値の線が点と点の真ん中にあることから、この考えが正しいことは分かる。

で、問題は第1四分位数と第3四分位数である。第1四分位数は下から数えて25%、第3四分位数は75%となるデータである。

データが偶数の場合は、半分に分けて、それぞれの真ん中の値が、それらに当たる。
この考えに基づくと、データが偶数で「箱ヒゲ図+プロット」を描くと、箱の下線と上線にプロットが重なるはずである、、、が重なっていない。

昨夜、投稿用の図を眺めていて、このことに気づいてヒヤっとした。何かミスがある〜、と。

ただ、summary()で四分位数を算出すると、図の値と一致している感じ。私の認識が間違っているのか?
> summary(ab)
       a              b            
 Min.   :  4.0   Length:10         
 1st Qu.: 14.0   Class :character  
 Median : 35.5   Mode  :character  
 Mean   : 46.6                     
 3rd Qu.: 84.0                     
 Max.   :100.0    
答えを教えてくれるサイトがありました。「四分位数の違いについて~みんなちがってみんないい~」さん、ありがとうございます。

四分位数には、複数の考え方があるそうです。

私の認識は「文部科学省推奨の定義」だった。文部科学省の指導要領の下で勉強してきた成果か。

で、サイトに「ExcelやRで使われる、“計算機タイプの四分位数”」と書かれていました。

概念は分かりにくいのだが、大量のデータに対して効率よく計算できる仕組みだそうです。計算方法は下記サイトに書かれています。計算方法自体は簡単です。

参考にしたサイト
四分位数の違いについて~みんなちがってみんないい~

2024年9月5日木曜日

draw.shape()で海面の色を指定

午前
・シーケンス準備
・PCR実験

午後
・論文書き(修論論文化)

死蔵になりそうな修論の論文化を強引に進める。

作図は随分前に終わっていたのだが、気になっていた点があった。

日本地図に円グラフを描くために、
ggplot() +
  geom_sf(地図データ) +
  geom_scatterpie(円グラフデータ)
を使っていたのだが、円グラフがやや縦長になってしまう。

緯度経度を同じ数値として扱ってしまうのかな〜と考えているのだが、改善方法が思いつかない。
このままで良いか、いやダメだ、と自問自答を繰り返した結果、作り直すことにした。

改善方法が思いつかないので、私が使っているもう一つの円グラフの作図方法「mapplots::draw.pie()」で対応した。

「mapplots::draw.pie()」は独特なデータ形式が使われるが、「mapplots::make.xyz()」で簡単に作成できる。

地図は、shpファイルを「mapplots::draw.shape()」で描く必要がある?

ここまでは他の論文のコピペなのですぐに終わったが、海面の色指定に苦戦してしまったのでメモ。海面をwhitesmoke、陸地をblackにする場合は下記のようにする。
jpn <- read.shapefile(PATH) #shpファイルの読み込み
xlim <- c(122.5, 147.5) #X軸の指定
ylim <- c(22.5, 47.5) #Y軸の指定

#ベースの指定;bgで海面の色を指定!
basemap(xlim, ylim, axes=F, xlab="", ylab="", bg="whitesmoke")

#地図=陸地の指定;colで陸地の色を指定!
draw.shape(jpn, col="black")
綺麗な円になった。
ただ、円グラフの位置が微妙ずれてしまい、微調整に2〜3時間かかった。

2023年12月8日金曜日

reveal.jsでMathJax

午前
・授業準備
・授業

午後
・卒論手伝い(SEM撮影)
・PCR実験
・授業準備
・電気泳動

授業のスライドは、Rmarkdownでreveal.jsを使って作成している。

Rmarkdownでの数式の表示には、MathJaxを使っているのだが、reveal.jsで作成したhtmlでは上手く表示できなかった。render_book( )では問題ない。

ローカルではJavaScriptが上手く動いていないのだと思ってサーバーに上げたがダメだった。

ネット検索の結果、MathJaxのコマンドをローカルにダウンロードして、それを指定する方法があった。

MathJaxのダウンロード方法は、ここで紹介されている。gitを使った。

2023年10月2日月曜日

円グラフと箱ひげ図

午前
・標本整理

午後
・授業準備

授業の準備で円グラフを作成したが、どハマり。Rで円グラフ作るのは辛い。でも、図を描く仕組みがなんとなく分かった。
ついでに論文用の図の修正も行なった。

箱ひげ図に色と模様をつけるだけなのだが、これも結構時間がかかった。

色と模様、それぞれで凡例が出てしまうので、その制御に手こずった。

2023年9月22日金曜日

ggRGB()

午前
・データ整理

午後
・論文書き

野外調査の予定だったが雨で延期。

ということで、昨日に引き続きデスクワークに励む。

論文、、、報告と言った方が正確か、、、を書く。一つは大体できたので、二つ目に手を出したところでハマった。

ggplot2を使って航空写真(geotif)を描く方法で、、、。

前にも同じようなことでハマったな〜、でも、こんな経験は何度もあるし、違うやつだったかなと、真面目に過去の記録を確認せずにネット検索したのが間違いだった。

ggolot2でgeotifをカラーで出力するにはRStoolbox::ggRGB()が便利なのだが、RStoolboxのinstallができなくなった。

で、同じことを7月に悩んでgithubのスクリプトを全て読み込めば実行できることを見つけていた、、、のだが、このチェックを怠った。

半日かけて色々と試して同じ結論に辿りついた。

ただ、本日は少しだけ成長して、「internalFunctions.R」と「ggRGB.R」のみ読み込めばOKであることが分かった。
ラグビーW杯、今週、日本は休みの週だったのか、、、。

1週間に1試合で、しかも、予選リーグの総当たり戦で(5チームなので)1回休みがあるという独特のシステム。

最初 or 最後に休みのチームと、中間に休みがあるチームではかなり不公平な気がする。

2023年7月26日水曜日

geom_sf

午前
・資料作成
・生態学ゼミ

午後
・生態学ゼミ
・卒論手伝い(野外調査)
・資料作成

航空写真(geotiff)をggplotで使用するために、RStoolbox::ggRGBを使っていたがエラーがでるようになった。

原因は良くわからない。version4.3が入っているmac mini(M1チップ)とsurfaceでは使えないが、version4.2が入っているmac book airでは動いているので、versionの影響を受けるのかも、、、。

terrainr::geom_spatial_rgbやplotRGBとかも検討したが、色々と問題が生じたので、再度、ggRGBを検討。

ひとまず、githubのスクリプトを全て読み込めばmac mini(M1チップ)でもggRGBは実行できることは分かった。

これで解決と思いきや、geom_sf(size=)がgeom_sf(linewidth=)に変更になっていることを知らず、結構な時間悩まされた。

昨日の、井上尚弥 vs フルトン戦は面白かった。

2023年5月9日火曜日

OTU名のみ色を変更

午前
・授業
・授業準備

午後
・授業準備
・データ解析

魔の火曜日だが、週末の準備を進めなければ。

以前から、絶対にやる方法があるはずだと思いつつ、これまで出来ていなかった、ggtree()を用いて系統樹を描く際にOTU名のみ色を変更する方法が、やっと分かった。
groupOTU()を使う方法や、もっとカッコいい方法もあったけど、理解できたのはこれ。
#読み込み方で形式が異なるの注意
tree <- read.beast("AAA.nex") 

#サンプル数が124の場合
#OTU名の列(label)と、blackだけの列(color)を作る;blackに意味はない
d <- data.frame(label=tree@phylo$tip.label, color=c(rep("black", 124)))

#変更したいOTU名(下記だと96、97)のcolor列のblackをmagentaに変更する
#厳密な意味で色指定ではなくグループ指定に過ぎない
dd <- d %>%
  mutate(color =case_when(
    label == 96 ~ "magenta",
    label == 97 ~ "magenta"))
    
 #上記のコマンドを実行すると指定しない行のcolor列にNAが入る
 #NAをblackに変更する
ddd <- dd %>% replace_na(list(color= "black"))

#作図
ggtree(tree, size=2.5) %<+% ddd +
  geom_tiplab(fontface='bold', size=8, aes(col=color), geom="text") +
  scale_color_manual(values=c("black", "magenta")) +  #実際の色はここで指定
  theme(legend.position = "none")
「%<+% 」でphyloデータにデータフレームを加えることができるみたい。

それを利用することで、ggplotと同じように扱えるようになる。

コマンドを流れを整理すると、、、

まず、OTU名とグループ名を含むデータフレームを作成する。

その後、case_whenを使って、変更したOTU名のグループ列の名称を変更する。上記は色を名称に使っているけど、グループにすぎないので適当な名前で良い。

(カッコ悪いやり方だけど)case_whenを使用すると、変更の指定をしてない行がNAになってしまう。そこで、NAを適当な名称に置換する。

geom_tiplab()のaes()でグループ名の列をcol=で指定して、scale_color_manual()で色を指定する。実際の色の指定はここで行う。グループ名のアルファベット順に色を指定する。

2023年3月13日月曜日

guides("none")

午前
・データ整理

午後
・ポスター作成

同じような日々が続く。

ggplot()でscale_***_***()を使用すると凡例の制御が難しく困っていたのだが、guides("none")で消せることが分かった。
scale_alpha_manual(values = plot_alpha) +
guides(alpha = "none") + #特定の凡例を消す
と、scale_***_***()の指定(最初の***)に合わせる必要がある。

2023年3月10日金曜日

矢印を描く

午前
・データ整理

午後
・ポスター作成

geom_segment()で矢印を描くとalphaがうまく機能しない?

下記のようにannotate()を使うと描くことができた。
annotate("segment",
	x=-30, y=-5, xend=-25, yend=3,
    arrow=arrow(length = unit(1, "cm")),
    lineend = 'round', size=10, color="blue", alpha=0.2) 

2023年3月9日木曜日

drposter

午前
・データ整理

午後
・ポスター作成

久しぶりのRどハマり、、、久しぶりでもないか。

今月末に使うポスターを作成するのに、一昨年の日本土壌動物学会で利用したdrposterを用いることにした。

コマンドの良いところで、コピペですぐ終わるはずだったのが、どハマりした。

一昨年、作成したポスターはこんな感じ。
同じコマンドをKnitすると、こんな感じになってしまう。この前にファイルの準備が不足して、色々とエラーも出ていた。
様式がおかしいのでcssがうまく機能していないみたいだな、とこの時考えたのだが、一昨年は実行できたので、その点について詳しく確認をせずにネットに頼ったのが良くなかった。

結局、解決方法がわからず、他のパッケージを試したりと時間を費やした。しかし、他パッケージはポスターの版組に自由度が少なかったりと良いものが見つからない。

で、やっぱりcssでは?と疑い始める。一昨年のRmdファイルをみると、
output:
  drposter::drposter_poster:
    self_contained: false
    css: custom.css
 
となっており、custom.css、が使われていることが分かる。

で、custom.cssを確認すると、こんな感じ。今思えば何も指定されていないのだが、この時は「以前はこれで実行できたし、、、」と思ってしまった。
しばらくウダウダしたのち、「dropster_files」というフォルダの中に、たくさんcssがあることに気付いた。
これがうまく指定できていないんだな、と分かったが、ここでも「以前はこれで実行できたし、、、」と思ってしまった。

さらに、ウダウダしたのち、ここにあるcssを全部読み込ませることを考えた。custom.cssに以下を記入。
@import url("drposter_files/drposter.css");
@import url("drposter_files/examples.css");
@import url("drposter_files/fonts.css");
@import url("drposter_files/grid.css");
@import url("drposter_files/header.css");
@import url("drposter_files/load_fonts.css");
@import url("drposter_files/page.css");
@import url("drposter_files/sections.css");
@import url("drposter_files/short_tags.css");
@import url("drposter_files/theme.css");
そうしたら無事に動いた!
PCを変更した影響なのか、一昨年は何か別に設定をしていたのか、なぜ、このような状態になったのかは良く分からないが、深刻な状況ではなくて良かった。

drposter.cssに他のcssが引用されているので、上記のようにcustom.cssは全てを引用する必要はなく、drposter.cssだけを@importすれば良いと思う。もしくは、Rmdファイルでdrposter.cssを指定すれば良いのかも。

2023年2月21日火曜日

labs(caption = )

午前
・修論手伝い(原稿添削)

午後
・標本整理
・論文下記

(Ma)、を書くだけに2時間ほど悪戦苦闘してしまった。
library(magick) を使って画像として処理すれば一発という思いだったのだが、エラーが出続け、結局、改善できなかった。

「geom_text」もY軸は負はだめだった。

「annotate」は負の値をとれるが、目盛りの下には描くことができなかった(目盛りがどんどん下がる)。

で、「labs(caption = )」と「 plot.caption = element_text(vjust=, hjust=), 」で無理やり作成した。

「labs()」で軸名以外を指定できることを始めて知った。

下記を参考にさせてもらいました。

【R】ggplotでの軸ラベルやタイトルの設定方法

2023年2月16日木曜日

超入門! Rでできるビジュアル統計学の感想。

午前
・覚えていない、、、

午後
・会議
・卒論手伝い(原稿添削)

3年生向けの統計ゼミで使用した「超入門! Rでできるビジュアル統計学」の感想。

コンセプトは面白く、単に「グラフを描く」と言っても色々なことを考えなければならない、ことを学ぶには有用に感じた。

個人的なクセなんだろうが、コマンドの書き方で2点気になるところがあった。

一つ目は、オブジェクトへの代入を最後に右向きで書く点。以前、見せてもらったコマンドでも、そのように書かれていることがあったので、よく使われているのかも知れないが、多くの本では最初に左向きに書いていると思う。
data %>%
	select(ID, BMI, SBP, DBP) %>%
	pivot_longer(-ID, names_to="variable", values_to="value") -> data_long
私なら
data_long <- data %>%
	select(ID, BMI, SBP, DBP) %>%
	pivot_longer(-ID, names_to="variable", values_to="value")
の方がしっくりくる。

上記以上に違和感を感じたのが、ggplotのコマンドでaesを単独で書く方法。
ggplot(data) +
	aes(x=Sex) +
	geom_bar(width=0.8, fill="black") +
	theme_classic()
私なら、このように書く。
ggplot(data, aes(x=Sex)) +
	geom_bar(width=0.8, fill="black") +
	theme_classic()
複数のgeom_にaesを使う時はどうするのだろう。そもそも、そのような書き方はダメということなのか。

また、質的変数、量的変数、そして、一変量、二変量の順に並んでいるので、散布図や折れ線グラフよりも先に、サンキー図やベン図などを学ぶことになる。

科学的データに慣れていない学生には、サンキー図の意義を説明することだけも一苦労する。

この本でRを勉強するよりは、多少、慣れてから頭の中を整理するのに使う本と考えた方が良いかも。それとも医学分野では、これがスタンダードなのかも。

2022年12月20日火曜日

サンキー図

午前
・標本整理

午後
・統計ゼミ
・データ整理

統計ゼミで勉強したサンキー図。

networkD3::sankeyNetwork( )で作成できるが、使用できるデータの様式に結構な癖があり、準備が難しい。
ゼミでは、超入門! Rでできるビジュアル統計学 学会・論文発表に役立つデータ可視化マニュアル 、を使って勉強した。
D3というJavaScriptが動いている?らしい。そのせいか、1/4名がエラーは出ないが、図が出力されなかった。

ネットではもっと綺麗な図のコマンドが紹介されている。また、networkD3で色々なネットワーク図が描ける。

networkD3:D3 JavaScript Network Graphs from R
35 フローチャート・サンキー図・タイムライン

2022年11月21日月曜日

面積を定量化

午前
・標本整理

午後
・データ解析
・会議

昨年からの悩みが、多少解決した。ただし、根本的な問題は解決できていない。

とある昆虫の砂丘内の分布データ。この分布の変化を定量的に評価したい、というのが昨年来の課題。
で例えば、geom_density_2d()を使ってカーネル密度推定を表示したり、
ggdensity::geom_hdr()を使ってhighest density regions (HDRs)を表示したり、
ggalt::geom_encircle()を使って輪郭を表示できることが分かった。
しかし、これらの面積を定量化する方法が分からなかった。

そこで、スマートな方法は諦めて、面積内を特定の色で塗りつぶして、そのピクセル数を計測する方法にたどり着いた。それはそれで面倒だったのだが、colorfindr::get_colors()というとても便利な方法を知った。

普通の散布図なら、グラフに使用していない色で塗りつぶせば良いのだが、上記のように地図上に描いている場合は、輪郭だけの図を作成した方が無難だろうと考えた。

その場合、地図の代わりに、空のラスタを準備して、そこに(例えば)黒色で塗りつぶした図を描く。

ちなみに、colorfindr::get_colors()は、一度、画像として保存して読み込まないと使えない(?)みたい。

プロットデータは、sf形式、平面直角座標のpointデータで、「Bn_2021_JGD2011」に入っているとする。
#枠となるラスタの作成
library(raster)
library(sf)
e <- raster::extent(-11000, -7000, -52500, -49500) #平面直角座標
sp_e <- as(e, "SpatialPolygons") #sp形式に変更
sf_e <- st_as_sf(sp_e) #sf形式に変更
st_crs(sf_e) <- "+proj=tmerc +lat_0=36 +lon_0=134.333333333333
		+k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0
		+units=m +no_defs +type=crs" #JGD2011;平面直角座標

#get_colors()は画像として一度保存する必要がある!
library(ggplot2)
png("Area_karnel.png", height=3000, width=3000)
ggplot() +
 	geom_density_2d_filled(aes(x=st_coordinates(Bn_2021_JGD2011)[,1],
  		y=st_coordinates(Bn_2021_JGD2011)[,2]), linewidth=0, col="white") +
	scale_fill_manual(values = c("white", "black", "black", "black", "black", "black", "black", "black")) +
	scale_x_continuous(limits = c(-11000, -7000)) +
 	scale_y_continuous(limits = c(-52500, -49500)) +
 	theme(panel.background = element_blank(), #背景を透明にする
		plot.background = element_blank(), #背景を透明にする
 		panel.grid.major = element_blank(), #グラフの線;邪魔なので白で消す
 		panel.border = element_blank(), #グラフの線;邪魔なので白で消す
		axis.title = element_blank(),
 		axis.text = element_blank(),
		axis.ticks = element_blank(), 
 		legend.position = 'none')
dev.off()
二値化の画像ができる。上記のカーネル密度推定は時計周りに回転させて、トリミングしているので一見異なるが、同じデータで作図している。
これを読み込んで、下記のようにget_colors()を実行すると、各色のピクセル数を計測できる。白と黒しかないのだが、他にもいくつか計測される、、、ごく僅か(ほとんど1ピクセル)なので無視して大丈夫だと思う。
library(colorfindr)
get_colors("Area_Bn_2021_karnel.png") #各色のピクセル数
これで算出できるのは、画像内の各色のピクセル数なので、比較する場合は、同じ大きさの画像で比較する必要がある。

(試していないが)面積そのものを計測したい場合は、例えば、100m2の方形区を作成して、そのピクセル数を基に計算すれば良さそう。

ただ、そもそも動物の分布を適切に表現できていないので、意味がないことも分かった。

2022年10月13日木曜日

イタリック

午前
・卒論手伝い(野外調査)

午後
・論文書き

ggplotで作成した軸名に学名が記されており、それをイタリックにするように指摘があった。

投稿時点で気づいていたが、やり方が分からず悩んでいるうちに修正することを忘れていた。

expression()を使うことでイタリックにはできるのだが、そもそも軸名をボールド表記するようにtheme()で指定してあるので、少し面倒だった。
ylab("No. of H. megaspinosa larvae")
とすると、ボールドの立字体となる。
そこで、次のようにする!
ylab(expression(paste(
    bold("No. of "),
    bolditalic("H. megaspinosa"),
    bold(" larvae"))))
イタリックかつボールドはbolditalic()で指定する。

expression()を実行するとtheme()の指定が消える?ので、bolditalic()以外のところはbold()で指定しないと普通の書体になってしまう。

あと、" "の中にスペースを入れる必要がある。

2022年9月14日水曜日

annotation_map_tile( )

午前
・卒論手伝い(野外調査)

午後
・資料作成

学生用の作図コマンド作成で、OpenSteetMapから国土地理院タイルに変更するのに思いっきりハマってしまった。
国土地理院タイルの全国最新写真(シームレス)を「annotation_map_tile()」で使う場合、「/{z}/{x}/{y}」を「/${z}/${x}/${y}」にする必要がある。
ggplot() +
  annotation_map_tile(zoomin = -1,
  type="https://cyberjapandata.gsi.go.jp/xyz/seamlessphoto/${z}/${x}/${y}.jpg")
エラーでずっと指摘されていたのが、意味が分からなかったのが、そのままだった。

2022年7月27日水曜日

ggtreeでNEXUSファイルの系統樹を描く

午前
・データ整理
・統計ゼミ

午後
・データ整理

ggtreeでNEXUSファイル(.nex)の系統樹を描くのにどハマりした。

そもそもの発端は、RAxMLで作成した系統樹を「ggtree」の「reroot()」を使うと、こんな感じでいまいちだったことに始まる。
「reroot()」のコマンドはこんな感じ。
#系統樹はML_tree.treeとして保管されている(nwk形式)
tree <- read.tree("ML_tree.tree")
tree_reroot <- root(tree, outgroup = "26") #外群のOTU名が26
str(tree_reroot)
ggtree(tree_reroot) +  geom_tiplab()
一応、「geom_rootedge()」も試したが改善できず。
ggtree(tree_reroot) +  geom_tiplab() + geom_rootedge(3)
ということで、Rでrerootをすることを諦めた。そこで、FigTree.appでrerootしてからRで作図する作戦に変更。

で、ここからどハマりした。

まず、FigTree.appでrerootをすると、「.nwk」だとboot strap値が保存できなくなる(?)ので「.nex」で保存しなければならない。

「.nex」で保存したので、「read.nexus()」を使えば良いだろうと思ったのだが、boot strap値が読み込めない(ファイルにはboot strap値はあることを確認)。

ここまでで半日以上を潰した。

試行錯誤の末、「read.beast()」を使えば良いことが判明!

普段、BEASTで作成したファイルを「read.beast()」で使っているので、これで終わった、、、と思ったら、別の意味で終わった。
#系統樹はML_tree_reroot_nex.nexとして保管されている
ML_TREE <- read.beast("ML_tree_reroot_nex.nex")
ggtree(ML_TREE) +  geom_tiplab()
これで作図できるはずなのだが、エラーがでる。「ggtree(ML_TREE)」だけだと実行できるので、OTUを描画する「geom_tiplab()」でエラーが出ているらしい。

悪戦苦闘の結果、原因が判明した。

まず、「read.beast()」で読み込んだ場合、下記のようにOTUは「tip.label」、 boot strap値は「label」として保存される。
#系統樹はML_tree_reroot_nex.nexとして保管されている
ML_TREE <- read.beast("ML_tree_reroot_nex.nex")
ML_TREE@phylo$tip.label #OTU
[1] "25" "24" "18" "16" "15" "13"
ML_TREE@data$label #boot strap値
label1  label2  label3  label4  label5 
    58      55      42       4      49
しかし、なぜか、「ggtree()」を経由すると、OTUは「label.x」、boot strap値は「label.y」になっていることが分かった!
ML_TREE <- read.beast("ML_tree_reroot_nex.nex")
ML_TREE_ggtree <- ggtree(ML_TREE)
ML_TREE_ggtree$data
どうやら、「geom_tiplab()」は、データ「label」を探しているが見つからないのでエラーが出ているようだ。

ということで、「label.x」を「label」に入れることで解決できた。
ML_TREE <- read.beast("ML_tree_reroot_nex.nex")
ML_TREE_ggtree <- ggtree(ML_TREE)
ML_TREE_ggtree_data <- ML_TREE_ggtree$data 
ML_TREE_ggtree_data$label <- ML_TREE_ggtree_data$label.x
ggtree(ML_TREE_ggtree_data) + geom_tiplab()
「.nex」を「read.beast()」で読むと発生する問題なのか、FigTree.appでrerootしたファイルだからなのか詳細な原因は良く分からないが、BEASTで作成したファイルではこのような問題は起こったことはない。 boot strap値についても同じような対応が必要な気がするが、私はboot strap値を一度取り出して別の変数に入れているので、こちらについては大きな問題となっていない。