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

解説

次世代シーケンサーデータの解析手法 第 15 回
RNA-seq 解析（その 3）
寺田 朋子 1、清水 謙多郎 1, 2、門田 幸二 1, 2*
12

東京大学 大学院農学生命科学研究科
東京大学 微生物科学イノベーション連携研究機構

Galaxy は、ウェブブラウザ上でマウスを操作して行う GUI ベースのデータ解析環境である。今回も前回に引き続いて、Galaxy 上で行う RNA-seq データ（GSE107337）の発現定量に関する解説を行う。まず、アノテーション情報を含む GFF ファイルの前処理（フィルタリング）を行い、遺伝子領域に対応するゲ
ノム中の塩基配列情報を抽出する。次に、得られた塩基配列群をリファレンス配列として Kallisto quant
プログラムを実行し、遺伝子ごとのカウント値や発現量に相当する TPM 値を得る。カウント値と周辺情報から CPM、CPK、FPKM、そして TPM といった様々な補正値を導き出す考え方について述べる。最後に、今回得られた結果を、GSE107337 の原著論文および第 14 回で得られたものと比較・検証する。ウェ
ブサイト（R で）塩基配列解析のサブ（URL: http://www.iu.a.u-tokyo.ac.jp/~kadota/r_seq2.html）中に本連載をまとめた項目（URL: http://www.iu.a.u-tokyo.ac.jp/~kadota/r_seq2.html#about_book_JSLAB）が存在する。ウェブ資料（以下、W）や関連ウェブサイトなどを効率的に活用してほしい。
Key words：NGS, Galaxy, RNA-seq, Lactobacillus rhamnosus GG

はじめに

今 回 の ス タ ー ト 地 点 は、 前 回（ 第 14 回 ）の W9-7 のGalaxy main 画面（ヒストリー名：GSE107337_3samples）

第 14 回 1）に 引 き 続 き、Lactobacillus rhamnosus GG 2）

である［W1］。しかし、必ずしも前回の内容を一通り行っの RNA-seq データ 3）を Galaxy 4）上で取り扱う。前回はておく必要はない。第 12 回 6）でも述べたように、ヒストマッピング時のリファレンスとしてゲノム配列を用いたリーの共有を我々にリクエストすればよいからである。こが、トランスクリプトーム配列にマップして高速にカウン

こでは新規ヒストリー（ヒストリー名：trans_map）を作

トデータを得るやり方も存在する。これは主にヒトやマウ成し［W2］、第 12 回の W13-3 でも解説したヒストリー間スなどの高等生物を対象とした戦略ではあるが、ゲノムサのデータコピーによって事前準備を行う。具体的には、ゲイズの小さい乳酸菌にも適用可能であることを示す。本稿

ノム配列ファイル（ASM2650v1.fa）、アノテーションファ

では、事前準備（ゲノム配列とアノテーションファイルをイル（ASM26501v1.gff3）、そして Trimmomatic 7）実行後用いたトランスクリプトーム配列取得）や後処理の解説をの計 6 個の FASTQ ファイル（SRR6322564, SRR6322567,

通じた Galaxy のスキルアップ、そして解析結果として得and SRR6322569 の計 3 サンプル分）をヒストリー trans_5）

map 上に配置して以降の作業を行う［W3］。一通りのコの導出を通じた RNA-seq データの正規化についても議論ピー作業が完了すれば、通常の 3 画面からなる状態に戻しする。

て解析準備完了である［W4］。

＊GFF ファイルの前処理られるカウントデータ・配列長情報・補正後の TPM 値

To whom correspondence should be addressed.
Phone ： +81-3-5841-2395
Fax ： +81-3-5841-1136
E-mail ： kadota@bi.a.u-tokyo.ac.jp

図 1 は、 解 析 準 備 が 完 了 し た、 ① ヒ ス ト リ ー trans_
map の Galaxy 画面である［W5］。中央パネルには、②ア25

Vol. 31, No. 1

Japanese Journal of Lactic Acid Bacteria

図 1．解析準備完了後の Galaxy 画面①ヒストリー名は trans_map。右側のヒストリーパネルには、ゲノム配列ファイル（ヒストリー 1）、アノテーションファイル（ヒストリー 2）、そして Trimmomatic 実行
後の paired-end RNA-seq リードファイル（計 6 ファイル；ヒストリー 3-8）が含まれている。中央パネルには、② GFF ファイルの中身が表示されている。③ Type 列（3 列目に相当）には、様々な種類の features が含まれていることがわかる。

ノテーションファイルの中身が表示されている。今回は、

リー 9（Extract features on data 2）に作成される［W6-

このファイル中の gene 領域を抽出してリファレンス配列10］。図 2 は、④ Extract features 実行結果を⑤中央パネル（マップされる側の配列）として利用する。理由について上に表示した Galaxy 画面である。図 1 では Type 列（3 列は次項で述べるが、Galaxy 上で最初に行う作業は、この目に相当）に様々な features が見られるが、図 2 では⑥

GFF ファイルを入力とした gene 行情報のみの抽出（他のType 列が gene のみになっていることがわかる［W6-12］。

features 行情報の除外）である。この作業は、Galaxy 画面

この gene 領域の座標情報のみからなるヒストリー 9 が、左側のツール選択パネルの Filter and Sort というカテゴ

ヒストリー 1 のゲノム配列ファイルと合わせてトランスク

リに含まれる、Extract features というプログラムを用いリプトーム配列取得に利用される。

て行う［W6-1］。
トランスクリプトーム配列取得

中央パネルの Extract features 操作画面上では、まずSelect GFF data で、入力として与える GFF ファイルを指定する［W6-3］。次に、どの列の情報を用いて目的の前項の前処理の必要性をどのようにして認識したのかにfeatures の行を抽出するかを指定すべく、「3 列目が gene

ついて、実際の作業と戦略立案の関係性を述べる。まずの行を抽出したい」ということを認識させる。ここではま

我々は、ゲノム上の座標情報を取り扱う有名なプログラず 3 列目がターゲットであることを指定すべく、From の

8）ムとして BEDTools が存在することを経験的に知っていところを、デフォルトの「Column 1 / Sequence name」

る。BEDTools は、ゲノム上の座標を表す BED（Browser

か ら、「Column 3 / Feature」に 変 更 す る［W6-5］。 次

Extensible Data）という形式ファイルを取り扱うためのプに、指定した列の中から抽出したい行情報を指定すべく、

ログラム群（なので Tools）を 1 つのパッケージとしてまExtract features のところで「gene」を選択する［W6-6］。

とめたものである。今手元にあるのは BED 形式ファイル

Execute ボタンを押すと、数分程度で実行結果がヒストではなく GFF 形式ファイルではあるものの、ファイル形
26

Jpn. J. Lactic Acid Bact.

Vol. 31, No. 1

図 2．Extract features 実行後の Galaxy 画面① Tools パネルの、② Filter and Sort カテゴリ内にある、③ GFF の、④ Extract
features プログラム実行後の Galaxy 画面。⑤で中央パネル上にプログラム実行結
果を表示させている。⑥ Type 列（3 列目に相当）には、gene のみの feature が含まれていることがわかる。

式の変換（BED ⇔ GFF）を行うプログラムは Galaxy 上でとの座標領域の塩基配列を抽出した multi-FASTA ファ

提供されているだろうという想定のもと、Galaxy のツーイルであり、ヒストリー 10 に作成される［W7-10］。図 3

ル選択パネル上で BEDTools を探すところから実際の作

は、bedtools GetFastaBed 実行結果を中央パネル上に表業を開始した。そして、ツール選択パネルで BED という示した Galaxy 画面である［W7-11］。この multi-FASTA

カテゴリを発見し［W7-1］、その中にリストアップされてファイル中の配列数は 2,949 個であり［W7-12］、第 14 回いるプログラム群を眺めることで、bedtools GetFastaBedの W9-8 で得られたカウント行列の行数と一致する。前回というプログラムが目的を達成する上で最も近いものであと同じ feature（GFF ファイル中の gene 情報）の塩基配列ると判断した［W7-2］。中央パネル上でこのプログラムの情報をリファレンスとして利用することで、次項で得られ操作画面を眺め、入力形式として BED 以外に GFF も取る数値行列との比較が可能となる。

り扱えると判断した［W7-4］。しかしながら、ヒストリー
発現定量

2 に相当するオリジナルの GFF ファイル（ASM26501v1.
gff3）を入力としても gene 領域の塩基配列情報を取得できないことを実体験し、試行錯誤を経て前項の Extract

トランスクリプトーム配列をリファレンスとして発現定量

features プログラムの利用に至った。以降では、前処理9）10）を行う代表的なプログラムとしては、Kallisto や Salmon

後のヒストリー 9 のデータを入力とするが、ヒストリーが挙げられる。これらは、大まかには「マッピングからカウ2 を指定してエラーが出ることを確認してもよいだろう

ント情報取得や RPKM 11）や TPM 5）のような補正後の発現［W7-5］。

量情報の取得」までを 1 つのプログラム内で行うものだと解釈すればよい。今回は、ツール選択パネル上で RNA-seq と前述したように、bedtools GetFastaBed プログラムの入力は、ヒストリー 9 の前処理後の GFF ファイルとヒスいうカテゴリ内に存在する Kallisto quant というプログラム

トリー 1 のゲノム配列ファイルである。出力は、gene ごの利用例を紹介する［W8-1］。
27

Vol. 31, No. 1

Japanese Journal of Lactic Acid Bacteria

図 3．bedtools GetFastaBed 実行後の Galaxy 画面① Tools パネルの、② BED カテゴリ内にある、③ bedtools GetFastaBed プログ
ラム実行後の Galaxy 画面。④で中央パネル上にプログラム実行結果（multi-FASTA
ファイルの中身）を表示させている。どの塩基配列も ATG から始まっており妥当と判断できる。

中央パネルの Kallisto quant 操作画面上では、まずリ

Kallisto はヒトやマウスを対象とした転写物ごとの発現定

ファレンスとして利用するトランスクリプトーム配列の量を目的としている。これらの生物種では、異なる転写物指定を行う。デフォルトでは、Galaxy 内で用意されて間で同じエクソンが共有されうる。このため、このようないるトランスクリプトーム配列を利用する（Use a built-

共有エクソン（shared exon）上にマップされたリード群

in transcriptome）オ プ シ ョ ン に な っ て い る た め、 こ

（正確には k-mer と呼ばれるリード由来の部分配列）をどれ を ヒ ス ト リ ー か ら 選 択 で き る（Use a transcriptomeの転写物に対してどれだけ分配するかという問題に取り組from history）オプションに変更する［W8-2］。変更後、む必要がある。Kallisto は、ブートストラップ（bootstrap）

FASTA reference transcriptome の候補としてヒストリーと呼ばれるテクニックを用いて信頼度の評価を行う枠組み

10（デフォルト）とヒストリー 1 がリストアップされるのを提供している。このテクニックは、データをリサンプリ

で、目的のトランスクリプトーム配列情報からなるヒスング（re-sampling）して同じ解析を行う作業を繰り返すた

トリー 10 を選択する。次に、Trimmomatic 実行結果であめ、一般に繰り返しの回数を増やすほど定量結果のばらつるリードデータが paired-end であることを認識させ［W8-きの推定精度は向上する。しかしながら、その分だけ計

3］、forward 側と reverse 側に相当するヒストリー情報を算時間が増えるだけでなく、試行ごとの計算結果の情報

それぞれ指定して Kallisto quant を実行する［W8-6］。

も増えていくことを正しく認識しておく必要がある。尚、後にこの解析結果を入力として sleuth 12）というプログ

W8-6 では Number of bootstrap samples オプションを変ラムを用いて発現変動解析まで行う際に重要になる事柄と更せずに 0 のままで実行したので、定量結果の推定精度にして、Kallisto quant 実行前に眺めたオプションのところ関する情報は持たない。

に存在する Number of bootstrap samples という項目（デ今回はリファレンス配列と paired-end の 3 サンプル

フォルトは 0）について少し触れておく［W8-6］。まず、分のリードデータを与えたので、入力は計 7 つであった28

Jpn. J. Lactic Acid Bact.

Vol. 31, No. 1

図 4．Kallisto quant 実行後の Galaxy 画面① Tools パネルの RNA-seq カテゴリ内にある、② Kallisto quant プログラム実
行後の Galaxy 画面。③ tabular の実行結果を中央パネル上に表示させている。④ 1列目（target_id）は配列名、⑤ 2 列目はリファレンス配列（この場合は gene）の長さ、⑥ 3 列目は有効配列長（effective length）、4 列目（est_counts）はカウントの推定値［W10-5］、そして 5 列目（tpm）は TPM と呼ばれる補正後の発現量［W11-1］である。

［W8-7］。このときは 10 分足らずで計算が終了し、計 6 個

回ブートストラップは行っていないので［W8-6］、sleuthの出力（ヒストリー 11 ～ 16）が得られた［W8-8］。サンの原著論文タイトルにもある read mapping uncertainty

プルごとに 2 種類の結果（tabular 形式と HDF5 形式）がの情報を現段階では持っていない。それゆえ、今回の場出力され、tabular はタブ区切りテキストファイル（拡張合は実質的に .h5 ファイルは無用であり、.tsv ファイルの子は .tabular）、そして HDF5（Hierarchical Data Format

み を 眺 め れ ば よ い と 判 断 し た［W9-2］。 尚、tsv は tab-5）は HDFView などの専用のビューワでしか見られない

separated values の略であり、タブ区切りテキストファ

バイナリファイル（拡張子は .h5）である［W8-10］。このイルであることを明示した拡張子である。第 8 回
14）

実行結果中にも出現しており、第 9 回

13）

の

15）

のとこ

HDF5 ファイルのみ、定量結果の推定精度に関する情報

DFAST（bootstrap estimates）を含んでいる。

ろでも簡単に触れてはいるものの、意外に知らないヒトが

もちろん我々は、Galaxy 上で得られる情報のみから多いので今回改めて解説した。従って、Galaxy からダウKallisto の全貌を把握しているわけではない。Kallisto のンロードしたファイルの拡張子は .tabular になってはいる

マニュアルも合わせて読むことで、注意すべき点や多少のものの、.tsv や .txt に拡張子を変更しても問題はない。

違いを気にしなくてもよい点を学習している［W9-1］。例
定量結果の解説

えば、我々は今回初めて HDF5 形式を知ったが、Kallisto
マニュアル中の「abundances.h5 is a HDF5 binary file …」を眺めることで、これがバイナリファイルだと認識した

Kallisto quant 実行結果の 1 つである、タブ区切りテキス［W9-2］。また、「This file can be read in by sleuth.」とトファイル（ヒストリー 12, 14, and 16）の中身を解説する

いう記述から、sleuth を用いた発現変動解析を行いたいと［W10-1］。図 4 は、ヒストリー 16 の中身を中央パネル上にきに入力として利用するものだと判断した。同時に、今表示させたものである。ダウンロードしたファイルを Excel
29

Vol. 31, No. 1

Japanese Journal of Lactic Acid Bacteria

で眺めるとより分かりやすいが、このデータは（ヘッダー

実際には、このデータの場合は Trimmomatic によって

行を除くと）2,949 行×5 列からなる。④ 1 列目（target_id

前処理済みであるため、リード長は一定ではない。また、列）は配列名である。ヒストリー 10 の中身を中央パネルで一定の割合で同一配列をもつリードが生成されうるため、表示させた、multi-FASTA ファイルの description 情報にユニークなリードの個数（つまり種類数）で考えると若干相当する［W7-11］。ヒストリー 10 は、ヒストリー 1 のゲノ減る。例えば、上記転写物配列の中から 3 塩基の長さ（μム配列ファイル（ASM2650v1.fa）とヒストリー 9 の gene 領

FLD

域の座標情報を入力として、gene 領域ごとの塩基配列を抜え た 場 合、l effective ＝l original －μFLD＋1＝10－3＋1＝8 個 生 成＝3）をもつリードを何個生成可能か」という問題で考き出したものである。このゲノムの Genbank accession 番

可能ということにはなる。しかし、ユニークなリードの号は FM179322 であり、ASM2650v1.fa 中の 1 行目にも記個数で考えた場合、7 種類（TCA, CAT, ATG, TGG, GGA,

載されている。一部の読者は、例えばヒストリー 16 の（ヘッGAA, and AAT）となり、このような事柄も実際には考慮

ダー行を除く）2 行目の配列名が「FM179322:1523-2663」でせねばならない。これは第 6 回 16）で述べた k -mer 解析にあり、ヒストリー 9 の対応する座標は終点が 2663 と一致しよるゲノムサイズ推定の戦略と基本的に同じである。第 6ているものの、始点が 1524 となっている点に違和感を持つ

回ではμFLD＝251 塩基の長さをもつリード群から、141 塩

かもしれない。しかし、これはゲノムの座標情報から配列基の長さをもつリードの部分配列（つまり k -mer with k

長を算出する際の利便性に関する些細な事柄である。これ＝141）を生成し、得られた k -mer の種類数を推定ゲノムの真の配列長は 1140 塩基である。そして、ヒストリー 9 のサイズとしている（実際には出現頻度が低いものを除いて始点情報を用いて計算すると 2663－1524＋1＝1140 とせねいる）。この例のように k -mer の長さが一定以上あれば、ばならないが、ヒストリー 16 の配列名のように予め始点側

k -mer の種類数をゲノムサイズの推定値として利用可能の数値を－1 しておけば、終点－始点＝2663－1523＝1140

である。

と単純化できる。数値の違いはほぼないに等しいため実害話を元に戻す。ヒストリー 16 の 4 列目（est_counts 列）を被ることは実質的にはないが、注意しておくに越したこ

はカウント数の推定値（estimated counts）である。これとはない。

は、single-end の場合は転写物配列上にマップされたリード数、paired-end の場合はマップされたフラグメント数

ヒ ス ト リ ー 16 の 2 列 目（length 列 ）は 転 写 物 配 列 の
長 さ l original 、 そ し て 3 列 目（eff_length 列 ）は そ の 有 効 長に相当する。“estimated” がついているのは、Kallisto が実（effective length; l effective ）である［W10-3］。配列によって際にはリードそのものではなく、k -mer に分割したもの

多少のばらつきはあるものの、概ね有効長のほうが 200 塩で疑似的なマッピングを行っていることに由来するためだ

基ほど短くなっている（つまり l original －l effective ≒ 200 だといと解釈すればよい。

うこと）ことがわかる。有効長の計算手順を完全に理解でCPM, CPK, FPKM, and TPM

きるほどの情報がマニュアルに記されているわけではな
いが、この約 200 塩基という数値はマップされた pairedend リードの平均フラグメント長（mean of the fragment

ヒ ス ト リ ー 16 の 5 列 目（tpm）は、 い わ ゆ る TPM

length distribution）、つまりマニュアル中のμFLD に相当（Transcripts per million）値 5） と 呼 ば れ る も の で あ る

する［W10-4］。フラグメント長という表現が難解な読者［W11-1］。 マ ニ ュ ア ル 中 の 記 述 に 従 う と、 解 析 サ ン プ

は、まず single-end リードで考えるとよい。長さが全てル で あ る RNA pool 中 の 転 写 物 の 割 合 を 測 定 し た 値（a一定の single-end リードで Kallisto quant を実行した場

measurement of the proportion of transcripts in your合で考えると、このリード長がμFLD に相当する。マニュ

pool of RNA）である［W11-2］。4 列目（est_counts）はカ

アル中の有効長を表現する式に対応させて書くと l effective ＝

ウント情報であり、発現情報ではない。5 列目の数値デー

l original －μFLD＋1 となる。これの意味するところは、「l originalタ（TPM 値）が、一般に発現データと呼ばれるものに相当という長さをもつ転写物配列の中から、μFLD という長する。2014 年出版の関連書籍 17）では（CPM/RPM, CPK/

さのリードを何個生成可能か考えたときに、その答えは

RPK, and FPKM/RPKM には言及しているものの）TPM

l original －μFLD＋1 という式で表すことができる」というこに言及していないが、2019 年出版の書籍 18）では頻繁に出とと本質的に同じである。例えば、「TCATGGAATG と力結果の一部として登場していることがわかる。以下でいう転写物配列（l original ＝10）の中から 5 塩基の長さ（μFLD

は、カウント値から周辺情報を用いて TPM 値を算出する＝5）をもつリードを何個生成可能か」という問題で考え手順を示す。具体的には、入力として 3-4 列目の情報のみ

た 場 合、l effective ＝l original －μFLD＋1＝10－5＋1＝6 が 答 え とを用い、5 列目の情報を出力として得る。

なる。つまり、6 個のリード（TCATG, CATGG, ATGGA,

例として、（ヘッダー行を除く）2 行目の転写物である

TGGAA, GGAAT, and GAATG）を生成可能ということでFM179322:1523-2663 の カウント 値（est_counts 列 の 1703

ある。自由度（degree of freedom）と似た概念だと思えばという数値）から、CPM・CPK・FPKM の算出を経て最よいだろう。

終的に 306.542 という TPM 値を導く［W11-3］。まず CPM30

Jpn. J. Lactic Acid Bact.

Vol. 31, No. 1（Counts Per Million）は、サンプル間で 総カウント数 のの FPKM 値は、以下のように計算する：

違いを補正することを目的としている。具体的には、Per
Million からも予想できるが「総カウント数を 100 万（onemillion）に揃える」だけである。このサンプル（ヒストリー
16 の SRR6322569）のカウント数の総和は 1,758,959 であり、4 列目の数値を全て足したものに相当する［W11-4］。CPMという補正後の値にする作業は、以下のように行う：

l effective の値は遺伝子ごとに異なるため、当然ながらカウ
ント値に掛ける係数も遺伝子ごとに異なる。これは FPKMの総和（または平均値）がサンプルごとに異なるという効
果 を も た ら す。 実 際、 ヒ ストリ ー 16（SRR6322569）、14
解釈としては、例えば「このサンプルの総カウント数は（SRR6322567）、12（SRR6322564）の FPKM 値の総和は、そ100 万より 1.759 倍ほど高いので、カウント値を 1.759 でれぞれ 3,435,847.75、4,005,980.92、5,275,181.61 となる［W11-6

割ってやれば、総カウント数が 100 万だったときのカウン～ W11-8］。この値（FPKM 値の総和）が 100 万に揃うようト数になる」のように理解してもよいだろう。全遺伝子にに補正を行ったものが TPM である。FM179322:1523-2663

対して同じ係数が掛かる（1.759 で割る）ので、このサンプの TPM 値は、以下のように計算する：

ルにおける CPM 値と元のカウント値の関係は以下のようになる：

もちろん係数はサンプルごとに異なる。Galaxy 上で実行済みの他のサンプルとして、ヒストリー 14（SRR6322567）

TPM 値を算出する他の考え方として、「CPK 値の総和とヒストリー 12（SRR6322564）の総カウント数は、それぞを 100 万に揃えた補正を行えばサンプル間で比較可能になれ 1,176,813 と 1,072,837 である。これは、両サンプルともに、る」でもよい。具体的には、FM179322:1523-2663 の CPK

CPM 値を得る際に掛ける係数がヒストリー 16 のものに比値（＝1,852.591）と、ヒストリー 16（SRR6322569）の CPK

べてより 1 に近いことを意味する。特に、ヒストリー 12 の値の総和（＝6,043,515.33［W11-6］）を用いて、以下のようサンプルにおいては、CPM 値が元のカウント値よりもほんに計算する：のわずかに小さくなる（1.073 で割る）程度だということである。ここまでが、連載第 13 回 19）の最後のほうで述べた「同一遺伝子の発現レベルの大小関係を異なるサンプル間で大まかに比較したい場合」に行う作業に相当する。
次に、「同一サンプル内で異なる遺伝子間の発現レベルの大小関係を大まかに比較したい場合」に行う CPK 補正について述べる。これは、遺伝子の長さ情報を用いた補補足情報として、l effective の値はサンプル間で一定ではな正に相当するものであり、Kallisto quant 上では 3 列目

く、若干異なる点にも注意してほしい。つまり、ヒストの eff_length 情 報（l effective ）を 用 い て TPM 値 が 算 出 さ れ

リー 16, 14, 12 の eff_length 列の数値が一定ではないといて い る［W11-5］。CPK は、Counts Per Kilobase の 略 でうことである。これはおそらく、リファレンストランスク

あり、作業的には「l effective を 1,000（one Kilobase）に揃え

リプトーム配列に RNA-seq リード（の部分配列）をマッる」だけである。CPK という補正後の値にする作業は、プした際に、複数配列にマップされる状況がサンプルごとFM179322:1523-2663 のカウント値である 1,703 と、l effectiveに若干異なるためであろう。

＝919.253 の情報を利用して以下のように行う：
原著論文と比較第 14 回で議論した 9 つの遺伝子群の倍率変化および
log2（倍率変化）との値と比較すべく、我々も改めて全サン
FPKM は、CPM と CPK の両方の係数を組み合わせる

プル（計 9 サンプル）について同様の手順で Kallisto quant

ことで、サンプル間補正と遺伝子間補正の両方を同時に行を実行した［W12-1］。表 1 は、②各群 3 反復のカウント

うことを目的とした補正値である。FM179322:1523-2663

データをマージし［W12-2］、③倍率変化と④ log2（倍率変31Vol. 31, No. 1

Japanese Journal of Lactic Acid Bacteria

表 1．原著論文提供情報との比較結果（a）原著論文提供のカウント情報、（b）第 14 回のパイプライン（bowtie2 → htseqcount）で得たカウント情報、（c）今回のパイプライン（Kallisto quant）で得たカウント情報。（a）と（b）は、第 14 回で述べた結果をまとめたものに相当する。①遺伝
子名、②カウント情報、③倍率変化、④ log2（倍率変化）。（c）の遺伝子名は、リファ
レンスゲノム上の遺伝子領域（始点 _ 終点）で示している。例えば、LGG_02372 の遺伝子領域は［2446836, 2449787］と解釈する。

化）を算出した結果である［W13-1］。（a）原著論文提供のて Kallisto quant を実行し、カウント情報および TPM 値

カウント情報と（b）第 14 回の手順で得られたカウント情を得た。Galaxy 上での操作感は、第 14 回の W7 で紹介し報に基づく結果は、第 14 回の原稿中で述べた事柄をまとた Bowtie2 によるマッピングとほぼ同じである。作業のめたものに相当する。（c）が今回得られた Kallisto quant

基本形を一度覚えておけば、Kallisto のような比較的最近

実行結果の中から、①対応する遺伝子領域と②そのカウン開発されたプログラムであっても、Galaxy 上で簡単に実ト情報を抽出した結果である。（c）のカウント値は全体的行できるのが非 Linux ユーザにとっての一番のメリットに大きいものの、その③倍率変化と④ log2（倍率変化）のと言えるだろう。

傾向は、（a）や（b）とよく似ていることがわかる。これは乳酸菌 RNA-seq データの原著論文 3）の Supplement でおそらく、乳酸菌を含むバクテリアは、ほとんどの場合 1

は、RPKM 値 の フ ァ イ ル（GSE107337_RPKM.csv）が 提つの遺伝子領域から転写されるものが 1 種類（遺伝子≒転供されている。なぜ TPM 情報を提供しないのかと疑問に思われるかもしれないが、そのこと自体は問題ではな

写物）だからであろう。

い。TPM 値が欲しければ、自分で RPKM 値の総和が 100
おわりに万になるように補正すればよいだけだからである。TPM
か RPKM/FPKM かの違いは、サンプル間比較をより正

今回は、Trimmomatic で前処理した RNA-seq リード

確に行いたいという思想に基づくものであり、TPM のほに対して、トランスクリプトーム配列をリファレンスとしうが戦略的にベターであることは間違いないだろう。し32

Jpn. J. Lactic Acid Bact.

Vol. 31, No. 1

かしながら（実験デザインにもよるが）サンプル間比較を

Analysis）を行う際に TPM 情報を利用するという考え方

より正確に行いたいのであれば、TPM のような単に数値もあるようだが、我々はカウントデータを利用している。

の総和を一定値に揃えるだけの正規化ではなく、おそら

長さを補正することで異なる遺伝子間の発現レベルの大小20）のような発現変動解析用プログラ

関係は変わりうるが、遺伝子数分だけの要素からなるサン

ムを直接実行するほうがよい。一般によく利用される倍プルベクトル間の類似度をスカラーで要約する際に、入力率変化情報であれば、大抵の解析結果に含まれているかが TPM かカウントかの違いは実質的に問題にはならない

らである。また、2019 年出版の書籍 18）でも言及しているというのが我々のスタンスである。

く edgeRや sleuth12）21）というアルゴリズムを実装した我々

それでもなお、定量 PCR のような他の手段で得られたの TCC パッケージ 22）、あるいはその GUI 版 23）を利用し結果と比較する際に TPM は有用だと思われるかもしれなてもよい。これは edgeR を内部的に繰り返し実行するこい。 し か し こ の 場 合 で も、CPK・RPKM/FPKM・TPMとで、正規化に悪影響を与えうる発現変動遺伝子の影響をの三者はともに比率尺度のものであり、直線性を評価する

できるだけ排除するという思想で開発されたアルゴリズム

際の傾きが多少異なるだけという結果になる。つまり、同である。

一サンプル内での評価という点では、これらは数学的に等が、DEGES 正規化

もちろん発現変動解析用プログラムは、同一サンプル内価である。

で異なる遺伝子間の発現レベルの大小関係を比較する目的
謝では利用できない。しかし少なくとも我々の周辺では、このような目的を設定したヒトを見たことがない。サンプ

辞本連載の一部は、JSPS 科研費 18K11521 の助成を受け

ル間クラスタリングのような探索的な解析（Exploratory

たものです。

		
参 考 文 献1）	寺田朋子，清水謙多郎，門田幸二（2019）次世代シーケンサー
データの解析手法：第 14 回 RNA-seq 解析（その 2）．日本
乳酸菌学会誌 30: 153-161．
2）	Kankainen M, Paulin L, Tynkkynen S, von Ossowski I,
Reunanen J, et al. (2009) Comparative genomic analysis
of Lactobacillus rhamnosus GG reveals pili containing ahuman- mucus binding protein. Proc Natl Acad Sci U S A
106: 17193-17198.
3）	B a n g M , Y o n g C C , K o H J , C h o i I G , O h S . ( 2 0 1 8 )Transcriptional Response and Enhanced Intestinal
Adhesion Ability of Lactobacillus rhamnosus GG after Acid
Stress. J Microbiol Biotechnol 28: 1604-1613.
4）	Afgan E, Baker D, van den Beek M, Blankenberg D,
Bouvier D, et al. (2016) The Galaxy platform for accessible,
reproducible and collaborative biomedical analyses: 2016
update. Nucleic Acids Res 44: W3-W10.
5）	Li B, Ruotti V, Stewart RM, Thomson JA, Dewey CN. (2010)
RNA-Seq gene expression estimation with read mapping
uncertainty. Bioinformatics 26: 493-500.
6）	寺田朋子，大田達郎，清水謙多郎，門田幸二（2018）次世代
シーケンサーデータの解析手法：第 12 回 Galaxy：ヒスト
リーとワークフロー．日本乳酸菌学会誌 29: 79-88．
7）	Bolger AM, Lohse M, Usadel B (2014) Trimmomatic: aflexible trimmer for Illumina sequence data. Bioinformatics
30: 2114-2120.
8）	Quinlan AR, Hall IM (2010) BEDTools: a flexible suite ofutilities for comparing genomic features. Bioinformatics 26:
841-842.
9）	Bray NL, Pimentel H, Melsted P, Pachter L (2016) Nearoptimal probabilistic RNA-seq quantification. Nat Biotechnol
34: 525-527.
10）	Patro R, Duggal G, Love MI, Irizarry RA, Kingsford C (2017)
Salmon provides fast and bias-aware quantification oftranscript expression. Nat Methods 14: 417-419.
11）	Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B(2008) Mapping and quantifying mammalian transcriptomes
by RNA-Seq. Nat Methods 5: 621-628.

12）	P imentel H, Bray NL, Puente S, Melsted P, Pachter L(2017) Differential analysis of RNA-seq incorporating
quantification uncertainty. Nat Methods 14: 687-690.
13）	谷 澤靖洋，神沼英里，中村保一，遠野雅徳，寺田朋子，清水謙多郎，門田幸二（2016）次世代シーケンサーデータの解
析手法：第 8 回アセンブリ後の解析．日本乳酸菌学会誌 27:187-195．
14）	T anizawa Y, Fujisawa T, Kaminuma E, Nakamura Y,
Arita M (2016) DFAST and DAGA: web-based integrated
genome annotation tools and resources. Biosci Microbiota
Food Health 35: 173-184.
15）	谷澤靖洋，真島 淳，藤澤貴智，李 慶範，中村保一，清水謙多郎，門田幸二（2017）次世代シーケンサーデータの解析手
法：第 9 回ゲノムアノテーションとその可視化、DDBJ への登録．日本乳酸菌学会誌 28: 3-11．
16）	谷 澤 靖 洋， 神 沼 英 里， 中 村 保 一， 清 水 謙 多 郎， 門 田 幸 二（2016）次世代シーケンサーデータの解析手法：第 6 回ゲノムアセンブリ．日本乳酸菌学会誌 27: 41-52．
17）	門田幸二（2014）シリーズ Useful R 第 7 巻 トランスクリプ
トーム解析，金明哲 編，共立出版，東京．
18）	坊 農秀雅 編（2019）RNA-Seq データ解析 WET ラボのための鉄板レシピ（実験医学別冊），羊土社，東京．
19）	寺 田朋子，坂本光央，清水謙多郎，門田幸二（2019）次世代
シーケンサーデータの解析手法：第 13 回 RNA-seq 解析（その 1）．日本乳酸菌学会誌 30: 38-45．
20）	Robinson MD, McCarthy DJ, Smyth GK (2010) edgeR: a
Bioconductor package for differential expression analysis ofdigital gene expression data. Bioinformatics 26: 139-140.
21）	Kadota K, Nishiyama T, Shimizu K (2012) A normalization
strategy for comparing tag count data. Algorithms Mol Biol
7: 5.
22）	S un J, Nishiyama T, Shimizu K, Kadota K (2013) TCC:
an R package for comparing tag count data with robust
normalization strategies. BMC Bioinformatics 14: 219.
23）	S u 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.

33

Vol. 31, No. 1

Japanese Journal of Lactic Acid Bacteria

Methods for analyzing next-generation sequencing data
XV. RNA-seq analysis (Part 3)
Tomoko Terada1, Kentaro Shimizu1, 2, and Koji Kadota1, 2
1

Graduate School of Agricultural and Life Sciences, The University of Tokyo.
2
Collaborative Research Institute for Innovative Microbiology,
The University of Tokyo.

		Abstract

Galaxy is an integrative data analysis environment run on the web browser. As before, we explain
a topic to quantify RNA expression from an RNA-seq dataset (GSE107337) which examined the acidstress response of Lactobacillus rhamnosus GG. First, we describe a preprocessing of annotation filewith GFF format and extract nucleotide sequence information of the gene regions on the genome.
Next, we obtain both count and TPM values for individual genes using the Kallisto quant program.
Lastly, we compare the current results with those obtained from the original study and our previous
ones. Supplementary materials are available online at: http://www.iu.a.u-tokyo.ac.jp/~kadota/r_seq2.
html#about_book_JSLAB.

34