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

総説

次世代シーケンサーデータの解析手法 第 20 回
RNA-seq カウントデータの性質と統計モデル
牧野 磨音 1、坂本 光央 2、清水 謙多郎 1, 3, 4、門田 幸二 1, 3, 4*
1

東京大学大学院 農学生命科学研究科
理化学研究所 バイオリソース研究センター
3
東京大学大学院 情報学環・学際情報学府
4
東京大学 微生物科学イノベーション連携研究機構
2RNA-seq 解析の目的の多くは、比較する状態または群間で発現の異なる遺伝子（DEG）の同定である。
ほとんどのプログラムは R のパッケージとして提供されており、その入力は、カウントデータとよばれる各行が遺伝子、各列がサンプルからなる数値行列である。本稿では、なぜ負の二項分布とよばれる統計モデルが DEG 検出目的でよく用いられるのかについて、数式を交えて解説する。また、このカウント
データの性質を説明する手段としてよく用いられる平均－分散プロットについて、データの前処理から
ggplot2 による描画まで述べる。
Key words：negative binomial distribution, RNA-seq, R Markdown, ggplot2

はじめに（W1）（pH4.5_24h）、そして 7 ～ 9 列目が対照群（pH7_CCG）である。各群につき 3 反復ずつ取得したこの実験デザインに

1）

前回 に引き続き、今回のウェブ資料（以下、W）も R

おいて、同一群内の反復データ間のばらつきがどのようにMarkdown（JSLAB20.Rmd）で作成している。レンダリン

なっているかを、要約統計量・数式・散布図などを用いてグによって作成した R スクリプトと実行結果から構成さ解説するのが本稿の主な内容である。

れる html ファイル（JSLAB20.html）も含めて、ウェブサ

DEG 同定は多くの場合、「比較する群間に差がない」とイト（R で）塩基配列解析のサブ（URL: http://www.iu.a.u-

いう帰無仮説を（形式的に）立て、統計的手法を用いて行

tokyo.ac.jp/kadota/r_seq2.html）で提供している。今回は特われる。この枠組みにおいて、我々は非常に低い p 値が得に、JSLAB20.html を併用してほしい。なお、R パッケー

られることを期待するが、これは「同一群内のばらつきのジのインストールやロード、RStudio の基本的な利用法に範囲内」という許容範囲を超えた結果が得られるのを望む

2）

ついては、第 18 回 でも解説している。

ことと同義である。本稿のトピックである平均 - 分散プ

入力として用いる RNA-seq カウントデータは、第 18

ロットは、同一群内のばらつき度合いを把握する描画手段

回から利用している Lactobacillus rhamnosus GG（LGG）のとしてよく用いられる散布図であり、これが帰無仮説に従酸ストレス応答を調べた 2,949 遺伝子×9 サンプルの数

う分布（つまり帰無分布）である。この分布は同一群内の3）値行列である 。第 18 回の図 5 でも示されているよう

反復の種類に大きく左右されるため、本稿は反復の種類のに、（行名の列を除く）最初の 1 ～ 3 列目が酸ストレス短説明からスタートする。

期暴露群（pH4.5_1h）、4 ～ 6 列目が酸ストレス長期暴露群
生物学的な反復と技術的な反復＊

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

反復データには、生物学的な反復（biological replicates）と技術的な反復（technical replicates）の 2 種類が存在する。
たとえばヒトの筋肉組織を塊として調べる bulk RNA-seq
21Vol.34,No.1JapaneseJournalofLacticAcidBacteriaの場合、「A さん、B さん、C さんという異なる個体からる一種の指標である。それゆえ、一般にこのばらつきの度

それぞれデータを取得して 3 反復としたものが生物学的な合いは、小さければ小さいほどよい。その一方で、生物学反復」に、そして「同一個体（たとえば A さん）の筋肉組

的な反復データ間のばらつき（biological variation）は、小織を 3 分割して 3 反復としたものが技術的な反復」に相当

さかろうが大きかろうがそれをそのまま受け止めるのが基する（図 1a）。また、筋肉組織を構成する細胞ごとの情報本である。重要なのは、どのような別個体のデータを取得を調べる scRNA-seq の場合、「A さん、B さん、C さんとしているのかである。

いう異なる個体からそれぞれ 1 つの細胞の発現データを取

具体的には、たとえば A さん、B さん、C さんがきょう

得して 3 反復としたものが生物学的な反復」に、そして「同だいの場合と赤の他人の場合で考えてみるとよいだろう。

一個体（たとえば A さん）の筋肉組織を構成する 3 つの筋

生物学的なばらつきは、後者よりも前者のほうがより小さ

肉細胞の発現データを取得して 3 反復としたものが技術的いのは明らかである。また、当該の 3 人が同じ赤の他人同な反復」に対応する。もちろん実際の scRNA-seq では、異士でも、たとえば同じ県内に住んでいる場合（例：全員沖なる個体から 1 つだけ、あるいは同一個体から数個の細胞縄県）と地理的に遠い別の都道府県（例：A さんが北海道、だけのデータを取得することはなく、ある個体の細胞塊にB さんが石川県、C さんが佐賀県）の場合でも、後者より含まれる数百～数千細胞のデータを取得するのが基本形でも前者のほうがより小さなばらつきとなるであろう。これある。

を突き詰めていくと、たとえば性別や人種の異なる別個体

ここまでで直感的に理解できることは、データのばらつ

で反復をとれば、ヒトという生物種により普遍的な現象をきは技術的な反復よりも生物学的な反復のほうが大きいと捉えうることを意味する。

いうことであろう。これは、発現データ取得手段がマイク

バクテリアの場合は、たとえば培地上の異なるコロニー

ロアレイの頃から議論されている実験デザインに関する事をそれぞれ異なる細胞ペレット（cell pellets）として回収し柄である 4）。技術的な反復データ間のばらつき（technical

たのち、ペレットごとに処理を行ってデータを取得したもvariation）は、実験機器や実験者の手技の安定性を計測すのが生物学的な反復に相当する。また、培地上の同一コロ
ニーを異なる細胞ペレットとなるように分割して回収し、ペレットごとにデータを取得したものが技術的な反復に相当する（図 1b）
。今回用いた乳酸菌データが実際にどのよ3）
うな手順で処理がなされたのかを原著論文 のみから正確に読み解くのは困難であるが、後述するデータの性質より、生物学的な反復に相当するものだと我々は判断している。
解析データの記号表現（W2）冒頭で言及した平均 - 分散プロットは、遺伝子ごとに同一群内の反復データの平均と分散の値を算出し、横軸に平均、縦軸に分散の値をプロットした散布図のことである。
RNA-seq カウントデータだけでなく、マイクロアレイデー
タの性質を説明する手段としても利用されている 5-7）。ただし、基礎情報である平均と分散の算出に先立ち、サンプル（反復または列）間の総カウント数の違いによる影響を排除する必要がある 7）。表 1 は、サンプルごとの総カウン
ト数を、他の基本情報と合わせて示したものである。第13 回 8）ではサンプル間のリード数が 10 倍以上異なることに言及したが、総カウント数についても同様の傾向になっていることがわかる。つまり、総カウント数の最大と最小の間には、2,641,162 / 233,195＝11.33 倍もの違いがある（W2.2）。
入力であるカウントデータ行列のうち、総カウント数
が最大のサンプルは 5 列目、最小のサンプルは 1 列目である。学術論文では、このような「総カウント数が最大ま

図 1．生物学的な反復と技術的な反復（a）ヒトの場合、（b）バクテリアの場合。TogoPictureGallery のイラスト（CC-BY-4.0）をそのままあるいは一部改変して利用。

たは最小のサンプル」を関数や記号でうまく表現せねばな
らない局面がある。この目的のためによく利用されるの
22

Jpn.J.LacticAcidBact.

Vol.34,No.1

表 1．解析データ（乳酸菌 RNA-seq カウントデータ）の基本情報
一番右の列がサンプルごとの総カウント数（W2.1 の実行結果と同じ）、それ以外の情報は第 13 回の表 1 と同じである。
サンプル名

SRRID

pH4.5_1h_1、そして s9＝pH7_CCG_3 だと解釈する。なお、
i や j は添え字またはインデックスとよばれるものである。
前述の「5 列目にあるカウント数が最大のサンプル」は s5 でリード数（片側のみ） 総カウント数

pH4.5_1h_1

SRR6322562

 301,126

 233,195

pH4.5_1h_2

SRR6322563

1,470,602

 840,602

pH4.5_1h_3

SRR6322564

1,760,461

1,072,851

pH4.5_24h_1

SRR6322565

1,375,368

 938,788

pH4.5_24h_2

SRR6322566

3,869,088

2,641,162

pH4.5_24h_3

SRR6322567

1,795,874

1,176,816

pH7_CCG_1

SRR6322568

3,095,834

2,381,981

pH7_CCG_2

SRR6322569

2,570,876

1,758,968

pH7_CCG_3

SRR6322570

 846,623

 569,947

あるが、この j＝5 というインデックス情報で当該サンプルを指し示すことができる。
前述の argmax や argmin は、総カウント数情報が格納されたサンプル数分の要素からなる数値ベクトル（233195,
840602, …, 569947）を入力として実行した場合、総カウン
ト数の最大値や最小値のインデックスを返す関数だという理解でもよい。以下に示すように、入力ベクトル中の要素
そのものを返す max や min の実行結果と見比べると、両者の違いがよくわかるだろう。

は argmax や argmin であるが、現実にはこの記号を見た時
点で思考停止してしまうヒトが多いのではないだろうか。

・max（233195, 840602, …, 569947）＝2641162

我々は、この原因はほとんど全て数式や記号のみで話が展・min（233195, 840602, …, 569947）＝233195

開されていくためであると認識している。このような経験・argmax（233195, 840602, …, 569947）＝5のある多くのヒトは、「たぶんこういうことだろうが…」

・argmin（233195, 840602, …, 569947）＝1という理解を確固たるものにするための分岐点（正しい理解のルートと曖昧な理解のルート）が読み進めていくうち

本 稿 の 表 記 法 と 完 全 一 致 し て い る わ け で は な い が、にどんどん増えていき、最終的に論点整理しきれなくなり

argmin については第 18 ～ 19 回で用いた MBCdeg 法の原撃沈という経過を辿ってきたであろう。しかし今一度、本著論文 9）で実際に利用しているので解読を試みてもよいだ

稿を足掛かりにして数式の克服（記号を用いた説明に慣れろう。なお、argmax や argmin に対応する R の関数としてること）に再チャレンジしてみてほしい。

は、which.max や which.min が用意されている（W2.4）。

たとえば、まず行列を A、各行を個々の遺伝子（ i＝1, 2, …,
総カウント数の計算（W3）G）
、各列を個々のサンプル（ j＝1, 2, …, S）、そして i 行 j 列目の要素 ai,j をカウント数と定義する。行列とその要素（元ともよばれる）などをまとめて、

さきほど定義した記号を用いると、j 列目のサンプル sjのように表してもよい。ここで、G は遺伝子数、そして S はの総カウント数 Tj は以下のように表すことができる。

サンプル数である。今我々が取り扱っているカウントデータの場合は、2949 遺伝子×9 サンプルの数値行列なので、G＝（1）

2949、S＝9 である。各要素の具体的な数値としては、たと
えば a3,1＝6579、a2947,2＝212、a2949,9＝6157 のように解釈する（図 2）。

具体例として、1 列目のサンプル s1（つまり pH4.5_1h_1）（g1, g2, …, gG）や遺伝子名やサンプル名についても、gi＝の総カウント数 T1＝233195 を得る手順を以下に示す。5 列

sj＝（s1, s2, …, sS）のように記号で表現することができる。こ目や 9 列目など、他の計算例を見比べるとより理解が深まの 場 合 は、g1＝EBG00001128470、g2949＝LGG_02944、s1＝

るであろう（W3.1）。

図 2．解析データの概要と記号表現

23Vol. 34, No. 1

Japanese Journal of Lactic Acid Bacteria

要 素 と す る 行 列 は、と表すことができる。但し、CPM の具体的な計算例を
含 め て 第 15 回 10）で も 述 べ て い る が、 正 規 化 後 の デ ー
タは群間比較時の入力として用いてはいけないので注意してほしい。ここではあくまでも同一群内の平均と分散を算出するための基礎情報として得ただけである。
それぞれのサンプルの総カウント数に対して Tj という平均と分散を算出（W5）

記号を割り当てることで、前述の「総カウント数情報が格納されたサンプル数分（つまり S＝9）の要素からなる数
値ベクトル（233195, 840602, …, 569947）」は、（T1, T2, …,

平均と分散は群ごとに算出するため、ここであらため

TS）または T（
j j＝1, 2, …, S）のように表すことができる（TSてどのサンプル（どの列）がどの群に属するかを示すクラは T9 でもよい）。記号を用いた表現手段に慣れてくると、スラベル L＝（l1, l2, …, lS）を作成する。一般には、たとえたとえば上述の数値ベクトルで表現した argmax（233195,

ば群数を K、群 k（＝1, …, K）の反復数を nk として、様々

840602, …, 569947）のほうが、むしろ違和感があるようなな実験デザインを表現可能としておく。今入力として用気がしてくるのではないだろうか。具体的な使用例としている行列

は、最初の 1 ～ 3 列目が酸ストレス短期は、u＝argmax（T1, T2, …, TS）とおくことで、総カウント数

暴露群（pH4.5_1h）、4 ～ 6 列目が酸ストレス長期暴露群が最大となるサンプルを su と表現することができる。

（pH4.5_24h）、そして 7 ～ 9 列目が対照群（pH7_CCG）である。このため、K＝3、k＝1 が pH4.5_1h 群でその反復数

CPM の計算（W4）が n1＝3、k＝2 が pH4.5_24h 群でその反復数が n2＝3、そして k＝3 が pH7_CCG 群でその反復数が n3＝3 だというこ

今は、平均 - 分散プロットの基礎情報である平均と分散とを意味する。この場合、クラスラベル L の実体は、計 9を得るための前処理段階にある。ここでは、サンプル（反つの要素からなる整数ベクトル（1, 1, 1, 2, 2, 2, 3, 3, 3）のよ

復または列）間の総カウント数の違いによる影響を排除す

うに解釈する（W5.1）。つまり、クラスラベル L の j 番目（ j

る手段の 1 つとして、一般によく用いられる CPM 正規化＝1, 2, …, S）の要素はそれぞれ、l1＝l2＝l3＝1、l4＝l5＝l6＝を実行する。CPM は Counts per million の略であり、正規

2、そして l7＝l8＝l9＝3 である。
クラスラベル L を用いることで、任意の群 k に属する

化後の総カウント数が 100 万に揃うような係数を掛ける正目の CPM 正規化後のカウント値のように表現すること

L 中の要素（元）を規化作業および正規化後の値を指す言葉である。i 行 j 列は、対応する正規ができる。たとえば、群 2（つまり k＝2）のクラスラベル

化前のカウント値 ai,j と j 列目のサンプル sj の総カウント

情報はと表現することができる。縦棒（|）

数 Tj を用いて以下のように表すことができる。

は条件を表す記号であり、縦棒より右側に書かれているのが具体的な条件である。L 自体は 9 つの要素からなるベ（2）という条件を満たす要素はそのクトルであるが、

うちの 3 つ、つまり { l4, l5, l6 } である。縦棒の左側にある
具体例として、2 列目のサンプル s2（つまり pH4.5_1h_2）

は、「の総カウント数（T2＝840602）を用いて、3 行 2 列目（ i＝3と j＝2）の CPM 正規化後のカウント値

は L に属する」という意味である。難解だと感じた読者は、たとえば L＝（l1, l2,…, lS）は

を得る手順

や

とも

を以下に示す。他のいくつかの計算例を見比べると、より表現できることを思い出せばよいだろう。この場合は、縦理解が深まるであろう（W4.1）。

棒の右側でオリジナルの 9 つの要素全てを満たす条件を記載しているのだと解釈すればよい。
実 用 上 は、のような特定の群 k に対応するクラスラベル L 中の要素だけではなく、そのイ ン デ ッ ク ス 情 報 を 利 用 し た い 局 面 も あ る。1 つ の 表式（2）の右辺の ai,j に掛かっているは、正現方法としては、クラスラベル L 自体に添え字をつけ規化係数とよばれるものに相当する。添え字（インデックて、 任意の群 k に属するクラスラベルのサブセットを、そのインデックスをjス）が（＝1,
2, …, S）であることからも明らかではあるが、この値はサンプル数分（つまり S＝9）の要素からなる数のように定義することである。たとえば群 2（つまり k＝

値ベクトルである（W4.2）。正規化係数の値が 1 より大き2）の 場 合 は、のインデックスはい（小さい）と、正規化後のカウント値は正規化前に比べて大きく（小さく）なる。正規化後のカウント値をことで、たとえば行列
24、そとなる。こうするの中から群 k に対応する列

Jpn. J. Lactic Acid Bact.

Vol. 34, No. 1のみ取り扱いたい場合に、その列のインデックス情報ををと指定できる。

と定義できる。たとえば

行列 V の 6 行 2 列の要素である v6,2 には、6 行目の遺伝子

ク ラ ス ラ ベ ル を 用 い る 一 番 の メ リ ッ ト は、群 k におの群 2 の分散が格納されている。分散の具体的な計算自ける i 行目の遺伝子の平均カウント値を mi,k、そして mi,k

体は、R の var 関数を用いて実行すればよい。一般的な不を要素とする行列をの偏分散の式を本稿の記号に置き換えれば、以下のようにように記号のみで過不足なく説明できる点である。たvi,k を定式化できる。さきほどの平均の具体的な計算同様、

かi に 6、j に 4, 5, 6、k に 2 を代入して確かめてみるとよいと え ば、 計 3 群（K＝3）の ラ ベ ル 情 報 を も つらの場合は、2949 行×3 列の平均カウントの行列（W5.3）。

が得られるのだと解釈する。
mi,k は以下のように定式化できる。
（4）（3）
平均と分散の関係性（W6）
たとえば行列 M の 6 行 2 列の要素である m6,2 には、6 行目の遺伝子の群 2 の平均が格納されている。つまり、g6＝

図 3 は、①平均 M および②分散 V の数値行列の概要で

LGG_00001 のカウントのうち、pH4.5_24h 群のみのカウン

ある（W6.1）。数値行列を入力として R の summary 関数をトの平均が格納されている。これを式（3）に近い形で表現

実行すると、列ごとの各種統計量（最小値、中央値、平均

すると、行列中の 6 行目の要素のうち、ク

値、最大値など）が返される。この場合は、V1 という列

ラスラベルが 2（）という条件を満たす j（l4＝l5＝l6＝が群 1（pH4.5_1h 群）に、V2 が群 2（pH4.5_24h 群）に、そ

2 なので j＝4, 5, 6）で構成される 3 つの値（つまりして V3 が群 3（pH7_CCG 群）の結果に相当する。①群ごと）の平均として 412.648 が得られている。

との平均カウント数の結果において、Mean の値が 339.1

以下に示すように、式（3）に具体的な数値を入れて確かめに揃っているのは、総カウント数を 100 万に揃える CPMとるとよいだろう（W5.2）。

正規化後のデータ（）に基づいているためである）。バクテリアのような数千程度（の遺伝子数であれば、ここで示されているような数百程度の Mean の値となる。その一方で数万程度の遺伝子数であれば、数十程度の Mean の値となる（たとえば遺伝子数 2万の場合は 1000000/20000＝50）。一見何気ない部分ではあるが、意図通りの前処理ができているかを真っ先に確認す
分散についても同様である。群 k における i 行目の遺伝

るチェックポイントである。

子のカウント値の分散を vi,k、そして vi,k を要素とする行列

次に最大値（③と④）を眺める。③平均の最大値は 5 桁、

図 3．summary 関数実行結果のスクリーンショット

25

Vol. 34, No. 1

Japanese Journal of Lactic Acid Bacteria④分散の最大値はそれよりもさらに 3 ～ 4 桁大きいことがれる。

わかる。一般にダイナミックレンジが広いデータは対数軸のような式を用

重要な点は、で全体像を把握することが多く、この場合も対数変換後のいれば、RNA-seq カウントデータの平均と分散の関係性

データで平均 - 分散プロットを描画する必要があるとこのをうまく捉えられるということである。今我々は、計 3 種

段階で判断できる。また、群ごとの平均と分散の関係性は、類の同一群内の反復データから算出した、2949 遺伝子×3

第一四分位数（1st Qu.）、中央値（Median）、そして第三四群＝8847 個の平均の分位数（3rd Qu.）を見比べることで大まかに把握可能であ類の情報をもつ。この式は、と dispersion のとという 2 種の情報があればる。たとえば群 1 はいずれも分散のほうが平均よりも大きがわかるということであり、RNA-seq カウントデーな値になっていることがわかる。群 3 も群 1 ほどではない

タが負の二項分布（negative binomial distribution; NB 分

が同様の傾向になっている。しかし群 2 については、平均布）に従うことと同義である。これは、統計学の分野でと分散が似たような分布になっている。具体的には、1st

は「データを特定の確率分布に従う確率変数とみなして

Qu. と Median で平均のほうが分散よりも大きく（44.0 > 17いる」ためである。この場合は、「RNA-seq カウントデーと 103.7 > 69）、3rd Qu. で逆転している（236.8 < 284）。しタを NB 分布に従う確率変数とみなしている」ということ

かし、いずれも他の群の違いに比べれば誤差範囲といえるになる。さらに具体的には、「RNA-seq カウントデータ

レベルである。

（）は、平均（）、dispersion（）という 2 つ

この summary 関数実行結果は、平均 - 分散プロットにのパラメータから構成される NB 分布に従う確率変数」だ

よって得られる散布図のイメージそのものである。たとえとみなせるということである。統計学の分野では、これをのように表記する 11）。

ば最も平均と分散が似ている群 2 の散布図では、2949 遺を i＝6 および k＝2 の具伝子分の点が y＝x の直線近辺に散らばっている（つまりvi,2＝mi,2）ことが容易に想像できる。その一方で、群 1 の体例で説明すると、k＝2（

散布図では、ほとんどの点が y＝x の直線の左上側に偏っ＝l5＝l6＝2 なので j＝4, 5, 6）で構成される 3 つの値（つまりて存在するであろう。これは vi,1 > mi,1 と表現してもよいが、
（

たとえば

j l4
）という条件を満たす （とは定数）のように

）は、

と
および

という 2 つの等号を用いて表現することができる。これは、一般的なパラメータから構成される NB 分布から得られたデータ（NBと似たような分布に従う確率変数）だとみなせるということである。確率多項式
ものである。2 次の多項式、て、定数部分を、そしてにおい変数という言葉がそれでも難解だと感じる読者は、これをとおけば、「確率的に変化する数値」に置き換えると同時に、k＝2 の群と見栄えが同じ（つまり酸ストレス長期暴露群（pH4.5_24h））における i＝6 行となる。

目の遺伝子に対して、4 つ目や 5 つ目の反復データを得た結

さきほど群 1 に限定して定式化したは、任意の群 k に対して果を想像してみるとよい。ほとんどのヒトは、既に得られてとといる 3 つの値（つまり表すことができる。一部の読者は、見慣れない （ふぁとい）の出現によって、理解が追い付いていないかもしれ

わず、これらの既出の値の範囲内か、範囲外だとしても既出ない。しかし今我々は、任意の群 k における i 行目の遺の値からそれほどかけ離れない 380 や 438 のような値が得ら

と分散

伝子の平均

）のいずれかと全く同じ値が出るとは思の具体的な数値情報が手元にれると予想するであろう。このように、手元にない仮想的なをデータを思い浮かべれば、確率変数とみなすほうがむしろ

あ る（W6.2）。 し た が っ て、
簡単な式変形によって

としてお

と

適切だと思えるのではないだろうか。

の値を代入すれば、任意の群 k

R では NB 分布に従う乱数を発生させる rnbinom 関数がにおける i 行目の遺伝子のを得ることができる。たと提供されているので、読者自身が実際に乱数を発生させ、とを 代 入 す れ ば、その分散が期待値と一致するか確認できる（W6.3）。たとき、具体的な
え ば、が得らを得る計算式の分子

れる。
解けるが、

は分散

と平均

え ば、

と

ば、
得られる。なお、この
れる

11）

。全ての

と同じサイズの行列を用いて 8からも読み個乱数を発生させると、（400, 392, 446, 430, 377, 423, 406,の差が小さいほど412）という結果が得られる。この分散は 484.786 であり、ゼロに近い値となり、マイナスの値も計算上は得られる。
た と え ば、

および

期待値（

）とかけ離れた結果だと思われるかを代入すれもしれない。しかし、たとえば 172000 個の乱数から計算がした分散は 716.481 となるように、発生させる乱数（試行は dispersion パラメータとよば回数）を増やせば期待値に収束していく。納得できるまでを計算することで、平均 M や分散 V

読者自身で検証してみるとよいだろう。

が得ら26

Jpn. J. Lactic Acid Bact.

Vol. 34, No. 1

平均 - 分散プロット（W7 ～ W10）値になっていることがわかる（W10.1）。第一四分位数（1stQu.）や第三四分位数（3rd Qu.）を眺めても、全体的にゼロ付

図 4 は、横軸に平均 M、縦軸に分散 V の値をプロット

近に位置していると判断できる。したがって、この群に限っした散布図である。群ごとにプロットした（a）～（c）につていえば、NB 分布というよりはむしろポアソン分布に従う

いては遺伝子数（＝2949）だけの点が、そして全データをと解釈するほうが自然といえよう。ポアソン分布は、NB 分）において、dispersion

まとめた（d）については 2949×3＝8847 個の点が存在す

布の数 式（る。左下から右上に伸びる y＝x の直線は、平均と分散がパラメータが全てゼロ（）の場合に相当する分布だ

完全一致する点に相当する。図の左上（右下）にある数値からである（つまり）。技術的な反復データはポ

は、y＝x の直線よりも左上（右下）側にある点の数、つまアソン分布に従うことも知られており 12）、このデータもそれり M < V（M > V）の遺伝子数である。この M < V（M > V）に匹敵するレベルで群内類似度が高いのだと解釈すればよは、より、dispersion パラメータい（W10.2）。

（ ）が正（負）の値になることと同義である。図 3 に基づ
統計モデルの意味

く議論と本質的に同じであるが、（a）と（c）では M < V となる遺伝子数が M > V に比べて 3 倍以上多いことがわかる。（b）では逆の傾向になっているものの、（d）全体とし一部の読者は、「我々が興味あるのは G1 群 vs. G2 群のを含む NB 分布で平均と分散の関ような群間比較であり、群内のばらつきではない」と思わてみると、明らかに係を定式化するほうがよいことがわかる。

れているかもしれない。そのような読者は、統計的検定の

図 4b の群 2（pH4.5_24h 群）については、遺伝子ごとのdispersion パラメータ（

意味や手続きを思い返してほしい。この場合の帰無仮説）の 平 均に ついても、（null hypothesis）は「比較する群間に差がない」であるが、とマイナスのどの程度であれば差がないとするかを具体的に定めないと

図 4．平均 - 分散プロット

27

Vol. 34, No. 1

Japanese Journal of Lactic Acid Bacteria

おわりにいけない。それが帰無仮説に従う分布、つまり帰無分布（null distribution）である。図 4 で我々が見ているのは帰無分布そのものである。

本稿では、RNA-seq カウントデータの性質を解説する

この場合の統計的検定は、たとえば比較する群を同一群という大枠の中で、バイオインフォマティクス分野の論由来だとみなした場合に、検定したい遺伝子の平均と分散

文を読み解く上で必要なリテラシーの向上（数式や記号の値が帰無分布からどれだけ離れているかを p 値として返表現に対する苦手意識の払しょく）を意識した構成にしす作業だと理解すればよい。もちろん p 値を得るためには、た。本文中ではほとんど触れなかったが、ggplot2 を用い「帰無分布をうまく数式で表現したもの」をあらかじめ構

た描画の基礎から応用については W08 ～ W09 で詳細に解築しておく必要がある。これがいわゆるモデル構築（model

説している。図 4 は、W9.8.3 で得た画像をベースとして、

construction）という作業に相当し、構築されたものは統計

PowerPoint で微調整しながら作成したものである。作図はモデル（statistical model）とよばれる。統計モデルの実体は凝りだすとキリがないうえ（作図沼に落ちる）、そもそも

ただの数式であり、この場合は「任意の平均と分散の値を

ggplot2 内で最初から最後まで描画せねばならないもので入力として与えると、それが帰無分布上のどのあたりに位もない。今回の図 4 作成手順のように、手作業の余地を残置するかを p 値で返してくれる数式」である。

すくらいの心構えでもよいだろう。

現実には、図 4 で見えている帰無分布を表現できる数

我々は今でも「反復なしデータで統計的手法を実行する

式（つまり統計モデル）は複数存在しうる。また、この帰とエラーになってしまうが、p 値を得るにはどうすればよ

無分布は全て 3 反復のみの限られたデータから得られたもいか？」という質問をたまに受ける。本稿の内容を正しくのである点も忘れてはならない。つまり、今手元にはない理解できた読者は、そもそも反復なしだと帰無分布を得よが 4 つめや 5 つめの反復データが将来得られると、帰無分

うがない（つまり統計モデルの構築ができない）ことがわ

布の形状は変わりうる。このような事柄も鑑みると、単にかるであろう。もちろん反復なしデータでもエラーを出さ

手持ちデータへのあてはまりのよさだけを唯一無二の評価

ずに結果を返すプログラムは存在するが、それらの結果は基準とすべきではないこともわかるであろう。このあたり

決して過信しないよう注意してほしい。

が統計学分野でよく行われている、複数の候補の中からど

次回は、今回得た乳酸菌 RNA-seq カウントデータの平の統計モデルを選択したほうがよいかという、モデル選択

均（M）と dispersion パラメータ（ ）を用いて、シミュレー（model selection）とよばれる議論である。

ションデータの生成を行うトピックを中心に解説する予定

最後のモデル云々で途端に難解だと感じた読者は、たとである。バイオインフォマティクス系の論文でよく行われえば居住地における過去 10 年の気温のデータから統計モているシミュレーションデータ解析の流れを知ることで、デルを構築し、3 日後・100 日後・1098 日後の気温をどう

表層的なノウハウだけでなく、バイオインフォマティシャ

やって予測するかという問題を考えてみるとよいだろう。

ンの思考回路にまで踏み込んでもらいたい。

多項式と正弦波（sin カーブ）の 2 択の場合、どちらの数式謝を選択するだろうか？著者らは当然後者を採用する。もち辞

ろん正弦波はあくまでも基本形であり、これをベースとし本内容の一部は、JSPS 科研費 21K12120 の助成を受けた改良版は複数存在しうる。この議論と同様、RNA-seq

たものです。図 1 のイラストは、TogoPictureGallery（CC-

カウントデータの統計モデルとしては NB 分布に従うモデ

BY-4.0）の図を（そのままあるいは一部改変して）利用さル（つまり NB モデル）が基本形だというのがこの分野のせていただきました。

コンセンサスだということである。
利益相反（COI）
牧野磨音、坂本光央、清水謙多郎、門田幸二：本論文発表の内容に関連して開示すべき COI 状態はない。

28

Jpn. J. Lactic Acid Bact.

Vol. 34, No. 1

		
参 考 文 献1） 牧 野磨音，清水謙多郎，門田 幸二（2022）次世代シーケン
サーデータの解析手法：第 19 回 R Markdown，日本乳酸菌学会誌 33: 195-205．
2） 牧野磨音，清水謙多郎，門田 幸二 (2022) 次世代シーケンサー
データの解析手法：第 18 回遺伝子発現データのクラスタリ
ング，日本乳酸菌学会誌 33: 87-94．
3） 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.
4） Zakharkin SO, Kim K, Mehta T, Chen L, Barnes S, Scheirer
KE, et al. (2005) Sources of variation in Affymetrix microarray experiments. BMC Bioinformatics 6: 214.
5） Archer KJ, Dumur CI, Ramakrishnan V (2004) Graphicaltechnique for identifying a monotonic variance stabilizing
transformation for absolute gene intensity signals. BMC
Bioinformatics 5: 60.
6） Anders S, Huber W (2010) Differential expression analysis

for sequence count data. Genome Biol 11: R106.
7） 門田幸二 著（金明哲 編）（2014）シリーズ Useful R ⑦ トラン
スクリプトーム解析，共立出版
8） 寺田朋子，坂本光央，清水謙多郎，門田 幸二（2019）次世代
シーケンサーデータの解析手法：第 13 回 RNA-seq 解析（その 1），日本乳酸菌学会誌 30: 38-45．
9） Osabe T, Shimizu K, Kadota K (2021) Differential expression
analysis using a model-based gene clustering algorithm forRNA-seq data. BMC Bioinformatics 22: 511.
10）寺田朋子，清水謙多郎，門田幸二（2020）次世代シーケンサー
データの解析手法：第 15 回 RNA-seq 解析（その 3），日本乳酸菌学会誌 31: 25-34．
11）Robinson MD, Smyth GK (2008) Small-sample estimationof negative binomial dispersion, with applications to SAGE
data. Biostatistics 9: 321-332.
12）Marioni JC, Mason CE, Mane SM, Stehpens M, Gilad Y (2008)
RNA-seq: an assessment of technical reproducibility andcomparison with gene expression arrays. Genome Res. 18:
1509-1517.

Methods for analyzing next-generation sequencing data
20. Properties and statistical models of RNA-seq count
data.
Manon Makino1, Mitsuo Sakamoto2, Kentaro Shimizu1, 3, 4,
Koji Kadota1, 3, 4
1

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

		Abstract

   RNA-seq is a tool for measuring gene expression and is commonly used for identifying differentiallyexpressed genes (DEGs) under different conditions or groups. Most of the programs for DEGs has been
provided as a free software environment called R. Users typically start the analysis with a numericalmatrix called “count data”, where each row a gene, each column a sample (a group’s replicate), and eachcell the number of counts. We describe the characteristics of the count data with a scatter plot (so-called“mean-variance plot”) displaying the relationship between the mean (x-axis) and variance (y-axis) within
replicates. We explain why a statistical model called the negative binomial distribution (NB model) has
been used for identifying DEGs.

29