初めに

近年の生物系の論文のグラフは生データを点で示し、検定結果とともに表示することが多くなっています。

このようなグラフはエクセルで作成するのは大変なのでRで作成したい、しかし記述量はできるだけ減らしたいというときは参考にしてください。

Rで最低限のグラフを出力し、パワーポイント等で整形する方法です。

ダイナマイト・プランジャー・プロット (平均値の棒グラフ＋エラーバー）は非推奨

追記:　コメントでご指摘いただきました。初版ではダイナマイト・プランジャー・プロットを使い記事を作成していましたが、このグラフは複数の測定値を集計して表示する際には適していない点があります。ダイナマイトプロットは以下のような図です。ダイナマイトプランジャーは、右の図のようなダイナマイトを発破する装置。

問題点としては、
・複数のデータを要約するグラフは、データが分布している領域を塗りつぶすべき。
・エラーバーは標準誤差や標準偏差、95％信頼区間などが使われるが、何が使われているのかわからず誤解を生む。正規分布をイメージしてしまうが実際わからない。
・平均値＋エラーバーでは情報を圧縮しすぎており、データの分布がわからない。
・元データを示しても平均値というただ一つの統計値が目立ちすぎ

参考

そこで今回はダイナマイトプロットの代替で書きます。

エクセルで読み込むファイルの作成

以下のような2列の表を作成します。

A1にはname, B1にはvalueと記述し、2行目以降のA列はサンプル名、B列は測定値を入れます。1列目をnameとvalueにしておけばこのページのコードをそのまま使えます。

今回はsample.xlsxというファイルのSheet1に入力しました。たくさん図を作りたい場合は、別シートで似たように作成してください。

サンプル名の並び順がそのままグラフの並び順になります。

name	value
WT_Cont	101.1
WT_Cont	80.7
WT_Cont	92.9
WT_Cont	106.8
Mutant_Cont	78.7
Mutant_Cont	73.7
Mutant_Cont	82.4
Mutant_Cont	76.1
WT_Treat	124.4
WT_Treat	129.5
WT_Treat	119.1
WT_Treat	125.7
Mutant_Treat	135.0
Mutant_Treat	121.6
Mutant_Treat	144.8
Mutant_Treat	120.0

作業フォルダーの変更。

読み込みたい エクセルファイルの存在する位置 に設定します。

Windowsの場合はRのウインドウの左上の　ファイル＞ディレクトリの変更 です。Macの場合は command+d です。出てきたフォルダ選択の画面で作成したエクセルファイルの存在するフォルダーを選んでください。この場合は以下のsetwd()コマンドはいらないです。

R studioやJupyterLab、Linuxを使っている場合はエクセルファイルのあるフォルダーを右クリックして、パスのコピー、貼り付けて以下のかっこの中を書き換えます。ダブルクォーテーションを忘れずに。
Windowsの場合はバックスラッシュまたは円マーク \ をスラッシュ / に書き換えます。

R: コマンドでしか作業フォルダーの変更ができない場合

setwd("C:/document/R")

ライブラリーのインストール

インストールされていないライブラリーのみインストールします。

以下の記事を参考にしました。

初回のみ以下の枠内を一気に貼り付けて実行します。

なんか色々出てくると思うので、終わるまで待ちます。結構長いです
">"が出てきたら終わりです

インストールに失敗したというエラーが出た場合は最新のRを入れなおしてください。

R: 初回のみ

libs <- c("multcompView", "tidyverse", "ggpubr", "openxlsx", "svglite","export")
for(lib in libs) { 
  if(!require(lib, character.only = T)) install.packages(lib) 
}

ライブラリーの読み込み

multcompViewはTukey検定で使い。tidyversはdataframeの操作に必要です。ggpubrはggplot2を使って便利に作図できます。openxlsxはエクセルファイルを開く際に使います。exportはパワーポイントを出力するのに使います。
Macではエラーになることがあります。OpenGLとrglに関連したエラーが出たら下のコマンドを試してみてください。

参考ページ:

library(multcompView)
library(tidyverse)
library(ggpubr)
library(openxlsx)
library(export)

R:MacでOpenGL関連のエラーになる場合

options(rgl.useNULL = TRUE)
library(rgl)
library(export)

ファイルの読み込み

作成したエクセルを読み込みます。sample.xlsxというファイル名でSheet1というシート名から読み込む場合は以下のコードです。
適せん自分の用意したファイル名とシート名に変更してください。

data1 = read.xlsx("sample.xlsx", sheet = "Sheet1")
data1

Rは並び順を勝手に変えるので、固定します。
サンプル名を抽出し、並び順を指定するための配列を作成。
data1[ ,1]で表の1列目を指定しています。

order1 = unique(unlist(data1[ ,1]))
data1[ ,1] = factor(data1[ ,1], level = order1)
order1

グラフのy軸の範囲や検定結果を表示するためにデータの最大値を取り出します。

ymax = max(data1[ ,2])

グラフの作成

ggpubrパッケージのggbarplotでグラフを作ります。

参考ページ:

サンプルデータのような小規模なデータ（n < 10）では散布図と平均値がよいでしょう。
3行目はy軸の範囲を指定しています。

2行目のx = "name", y = "value"は作成した表の列名に書き換えます。

R: 散布図と平均値

fig = ggstripchart(data1, x = "name", y = "value", size = 3, 
        add = c("mean") , add.params = list(color = "red") )
fig = ggpar(fig, ylim = c(0, ymax * 1.1)) + scale_y_continuous(expand = c(0, 0, 0.1, 0)) 
fig

比較用として棒グラフ。できるだけ上の散布図を使います。add=c()の"mean"を"mean_sd"や"mean_se"にすると上で紹介したダイナマイトプロットになってしまいます。

R: 平均値の棒グラフ

fig = ggbarplot(data1, x = "name", y = "value", add = c("mean","jitter"))
fig = scale_y_continuous(expand = c(0, 0, 0.1, 0)) 
fig

箱ひげ図、サンプルサイズ n > 10 程度。あまり少ないデータだと四分位数の意味をなさないのと、一つのデータによってグラフの概形が大きく変わってしまいます。

R: 箱ひげ図

fig = ggboxplot(data1, x = "name", y = "value", add = c("jitter"))
fig = ggpar(fig, ylim = c(0, ymax * 1.1)) + scale_y_continuous(expand = c(0, 0, 0.1, 0)) 
fig

さらにn数が多ければバイオリンプロット。add の中身の "jitter"を消せば、散布図は消えます。

R: バイオリンプロット

fig = ggviolin(data1, x = "name", y = "value", add = c("jitter", "boxplot"))
fig = ggpar(fig, ylim = c(0, ymax * 1.1)) + scale_y_continuous(expand = c(0, 0, 0.1, 0)) 
fig

検定結果を追加

t検定か tukey検定の結果を追加できます。

1. T検定とともに出力

一行目は、検定したい組み合わせを書く c("A", "B")の形式で書きます。増やす場合はカンマでつなぎ　c("A", "B"), c("C", "B"), c("A", "D")のようにします。

5行目はt検定。 var.equal = TRUE でstudentのt-tesになり、var.equal = FALSE でWelchのt-testになります。

comp = list( 
    c("WT_Cont", "Mutant_Cont"),
    c("WT_Treat", "Mutant_Treat") ) 
    
fig = fig + stat_compare_means(
    comparisons = comp, 
    method = "t.test",
    method.args = list(var.equal = TRUE) )  
fig

2. Tukey test　ともに出力

1行目の(value ~ name)は作成した表の列名に書き換えます。nameとvalueにしておけばそのままでよいです

anova = aov(value ~ name, data = data1)
tukey = TukeyHSD(anova)
tukey_res = multcompLetters4(anova, tukey, threshold = 0.05)[[1]]$Letters
tukey_res = data.frame(tukey_res)
tukey_res$name = factor(row.names(tukey_res),levels = order1)
tukey_res = arrange(tukey_res, name) #並び替え
tukey_res$group1 = tukey_res$group2 = 1
tukey_res

A data.frame: 4 × 4
	tukey_res	name	group1	group2
	<chr>	<fct>	<dbl>	<dbl>
WT_Cont	b	WT_Cont	1	1
Mutant_Cont	b	Mutant_Cont	1	1
WT_Treat	a	WT_Treat	1	1
Mutant_Treat	a	Mutant_Treat	1	1

fig = fig + stat_pvalue_manual(
  tukey_res, x = "name", y.position = ymax, label = "tukey_res",
  position = position_dodge(0.8), size = 8
)
fig

テキスト部分の編集

パワポで編集もできますが、ここで編集する方法もあります。

Font sizeの変更

xy.textは軸の数値のフォントサイズ、xyは軸ラベルのフォントサイズ

fig  = fig + font("xy.text", size = 20) + font("xy", size = 20)
fig

x軸ラベルを傾ける

fig = ggpar(fig, x.text.angle = 45)
fig

パワーポイントで保存

ファイル名は適当なものに変更してください。

別窓で出ているグラフのウインドウの大きさを変更すると保存されるファイルにも図のサイズが反映されます。保存したいサイズにマウスで合わせてもよいです。後にパワーポイント上で複数枚並べたい時などは、サイズを指定して揃えられます。

どちらか一方を実行してください

graph2ppt(fig, "sample.pptx") #表示されているサイズで保存
graph2ppt(fig, "sample.pptx", width = 8, height = 8) #縦横の長さを指定して保存

エクセルファイルの存在するフォルダに出力されるはずです。

SVGで保存

上のコマンドが動かない時や、パワーポイント以外（Inkscapeやイラストレーター等）で編集したいときはSVGで保存します。

ファイル名は自由に変更してください。(拡張子が.svgならSVG画像。.tiffならTIFF画像。.png にすればPNG画像として保存できます）最初に選択したフォルダー内にSVGファイルはできています。

ggsave("sample.svg", fig) #表示されているサイズで保存
ggsave("sample.svg", fig, width = 8, height = 8) #縦横の長さを指定して保存

パワーポイントで編集可能です。ドラックアンドドロップで図をパワーポイントに貼り付け、図を一度左クリックし、"Ctr+Shift+g"でグループ化の解除（windowsの場合）をして、ポップアップがでたら[はい]を選択します。Macの場合はCommand+option+Shift+gです。

その後2回ほどグループ化を解除すれば、ばらばらになるので、グラフの幅やフォントサイズ、色を変更可能です。

パワーポイントで加工

生データの点は枠線を太くすれば大きくなります。透明な四角形が重なってきているので、消すと、下のオブジェクトがクリックできるようになります。
散布図の平均値を横棒にパワポと手作業で置き換えて色付けするとこうなります。

Rで生データの点と検定結果をつけたグラフを出力し、パワーポイントで整形する方法

初めに