【R】遺伝子発現をggplot2で論文体裁な棒グラフ作成 - ②

35708 ワード

テキストリンク

はじめに

前回の記事の続き的な記事になると思います。

今回のR script作成にあたって、Aesthetic　errorsなるものに相当苦しめられました。これを回避するために、何度もdataframeを定義し直すことになってしまって、少しコードが汚いです。

どなたか上手いやり方を教えてください。

目標

記事のコンセプトとしては、論文体裁な遺伝子発現量のグラフを、Rだけで作成できるようになることです。

より具体的には、遺伝子発現（Real-Time qPCR）のcsvデータをRStudioに読み込ませて、よく論文で見るようなシンプルな見た目かつ、エラーバーや有意差検定（多重比較）の結果まで載っているfigureの作成ができるコードを書いていきます。

今回の記事では、こんな感じのグラフを描画できるようになります。
前回のものと比べるとマイナーチェンジに見えるかもしれませんが、論文図表感がめちゃくちゃ増してませんか？？

実験条件

（前回の記事と全く同じなので、読んだ方は読み飛ばしてください。）
#使用する実験データ　←ページ内リンク

今回のfigは、植物が特定のストレスを受けた際に遺伝子発現をどのように変化させているのか、棒グラフで示します。

本実験は、モデル植物（シロイヌナズナ）のストレス応答について調べます。

植物は動物と違い、環境ストレス(乾燥、温度、塩、光など)に晒されても、ストレスのない快適な場所へ移動することができません。そのため植物はストレスを感知して、植物体内の生理変化（植物ホルモンや遺伝子発現の調節）を引き起こしてその場で耐えます。

BioRenderさんで作成

その遺伝子発現量がストレス処理前後でどう変化しているのかを調べるため、特定のストレス処理を行った植物体からRNAを抽出、cDNA合成、そしてストレス誘導性遺伝子をマーカーとしてRT-qPCRを行いました。

使用する実験データ

以下、今回のグラフ作成に使用するcsvデータになります。

前回のものと少し違うので注意してください！

sample-time-gene.csv

sample,time,gene
A,0,0.043122871
A,0,0.063175676
A,0,0.032600745
A,3,0.992119453
A,3,0.953639314
A,3,1.562149158
B,0,0.050309051
B,0,0.040279368
B,0,0.038257683
B,3,2.883036936
B,3,2.791183295
B,3,3.036106124

データの見方

これも前回の記事とほぼ同じですが、念の為に書いておきます。

sample-time-geneのcolumn（列）は、ハイフン（-）を挟んで左側がサンプル名（ecotype名）、真ん中がストレスの処理時間、右側がストレス応答性遺伝子（gene_Pとしましょう）の発現量です。厳密に言うと、ハウスキーピング遺伝子（ACT2）の発現量でNormalizeしたrelative expressioionです。

一番下のrow（行）を例に見てみると…

B,3,3.036106124

これは、エコタイプBの植物にストレスを3日間処理した結果、gene_Pの発現量が、ACT2に比べて約3倍だったということです。

要求package

ggplot2
dplyr
multcomp

動作環境

以下の環境で動作確認済みです。

MacBook Air(M1, 2020), macOS Monterey 12.0.1
RStudio, "Ghost Orchid" Release for macOS Mozilla/5.0 (Macintosh; Intel Mac OS X 12_0_1)

1. ワーキングディレクトリの設定とパッケージのロード

例によって、忘れないうちにやっておきます。

ワーキングディレクトリの設定.R

setwd("~/Desktop/Zenn/bar-graph_2")
#グラフ化したいデータのcsvファイルが入っているディレクトリをワーキングディレクトリに設定

僕はDesktopにおいてある、Zennフォルダ下のbar-graph_2フォルダに、sample-time-gene.csvを置いていたので、そこにsetwd()で仕事場を指定。

続いて、パッケージのロードのためにlibrary()を実行します。今回ロードするパッケージは以下3つ。

ggplot2
dplyr
multcomp

未インストールの方はinstall.package()であらかじめインストールしておきましょう。（インストールは初回のみ）

パッケージのロード.R

library(ggplot2) #グラフの描画に使う
library(dplyr) #データの加工に使う
library(multcomp) #検定に使う

2. csvファイルの読み込みと加工

次に実験データの入っているcsvファイルをRStudioに読み込ませます。

csvファイルの読み込み.R

d <- read.csv("sample-time-gene.csv", header = T)

上の操作では、read.csv()を用いて読み込んだcsvファイルを、dという文字に代入する操作をしています。（dはなんの文字を指定しても良いです。僕の場合は、短い方が楽なのでdataの頭文字を指定しています。）

header = Tはデータの１行目がサンプル名とか、数値データじゃない時に指定してあげます。１行目から数値データの時にはheader = Fを指定します。

3. 検定用のデータフレーム作成〜Tukeyの検定

まず、グラフに描画する検定結果を先に計算します。

グラフの描画を、ストレス無処理/処理群で色分けしたいのでわざわざsample（ecotype名）とtime（ストレスの処理時間）でcolumnを分けてデータを作っています。しかし実のところは、統計解析をするにはまとまっていた方が都合が良いです。

というわけで、dplyrのmutate()関数をつかってsampleとtimeをくっつけた新しいcolumnを作ります。

検定用データフレーム作成

統計検定用のデータフレーム作成.R

d1 <- mutate(d, sample_time = paste(d$sample, d$time)) #1行目
d_tukey <- data.frame(sample_time = d1$sample_time, ex = d1$gene)　#2行目
d_tukey$sample_time <- factor(d_tukey$sample_time)　#3行目

1行目のmutate()関数のイメージ

2行目: d_tukeyに、新しく作成したsample_timeと、遺伝子発現の定量値をexとして格納しました。

3行目: 後に使用する関数がsample_timeを正しく認識できるようにデータの持つ属性を編集しています。

d_tukeyを実行して、データを確認してみます。↓

ここまでで、念のためにTukey検定ができるかも確認しておきましょう。

一応確認しておく.R

TukeyHSD(aov(ex ~ sample_time, data = d_tukey))

ちゃんとできていました〜

Tukeyの検定

あとは、【R】multcompで多重比較結果をアルファベット表示させるでやったように、検定結果をアルファベットで出力します。

アルファベット.R

res <- aov(ex ~ sample_time, data = d_tukey)
tuk <- glht(res, linfct = mcp(sample_time = "Tukey"))
mltv <- cld(tuk, decreasing = F)
abc <- mltv[["mcletters"]][["Letters"]]

以前の記事では書いていませんでしたが、実はこの検定結果（mltv）は、["mcletters"]][["Letters"]という場所にアドレスされています。

なので、abc <- mltv[["mcletters"]][["Letters"]]これを実行してあげることで、検定の結果のアルファベットを、abcという函に格納することができます！（これをグラフ描画に使います。）

4. グラフ描画準備

続いて、グラフ描画用のデータフレームを作っていきます。

僕たちが作りたいグラフの作成に必要な情報（Matrixで用意したいもの）は、

sample名（ecotype AとかBとか）
ストレス無処理/処理の条件（0 dayとか3 daysとか）
遺伝子発現の平均値
標準偏差（エラーバーつけるため）
有意差を示すアルファベット

dplyrのgroup_by関数とsummarise_all関数

まずはグラフとエラーバーの描画のために、遺伝子発現の定量値について、平均値と標準偏差を算出しておきます。

dplyrの機能、group_by()関数と、summarise_all()関数を使いましょう。

statistics.R

d_stat <- d %>% 
  group_by(sample, time) %>% 
  summarise_all(list(mean = mean,
                     sd = sd))

上のコードでは、d1のデータをgroup_by()関数でsample、timeごとにまとめてsummarise_all()関数に渡し、サンプルごとの平均値（mean）と標準偏差（standard deviation）を算出したものをd_statに格納しています。

d_statの実行結果
こんな感じのテーブルになってます。

グラフ描画のためのデータフレーム作成

d_statをもう少し加工して、グラフ描画に使える一つのデータフレームを作ります。
（美しくないコードを紹介するのは恥ずかしいですが、上手いやり方が見つからないのでこんな感じになってます。どなたかシンプルにできる方いませんか…）

graph_matrix.R

d_graph <- data.frame(d_stat$sample, d_stat$time, abc)
d_graph <- mutate(d_graph, sample_time = paste(d_stat.sample, d_stat.time)) 
d_graph <- data.frame(sample = d_graph$sample_time,
                     time = as.character(d_stat$time),
                     mean.ex = d_stat$mean,
                     sd.ex = d_stat$sd,
                     abc = d_graph$abc)

この時timeには、as.character()関数を使って文字属性を付与しておきます。（数値としてグラフ描画を行うと、少しおかしなことになります。）

d_graphを確認してみます。

ちゃんとできています！

エラーバーの準備

グラフにエラーバーを描画させたいので、aes()関数で、エラーバーに“長さ”という属性を付して、さらにこれをerrorという函に格納しておきます。

error.R

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

これでようやく、グラフ描画に盛り込みたい要素を全て用意することができました〜！

5. グラフ作成

早速ggplotでグラフを作成していきます。

bar-graph.R

g <- ggplot(d_graph, aes(x = sample, y = mean.ex, fill = time)) + 
  geom_bar(stat = "identity", width = 0.6) +
  geom_errorbar(errors, width = 0.1) +
  scale_fill_brewer(palette = "Greys") +
  stat_summary(geom = 'text', label = abc, fun = max, vjust = -5) +
  ggtitle("gene P") +
  scale_y_continuous(limits = c(0, 3.5), name = "Relative expression level") +
  theme_classic() +
  theme(axis.title.x = element_text(size = 10, family = "Arial"), 
        axis.title.y = element_text(size = 10, family = "Arial"), 
        axis.text.x = element_text(size = 10, colour = 1, family = "Arial"), 
        axis.text.y = element_text(size = 10, colour = 1, family = "Arial"))

plot(g)

plot(g)の出力結果

ggplot()
どのデータセット（データフレーム?）を使うのか、縦、横軸に何をセットするかを指定します。
geom_bar
グラフの種類によって変更する。geom_lineとかgeom_histogramとか色々あります。
geom_errorbar
エラーバーの描画に使いました。
scale_fill_brewer()
グラフの色を指定することができます。今回はGreysでシンプルにしてみました。
RNA-seqなどのヒートマップにはpalette = "Spectral"で青〜赤のグラデーションにしてあげるとよいと思います。
stat_summary()
今回は、検定結果を表示させるのに使用しました。エラーバーとアルファベットが重ならないように、vjustで位置の指定をします。
ggtitle
グラフタイトルをつけられます。
scale_y_continuous()
デフォルトだとバーが一番高いところまでグラフ描画をしてくれます。今回は検定結果のアルファベットと、エラーバーが見切れてしまったので、見切れない高さまでlimit()で描画してもらいます。
theme()
フォント、フォントサイズの指定など、ディティールを指定します。論文体裁にしたかったので、フォントはArialで、あとは適当にやりました。

おわりに

今回は、Rを使ってデータの読み込み〜検定〜グラフ出力を行いました。

全ての工程において、なるべく掘り下げて易しく書きました。理解の助けになれば幸いです。もっと説明、プレゼン上手になれるように頑張ります。

実験結果

一応、実験条件についてお話したので、今回の個人的な出力データの解釈についても書いておきます。

ecotype AもBも、ストレス無処理のサンプルに比べて、ストレス処理3日後のサンプルの方でgene Pの発現レベルが有意に増加しており、その増加のレベルは、ecotype AよりもBの方が大きいです。

gene Pは、シロイヌナズナのある特定のストレス応答のマーカー遺伝子として用いられる、ストレス応答シグナリングの下流因子で、この遺伝子の発現上昇はストレスに耐えようとしたかどうかの結果と解釈できます。（遺伝子発現は秒単位で起こるのに、ストレス処理後3日間もおいているのはgene Pがシグナル伝達経路のほぼ最下流だから）

つまり、ecotype A, Bともにストレス応答をしていますが、その程度はBの方がより高いことが示唆される結果となりました。

今回の全R script

シンプルにできる方、どなたかその方法についてご指摘、ご教授くださると嬉しいです。

bar-graph.R

setwd("~/Desktop/Zenn/bar-graph_2")

library(ggplot2)
library(dplyr)
library(multcomp)

d <- read.csv("sample-time-gene.csv", header = T)

d1 <- mutate(d, sample_time = paste(d$sample, d$time)) 
d_tukey <- data.frame(sample_time = d1$sample_time, ex = d1$gene)
d_tukey$sample_time <- factor(d_tukey$sample_time)

res <- aov(ex ~ sample_time, data = d_tukey)
tuk <- glht(res, linfct = mcp(sample_time = "Tukey"))
mltv <- cld(tuk, decreasing = F)
abc <- mltv[["mcletters"]][["Letters"]]

d_stat <- d %>% 
  group_by(sample, time) %>% 
  summarise_all(list(mean = mean,
                     sd = sd))

d_graph <- data.frame(d_stat$sample, d_stat$time, abc)
d_graph <- mutate(d_graph, sample_time = paste(d_stat.sample, d_stat.time)) 
d_graph <- data.frame(sample = d_graph$sample_time,
                     time = as.character(d_stat$time),
                     mean.ex = d_stat$mean,
                     sd.ex = d_stat$sd,
                     abc = d_graph$abc)

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

g <- ggplot(d_graph, aes(x = sample, y = mean.ex, fill = time)) + 
  geom_bar(stat = "identity", width = 0.6) +
  geom_errorbar(errors, width = 0.1) +
  scale_fill_brewer(palette = "Greys") +
  stat_summary(geom = 'text', label = abc, fun = max, vjust = -5) +
  ggtitle("gene P") +
  scale_y_continuous(limits = c(0, 3.5), name = "Relative expression level") +
  theme_classic() +
  theme(axis.title.x = element_text(size = 10, family = "Arial"), 
        axis.title.y = element_text(size = 10, family = "Arial"), 
        axis.text.x = element_text(size = 10, colour = 1, family = "Arial"), 
        axis.text.y = element_text(size = 10, colour = 1, family = "Arial"))
plot(g)

References

たくさん参考にさせていただきました。ありがとうございます。

Author And Source

この問題について(【R】遺伝子発現をggplot2で論文体裁な棒グラフ作成 - ②), 我々は、より多くの情報をここで見つけました https://zenn.dev/taki_san/articles/410c66eab338fa

著者帰属：元の著者の情報は、元のURLに含まれています。著作権は原作者に属する。

Collection and Share based on the CC protocol

Ubuntu22.04にEclipseを入れて日本語化する

Flutterで動画収録時にiOS不具合を回避する方法