ggplot2でお絵かき:棒グラフ

バイオ系博士学生のGoro (@BioDr_goro) です。

実験データを図としてまとめる際に，"平均値" と "標準偏差 (誤差)" しか示さないのはよくないので生データも載せよう，という風潮が高まっている。

journals.plos.org

上記論文のFig. 1では，パネルB–EすべてがパネルAの棒グラフで図示できることの危険性を説いている。
もちろんp値はB–Eで異なる。
生データを記載することで読者のミスリードを防ぎ，著者が適切に統計的処理できているかどうかが判断できる。

実際に，データのプロット数が3つ以下の場合はすべてのデータをプロットすることがNature Cell Biology誌の統計ガイドラインで決まっているらしい。
https://www.mbsj.jp/admins/ethics_and_edu/PNE/6_article.pdf

例えば，以下の論文では生データ＋平均＋標準誤差を図示している。

www.ncbi.nlm.nih.gov

ここからは私が作成した模擬データを使って説明する。
最初に紹介したPerspectiveでは，Aのような図ではなくBのような図を作ろう，ということを訴えている。

f:id:LPasteur:20190415014139p:plain — 遺伝子AとBの実験開始時および12時間後における発現量. (A) 棒グラフ, (B) 散布図.

上記のデータはエクセルで正規分布に従う乱数を発生させ，エクセル上で棒グラフと散布図を作った。

実験条件は，
「試験管を10本用意し，大腸菌を培養した。実験開始時 (0 h) に増殖阻害剤を添加した。実験開始時と12時間後 (12 h) にRNAを抽出し，逆転写後に遺伝子AとBのcDNAを標的として定量 PCRをした。」
とする。

元データはこちら。
demo_data - Google ドライブ

棒グラフや散布図を描くだけならエクセルで十分であるが，生データ，平均，標準 (偏差) 誤差を1つの図にまとめるのは至難の技。
またエクセルでは図示できるデータ数に限界があるため，ビッグデータを扱うには不向きである。
そこで，描画システムとしての機能も持ち合わせているRの出番だ。
ggplot2というパッケージは非常に優れもので，使いこなすことができれば論文で見かけるほとんどの図表は簡単に作れるようになる。

では図を作っていく。この記事ではタイトルにもあるようにパネルAのような棒グラフを作ることを目的とする。
結局全部のプロットを載せないんかい！というツッコミはなしで…。次の記事でまとめる予定だ。

なお，今回の環境は以下の通り。

macOS Mojave version 10.14.4
R version 3.5.0 GUI 1.70 El Capitan build
RStudio version 1.2.1335
ggplot2 version 3.1.1

1. 作業ディレクトリを変更する
2. 今回使うパッケージを呼び出す
3. データをRに読み込ませる
4. 遺伝子ごと，時間ごとに平均値と標準偏差を求める
5. エラーバーをつくる
6. とりあえず棒グラフを作ってみる
7. 体裁を整える
8. 有意差バーをつける
参照したサイト

1. 作業ディレクトリを変更する

RStudioの場合，ツールバーからディレクトリを変更できる。
f:id:LPasteur:20190416003921p:plain

2. 今回使うパッケージを呼び出す

library(ggplot2)
library(ggthemes)
library(dplyr)
library(ggsignif)

そんなパッケージないよ！って怒られたらエラーコードをググれば解決策がすぐにヒットする。

3. データをRに読み込ませる

もし表全体をコピーして，クリップボードから読み込ませたい場合は，

#Windows
demo_1 <- read.table("clipborad", header = T, sep = "\t")

#Mac
demo_1 <- read.table(pipe("pbpaste"), sep = "\t", header = T)

となる。

もしcsvファイルを読み込ませたいのであれば，

demo_1 <- read.csv("demo_data.csv", header = T, sep = ",")

とする。
このとき作業ディレクトリ内にファイルがないとエラーを起こす。
右上のウィンドウにdemo_1が表示されたら問題ない。
demo_1をクリックするとRStudio上でデータを閲覧できる。

f:id:LPasteur:20190416004438p:plain

4. 遺伝子ごと，時間ごとに平均値と標準偏差を求める

demo_mean_sd <- demo_1 %>% group_by(group, time) %>% summarize(mean = mean(value), sd = sd(value))

Macの場合，%>% (パイプと読む) はコマンド＋シフト＋Mで出せる。
次のようなオブジェクトができているはずだ。

f:id:LPasteur:20190416005933p:plain

5. エラーバーをつくる

errors <- aes(ymax = mean + sd, ymin = mean - sd)

6. とりあえず棒グラフを作ってみる

plot <- ggplot(demo_mean_sd, aes(x = group, y = mean, fill = time)) +
geom_errorbar(errors, position = position_dodge(0.9), width = 0.1) +
geom_bar(stat = "identity", position = "dodge")

plot

ggplotはレイヤーを重ねていくように描画していく。

ggplot(demo_mean_sd, aes(x = group, y = mean, fill = time))
→オブジェクトを読み込ませ，x軸にgroup列 (遺伝子名) を，y軸にmean (平均値) をとり，色をtime列 (時間) で識別

geom_errorbar(errors, position = position_dodge(0.9), width = 0.1)
→エラーバーをつける。positionをいじらないとエラーバーが棒グラフの中心にこないので設定した。widthではエラーバーのTの横線の長さを変えられる。

geom_bar(stat = "identity", position = "dodge")
→stat = "identity"はx，y軸の両方に値を指定していることを意味する。position = "dodge" は積み上げではなく横にグラフを配置することを指定している。

RStudioは図の出力も簡単にできる。

f:id:LPasteur:20190416011240j:plain

7. 体裁を整える

plot <- ggplot(demo_mean_sd, aes(x = group, y = mean, fill = time)) +
geom_errorbar(errors, position = position_dodge(0.9), width = 0.1) +
geom_bar(stat = "identity", position = "dodge", colour = "black") +
scale_fill_manual(values = c("gray", "white")) +
labs(x = "", y = "Expression level [copies/ng]") +
theme_classic() +
theme(axis.title = element_text(size = 15, family = "Helvetica", colour = "black"), axis.text = element_text(size = 15, family = "Helvetica",  colour = "black")) +
scale_y_continuous(expand = c(0,0), limits = c(0,8))

plot

fillは塗りつぶしの色を，colourは枠線や文字の色を指定している。

scale_fill_manual(values = c("gray", "white"))
→ggplotで指定した fill = time の色をいい感じに。

labs(x = "", y = "Expression level [copies/ng]")
→軸の名前を設定。

theme_classic()
→見た目の印象をクラッシックな感じに。

theme(axis.title = element_text(size = 15, family = "Helvetica", colour = "black"), axis.text = element_text(size = 15, family = "Helvetica", colour = "black"))
→軸ラベルの文字ポイントやフォントを決定。

scale_y_continuous(expand = c(0,0), limits = c(0,8))
→y軸のスケールを決定。デフォルトのままだと棒グラフがx軸に浮いているのが気になる。

他にも，テーマを変えたり，表を分割したりできる。

plot <- ggplot(demo_mean_sd, aes(x = group, y = mean, fill = time)) +
geom_errorbar(errors, position = position_dodge(0.9), width = 0.1) +
geom_bar(stat = "identity", position = "dodge", colour = "black") +
scale_fill_manual(values = c("gray", "white")) +
labs(x = "", y = "Expression level [copies/ng]") +
theme_few() +
theme(axis.title = element_text(size = 15, family = "Helvetica", colour = "black"), axis.text = element_text(size = 15, family = "Helvetica",  colour = "black")) +
scale_y_continuous(expand = c(0,0), limits = c(0,8))

plot

f:id:LPasteur:20190416014310j:plain

plot <- ggplot(demo_mean_sd, aes(x = group, y = mean, fill = time)) +
geom_errorbar(errors, position = position_dodge(0.9), width = 0.1) +
geom_bar(stat = "identity", position = "dodge", colour = "black") +
scale_fill_manual(values = c("gray", "white")) +
labs(x = "", y = "Expression level [copies/ng]") +
theme_bw() +
theme(axis.title = element_text(size = 15, family = "Helvetica", colour = "black"), axis.text = element_text(size = 15, family = "Helvetica",  colour = "black")) +
scale_y_continuous(expand = c(0,0), limits = c(0,8))

plot

f:id:LPasteur:20190416014339j:plain

plot <- ggplot(demo_mean_sd, aes(x = group, y = mean, fill = time)) +
geom_errorbar(errors, position = position_dodge(0.9), width = 0.1) +
geom_bar(stat = "identity", position = "dodge", colour = "black") +
scale_fill_manual(values = c("gray", "white")) +
labs(x = "", y = "Expression level [copies/ng]") +
theme_few() +
theme(axis.title = element_text(size = 15, family = "Helvetica", colour = "black"), axis.text = element_text(size = 15, family = "Helvetica",  colour = "black")) +
scale_y_continuous(expand = c(0,0), limits = c(0,8)) +
facet_wrap(~group, scale = "free")

plot

f:id:LPasteur:20190416014424j:plain

8. 有意差バーをつける

Rで有意差を求めることは当然できるのだが，同じ遺伝子内で0 hと12 h間でt検定を指定することが私にはできなかったので諦めてエクセルでt検定をした。

ただし，以下のようなデータ型にすればRを使ってt検定できた。

f:id:LPasteur:20190416020653p:plain

今回は同じサンプルを0 hと12 hで測定しているのでデータは対応していることになる。
また阻害剤を添加したので，遺伝子の発現量は減少するかどうかを判断したいので片側検定だ。

t.test(demo_3$X0h, demo_3$X12h, alternative="greater", paired = T)

f:id:LPasteur:20190416021127p:plain

p値は2.852e-07と求められた。エクセルで求められた値と一致。

両側検定や2標本が不等分散するサンプルの検定をしたい場合は，以下のサイトが参考になる。
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/62.html

求めたp値の情報をグラフに盛り込む。

plot <- ggplot(demo_mean_sd, aes(x = group, y = mean, fill = time)) +
geom_errorbar(errors, position = position_dodge(0.9), width = 0.1) +
geom_bar(stat = "identity", position = "dodge", colour = 1) +
scale_fill_manual(values = c("gray", "white")) +
labs(x = "", y = "Expression level [copies/ng]") +
theme_few() +
theme(axis.title = element_text(size = 14, family = "Helvetica"), axis.text = element_text(size = 14, colour = 1, family = "Helvetica")) +
scale_y_continuous(expand = c(0,0), limits = c(0,8)) +
geom_signif(y_position=c(7, 7), xmin=c(0.8, 1.8), xmax=c(1.2, 2.2), annotation=c("p < 0.001", "p < 0.05"), tip_length=0)

plot