ggplot の書式、特殊文字
ggplot2 内で、モデル式やR2 といった累乗の記号、ギリシャ文字を書きたいことが多々あると思います。
その場合は、expression 関数の例が多いと思うのですが(ggplot2でギリシャ文字や数式を表示したい - Qiita)、markdown や、文字コードを使っても色々かけます。
ggtext::element_markdown
markdown 形式で書くためのパッケージ・関数が、ggtext パッケージの element_markdown です。
Markdown theme elements • ggtext
使い方は、文字列を markdown 形式で書いて、theme で指定してやるだけです。
プロットエリア内は、geom_richtext を使います。
Improved Text Rendering Support for ggplot2 • ggtext
library(ggplot2)
library(ggtext)
g = ggplot()+
geom_richtext(aes(x=1, y=1, label="*Hello*<br>world</br>"))+
labs(x="*italic*", y="**bold**",
title="square<sup>2</sup> subscript<sub>2</sub>")+
theme(
plot.title = element_markdown(),
axis.title.x = element_markdown(),
axis.title.y = element_markdown()
)
plot(g)
自分が markdown に詳しくないのですが、どんな markdown でも対応しているわけではなさそうです。
書き方を調べるうえでは、Rmarkdownの書き方を参考にすれば、うまくいくことが多いです。
以下、element_markdown での書き方の例です。
順次更新していこうと思います。
| 斜体 | *word* <en>word<en/> |
| 太字 | **word** |
| 上付き | <sup>word<sup/> |
| 下付き | <sub>word<sub/> |
| 下付き | <sub>word<sub/> |
| アルファ(α) 斜体アルファ(α) |
α <i>α</i> |
| ベータ(β) | \β <i>β</i> |
| シータ(θ) | \θ <i>θ</i> |
RStanインストールのメモ
Windows、R 4.2.3、Rtools42の環境下では以下のページに従えばうまくいった。
discourse.mc-stan.org
remove.packages(c("StanHeaders", "rstan"))
install.packages("StanHeaders", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
install.packages("rstan", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
ちなみに失敗例は
- 何も考えずにrstanをインストールして、stan関数を動かしたときは「collect2.exe: error: ld returned 1 exit status no DLL was created」とエラーが出た
- 1の解決のために、古いバージョン（rstan 2.19.xx)をインストールしても、同様のエラーが発生
- 公式にしたがって、以下のコマンドを実行
Sys.setenv(MAKEFLAGS = paste0("-j",parallel::detectCores()))
install.packages(c("StanHeaders","rstan"),type="source")消えるはずの「Error in compileCode(f, code, language = language, verbose = verbose) : 」のエラーが出続ける。
最終的に、冒頭のページに従えば解決した。
開発中のパッケージバージョン（2.26.22）をインストールしたっぽい。
Visual studio codeのテーマを作ってみた
仕事の上でVisual studio codeをみんな使っていたので、自分も使ってみました。
しかし、テーマカラーの色で気に入ったものがなかったため、自分で作ってみました。
作ったからには人に見せたいと思ったので、テーマを置いてあるGithubのリンクを張っておきます。
GitHub - hiroakif93/VScode_theme: My visual studio code theme
使い方は、Githubからおとした「snow」フォルダを、「ホーム」ディレクトリにある「.vscode」フォルダ内の「extensions」フォルダに置くだけ（「.vscode」は隠しフォルダになってたので、コマンドラインからコピー）。
そしたら、VS code内の設定から変更できるようになります。
エディターは白のほうが好き（慣れてしまった）なのですが、他のテーマ見てたら白が強くてさすがに目が疲れたので、少しぼかしてるのが特徴です。
あと、白の要素を減らすためにターミナルは黒にしてみました。
名前の由来は、以下のGoogle Chromeのテーマを真似したからです。
https://chrome.google.com/webstore/detail/red-fox-snow-theme/cgaadipmojdihomphfmjphmelinpdalg
作るにあたっては以下のサイトを参考にさせていただきました。
atmarkit.itmedia.co.jp
構文ハイライトはこれからいじっていきたいと思います。
追記
2023/5/15
別のパソコンに落としてもうまくいかなかった。。。
ggplot2 で追加のscale_color / scale_fillを指定する方法
とても便利なのに、日本語サイトで見ないのでメモを兼ねて書きます。
よく ggplot2 を使っていると、geom_point と geom_bar で、違うcolor_scale を指定したいことがあると思います。また、色の指定の際に、連続値(continuous scale) と離散値（discrete scale）の両方を使うとエラーが出てしまいます。
それを解決してくれるのが以下のパッケージです。
使い方はとても簡単で、scale_color_,,, の後に、new_scale() を追加するだけです。
使用するパッケージを読みこみます
library(ggplot2) library(ggnewscale) ## 新しくscaleを追加するためのパッケージ library(tidyr)
テストデータを作成。
site=10
sp=5
# サンプル×種の表
mat = matrix(runif(site*sp), nrow=site, ncol=sp,
dimnames=list(paste('S', formatC(1:site,width=2,flag="0"), sep='_'),
paste('Sp', formatC(1:sp,width=2,flag="0"), sep='_')))
# 環境データ
env = data.frame('Sample' = rownames(mat),
'Env1' = as.integer(runif(site, 0, 2)),
'Env2' = as.integer(runif(site, 0, 4)),
row.names = rownames(mat))
# 2つの表を引っ付ける
df <- cbind(env, mat)
# ggplot2 を使うため、ロングフォーマットになおす。
lf <- gather(df, key, value, -c(1:ncol(env)))
バープロットの作成
bar_color <- c("firebrick", "skyblue3", "forestgreen", "darkorange", "aliceblue")
bar1 = ggplot(lf, aes(x=Sample))+
geom_bar(aes(y=value, fill=key), stat='identity')+
geom_point(aes(y=Env1, fill=Env1),
shape=21, # fill を使いためpoint の形を変更
stroke=0.3# 枠線の太さ
)+
scale_fill_manual(values=bar_color)
ここで、geom_point を環境の値によって色を変えたいのですが、連続値(環境変数)と離散値(種名)を入れているため、以下のようなエラーが出ます。
Error: Continuous value supplied to discrete scale
次に、new_scale 関数を使用します。
bar_color <- c("firebrick", "skyblue3", "forestgreen", "darkorange", "aliceblue")
bar1 = ggplot(lf, aes(x=Sample))+
geom_bar(aes(y=value, fill=key), stat='identity')+
scale_fill_manual(values=bar_color)+ ## <- ここでバープロット用のscale_fill を指定
ggnewscale::new_scale_fill()+ ## <- geom_pointの前にnew_scale_fill を挟む
geom_point(aes(y=Env1, fill=Env1),
shape=21,
stroke=0.3
)+
scale_fill_gradient2() ## <- geom_pointの用のscale_fill を指定
以下の図が作成できます。
離散値と連続値の色を指定できるようになるのは、とても便利ですね。
また、色以外の scale (scale_size、scale_shape)にも対応していそうです。
よくある分類群ごとの積み上げグラフ（バープロット）の作り方の例
研究室内で、「分類群ごとの積み上げグラフをどうやって作るか」という質問があったので簡単な作り方の例を紹介します。
（こういうやつ↓）
Networks Depicting the Fine-Scale Co-Occurrences of Fungi in Soil Horizons | PLOS ONE
用意するファイルは以下の3つです。
- サンプル×種の表 (行がサンプル、列が種、中身は個体数情報など）
- 分類群の情報 ( 行が分類群の最小単位、種でもいいしNGSデータならOTU・ASVです。列が最小単位の分類群の上位階層）
- あれば環境データ(行がサンプル、列が環境データ）
テストデータです
site=10
sp=5
# サンプル×種の表
mat = matrix(runif(site*sp), nrow=site, ncol=sp,
dimnames=list(paste('S', formatC(1:site,width=2,flag="0"), sep='_'),
paste('Sp', formatC(1:sp,width=2,flag="0"), sep='_')))
# 分類群の情報
taxa = data.frame('Phylum' = formatC(as.integer(runif(sp, 0, 2)),width=2,flag="0"),
'Order' = formatC(as.integer(runif(sp, 0, 4)),width=2,flag="0"),
'Genus' = formatC(as.integer(runif(sp, 0, 8)),width=2,flag="0"),
row.names = colnames(mat))
for(i in 1:ncol(taxa)){ taxa[,i] <- paste(colnames(taxa)[i], taxa[,i], sep='_')}
# 環境データ
env = data.frame('Sample' = rownames(mat),
'Env1' = as.integer(runif(site, 0, 2)),
'Env2' = as.integer(runif(site, 0, 4)),
row.names = rownames(mat))サンプル×種の表の列が、最小分類群単位になっているので、これを見てみたい分類群の階層で、足し合わせます。(例えば、Sp_01 から Sp_05 は、同じ Order_03 という目に入るので、それらを足してあげます）
ただし、これを行うのは、種数が多いと大変です。
なので、関数を用意しました。
# x がサンプル×種の表
# y が分類軍の情報
# taxaLabel がまとめたい分類群の名前
# function は足しわせるのに都合が悪い時用に、用意してあります。
# 注意点は x のcolnames と y のrownamesが対応してある必要があります
Taxa.mat <- function(x, y, taxaLabel, func=function(x){sum(x)}){
colnames(x) <- y[colnames(x), taxaLabel]
summary <- do.call(cbind,
lapply(unique(colnames(x)),
function(a){ num <- which(colnames(x)==a)
apply(as.matrix(x[,num]), 1, func)}) )
colnames(summary) <- unique(colnames(x))
summary
}
あとは、ggplot の geom_bar にかけてあげるだけできます。
library(ggplot2)
library(tidyr)
merge.mat <- Taxa.mat(mat, taxa, 'Genus')
# ロングフォーマットになおす。その際に環境データをくっつける
lf <- gather(cbind(env, merge.mat), key, value, -c(1:ncol(env)))
head(lf)
# bar1 がただの積み上げ。geom_bar の引数のposition = 'fill を指定すると、100 % 積み上げグラフになります
bar1 = ggplot(lf)+
geom_bar(aes(x=Sample, y=value, fill=key), stat='identity')
bar2 = ggplot(lf)+
geom_bar(aes(x=Sample, y=value, fill=key), stat='identity', position='fill')
# ｘ軸を環境データ(Env2 などにしたら、さらに環境でまとめられます）
bar3 = ggplot(lf)+
geom_bar(aes(x=Env2, y=value, fill=key), stat='identity', position='fill')
x軸ラベルの重なりなどは、以前の記事を参考にしてもらえればと思います。
ggplot2 で、軸ラベルの重なりを防ぐ方法 - 微生物屋のノート
以上が一例になります。
この記事を書きながら、ggplot の段階で色分けする方法もあるなと気づいていしまったので、今回の記事はやりやすい方法を見つける参考にしていただければと思います。
anova や glm を、for 関数で処理する方法
色々な処理区があり、各処理区の群集を特徴付ける指標がたくさん得られてきた。
そして、その処理区ごとの各指標の平均値の違いは、処理区ごとに説明できるのかを統計的に示す時に、分散分析を行うと思います。
Rプログラムの中では、anova 関数や aov 関数が分散分析を行う関数になり、この関数は、formula 型の Input ( y ~ x, data=df )を要求してきます。
この aov 関数と for 関数を使って、いろんな指標の分散分析を一括で処理する時、formula 関数が便利です。
次のような data frame (省略版) があるとします。
A~Gの処理区があって、3つの指標を測ったマトリックスになっています。各処理区にはに100個の値（１００サンプル）が、ほんとはあります。
treatment Richness Abundance Interaction 1 A 63.03062 99.183025 61.779550401 2 A 59.28171 92.629037 62.014830418 3 A 60.02215 99.844330 59.648745452 4 B 77.90848 46.798696 54.998186842 5 B 71.64935 41.451237 61.789154721 6 B 71.02360 52.430302 56.806729181 7 C 32.37181 79.320784 44.909000587 8 C 40.65879 65.529647 55.526555285 9 C 31.65834 65.384726 44.569884241 10 D 93.80565 62.985476 86.305609523 11 D 73.20372 62.672580 83.021430326 12 D 76.77968 69.194962 81.886503028 13 E 56.08483 -12.653348 4.713663180 14 E 61.87220 -5.530885 0.005947824 15 E 50.68268 5.953438 -2.610380320 16 F 72.80016 3.270578 57.132157097 17 F 76.54925 1.466419 85.587994020 18 F 70.88253 3.530791 79.221027188 19 G 97.02695 12.273005 7.767906549 20 G 86.16959 -1.815169 0.361944723 21 G 77.37816 26.637699 8.648186671
図にすると以下の通りです。
それでは、各指標、１列ごとにfor 関数を回して、分散分析を行います。
この時、paste 関数で式の形をした文字列を作り、それを formula 関数で type : formula にしてあげます。
result <- list()
for(i in colnames(summ)[-1]){
formula.tmp <- paste(i, 'treatment', sep='~') # 中身としては"Interaction~treatment"という形
formula.object <- formula(formula.tmp) # object の型をformulaにする
result[[i]] <- aov(formula.object, data=summ) # aov 関数にかける
}これで、result の中に各指標の結果を放り込むことができました。
おまけ
分散分析で有意差が出た後、Post-hoc検定、処理区間のペアワイズで何らかの検定を行うそうです。
処理区ごとに、有意差があると見せる時はgeomsignif package が便利そうです。
Rで解析：ggplot2のプロットに有意差バーを追加！「ggsignifr」パッケージ
しかし、処理区間で有意差が出まくると、線がたくさん引かれて微妙な気持ちになったので、別の見せ方のコードを提示しておきます。
線ではなく、有意差があった処理区の名前を表示するようにして見ました。
ここでのポイントは、geom_textの高さをvjust で揃えているところです。
実データではきれいにできたのですが、テストデータだと微妙でしたね...
＊TukeyHSDで２群間比較を行なっていますが、どの方法でPost-hoc検定を行うかは使うデータによって判断してください
summ <- c()
for(l in LETTERS[1:7]) {
df <- matrix(NA, nrow=3, ncol=3, dimnames=list(NULL, c('Richness', 'Abundance', 'Interaction') ))
for(i in 1:ncol(df)){
df[,i] <- rnorm( n= 3, mean=sample(1:100, 1), sd=sample(1:10,1))
}
summ <- rbind(summ, cbind(treatment=l, as.data.frame(df)))
}
library(tidyr)
library(ggplot2)
library(RColorBrewer)
lf.sum <- gather(summ, key=index, value, -1)
head(lf.sum)
ggplot(lf.sum) +
geom_boxplot(aes(x= treatment, y=value)) +
facet_wrap(~key)
result <- list()
for(i in colnames(summ)[-1]){
p.mat <- matrix(NA, ncol=7, nrow=7, dimnames=list(unique(summ$treatment), unique(summ$treatment)))
formula.tmp <- paste(i, 'treatment', sep='~')
formula.object <- formula(formula.tmp)
result[[i]] <- aov(formula.object, data=summ)
if( summary(result[[i]])[[1]][[5]][1] < 0.05 ){
test <- TukeyHSD(aov(formula.object, data=summ))$treat
for(l in 1:nrow(test)){
a <- strsplit(rownames(test), '-')[[l]]
p.mat[a, a] <- test[l,4]
}
diag(p.mat) <- NA
}
for(l in 1:nrow(p.mat)){
sig <- names(which(p.mat[l,] < 0.05))
lf.sum[which( lf.sum$treatment==rownames(p.mat)[l] & lf.sum$index==i ), 'pval'] <- paste(sig, collapse='\n')
lf.sum[which( lf.sum$treatment==rownames(p.mat)[l] & lf.sum$index==i ), 'y'] <-
quantile(lf.sum[which( lf.sum$index==i ), 'value'], probs=seq(0,1,0.1))[11]
}
}
gg <- ggplot(lf.sum)+
geom_boxplot(outlier.size = 0,aes(x=treatment, y=value, color=treatment, fill=after_scale(alpha(color, 0.5) )), width=0.5)+
geom_text( aes(x=treatment, y=y, label=pval),vjust="inward" )+
facet_wrap(~index, scales='free', ncol=2) +
scale_color_manual(values=brewer.pal(8,'Accent')[-c(6)])+
scale_shape_manual(values=c(21:25,6:8))+
theme_minimal(base_size=18) + theme(panel.border=element_rect(fill=NA, size=1),)+
scale_x_discrete(guide = guide_axis(n.dodge = 2))
ggplot2 で、軸ラベルの重なりを防ぐ方法
自分が忘れてたので、使い方と感動を忘れないために。
たまたまggplot2 のバージョンが更新されたのを見て、何が変わったのか見てみたら、すごく感動する仕様になっていた。
その一つが、軸ラベルの重なりを防ぐ方法である。
Test data を作成し、ggplotでプロットしてみる
test <- data.frame(x=c(rep('whole genome seq', 4), rep('whole metagenome seq', 4), rep('whole genome shotgun seq', 4), rep('whole metagenome shotgun seq', 4) ),
value=rnorm(4*4))
gg <- ggplot(test)+
geom_jitter(aes(x=x, y=value), width=0.2) +
geom_boxplot(aes(x=x, y=value), width=0.5)
plot(gg)
見てわかる通り、ｘ軸の名前がわからない。
これを防ぐために、scale_x_discrete(guide = guide_axis(n.dodge = 2)) を追加する。
test <- data.frame(x=c(rep('whole genome seq', 4), rep('whole metagenome seq', 4), rep('whole genome shotgun seq', 4), rep('whole metagenome shotgun seq', 4) ),
value=rnorm(4*4))
gg <- ggplot(test)+
geom_jitter(aes(x=x, y=value), width=0.2) +
geom_boxplot(aes(x=x, y=value), width=0.5)+
scale_x_discrete(guide = guide_axis(n.dodge = 2))
plot(gg)
scale_x_discrete(guide = guide_axis(n.dodge = 2)) を追加することで、重なったラベルを縦にずらすことが可能となった。
また、よく使われる方法としては、軸の角度を変えたり、そもそも字の大きさを小さくする方法がある。
・軸ラベルの角度を変える
gg <- ggplot(test)+
geom_jitter(aes(x=x, y=value), width=0.2) +
geom_boxplot(aes(x=x, y=value), width=0.5)+
theme(axis.text.x = element_text(angle=60, vjust=1, hjust=1))
plot(gg)
さらに思いつくのは、少々めんどくさいが、ラベルを改行する方法である。
やり方は、正規表現の改行を意味する \n を、名前にはさんでやる。
・軸ラベルを改行する
test <- data.frame(x=c(rep('whole\ngenome\nseq', 4), rep('whole\nmetagenome\nseq', 4), rep('whole\ngenome\nshotgun\nseq', 4), rep('whole metagenome\nshotgun\nseq', 4) ),
value=rnorm(4*4))
# または、gsub 関数で列単位で処理。今回は、スペースを \n に置き換えた。
#test <- data.frame(x=c(rep('whole genome seq', 4), rep('whole metagenome seq', 4), rep('whole genome #shotgun seq', 4), rep('whole metagenome shotgun seq', 4) ),
# value=rnorm(4*4))
#
# test[,'x'] <- gsub(' ', '\n', test[,'x'])
gg <- ggplot(test)+
geom_jitter(aes(x=x, y=value), width=0.2) +
geom_boxplot(aes(x=x, y=value), width=0.5)
plot(gg)
以上、軸ラベルの重なりを防ぐ方法を紹介した。
軸ラベルを斜めにすると場所を取るので、ずらして表示することができるようになって、すごく感動した。
ただし、今日のように、軸ラベルが長すぎるとプロット枠の外に出てしまうので、余白を設定する必要があるのかもしれない。
そういう点では、どれか一つの方法にこだわらずに、いくつかをラベルずらしと改行を組み合わせてやるなどをしてもいいのかもしれない。
