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

解説

次世代シーケンサーデータの解析手法 第 18 回
遺伝子発現データのクラスタリング
牧野 磨音 1、清水 謙多郎 1, 2, 3、門田 幸二 1, 2, 3*
1

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

遺伝子発現データは、各行に遺伝子、各列にサンプルが配置された数値行列のことを指し、各要素には対応するサンプル中の遺伝子がどれだけ働いているかを表す発現量の数値情報が格納されている。本稿では、第 15 回で得た 2,949 遺伝子×9 サンプルからなる Lactobacillus rhamnosus GG の酸ストレス応
答を調べた数値行列データを用いる。まず、データ解析環境 RStudio の基本的な利用法として、パッケー
ジのインストールやロードといった基礎的な事柄を述べる。次に、似た発現パターンを示すサンプルのクラスタリングについて、その意義や目的、そして結果の解釈について述べる。遺伝子のクラスタリン
グについては、RNA-seq 用の代表的なパッケージである MBCluster.Seq について、その概要と入出力形式を述べる。最後に、我々が最近開発した MBCluster.Seq の改良版という位置づけの MBCdeg 法について紹介する。ウェブサイト（R で）塩基配列解析のサブ（URL: http://www.iu.a.u-tokyo.ac.jp/kadota/r_
seq2.html）中のウェブ資料（以下、W）を併用してほしい。
Key words：clustering, RNA-seq, Lactobacillus rhamnosus GG, RStudio

はじめにの RStudio があればそれで実行してみてもよいが、R ver.
4.1.0 以前のもので不具合が生じた場合は、最新版をイン

RStudio は、フリーソフトウェア R の統合開発環境であストールして再度試してほしい。

る。大まかには、内部的に R を動かして、R 単体での利今回用いるのは、数値行列データ（JSLAB18.xlsx）と用よりもさらに便利で使いやすくしたものという理解でよR ス ク リ プ ト（JSLAB18.R）の 2 つ の フ ァ イ ル で あ るい。連載第 1 回 1）の図 2 では、R 単体の実行結果のみを示（W02）。本稿ではデスクトップ上に作成した hoge フォルしている。第 17 回 2）では、アグリバイオインフォマティ

ダ内にこれらのファイルが存在するという前提で話を進め

クス教育研究プログラムで R（RStudio を含む）を学べる

るが、RStudio に慣れている読者は適宜変更してもらって科目や R の位置づけについて述べている。しかし RStudioも構わない。

を用いた実際のデータ解析については、本連載の枠組み
作業ディレクトリの変更

ではこれまで解説してこなかった。従って本稿では、今回初めて R と RStudio をインストールし終えたばかりの RStudio 初心者が、管理者として RStudio を起動した

RStudio 起動後に、通常最初に行うのは作業ディレク状態を想定してスタートする（W01）。インストール済み

トリの変更である（W03）。この作業は、JSLAB18.xlsx といったような名前のみで指定したファイルを「デスクトッ＊プ上にある hoge フォルダ」内に限定して探索せよという

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

宣言を RStudio に対して行っているのだと解釈すればよい。こうすることで、中身が異なる同じ名前の複数のファ
イルが PC 内の他のフォルダ内に存在する場合でも、問題
94Jpn. J. Lactic Acid Bact.

Vol. 33, No. 2

図 1．作業ディレクトリ変更後の RStudio 画面

なく目的のファイルを一意に定めることができる。

は全て無視されるので、例えば # から始まる行は実行し

図 1 は、作業ディレクトリ変更後の RStudio の画面でてもしなくてもよい。スクリプトはソースコードとも呼ばある。① Files タブ上で見えている② 2 つのファイルが今れ、拡張子として .R がつけられる場合がほとんどである。

回用いるものである。作業ディレクトリの変更で行った一

中身自体はテキストファイルであるため拡張子を .txt にし連の操作は、③のコマンド実行と同義である。③の実行とても特に問題はない。しかし RStudio 上で中身が R スク

同時に④の部分も変更され、④の場所が作業ディレクトリ

リプトであることが明示されている .R ファイルとして開である。④と⑤は実質的に同じ場所であるが、⑤は単にくことで、文法的に問題がある箇所を警告してもらえるな

RStudio 経由で任意のフォルダを指定して表示させている

どのメリットがある。従って、R スクリプトファイルはだけである。しかし、初心者が遭遇するエラーは、「作業RStudio などの R 専用のエディタ上で開くことを強く推ディレクトリの変更忘れ」と「解析したいファイルが作業奨する。

ディレクトリ内にそもそも存在しない（JSLAB18.xlsx.txt

図 2 は、R スクリプトファイル（JSLAB18.R）を RStudioのような感じでファイルのダウンロード時に余分な拡張子

上で開いた状態である（W04）。作業自体は① Files タブ上がついてしまう不幸な事例を含む）」が原因であることがで見えている②当該ファイルをクリックするだけである。

多い。そのため、④と⑤が一致しているか、そして解析し③ JSLAB18.R という名前のタブが開かれ、ファイルの中身

たいファイルが②の場所にリストアップされているかに多が左側の領域の上半分を使って表示されていることがわかくの注意を払うことを勧める。

る。これが R のエディタ画面である。一方、図 1 で左側の領域全体を占めていて、図 2 で左側の下半分になった部分

エディタの起動は④コンソール（Console）と呼ばれるもの（正確にはペインとよぶ）である。③スクリプトファイル中の実行したい

R スクリプトとは、R 上で実行可能な一連のコマンドを行または領域を選択して、④ Console 画面上で実行し結果

1 つ 1 つ逐次に実行できる形でまとめたものである。多くを眺めるというのが基本形である。新規に出力ファイルをの場合、行ごとに実行可能であり、一定のまとまりごとに、生成するスクリプトを実行した場合は、① Files タブ上でそのまとまりがどのような役割を果たすのかを説明するコ正しく生成されたかどうかを確認することもできる。

メントがつけられている。シャープ（#）より右側の文字
95

Vol. 33, No. 2

Japanese Journal of Lactic Acid Bacteria

図 2．スクリプトファイル（JSLAB18.R）を開いた状態の RStudio 画面

パッケージのインストール

図 2 の左上に見えているスクリプトファイル（JSLAB18.
R）中 の 2 ～ 9 行 目 の コ マ ン ド は、 計 7 個 の パ ッ ケ ー

第 1 回 1）のおさらいを一部含むが、R と RStudio をイジ（BiocManager, TCC 4）, MBCluster.Seq 5）, openxlsx,

ンストールしただけの状態では機能が限られている。こggplot2, ggsci, and reshape2）のインストール作業に相当

れは、新しい PC を購入しただけの状態と同じだと考えする。前述のとおり、これらのコマンドの大部分、具体的ればよい。PC に Zoom や Slack などをインストールしてには Bioconductor から提供されている TCC 以外の 6 パッ

利用するのと同様、RStudio の機能を有効利用するためにケージのインストール作業は、⑤のポップアップメッセー

は、CRAN や Bioconductor 3）などから配布されているパッジに対して Install を選択することでまかなえる。しかし、ケージと呼ばれるものをインストールする必要がある。具

初心者の視点で見ると、どのパッケージがどのリポジトリ

体的には、今回入力ファイルとして用いる xlsx 形式のファ（CRAN や Bioconductor のようなパッケージ提供サイト

イルを RStudio で読み込むためには、例えば CRAN からのこと）から提供されているかを把握するのは難しい。ま

提供されている openxlsx というパッケージを予めインスた、この 2 ～ 9 行目のコマンドは、リポジトリに依存せずトールしておかねばならない。

統一的にパッケージをインストール可能なやり方であるた

一部の読者は、PDF ファイルを開くためのソフトウェ

め、今回はこのやり方で行う（⑤のポップアップは Don’t

アがインストールされていない状態で、手元にある PDF

Show Again をクリックする）。

ファイルを開こうとした経験があるのではないだろう

図 3 は、パッケージのインストール作業終了後の状態でか。そのような場合、Acrobat Reader をインストールしある（W05）。① JSLAB18.R 中の 2 ～ 9 行目を反転させてますか？のようなメッセージが出て、ほぼ無意識に「は② Run ボタンを押す作業は、③コンソール画面上でのこれい」を選択するようなステップを経たはずである。最近のら一連のコマンド実行に相当する。終了したかどうかは、

RStudio は高機能化が進んでおり、開いたスクリプトファ④の部分がコマンド入力待ち状態（コマンドプロンプトとイルを調べてインストールが必要なパッケージを提示し、呼ばれる > が出ている）かどうかで判断できる。⑤のフレー

インストールするかどうかを問い合わせてくれるようにズ（「reshape2 は無事に…」）から、⑥の最後のコマンドまなっている。それが、図 2 の⑤の部分に出現しているメッで成功裏に終わったことが示唆される。⑥のコマンドは、セージである。

「BiocManager パッケージが提供する install という関数を96

Jpn. J. Lactic Acid Bact.

Vol. 33, No. 2

図 3．パッケージのインストール完了後の RStudio 画面

用いて、reshape2 というパッケージを update＝F というオある点が挙げられる。特に後者はそのまま進めて問題ない

プションで実行せよ」という意味である。尚、⑤の MD5 サ

かどうかの判断自体が初心者には難しいが、② 2 回目（20ムは、第 3 回 6）でも解説しているが、ダウンロードしたファ～ 25 行目）の実行結果として、③のように 2 回目で何もメッ

イル（この場合は reshape2）が提供元と完全に一致していセージが表示されていなければ基本的に問題ないと判断しるかどうかを確認するためのものである。

てよい。ここまでで RStudio の解析環境構築は完了である。

パッケージのロード

データの概要パッケージのロードは、インストールしたパッケージを今回入力データとして用いるファイル（JSLAB18.xlsx）

読み込んで、使える状態にする作業である（W07）。大まは、第 13 ～ 15 回で解説した L. rhamnosus GG（LGG）の

かには、PC にインストールしたソフトウェアの起動のよ

酸ストレス応答を調べたカウントデータである 7-9）。この研うなものだという理解でよい。R の場合は、複数個のパッ究では、以下に示す計 3 状態の RNA-seq データを各状態

ケージを一度にインストールすることが多いため、図 3 左につき 3 反復ずつ取得して比較している 10）。他分野の読者下のコンソール画面の結果のみでは全てのパッケージが無は、例えば「低濃度抗がん剤投与群、高濃度抗がん剤投与

事インストールできているか確認しきれない。それゆえ、群、対照群」のように解釈してもよいし、よりシンプルに

我々は通常、当該パッケージ群のロードを通じてインス「G1 群、G2 群、G3 群」の 3 群間比較だとみなしてもよい。
・酸ストレス短期暴露群（pH4.5_1h）

トールの成否を確認する。
図 4 は、スクリプトファイル（JSLAB18.R）中の 12 ～ 25・酸ストレス長期暴露群（pH4.5_24h）・対照群（pH7_CCG）

行目のコマンド実行結果である。① 11 行目と② 19 行目のコメントからも想像できるが、同じ作業を 2 回行っている
ことがわかる。この理由は、① 1 回目（12 ～ 17 行目）の実

図 5 は、JSLAB18.xlsx の概要である。1 行目は列名の行のみでは、大量のメッセージが表示されて判読しづらい

情報からなるヘッダー行、① 1 列目は行名（遺伝子 ID）情こと、実用上問題ない警告メッセージが表示されることが報から構成されていることがわかる。これが典型的な入力
97

Vol. 33, No. 2

Japanese Journal of Lactic Acid Bacteria

図 4．パッケージのロード後の RStudio 画面

図 5．入力データ（JSLAB18.xlsx）の概要

ファイルの形式であり、実際のデータは 2,949 行×9 列のることから妥当である。重要な点は、本格的なデータ解析

数値行列である。② pH4.5_1h_1 というサンプル名のデーを進めていく際に、これらの事象が何らかの悪影響を与え

タが他に比べて全体的に数値が小さいことに気づくが、こている可能性を排除しないことであろう。

れは第 13 回の表 1 を見れば納得できる。つまり、このデー
サンプルのクラスタリング

タは、生のリード数の点で他よりも明らかに少ないためであろう。逆に、③ pH4.5_24h_2 が全体的に大きな数値となっているのは、リード数やファイルサイズが最大値であ

我々は通常、探索的データ解析の一環として、まずサ
98

Jpn. J. Lactic Acid Bact.

Vol. 33, No. 2

図 6．サンプルクラスタリング実行後の RStudio 画面

ンプルのクラスタリングを実行する（W08）。これは、全の状態で見れば容易に理解できるだろう。

体的なサンプル間の類似度を評価する枠組みであり、前図 7 は、サンプルのクラスタリング実行結果ファイル

述のような全体的な数値の大小の違いが同一群内の反復（fig7.png）の中身である。大きさがわかるように外枠を示

データ間の類似度に悪影響を与えていないかなどの観点しているが、図 6 の②で指定した縦横比（500×400 ピクセ

で結果を眺める。ここでは、TCC パッケージが提供する

ル）になっていることがわかる。①酸ストレス長期暴露群

clusterSample 関数を用いてサンプルのクラスタリングを（pH4.5_24h）は 3 つの反復データのみで 1 つのクラスタを行う。サンプル間の類似度を順位尺度である Spearman 相形成している一方で、②対照群（pH7_CCG）のクラスタ中

関係数で定めているため、全体的な数値の大小に影響されに③ pH4.5_1h_1 が紛れ込んでいるのがわかる。我々は、ないのが特徴である。

もし③のデータが他と比べて特に変わったことがなけれ

図 6 は、28 ～ 36 行目のコマンド実行結果である。コメば、図 7 の結果を素直に受け止めていたであろう。

ント行（28, 30, 33 行目）の記載内容からもわかるが、① 29

しかしながら実際には、③ pH4.5_1h_1 は他と比べて極行目で入力ファイル（JSLAB18.xlsx）を読み込み、② 34 行

端に全体的に数値が小さいことや、リード数が少ないとい目で出力ファイル名（fig7.png）を定めていることがわかる。

うことを確認している。また、クラスタリングアルゴリズ③実行時にエラーメッセージが出ていない（null device とムとしても、サンプル間の類似度を定義する数式として、いうのは特に気にしなくてもよい）こと、④作業ディレクト

ユークリッド距離のような数値ベクトル中の対応する要素

リに fig7.png が作成されていることもわかる。

間の差（の二乗や絶対値）を足していくような類のもので入力ファイル読み込み時に利用している read.xlsx 関数はないため、手法選択が間違っているとも考えにくい。まは、図 4 でロードした openxlsx パッケージが提供していた、もし③を除いて考えると、大きく「対照群 vs. 酸スト

るものである（から利用できる）。①で読み込んだ情報は、レス暴露群」の 2 つのクラスタとして分かれることや、酸

左向きの矢印の先にある data という名前のもの（これをストレス群の中でもさらに「④短期暴露群 vs. ①長期暴露

オブジェクトという）に格納されている。この中身の概要

群」に分かれている。これらの状況を鑑み、我々は③を外は、⑤ Environment タブで見えている⑥ data のところにれサンプル示されている（W09）。「2949 obs. of 9 varia…」の意味は、
このデータが 2,949 行×9 列の数値行列であることが既知99であると判断した。

11）Vol. 33, No. 2

Japanese Journal of Lactic Acid Bacteria

図 7．サンプルクラスタリング実行結果ファイル（fig7.png）

遺伝子のクラスタリング（MBCluster.Seq）

群 vs. G2 群）の結果として得られた 100 個の DEG のカウ
ントデータを入力として、K ＝3 として MBCluster.Seq を遺伝子のクラスタリングとは、似た発現パターンを示す

実行すると、出力結果として「計 3 つのクラスタ中心の発グループに遺伝子を分類する作業のことである。類似した

現パターン（3 行×2 列のデータ）」と「100 遺伝子×3 列の発現パターンを持つ遺伝子群は、機能的にも類似している

事後確率の数値行列」が得られる。

ことが期待される

12）。例えばある遺伝子の機能が未知でDEG 検出を兼ねた遺伝子のクラスタリング（MBCdeg）

ある場合、同じクラスタに分類された機能既知の遺伝子群と類似した機能を持つのではないかという推測に利用する

もちろん遺伝子クラスタリングの考え方を DEG 検出ことができる。
遺伝子発現解析の多くの実験デザインは 2 群間比較（G1そのものに適用するアイデアもごく少数ではあるが存在

群 vs. G2 群）や 3 群間比較（G1 群 vs. G2 群 vs. G3 群）でする 13, 14）。そのうちの 1 つが、我々が 2021 年に提唱しあり、その主目的は発現変動遺伝子（DEG）の検出であた、MBCluster.Seq をベースとして DEG 検出を行う戦略

る。特に 2 群間比較の場合は、DEG のとりうるパターン（MBCdeg と命名）である 13）。簡単にいえば、（DEG 検出

は、基本的に G1 群で高発現パターン（DEG_G1）か G2 群も目的に含まれるので）「全遺伝子のカウントデータ」とで高発現パターン（DEG_G2）のいずれかであるため、遺「想定する発現パターン数に相当するクラスタ数 K 」を与伝子のクラスタリングを行う概念自体がない。それゆええて MBCluster.Seq を実行するものである。

遺伝子クラスタリングは、3 群間比較や時系列データのよ

前述のとおり、MBCluster.Seq 自体は本来 DEG 検出を

うな、とりうるパターンが複雑になりうる実験デザインの想定しておらず、DEG 検出後のパターン分類目的で用いDEG 検出結果のパターン分類目的でしか利用されてこな

るものである。しかし全遺伝子のカウントデータを入力と

かった

13）

して MBCluster.Seq を実行すること自体は手続き的に可。

能である。この場合、直感的には「DEG ではない遺伝子

代表的な RNA-seq カウントデータの遺伝子クラスタリ5）である。このパッ群（non-DEG）から構成される大きなクラスタ」が形成さ

ケージは、入力として「DEG のカウントデータ」と「想定れ、「そのクラスタ中心の発現パターンはフラット（比較する発現パターン数に相当するクラスタ数 K 」を与えてクする群間でほぼ同じ値）」であることが期待される。我々

ラスタリングを実行する。出力は、「クラスタ中心の発現の MBCdeg は、この性質を利用して non-DEG クラスタ

ング用パッケージは MBCluster.Seq

パターン」と「遺伝子ごとの各クラスタへの属しやすさをを同定し、そのクラスタに属する事後確率が低い順に遺表す事後確率の数値行列」である。例えば 2 群間比較（G1

伝子をランキングするというのが基本戦略である。non100

Jpn. J. Lactic Acid Bact.

Vol. 33, No. 2

図 8．MBCdeg 実行結果
スクリプト実行結果として自動生成される fig8.png をベースとして、右上部分の数値情報を後から追加している。

MBCdeg の実行

DEG クラスタの事後確率は、統計的検定を行った結果として得られる p 値と似たようなものである。つまり事後

MBCdeg には、データ正規化の部分で MBCluster.Seq

確率の値が低い遺伝子ほど、non-DEG とみなすには無理の デ フ ォ ル ト の 方 法 を 実 装 し た MBCdeg1 と、 頑 健 ながある（つまり DEG である）と解釈すればよい。
MBCluster.Seq は、遺伝子を事後確率が最も高いクラRNA-seq デ ー タ 正 規 化 法 で あ る DEGES 18）を 実 装 し た

スタに割り当てることで発現パターン分類を行うのが本

MBCdeg2 の 2 種類が存在する。ここでは、単純化のため来の利用法である。non-DEG 以外のクラスタは、当然に「酸ストレス長期暴露群（pH4.5_24h）vs. 対照群（pH7_

フラットではない何らかの発現変動パターンをもつため、CCG）」という 2 群間比較を MBCdeg2 で実行する。クラ

MBCdeg を使うことで DEG 同定と発現パターン分類がスタ数については、I 群間比較の場合は、K ＞I ＋1 を推奨同時にできるというメリットもある。例えば、2 群間比較している。つまりこの場合は K ＞ 3 が推奨ということにな

データを K ＝3 で実行すると、直感的には得られた 3 つのるが、ここでは挙動確認の意味も込めて K ＝3, 4, 5 を独立

クラスタが「non-DEG、DEG_G1、DEG_G2」のいずれかに実行する（W10）。

のパターンを示すことが期待される。それゆえ、MBCdeg

図 8 は、クラスタごとのクラスタ中心の発現パターンをの実行結果の段階で、得られた DEG が DEG_G1 と DEG_

計 3 種類の K 値で実行した結果である。縦軸の値は log 比、G2 のどちらのパターンに属するかまで得られることに

横軸は左側が pH4.5_24h 群、右側が pH7_CCG 群である。

なる。

右上の数値は、上が pH7_CCG 群側の log 比の値、下がそ

MBCdeg は、DEG 検出専用の有名なパッケージであのクラスタに属する遺伝子数を表している。真ん中の K ＝る edgeR 15）や DESeq2 16）はもちろんのこと、我々が提供

4 の結果を例にとると、1 つめのクラスタ（c1）には 94 個

13）

。あくまで

の遺伝子が属しており、その代表パターンは pH7_CCG 群

も TCC が生成するシミュレーションデータに基づく結が pH4.5_24h 群に比べて e2×（-2.231）＝0.0115 倍高発現（86.66

論ではあるものの、データ中に占める DEG の割合や偏り

倍低発現）だと解釈する。2 つめのクラスタ（c2）には 684

といった様々なシナリオを用いた性能評価の結果として、個の遺伝子が属しており、その代表パターンは e2×（0.398）＝

TCC を上回る性能であったという点が重要である 17）。

2.217 倍高発現だと解釈する。3 つめのクラスタ（c3）にする TCC をも圧倒する性能を示しているは 1,270 個の遺伝子が属しており、その代表パターンは
e2×（0.068）＝1.146 倍高発現だと解釈する。4 つめのクラスタ（c4）には 901 個の遺伝子が属しており、その代表パター
101Vol. 33, No. 2

Japanese Journal of Lactic Acid Bacteria

ンは e2×（-0.459）＝0.399 倍高発現だと解釈する。この場合は、現）という結果と比較すると、満足のいく結果とは言い難

c3 が最も変動の少ない non-DEG クラスタと定義される。

い。しかしそもそも 2 群間比較の場合は K ＞ 3 を推奨して

K ＝4 の結果の場合、pH7_CCG 群で低発現パターンをおり、この意味において原著論文のガイドラインは妥当と示すクラスタは c1 と c4、高発現パターンのクラスタは c2いえる。

だと解釈する。K ＝5 の場合は、non-DEG クラスタは 984
おわりに個の遺伝子から構成される c2、pH7_CCG 群で低発現パターンを示すクラスタは 84 個の c1 と 763 個の c3、高発

本稿の前半は、RStudio の基本的な利用法、パッケージ

現パターンのクラスタは 263 個の c4 と 855 個の c5 だと解のインストールやロード、そして解析データの全体像を述釈する。
一部の読者は non-DEG だと判定されたクラスタを構べた。後半は、RNA-seq カウントデータ専用のパッケー

成する遺伝子数が、K ＝4 の c3 が 1,270 個、K ＝5 の c2 がジを用いたクラスタリングについて述べた。最初のサンプ

984 個と大きく異なっていると思われたかもしれない。主

ルクラスタリングでは、パッケージと関数の関係性、あな理由は、K ＝4 の c3 は pH7_CCG 群で若干の高発現を示りがちなミスとその対処法、そして結果の解釈のしかたす遺伝子を含んでいるからだと解釈すればよい。実際、Kについて解説した。次に代表的な遺伝子クラスタリング＝5 では大きめの高発現パターンを示す c4 が 263 個、小用パッケージである MBCluster.Seq の概要を述べるととさめの高発現パターンを示す c5 が 855 個となっており、

もに、このパッケージが本来想定していなかった DEG 検

K ＝4 で高発現パターンを示す c2 の 684 個よりも全体とし出にも拡張して応用可能であり、かつその精度が DEG 検て多くなっている。K ＝4 の c3 にカテゴライズされていた出専用のパッケージをも凌駕することを述べた。最後に、一部の若干の高発現を示していた遺伝子群が、K ＝5 とク

MBCluster.Seq に 基 づ く DEG 検 出 ア ル ゴ リ ズ ム で あ る

ラスタ数を増やしたことによって、c5 に含められるよう

MBCdeg を実行し、指定したクラスタ数の違いによる結になったのだと解釈すればよい。

果の違いや解釈の仕方について述べた。

ただしこれらの解釈は、あくまでも図 8 の結果を見比べ指定した K 値ごとの遺伝子の事後確率情報は、今回て導き出しただけの想像にすぎない。原著論文レベルのクの W10 のスクリプト実行結果ファイルとして生成されるオリティに仕上げるためには、「In fact, …」で書けるよう

MBCdeg_K ○ .xlsx（○の部分には 3, 4, 5 の数値が入る）

な、実際のところどうだったかを丁寧に調べる必要がある。

中に含まれている。このスクリプトについては、どこにど

この場合は、例えば K ＝4 の c3 に属していた遺伝子のどれのような出力結果が格納されているかや、それらを用いてだけが K ＝5 の c2 や c5 に含まれているのか？

どのように集合演算を行うのかなどを含め、次回以降丁寧などは確認しておきたいところである。この種のちょっとした解析に解説していく予定である。

を自在に行うためには、出力結果のどこにどのような情報謝が格納されているかを把握したり、それらの情報を用いて集合演算を行うスキルが必要となる。

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

な お、K ＝3 の 結 果 は、 最 も 変 動 の 少 な い c3 が non-

たものです。

DEG クラスタとみなされるものの、その代表パターンは
e2×（-0.180）＝0.698 倍 pH7_CCG 群で高発現（1.433 倍低発現）となっている。K ＝4 の c3 が e2×（0.068）＝1.146 倍高発現、そして K ＝5 の c2 が e

利益相反（COI）
牧野磨音、清水謙多郎、門田幸二：本論文発表の内容に2×（-0.034）＝0.934 倍高発現（1.070 倍低発

関連して開示すべき COI 状態はない。

		参 考 文 献1）	門田幸二，孫建強，湯敏，西岡輔，清水謙多郎（2014）次世
代シーケンサーデータの解析手法：第 1 回イントロダクション，日本乳酸菌学会誌 25: 87-94．
2）	門田幸二，大森良弘，寺田朋子，三浦文，寺田透，清水謙多郎（2021）次世代シーケンサーデータの解析手法：第 17 回
バイオインフォマティクス教育の今後，日本乳酸菌学会誌 32:129-37．
3）	Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M,
et al. (2004) Bioconductor: open software development forcomputational biology and bioinformatics. Genome Biol 5:
R80.
4）	Sun J, Nishiyama T, Shimizu K, Kadota K (2013) TCC: an R
package for comparing tag count data with robust normal-

ization strategies. BMC Bioinformatics 14: 2195）	Si Y, Liu P, Li P, Brutnell TP (2014) Model-based clusteringfor RNA-seq data. Bioinformatics 30: 197-205.
6）	孫建強，三浦文，清水謙多郎，門田幸二（2015）次世代シーケ ン サ ー デ ー タ の 解 析 手 法： 第 3 回 Linux 環 境 構 築 か ら
NGS データ取得まで，日本乳酸菌学会誌 26: 32-41．
7）	寺田朋子，坂本光央，清水謙多郎，門田幸二（2019）次世代
シーケンサーデータの解析手法：第 13 回 RNA-seq 解析（その 1），日本乳酸菌学会誌 30: 38-45．
8）	寺田朋子，清水謙多郎，門田幸二（2019）次世代シーケンサー
データの解析手法：第 14 回 RNA-seq 解析（その 2），日本乳酸菌学会誌 30: 153-61．
9）	寺田朋子，清水謙多郎，門田幸二（2020）次世代シーケンサー
データの解析手法：第 15 回 RNA-seq 解析（その 3），日本

102

Jpn. J. Lactic Acid Bact.

Vol. 33, No. 2

乳酸菌学会誌 31: 25-34．
10）	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.
11）	Kadota K, Tominaga D, Akiyama Y, Takahashi K (2003)
Detecting outlying samples in microarray data: A criticalassessment of the effect of outliers on sample classification.
Chem-Bio Informatics J 3: 30-45.
12）	Eisen MB, Spellman PT, Brown PO, Botstein D (1998) Cluster analysis and display of genome-wide expression patterns. Proc Natl Acad Sci USA 95: 14863-8.
13）	Osabe T, Shimizu K, Kadota K (2021) Differential expression
analysis using a model-based gene clustering algorithm forRNA-seq data. BMC Bioinformatics 22: 511.

14）	Vavoulis DV, Francescatto M, Heutink P, Gough J (2015)
DGEclust: differential expression analysis of clustered countdata. Genome Biol 16: 39.
15）	Robinson MD, McCarthy DJ, Smyth GK (2010) edgeR: a
Bioconductor package for differential expression analysis ofdigital gene expression data. Bioinformatics 26: 139-40.
16）	Love MI, Huber W, Anders S (2014) Moderated estimation of fold change and dispersion for RNA-seq data with
DESeq2. Genome Biol 15: 550.
17）	門田幸二 , 清水謙多郎（2021）次世代シーケンサーデータの解析手法：第 16 回なぜ次から次へと新規手法が開発されるのか ?, 日本乳酸菌学会誌 32: 123-8.
18）	Kadota K, Nishiyama T, Shimizu K (2012) A normalization
strategy for comparing tag count data. Algorithms Mol Biol
7: 5.

Methods for analyzing next-generation sequencing data
18. Clustering for gene expression data.
Manon Makino1, Kentaro Shimizu1, 2, 3, Koji Kadota1, 2, 3
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.

1		

Abstract

Gene expression data refers to a numerical matrix in which genes are arranged in each row andsamples in each column, and each element stores numerical information on the amount of expression.
In this paper, we use the numerical matrix which measured the acid stress response of Lactobacillus

rhamnosus GG, consisting of 2,949 genes × 9 samples. We first describe the basic usage of RStudio,
an integrated development environment for R. Next, we discuss the significance of clustering samples
having similar expression patterns, and the interpretation of the results. For gene clustering of RNAseq data, we outline MBCluster.Seq, a representative package for the purpose. Finally, we introduce
a modified version of the package (called MBCdeg) that can also be used for detecting differentiallyexpressed genes (DEGs). Supplementary materials are available online at https://www.iu.a.u-tokyo.
ac.jp/kadota/r_seq2.html.

103