Japanese Journal of Lactic Acid Bacteria
Copyright © 2022, Japan Society for Lactic Acid Bacteria

解説

次世代シーケンサーデータの解析手法 第 19 回
R Markdown
牧野 磨音 1、清水 謙多郎 1, 2, 3、門田 幸二 1, 2, 3*
1

東京大学 大学院 農学生命科学研究科
東京大学 大学院 情報学環・学際情報学府
3
東京大学 微生物科学イノベーション連携研究機構
2

マークダウン（Markdown）は、HTML を手軽に生成するために開発された軽量マークアップ言語である。R Markdown は、R のコマンドと Markdown を組み合わせたものであり、近年主流となっている R の実行手段である。多くのユーザにとっては見慣れないファイルの拡張子（.Rmd）かもしれないが、アグリバイオインフォマティクス教育研究プログラムの中でも利用されている。本稿ではまず、R
Markdown の基本的な利用法として、チャンクの概念や HTML の生成などを述べる。次に、前回作成した R スクリプトファイルの内容を R Markdown 化し、相違点について述べる。最後に、個々の R コマ
ンドや実行結果として得られるオブジェクトについて解説する。ウェブサイト（R で）塩基配列解析のサブ（URL: http://www.iu.a.u-tokyo.ac.jp/kadota/r_seq2.html）中のウェブ資料（以下、W）を併用してほしい。
Key words：R Markdown, normalization, RNA-seq, RStudio

はじめにない。HTML ファイルをテキストエディタで眺めると、
例えば一番大きな見出しは <h1></h1>、イタリックは

今回も、前回に引き続き管理者として RStudio を起動し、<i></i> のタグで囲まれていることがわかる。このようなデスクトップ上に作成した hoge フォルダを作業ディレク

タグを利用して見出しなどの全体的な構造を表現するのがトリにした状態からスタートする（W01）。R Markdown

マークアップである。

の使用感は、ウェブブラウザでの Python 実行環境である

マークダウンは、マークアップの表記法を簡単にした

Google Colaboratory（ 以 下、Colab）や Jupyter Notebook

ものであり、軽量マークアップ言語というカテゴリに属すとよく似ている。それゆえ、これらを知っている読者は、る。マークダウンの記法で書かれた文書は、HTML（や前回の内容よりも今回の R Markdown のほうがわかりやす

PowerPoint や Word）に簡単に変換することができる。もいかもしれない。また、Colab などを使ったことがない（がちろんマークダウン特有の記法を覚える必要はあるものの、興味はある）読者も、R Markdown を学ぶことによって、

GitHub や Qiita などのウェブサービスでも広く普及してい利用上のハードルが低くなることが期待される。

る。例えば我々が提供している TCC-GUI

1）

や MBCdeg 2）の GitHub サイトにアクセスすると、README.md というマークダウンという言葉に馴染みがない読者は、HTML

マークダウンファイルの中身が表示される。

の 正 式 名 称（HyperText Markup Language）に 含 ま れ

R Markdown を利用する一番のメリットは、R コードだ

るマークアップという言葉と対比させるとよいかもしれ

けでなくその実行結果も含めて HTML で保存できる点に＊ある。例えば、前回 3）の R スクリプトファイル（JSLAB18.

To whom correspondence should be addressed.
Phone ： +81-3-5841-8155
Fax ： +81-3-5841-1136
E-mail ： koji.kadota@gmail.com

R）は、実行に必要な情報しか含んでいない。我々の経験上、このスクリプトを数年後に再度実行すると、何らかの不具合が発生することが予想される。しかし、HTML で195

Vol. 33, No. 3

Japanese Journal of Lactic Acid Bacteria

保持しておくと、たとえ数年後に作成当時のコードがうま配慮しきれないのが現実である。今回インストールした

く動かなくとも、当時どのようなコマンドでどのような結

rmarkdown ver. 2.14 は、要件として R 本体が ver. 3.0 以果が得られていたかがわかる。

上、pandoc が ver. 1.14 以上だと記されている（W03）。この要件自体は、本稿執筆時点の最新版（2022 年 5 月末時点

パッケージやソフトウェアの依存関係

で R ver. 4.2.0 および pandoc ver. 2.18）との比較から、全く厳しいものではないことがわかる。つまり、インストー

RStudio では、様々なファイル形式を変換するためのソフ

ルで躓くことはほぼない。

トウェアである pandoc を内部的に用いて、R Markdown のR パッケージ同士の依存関係で実際に躓く例としては、

HTML などへの変換を実現している。RStudio 上で pandoc

R ver. 4.0.5 での TCC

4）6）が挙げられる。TCC は内部的に7）を実際に利用しているのは、knitr という R パッケージでDESeq2ある。R Markdown 本体に相当するのが rmarkdown 5）といいている locfit の要件が R ver. 4.1.0 以上となっているたを用いているが、DESeq2 がさらに内部的に用う R パッケージであり、rmarkdown 自体は knitr の機能をめである。これは 2022 年初頭に発生していたエラーであ

内部的に用いている（つまり依存関係がある）。このような場り、バイオインフォマティクス関係の Q&A サイトである合、rmarkdown のみをインストール対象として指定すれば、Biostar 8）でも 2022 年 4 月末にこのエラーに関するやりと依存関係にあるパッケージ（この場合は knitr）も自動的にイりが行われている。前回 3）の冒頭部分で「R ver. 4.1.0 以前

ンストールされる（W02）。なお、管理者として RStudio をのもので不具合が生じた場合は、最新版をインストールし起動するよう強調していたのは、一般ユーザとしてパッケーて再度試してほしい。」と書いたが、その理由は TCC のイジのインストールを行った場合に、一部の PC 環境で不具合

ンストールエラーを想定したものである。

が生じるためである。
R Markdown の基本的な利用法（W04）

Pandoc 自体は独立したソフトウェアであるため、インストールも独立して行わねばならない。バージョンアップも独自に行われるため、本稿作成時のバージョン（ver.

図 1 は、R Markdown ファイル の 新 規 作 成 例 で ある。

2.18）と読者が実際に試すバージョンは異なりうる。そRStudio 上では、「① File →② New File →③ R Markdown」

して依存される側（pandoc）は、依存する側（knitr）にはという手順で行う。④最初にこのドキュメントの Title（デ図 1．R Markdown の新規作成

196

Jpn. J. Lactic Acid Bact.

Vol. 33, No. 3

フォルトは Untitled）と Author（デフォルトは空欄）を入

RStudio は R だけでなく Python も実行することができる。

力するように促される。R Markdown ファイルの変換はこれらを区別できるようにするため、チャンクの開始行

HTML と PDF と Word が選択可能であるが、基本的にデ（この場合は 8 行目）で明示的に宣言しているのだと解釈

フォルトのままでよい（⑤ HTML で⑥ OK ボタンを押す）。

すればよい。

R Markdown は、R の コ ー ド と Markdown 形 式 の 文

図 3 は、①図 2 で作成した新規 R チャンク内に簡単な書 を 組 み 合 わ せ た も の で あ る。R コ ー ド は、 チ ャ ン ク

数値計算を行うコマンド（10 ～ 11 行目）を書き込んで実（chunk）と呼ばれるブロック単位で実行される。例えば、

行（②の部分をクリック）し、③チャンク外に Markdown

第 18 回の図 3 では、「複数のパッケージのインストールをを 2 行分追加した R Markdown である。④ test19.Rmd と行う 2 ～ 9 行目のコマンドを選択範囲として反転させたのいう名前で保存したのち、⑤ Knit（ニットと読む）ボタン

ち、Run ボタンを押して実行」した。マウス操作の場合はを押すと（これをレンダリングという）、⑥ test19.html が誤った選択範囲で実行してしまうことがあるが、これらの生成されると同時に、その中身も表示される（中央下部のコマンド群を最初から 1 つのブロックとして取り扱える枠内）。③の 7 行目が⑦に対応することからもわかるようようにするのがチャンクである。このやり方は Colab やに、行頭の # は Markdown では見出しを意味する。

Jupyter Notebook でも採用されている。

HTML ファイルの表示結果には、著者や日付、R スク

図 2 は、図 1 で新規作成した R Markdown の 8 行目以

リプトとその実行結果が含まれる。③チャンク外に記載

降の記述内容を一旦削除したのち、8 行目にカーソルがあする Markdown として⑧のような結果の考察なども書きる状態で①および②をクリックして新規 R チャンクを挿こむことができるため、これで 1 つのレポートとして完入し終えた状態である。③作成された新規 R チャンクは 8

結させることができる。教員側（レポートを評価する側）～ 10 行目の部分に相当し、8 行目が開始行、10 行目が終も、R スクリプトファイルと実行結果を別々に提出され了行である。9 行目の部分が実際に実行させたいコードをるよりハンドリングが楽である。これがレポート作成と書き込んでいく部分である。8 行目に {r} という記述があ絡めて R Markdown がよく紹介されるゆえんである。なるが、これは 9 行目の記述内容が R コードだと認識させお、第 18 回の R スクリプトファイルでは、# はコメント

るためのおまじないのようなものだと解釈すればよい。④を意味していた。しかし R Markdown では、チャンク外で Python が選択肢にあることからも予想できるように、の文書は基本的に全てコメントと同様の扱いになるため、

図 2．新規 R チャンク挿入後の状態

197

Vol. 33, No. 3

Japanese Journal of Lactic Acid Bacteria

Markdown 中の # の有無は、見出しかそうでないかの違

当該チャンク内に 4 行分のコードが書かれていて、
「実行

いを表す。

結果は全て含めるが 1 行目と 4 行目のみコードは含めない」
場合は、echo＝2:3 とすればよい。これは、コード部分のJSLAB18.R の一部を Rmd 化（JSLAB18.Rmd の作成）

1 行目が FALSE、2 ～ 3 行目が TRUE、そして 4 行目がFALSE という echo オプションの指定だと解釈してもよい。

より実践的な R Markdown を作成すべく、ここでは第第 18 回のスクリプトファイル（JSLAB18.R）では、「サン

18 回で実行した「パッケージのロード」と「サンプルのク

プルのクラスタリング」を 3 ステップ（データの読み込み、ラスタリング」に焦点を絞って解説する（W05）。まず前者

クラスタリング本番、作図）に分けている。Markdown のについては、R チャンク内にどのパッケージをロードする

見出しの観点では、これらは最上位のレベル 1 より下のレかが明記されている。また、パッケージロード時に表示さベル 2 に相当する。HTML ファイルでは文字の大きさが 1

れるメッセージが HTML レポートに含まれると、全体と段階小さくなるのみであるが、見出しのレベル分けは目次

して冗長になってしまうという側面もある。R Markdown（Table of Contents；TOC）の構成に直結するため、注意

では、いくつかの R チャンクオプションが用意されてい深く行ったほうがよい。

る。デフォルトでは R チャンク内のコードと実行結果がR Markdown では、Rmd ファイル上部の YAML ヘッ

全て HTML レポートに反映されるが、例えば include＝ダと呼ばれる部分で「toc: TRUE」とすれば、目次を表示

FALSE というオプションをつけることで、当該 R チャン

させることができる（W06）。デフォルトはウェブページ

ク内のコードと実行結果を HTML レポートから除外するの冒頭部分であるが、さらに「toc_ﬂoat: TRUE」を追加すことができる。

ることで、ページの左側にサイドメニューとして表示させ

他の R チャンクオプションとしては、「実行結果は含める

ることができる。図 4 は、これまでの内容を反映した①がチャンク内のコードは含めない」場合は echo＝FALSE、

JSLAB18.Rmd から② JSLAB18.html を生成した結果であ「実行時にチャンク内を R コードとして実行しない（評価しる。この Rmd ファイルの場合は、YAML ヘッダは 1 ～ない）」場合は eval＝FALSE とすることで、当該チャンク

9 行目に相当する。③サイドメニューの目次は、図 3 の全体を制御することができる。さらに細かい設定として、

5 行目を図 4 の 5 ～ 8 行目のようにすれば追加できる。RR チャンク内の一部のみを制御することもできる。例えば

Markdown は、通常の Markdown に「YAML ヘッダと R

図 3．単純な数値計算を含む R Markdown および HTML レポート生成例

198Jpn. J. Lactic Acid Bact.

Vol. 33, No. 3

チャンク」を組み込んだものという理解でもよい。

JSLAB19.Rmd 中の 29 行目以降が、MBCdeg2（TCC パッ
ケージ中の頑健な RNA-seq データ正規化法 DEGES 10）をJSLAB18.R の全部を Rmd 化（JSLAB19.Rmd の作成）

MBCdeg に組み込んだ方法）を実行する部分である。全部で約 100 行からなるため、MBCdeg2 全体でレベル 1 のさらに実践的な R Markdown として、JSLAB18.R の内見出しとし、それをさらに計 8 つに細分化して説明すべ

容を全て Rmd 化したものが JSLAB19.Rmd である。これ

くレベル 2 の見出しとした。こうすることで、R エディ

は、前節で作成した JSLAB18.Rmd の内容に加えて、遺伝

タ上で任意の見出しやチャンクにジャンプすることがで

9）子クラスタリング用パッケージである MBCluster.Seq をきる（W09）。実用上も、特定の場所に任意のチャンクや発現変動遺伝子（DEG）検出に転用した MBCdeg 2）の実行

Markdown を挿入したりするため、何行目から何行目とまで含めたものである。JSLAB18.Rmd との共通部分には、いった表現はせず、見出し名などで場所を特定するのが一

多少変更を加えている。例えば、パッケージのロード部分

般的である。

の R チャンクは、JSLAB19.Rmd でも include＝FALSE オ
ファイルの読込と subsetting（W10）

プションをつけているので HTML レポートには反映されない。このため、JSLAB18.Rmd の 10 行目に存在していた「レベル 1 の見出し」は削除した（W07）。

ここからは、MBCdeg2 実行におけるレベル 2 の見出し

JSLAB18.Rmd のサンプルのクラスタリング部分では、（R チャンク）単位で解説する。図 5 は、ファイルの読込と見出しのレベルの違いと HTML レポート中の目次の見え

subsetting（サブセットを得る作業のこと）を行うチャンク方の違いがわかるように、レベル 2 の見出しをつけた 3 つの実行結果である。① JSLAB19.Rmd 中の、②のチャンクのチャンクに分けて実行した（W08）。実際には、この程では、JSLAB18.xlsx 11）を読み込んだ結果を data_all という度の分量のものは 1 つのチャンクにまとめて実行する。こオブジェクトに格納している。読み込んだ直後の③ data_の R チャンクでは、オプションとして「fig.width＝5, fig.

all は、2,949 遺伝子×9 サンプルの数値行列である（第 18height＝4」を指定している。この数値の単位はインチで回の図 5）。但し、単純化のために「酸ストレス長期暴露群あり、クラスタリング結果の樹形図を plot 関数で描画す（pH4.5_24h）vs. 対照群（pH7_CCG）」の 2 群間比較を行いる際に利用している。

たいので、元の 9 列分のデータから、4 ～ 9 列目のサブセッ

図 4．JSLAB18.R の一部を Rmd 化し HTML を生成した結果

199

Vol. 33, No. 3

Japanese Journal of Lactic Acid Bacteria

トを抽出した結果を data というオブジェクトに格納してるかを示す情報のことである。バイオインフォマティクス

いる。

業界ではクラスラベルと呼ぶのが一般的であるが、クラス

この種の行列データの抽出・加工・操作を行うモダンという表現が難解だと感じる読者に配慮して、ここでは群な手段は、tidyverse と呼ばれる R パッケージ群の利用でラベルとした。TCC 正規化とは、ここでは TCC パッケーある。実際我々も、TCC のウェブアプリ版である TCC-

ジが提供する頑健な正規化法である DEGES のことを指に、（R Markdown はもちろんのこと）tidyverse をす。実用上は、「TCC パッケージ（ver. X .Y .Z ）のデフォル構成するパッケージ群（dplyr, tidyr, and plotly）を利用しトオプションを用いて実行した。」のように書けば十分でた行列データの操作やインタラクティブな描画を実装してあるため、本稿ではアルゴリズムの詳細は省く。

GUI

1）

いる。しかし、本稿であまり沢山の事柄を盛り込みすぎて

図 6 は、群ラベル情報の作成と TCC 正規化を行うチャも理解が追いつかないため、このチャンクではベーシック

ンクの実行結果である。①群ラベル情報の作成は 39 行目

な subsetting を行っている。

に相当し、酸ストレス長期暴露群（pH4.5_24h）に 1、対照

34 行目では、2,949 行×9 列からなる data_all オブジェ群（pH7_CCG）に 2 というラベルを割り当てている。このクトの中から、4 ～ 9 列目のみ抜き出した情報を data とチャンクの目的は、MBCluster.Seq の入力として与えるたいう名前で保存している。任意の列を抽出したい場合は、めのサイズファクターと呼ばれるデータ正規化に用いる情④コンマ（,）の右側に抽出したい列情報を指定する。同様

報の取得である。但し、TCC や edgeR 12）は正規化係数とに、任意の行を抽出したい場合は、④コンマの左側に抽出呼ばれる値で取り扱うのが基本形である（43 行目の norm.

したい行情報を指定する。目的の 2 群のサブセットを抽出

factors）。それゆえ、得られた正規化係数をサンプルごとできているかどうかは、⑤ environment タブ上で見えての総カウント数（44 行目の colSums（data）に相当）に掛いる⑥ data オブジェクトを眺めることで確認できる。

けて有効ライブラリサイズ（44 行目の ef.libsizes）とよばれるものをまず作成し（44 行目全体）、それを有効ライブ

群ラベル情報の作成と TCC 正規化（W11）

ラリサイズの平均値で割ることによってサイズファクター（45 行目の size.factors）情報を得ている。

群ラベルとは、2,949 遺伝子×6 サンプルの数値行列でサイズファクターという概念は、edgeR と双璧をなすある data オブジェクト中のどのサンプルがどの群に属す

もう 1 つの有名な R パッケージである DESeq2 が使って

図 5．ファイルの読込と subsetting（チャンク 3）の実行結果

200

Jpn. J. Lactic Acid Bact.

Vol. 33, No. 3いるものであり、MBCluster.Seq もそれに準じている。正

させたい場合に利用する関数である。括弧の中の数値は、規化係数とサイズファクターは、ともにデータ正規化時に発生させる乱数種の番号だという理解でよい。たまたま今利用される情報という意味で似た概念である。しかし、②年が 2022 年なので 2022 という数値にしているが、999 や

43 ～ 45 行目で全体的な関係性が明確に示されているよう

27 など任意の整数であればなんでもよい。

に、両者が別物であるという点は正しく認識しておかねば
MBCluster.Seq 本番（W13）ならない。

図 8 は、K ＝3 での MBCluster.Seq の実行結果である。

MBCluster.Seq 実行の共通部分（W12）① 57 ～ 58 行目で、K の値と「遺伝子ごとの各クラスタへ
ここからは、MBCluster.Seq パッケージ内での解析となの属しやすさを表す事後確率の数値行列」を格納するためる。ここでは、① RNASeq.Data という関数を用いて、②の出力ファイル名（MBCdeg_K3.xlsx）情報を与えている。

data オブジェクト、③群ラベル情報である data.cl オブジェ

59 ～ 62 行目が MBCluster.Seq 実行のメイン部分であり、クト、そして先ほど算出したサイズファクター情報である

cls というオブジェクトに実行結果を格納している。63 ～④ size.factors オブジェクトの対数を入力として、⑤ hoge

66 行目では、cls の中から事後確率情報部分のみを PP とというオブジェクトを作成している（図 7）。hoge の中身いう名前で抽出し（63 行目）、その列名を変更して（64 行

は、MBCluster.Seq 内でその後の解析に利用していくため目）、その他の出力したい情報を合わせたものを（65 行の特有の情報格納形式となっている。入力として与えた以目）、XLSX 形式ファイルに書き出している（66 行目）。

上の情報にはなりようがないため、必要に迫られたときのもう 1 つの主要な結果である「クラスタ中心の発現パ

み内部構造を理解しにいくというスタンスでよいだろう。

ターン（代表パターン）」情報は、67 行目の実行結果とし

MBCluster.Seq は、似た発現パターンを示す K 個のグて② Console 画面上に表示される（但し②の画面上では隠ループに遺伝子を分類する目的で利用するのが基本形であれている）。ここでは③で当該部分のみ独立に示している

る。内部的には K -means クラスタリングという乱数を発ように、計 3 つの代表パターン（3 行×2 列の数値データ）生させるアルゴリズムを採用しているため、試行ごとに結情報が得られる。1 列目と 2 列目の数値の絶対値は同じで果が異なりうる。⑥ 50 行目の set.seed は同じ乱数を発生あり、2 列目の列名（つまり 2）が「対照群（pH7_CCG）」

図 6．群ラベル情報の作成と TCC 正規化（チャンク 4）の実行結果

201

Vol. 33, No. 3

Japanese Journal of Lactic Acid Bacteria

図 7．MBCluster.Seq 実行の共通部分（チャンク 5）の実行結果

図 8．K =3 での MBCluster.Seq（チャンク 6）の実行結果

202

Jpn. J. Lactic Acid Bact.

Vol. 33, No. 3

に対して割り当てたラベル情報である。行数は指定したクラスタに割り振られたかという数値ベクトル」であり、こラスタ数 K に対応し、例えば 1 行目がクラスタ 1（c1）のの元情報が図 8 の④で見えている cls$cluster というオブ代表パターンである。

ジェクトである。R では、このようなベクトルを入力とし2 列目を基準として考えると、1 行目（c1）の代表パターて実行し、⑤ベクトル中の要素の種類ごとにその出現回数

ンは、pH7_CCG 群が pH4.5_24h 群に比べて e2×（0.2996681）＝

を返す table という関数が用意されている。⑥で見えてい1.82091 倍高発現だと解釈する。同様にして、c2 は pH7_

るものが実行結果であり、例えば c3 に属する遺伝子数はCCG 群 が pH4.5_24h 群 に 比 べ て e2×（-1.7270041）＝0.03161865

1,424 個である。

倍高発現（31.626910 倍低発現）、そして c3 は e2×（-0.1801454）

K ＝4 および 5 の MBCluster.Seq の実行結果についても、＝0.6974735 倍高発現（1.433746 倍低発現）だと解釈する。

K ＝3 と同様の手順で実行することができる（W14）。重要

このような計算をする理由は、MBCluster.Seq が遺伝子ごな点は、図 8 の 59 ～ 68 行目に相当する部分が共通して利とのカウント値をその平均値で割り、さらにその対数を用されているということである。①K 値と出力ファイル

とった状態で取り扱っているためである。

名、そして⑦主要な実行結果である cls オブジェクトをそ

図 8 の③で見えている数値行列において、1 列目と 2 列れぞれの K 値由来であることがわかるように変更するこ

目の数値の絶対値が同じ理由は、平均値で割った値をそれとで、コードの大部分を使いまわすことができる。もちろ

ぞれの群に割り振っているからである。それゆえ、例えばんこれらをより効率的なコードにすることは可能である

c1 のクラスタ中心の群間の発現変動の度合いとして示したが、まずはこのような書き方からスタートしていくとよい

e2×（0.2996681）が難解だと感じた読者は、e|-0.2996681|＋|0.2996681| と読みかもしれない。

替えてもよい。なお、これらの実際の発現パターンは、第作図（元データの作成）（W15）

18 回の図 8 に示されているものと同じである。最も変動の少ない non-DEG クラスタは、視覚的にも c3 と判定できる。
MBCluster.Seq のような発現パターン分類の解析結果で

図 9 は、K ＝3 ～ 5 で得られたクラスタごとの代表パター

は、通常どのクラスタにいくつの遺伝子が割り振られたかン情報をまとめたオブジェクト matome を作成するチャが論じられる。この議論の基礎となる情報は、出力ファインクの実行結果である。この① matome オブジェクトは、

ル中の一番右側の列に示されている「どの遺伝子がどのク

30 行×4 列の行列データであり、ggplot2 というパッケー

図 9．
「作図（元データの作成）」チャンクの実行結果

203Vol. 33, No. 3

Japanese Journal of Lactic Acid Bacteria

ジを用いて描画する際の入力として用いられる。中央下部

×5 列からなり、②列がクラスタ番号、③行が K 値であの黒枠内が matome オブジェクトの最初の 12 行分に相当る。行と列を入れ替えたい（行列の転置をしたい）場合は、する情報である。② 1 列目（variable 列）はクラスタ番号④「k ~ variable」を「variable ~ k」とすればよい。122 ～情報であり、クラスタあたり 6 行分×計 5 クラスタ（K ＝5128 行目の記述の仕方が独特だと感じるかもしれないが、が最大値なので c1, c2, …, c5）なので、30 行分の情報とな

例えば 123 ～ 127 行目を削除した状態でも実行可能であるる。③ 2 列目（value 列）は代表パターンの情報であるが、ため、いろいろと試してみるとよいだろう。

④ 3 列目（k 列）と⑤ 4 列目（u 列）の情報を合わせると理おわりに解しやすいであろう。
例 え ば、 ⑥ 1 ～ 6 行 目 は、 ④K ＝3 ～ 5 の MBCluster.
Seq 実行結果の、②クラスタ 1（c1）の③代表パターン情

本稿では、R Markdown の基礎的な事柄からスタート報が格納されている。④からも想像できるが、1 ～ 2 行

し、前回の R スクリプトファイルの Rmd 化、そして前目が K ＝3 の結果、3 ～ 4 行目が K ＝4 の結果である。⑤回より踏み込んだ形で具体的な R コードの解説を行った。

は 群 ラ ベ ル 情 報 で あ り（ 図 6）、1 が pH4.5_24h 群、2 がマークダウン・チャンク・レンダリングなど聞きなれない

pH7_CCG 群である。重要な点は、フォーマットの詳細を用語に最初は戸惑うが、W04 以降を参考にしながら進め理解するというより、「確かにこれだけの情報が含まれてていけばすぐに慣れるであろう。

いれば、クラスタごと、K 値ごとの代表パターンを描画できるはず」と納得することであろう。

謝辞

作図（本番）（W16）

たものです。

本内容の一部は、JSPS 科研費 21K12120 の助成を受け

図 10 は、クラスタごと、K 値ごとの代表パターンを描

利益相反（COI）

画するチャンクの実行結果である。121 行目で matome オ

牧野磨音、清水謙多郎、門田幸二：本論文発表の内容にブジェクトを入力として与え、ggplot 関数を実行してい関連して開示すべき COI 状態はない。

るのがわかる。①右下の黒枠が得られる図である。3 行

図 10．
「作図（本番）」の実行結果

204

Jpn. J. Lactic Acid Bact.

Vol. 33, No. 3

		
参 考 文 献1）	Su W, Sun J, Shimizu K, Kadota K (2019) TCC-GUI: a
Shiny-based application for differential expression analysisof RNA-Seq count data. BMC Res Notes 12: 133
2）
	Osabe T, Shimizu K, Kadota K (2021) Differential expression
analysis using a model-based gene clustering algorithm forRNA-seq data. BMC Bioinformatics 22: 511.
3）	牧野磨音，清水謙多郎，門田幸二（2022）次世代シーケンサー
データの解析手法：第 18 回遺伝子発現データのクラスタリ
ング，日本乳酸菌学会誌 33: 87-94．
4）	Xie Y. (2022) knitr: A General-Purpose Package for Dynamic Report Generation in R. R package version 1.39, https://yihui.org/knitr/.
5）	Allaire J, Xie Y, McPherson J, Luraschi J, Ushey K, et al.
(2022) rmarkdown: Dynamic Documents for R. R package
version 2.14, https://github.com/rstudio/rmarkdown.
6）	Sun J, Nishiyama T, Shimizu K, Kadota K (2013) TCC: an R
package for comparing tag count data with robust normalization strategies. BMC Bioinformatics 14: 2197）	Love MI, Huber W, Anders S (2014) Moderated estima-tion of fold change and dispersion for RNA-seq data with
DESeq2. Genome Biol 15: 550.
8）	Parnell LD, Lindenbaum P, Shameer K, Dall’Olio GM, SwanDC, et al. (2011) BioStar: an online question & answer resource for the bioinformatics community. PLoS Comput
Biol. 7: e1002216.
9）	Si Y, Liu P, Li P, Brutnell TP (2014) Model-based clusteringfor RNA-seq data. Bioinformatics 30: 197-205.
10）	Kadota K, Nishiyama T, Shimizu K (2012) A normalization
strategy for comparing tag count data. Algorithms Mol Biol
7: 5.
11）	Bang M, Yong CC, Ko HJ, Choi IG, Oh S (2018) Transcriptional Response and Enhanced Intestinal Adhesion Abilityof Lactobacillus rhamnosus GG after Acid Stress. J Microbiol Biotechnol 28: 1604-13.
12）	Robinson MD, McCarthy DJ, Smyth GK (2010) edgeR: a
Bioconductor package for differential expression analysis ofdigital gene expression data. Bioinformatics 26: 139-40.

Methods for analyzing next-generation sequencing data
19. R Markdown.
Manon Makino1, Kentaro Shimizu1, 2, 3, Koji Kadota1, 2, 3
1

Graduate School of Agricultural and Life Sciences, The University of Tokyo.
2
Interfaculty Initiative in Information Studies, The University of Tokyo.
3
Collaborative Research Institute for Innovative Microbiology,
The University of Tokyo.

		Abstract

R Markdown is a recent mainstream way of executing R, consisting of R commands and a uniquenotation called Markdown. Its file extension is Rmd, which may not be familiar with many researchers,
but is also used in our graduate educational and research program. We first describe the basic usageof R Markdown, including the concept of chunks and HTML generation. We next convert the contentsof the previously created R script file into R Markdown and discuss the differences. Finally, we
describe individual R commands and objects resulting from their execution in detail. Supplementarymaterials are available online at https://www.iu.a.u-tokyo.ac.jp/kadota/r_seq2.html.

205