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

解説

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

東京大学 大学院農学生命科学研究科
理化学研究所 バイオリソース研究センター 微生物材料開発室
3
東京大学 微生物科学イノベーション連携研究機構

Lactobacillus rhamnosus は、健康なヒトの胃腸粘膜からしばしば単離されるプロバイオティクス乳酸菌である。今回は、まず L. rhamnosus GG に関する 2 つのゲノム配列決定論文の内容を簡単に紹介する。
そして、Gepard によるゲノムスケールのドットプロットを用いて、一部領域が反転している状況を確認する。次に公共データベース上で検索する際の Tips を述べ、この菌株の酸ストレス応答を調べた RNAseq データ（SRP125628 or GSE107337）を取得したのち、実験デザインの解説を行う。今回の内容は、次回以降解説予定のマッピング・カウントデータ取得・発現変動解析を含む一連の RNA-seq データ解析を行うために必要なファイルの取得や全体像を把握する準備編という位置づけである。ウェブサイト（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

はじめにていた（forward 側が 107 bp、reverse 側が 93 bp）。これらの状況証拠を鑑み、SRR616268 は例題として不適切だと

本連載では、乳酸菌の RNA-seq データを第 3 ～ 5 回で判断した。結果的に第 3 ～ 5 回まで適切とはいえないデー

取り扱っている 1-3）。具体的には、Illumina HiSeq 2000 でタを用いたことに対し、率直に読者にお詫び申し上げる。

取 得 し た L. casei 12A 株 4）の stranded paired-end デ ー

本連載開始（2014 年）から 5 年が経過し、公共データ

タ（SRR616268）である。これが選ばれた主な理由は、第ベース（以下、DB）に登録された乳酸菌の RNA-seq デー

1 ～ 2 回 5-6）で取り扱ったゲノム配列と同じ株由来のデータも増加している。本稿では、L. rhamnosus GG 7）の酸タであるため、RNA-seq リードのリファレンスゲノムへ

ス ト レ ス 応 答 を 調 べ た RNA-seq デ ー タ（SRP125628 orの高いマップ率が期待されたからである。しかしながら

GSE107337）のダウンロードや実験デザインを中心に解説

実際には、第 5 回の図 2 で示されたように、forward 側する。前回同様、予期しない不具合を避けるため、ウェの 100-107 bp 領域が原因で全リードの約 0.4% しかマップ

ブ ブ ラ ウ ザ は Google Chrome ま た は Firefox（Internet

されなかった 3）。また、Illumina のリード長は forward とExplorer は非推奨）を用いてほしい。

reverse ともに同じはずだが、このデータはなぜか異なっ
リファレンスゲノム配列＊

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

リファレンスゲノムとして用いる L. rhamnosus GG は、
健康なヒトの胃腸粘膜からしばしば単離されるヘテロ発酵型の乳酸菌である。このゲノム配列決定に関する原著
38

Jpn. J. Lactic Acid Bact.

Vol. 30, No. 1

論文によると、1 本の環状染色体のみから構成され（ゲFM179322.fasta、y 軸に AP011548.fasta を配置するように

7）

ノムサイズは 3.01MB）、プラスミドはもたない ［W1］。

入力ファイルを指定しているので、配列長を示すドットプ

遺 伝 子 数 は 2,944 個 で あ り、DDBJ/GenBank/EMBL のロットの②右上端が 3010110、③左下端が 3005050 となっ

Accession 番号は FM179322 である。原著論文から辿れるているのは妥当である。実際の塩基数よりも 1 少ないのは、

FM179322 のリンク先は「GenBank 形式」で表示されてい始点が 0 だからだと解釈すればよい。このドットプロット

るが、例えば「FASTA（text）形式」に切り替えることで、から読み取れることは、対角線上に位置する塩基が同じだ

リファ レ ン ス 配 列の 一 般 的 な 形 式 で あ る FASTA 形式という点である。つまり、同一配列間のドットプロットのファイル（FM179322.fasta）として保存することができる

特徴と酷似しているということである。

［W2］。また、GenBank 形式の情報を眺めることで、正しかしながら、全ゲノム同士のドットプロットでは、確なゲノムサイズ（3,010,111 bp）やこの菌株が American

AP011548 のゲノム配列のほうが約 5 KB 短く、約 8.9 KB

Type Culture Collection の 53103 に相当するもの（ATCCの領域が反転しているという状況を読み取ることは当然で53103）であることなどがわかる［W3］。Taxonomy ID がきない。理由は、全ゲノム中に占める当該領域の割合が568703 であることもわかる。

（5,060＋8,880）/ 3,005,051＝0.46% と極めて小さいからであ

特記事項として、この Taxonomy ID のリンク先を辿っる。もちろん Gepard は、注目する領域間のドットプロッていくと、別グループの論文 8）が Reference として提示

トを描画する機能がある。図 1b は、比較する 2 つ配列ともされる点を挙げておく［W4］。このゲノム配列決定論文でに、領域［610000, 640000 bp］に限定して再描画した結果

は、ゲノムサイズは 3,005,051 bp と報告されている［W5］。

である。この領域は反転領域［618415, 627294 bp］を約 1/3

また、2,834 のタンパク質コード遺伝子が予測されておの割合で中央付近に含んでいる。④若干右上側にずれてはり、その内訳は 1,939 個（68%）の既知機能遺伝子、およびいるものの、左上端から右下端にプロットされている直線

895 個（32%）の仮想遺伝子となっている。Accession 番号部分が同一領域に相当する。そして、⑤左下方向から右上は AP011548 と記載されており、前述の FM179322 とは方向にプロットされている直線部分が、反転領域に相当す異なる。この理由は AP011548 の論文で言及されており、る［W9］。このように、メインの直線と垂直に交わるような

FM179322 のほうが登録および論文公開が早かったという形でプロットされるのが反転領域の典型的な特徴である。
対比として、第 7 回の図 1b と比較してもよいだろう。

ことである。AP011548 論文では、両者のゲノム配列の違いについても言及されている。AP011548 のゲノム配列は、

このときに着目していた重複配列領域は、メインの左上

FM179322 中のものよりも約 5 KB（3,010,111－3,005,051＝

端から右下端の直線と平行にプロットされていた。尚、

5,060 bp）短く、約 8.9 KB の領域［618415, 627294 bp］が

図 1b の領域［610000, 640000 bp］に限定したドットプロッ

反転している。後で利用するため、AP011548 の FASTA

トにおいて、④メインの直線が若干右上側にずれている

形式ファイル（AP011548.fasta）も作成しておく［W6］。

理由は、FM179322 にはあるが AP011548 にはない領域が存在するためだと解釈すればよい。例えば、領域［238000,

Gepard でドットプロット

244000 bp］のドットプロットを眺めると理解できるであろう［W10］。

上述した一部領域の反転状況は、第 7 ～ 8 回 9-10）でも解
seqinr でドットプロット

説したドットプロット 11）で確認することができる。この場合は、FM179322 と AP011548 とのゲノム配列を x 軸とy 軸にそれぞれ並べ、同一塩基部分をハイライトさせるこ

ここでは、反転の具体的なイメージを掴んでもらうべとに相当する（図 1）。第 7 回 9）では、数万塩基程度までのく、以下の 2 つの仮想塩基配列を用いて示す。

比較的短い塩基配列同士のドットプロット作成手段としてdotter プログラム 12）を用いた。そして第 8 回 10）では、約配列 k ：ACTCGTAGTCTATCATACGA

230 万塩基（2.3MB）の配列を入力とした場合に、dotter

配列 l ：ACTCGACTATCTGATTACGAが実質的に実行不可能であることを示した。今回の約 300
万塩基の入力ファイルに対しても dotter は適用できない

l の下線部分は、k の下線部分の左右を入れ替えたもの

ため、ここでは平成 29 年度 NGS ハンズオン講習会でも紹に相当し、これが反転された状態である。それゆえ、配列介した Gepard 13）を利用する［W7］。

k と l は下線部分が異なっている。図 2a は、配列 k 同士のGepard は Java プ ロ グ ラ ム で あ り、FM179322 とドットプロットである［W11］。この図は第 7 回でも利用しAP011548 の全ゲノム同士の比較でも、10 秒程度でドットた seqinr パッケージ 14）中の dotPlot 関数を用いて作成してプロットを描画できる（図 1a）［W8］。プログラムによっいるため、図 1 と異なり原点の位置が①左下端になっていて原点の位置は異なるが、① Gepard の場合も dotter とる点に注意してほしい。同一配列同士の比較であるため、同じく左上となっていることが分かる。ここでは x 軸に左下端から右上端にかけて、キレイに直線状のドットが描39

Vol. 30, No. 1

Japanese Journal of Lactic Acid Bacteria

図 1．Gepard による FM179322 と AP011548 のゲノムスケールのドットプロット。
（a）全ゲノム同士の比較。 （b）領域［610000, 640000 bp］の比較。

かれていることがわかる。図 2b は、配列 k と l のドットプの紹介を兼ねて解説した。どちらのゲノム配列も間違

ロットである［W12］。枠で囲った反転領域のみ、メインのい で は な い と 思 わ れ る が、Ensembl Bacteria 15）中 の L.

直線と垂直に交わるような形で、②左上方向から右下方向

rhamnosus GG の 情 報 は FM179322 と 同 じ で あ る。 まにかけて直線状のドットが描かれていることがわかる。

た、FM179322 の原著論文のほうが先に公開されている

ここまでは、RNA-seq 解析を行う際に用いるリファ

ことから、以後は Ensembl Bacteria から提供されている

レンスゲノム配列について述べた。今回は FM179322 とASM2650v1 の ゲ ノ ム（ASM2650v1.fa）お よ び ア ノ テ ー

AP011548 の両方を偶然発見したため、疑問点の解消やション情報（ASM2650v1.gff3）を利用する［W13］。

違いを確認する目的で、ドットプロットの復習や Gepard
40Jpn. J. Lactic Acid Bact.

Vol. 30, No. 1

図 2．seqinr による仮想塩基配列のドットプロット。
（a）同一配列（配列 k vs. k）の場合。 （b）反転領域を含む場合（配列 k vs. l）。

RNA-seq データの取得者らのうち 2 名が試した限りでは、1 名はすんなり成功したが［W20］、後程動作確認したもう 1 名は一度失敗した

こ こ で は、 公 共 DB を 眺 め て L. rhamnosus GG の［W21］。失敗した場合の対応策としては、ENA 上でファRNA-seq データ（SRP125628 or GSE107337）の発見、そイルを②個別にダウンロードするか、R の SRAdb パッの同定、そしてデータ取得までを述べる。

ケージ 21）を利用した一括ダウンロードを試してもよいだ

公共 DB の基礎知識や基本的な利用法については第 3 回でろう［W22］。但し SRAdb を利用する際は、原著論文中に述べている。SRR・DRS・PRJ といった様々な ID の対応記載されている GSE107337 ではなく、その ID から辿れ関係について復習しておくとよい。

る SRP125628 を指定しなければならない。著者らの経験

の原著論文

16）

ま ず NCBI の SRA

17）

で「Lactobacillus RNA-seq」で上、失敗するときは何をやっても失敗するが、翌日のリト

キーワード検索を行うと、14 件のプロジェクトがヒット

ライで何事もなかったかのように成功する場合が多い。

16）

にた

このデータセットは、全 18 ファイル（paired-end なのどり着いた際に確認に用いた GSE107337 は、この 14 件ので 1 サンプルにつき 2 ファイル）合わせても約 6GB と比中から「実験デザインが分かりやすく、それほどサンプル

較的サイズが小さいため、ダウンロードできないという事

数が多くない」という条件を最初に満たした Accession 番

態は想像し難い。しかし万が一そういう事態に遭遇した

した（2019 年 2 月 1 日調べ）
［W14］。後に原著論文号である［W15］。このデータセットは NGS 機器の 1 つでとしても、自分の PC に一旦ダウンロードすることなく、ある Illumina MiSeq を用いて得られており、全部で 9 サ

Galaxy 22）で解析するやり方もある。それが図 3 の右端にンプルからなる［W16］。このデータの原著論文は一見す見えている③ FASTQ files（Galaxy）という名前の列を利るとまだないように見えるが、Registration date が 25-用するやり方である。ENA 上の FASTQ ファイルを直接Nov-2017 であることから、うまく探せば見つかるのでは

Galaxy 上にアップロードする手段と理解すればよく、具という視点で探す。我々は試行錯誤の末、NBCI が運営す体的なやり方については次回以降で述べる予定である。

る発現 DB である GEO 18）の情報（正確には GSE107337 の尚、DRA 上で検索を行う場合には、検索する場所に

著者情報）を用いることで、2018 年 10 月 28 日に公開されも注意してほしい。エンドユーザから見れば GSE107337ている GSE107337 の原著論文 16）にたどり着くことができ

は Accession 番 号 で あ る が、DRASearch の Accession

た［W17］。

という場所で GSE107337 を検索しても何も見つからな第 3 回 で も 述 べ た よ う に、FASTQ フ ァ イ ル の ダ ウ

い［W23］。 し か し、Keyword と い う 場 所 で GSE107337

ン ロ ー ド は ENA 19）ま た は DRA 20）で 行 う。ENA 上 でを検索すると 3 件ヒットする［W24］。検索結果から辿GSE107337 を 検 索 す る と、PRJNA419802 や SRP125628

れ る SRP125628 を 眺 め る こ と で、 確 か に SRA・SRX・などのデータセット全体を指し示す ID や、全 9 サンプ

SRS などのメタデータ情報の存在を確認することができ

ルに付随する様々な ID 情報を俯瞰できる（図 3）［W18］。

るのである。この事実から予想できるように、DRA 上ENA は①一括ダウンロードにも対応している［W19］。著の SRP125628 の URL（http://ddbj.nig.ac.jp/DRASearch/
41Vol. 30, No. 1

Japanese Journal of Lactic Acid Bacteria

図 3．ENA 上の GSE107337 のページ

study?acc＝SRP125628）は、DRASearch の Accession とcomponents in R packages with LGG as the reference

いう場所で SRP125628 を検索した結果としても得られ

genome.” の み で あ る（LGG は L. rhamnosus GG の こる［W25］。原著論文から得られる情報は GSE107337 でと）。LGG の酸ストレスの影響を調べたこの研究では、計はあるものの、現実にはそこから辿れる別の ID である

3 状態の RNA-seq データを各状態につき 3 反復ずつ取SRP125628 を利用するほうが有意義なのである。また、得して比較している［W26］。具体的には、酸ストレス短第 3 回では DRA 上で FASTQ ファイルの提供が遅れて期暴露条件群（pH4.5_1h）と酸ストレス長期暴露条件群いる例を示したが、今回の SRP125628 は FASTQ と SRA（pH4.5_24h）を、対照群（pH7_CCG）と比較する実験デザ

ファイルの両方の提供がまだなされていない例でもある

インである。表 1 は、サンプルごとの SRR ID、リード数、（2019 年 2 月 1 日調べ）。

そして gzip 圧縮状態でのファイルサイズ（単位はバイト）を示したものである［W27］。リード数とファイルサイズ

RNA-seq データ（GSE107337）の詳細

はほぼ比例関係となっており、リード数の最大と最小の間には 10 倍以上の違い（3,869,088 / 301,126＝12.85）がある

GSE107337 の 原 著 論 文

16）で は、 残 念 な が ら Illumina

こともわかる。

MiSeq を用いた塩基配列データ取得以降の詳細な解析

原著論文 16）中では示されていないものの、論文著者ら手順が示されていない。論文中に記載されていることはは GEO の GSE107337 において、データ解析手順に関す“Data analysis was performed using EdgeR, Bioconductor

るより詳細な情報を記載している。例えば、GEO 内で
42

Jpn. J. Lactic Acid Bact.

Vol. 30, No. 1

表1
サンプル名

SRR ID

pH4.5_1h_rep1
pH4.5_1h_rep2
pH4.5_1h_rep3
pH4.5_24h_rep1
pH4.5_24h_rep2
pH4.5_24h_rep3
pH7_CCG_rep1
pH7_CCG_rep2
pH7_CCG_rep3

SRR6322562
SRR6322563
SRR6322564
SRR6322565
SRR6322566
SRR6322567
SRR6322568
SRR6322569
SRR6322570

リード数（片側のみ）  301,126
1,470,602
1,760,461
1,375,368
3,869,088
1,795,874
3,095,834
2,570,876  846,623

ファイルサイズ（forward 側）53,833,728
265,838,592
319,807,488
241,553,408
682,504,192
315,613,184
551,383,040
451,117,056
148,500,480

ファイルサイズ（reverse 側）
59,662,336
296,910,848
354,557,952
268,341,248
759,697,408
353,099,776
615,870,464
510,308,352
167,350,272

pH4.5_1h_rep1 のサンプル ID に相当する GSM2864941 をrep1」と「 リ ー ド 数 最 大 の pH4.5_24h_rep2」間 で 行 う 場

眺めると、リファレンスゲノム配列として NC_013198.1

合を考えるとよい。もしマップされたカウント数のみかを用いていることがわかる［W28］。また、リファレンスらなるカウントデータを補正なしで用いると、おそらくゲノムへのリードのマッピング（アラインメント）手段とほとんど全ての遺伝子のカウント数が pH4.5_24h_rep2 ＞

して Bowtie 23）が用いられており、RPKM 値 24）の作成に

pH4.5_1h_rep1 となるであろう。また、RNA-seq は転写

は edgeR パッケージ 25）が利用されていることもわかる物の断片配列をシークエンスしたものであるため、原理的［W29］。発現変動遺伝子（以下、DEG）同定の入力データに長い転写物ほどカウント数が多くなる傾向にある。それとして用いられるカウントデータは、マッピング結果ファ

ゆえ、遺伝子の長さ情報を用いた補正なしに遺伝子間のカ

イルとアノテーション情報を入力として作成される。一般

ウント数を比較すると、長い転写物ほど多く発現している的なカウントデータの形式は、各行が遺伝子、各列がサンという誤った結果を導いてしまう。このような基本的な補プルからなる数値行列である。カウントデータファイル中

正の考え方を具体化したものが RPKM 値であり、広く普の行数は、アノテーションファイル中の遺伝子数に応じ及している。

て変化する。例えば、遺伝子数が 2,944 個の場合は、2,944
今後の予定

行分の数値情報を含むカウントデータファイルとなる。アノテーションファイル中の遺伝子領域内にマップされた
リード数をカウントしたデータであることが、カウント

GSE107337 論文

16）の著者らは、GEO のサイト上でカウントデータと RPKM のファイルを提供している［W30］。

データと呼ばれる所以である。
当然ながら、アノテーションファイル中の情報は、マッ両者ともに反復データの平均値しか示されていないもの

ピング時に用いたリファレンスゲノムと完全に対応していの［W31］、論文中の記載内容と一致した計算結果になる

なければならない。最初のほうでゲノム配列の違いを議論ことを確認済みである［W32］。次回以降は、この乳酸菌

し た FM179322 と AP011548 を 例 に と る と、「FM179322

RNA-seq データ（GSE107337）を用いて、Galaxy 上でク由来のリファレンスゲノム配列にマッピングした結果ファオリティコントロール→リファレンスゲノム配列へのマッイル」と「AP011548 由来のアノテーションファイル」ではピング→カウントデータ取得までを行う。また、RPKM

同一遺伝子の座標情報が（一部については同じかもしれな値の算出や DEG 同定も行い、原著論文の結果と比較検討いが）異なる。そのため、当然ながら正しいカウント情報する。さらに発展的な解析として、発現パターン分類 26）を取得することができない。大抵の場合は、カウントデーやシルエットスコアの計算 27）なども行う予定である。

タ取得時にアノテーションファイル中で示された遺伝子領
おわりに域の座標情報がリファレンスゲノム中に存在しないなどの理由とともに正しくエラーを出してくれるものの、気を付

今回のストーリー展開は、実際に行った作業の流れ

けるべき大事なポイントである。
RPKM 値は、「カウント情報」と「遺伝子の長さ情報」を（RNA-seq データを取得したのち、同じ菌株のゲノム配列

入力として、「同一遺伝子の発現レベルの大小関係を異なおよびアノテーションファイルを取得）とは異なる。また、

るサンプル間で大まかに比較したい場合」と「同一サンプ

ターゲットとした L. rhamnosus GG のゲノム配列を取得

ル内で異なる遺伝子間の発現レベルの大小関係を大まかする際に、我々は実際には先に AP011548 の原著論文 8）をに比較したい場合」の両方の目的を達成するために補正さ発見した。そして論文中に記載されていたゲノムサイズ

れた数値である。これは、例えば前者の同一遺伝子のサンや遺伝子数が Ensembl Bacteria 15）中の情報と明らかに異プル間比較を、表 1 で示した「リード数最小の pH4.5_1h_なっている点に疑問をもち、試行錯誤した末に FM179322
43Vol. 30, No. 1

Japanese Journal of Lactic Acid Bacteriaの原著論文 7）を発見した。ある程度の全体像を把握した後から振り返ってみれば、AP011548 の原著論文

比べて 1 桁少ないが、今後新たな問題に遭遇する可能性が

8）

を通読す

ゼロでないことは明記しておきたい。

ればよかっただけである。本連載初期（第 3 ～ 5 回）に用重 要 な お 知 ら せ と し て、 第 7 回 で 紹 介 し た DDBJ

いていた RNA-seq データ（SRR616268）が不適切である

Pipeline 28）のサービスが 2019 年 2 月中旬をもって終了す

可能性については、第 5 回原稿執筆時に認識はしていた。

る。DDBJ Pipeline ユーザは、メールによるお知らせを 2

しかし、そうだと確信したのはさらに数年後である。

月初めに受け取っているであろう。また、本連載のウェブ

次回以降本格的に利用予定の乳酸菌 RNA-seq データ

ページを（R で）塩基配列解析から（R で）塩基配列解析の（GSE107337）についても、一抹の不安は残っている。この

サブに変更している。もちろん主要な項目についてはこれ

データは Illumina MiSeq で取得されているが、著者らの経までのリンク先からも辿れるようにしているが、本連載開験上 MiSeq は配列解析用がほとんどであり、発現解析に用始（2014 年）から 5 年が経過し、これまで紹介してきた膨いられるのは稀だからである。著者らは既に Galaxy 上で 1

大な情報の中にはリンク切れなど様々な不具合も生じはじ

サンプル分（SRR6322562）のカウントデータ取得まで終えめているであろう。我々も継続的に注意深くチェックしてており、GEO のサイト上で提供されているカウントデータいるものの、もし発見したら些細なことでもよいのでぜひ（GSE107337_RawCounts.csv）と似た傾向を示すところま指摘してほしい。

では確認した。しかし、表 1 からもわかるように、これま謝でよく発現解析に用いられてきた Illumina HiSeq シリーズ

辞本 連 載 の 一 部 は、JSPS 科 研 費 JP15K06919 お よ びの一般的なリード数（数千万～数億リード）よりも 1 ～ 2 桁少ない。もちろん遺伝子数（～ 3000 遺伝子）もヒトなどと18K11521 の助成を受けたものです。

		
参 考 文 献1）	孫建強，三浦文，清水謙多郎，門田幸二（2015）次世代シーケ ン サ ー デ ー タ の 解 析 手 法： 第 3 回 Linux 環 境 構 築 か ら
NGS データ取得まで．日本乳酸菌学会誌 26: 32-41．
2）	孫建強，湯敏，清水謙多郎，門田幸二（2015）次世代シーケ
ンサーデータの解析手法：第 4 回クオリティコントロールとプログラムのインストール．日本乳酸菌学会誌 26: 124-132．
3）	孫建強，清水謙多郎，門田幸二（2015）次世代シーケンサー
データの解析手法：第 5 回アセンブル、マッピング、そしてQC．日本乳酸菌学会誌 26: 193-201．
4）	Broadbent JR, Neeno-Eckwall EC, Stahl B, Tandee K, Cai H,
et al. (2012) Analysis of the Lactobacillus casei supragenome
and its influence in species evolution and lifestyle
adaptation. BMC Genomics 13: 533.
5）	門田幸二，孫建強，湯敏，西岡輔，清水謙多郎（2014）次世
代シーケンサーデータの解析手法：第 1 回イントロダクショ
ン．日本乳酸菌学会誌 25: 87-94．
6）	孫建強，湯敏，西岡輔，清水謙多郎，門田幸二（2014）次世
代シーケンサーデータの解析手法：第 2 回 GUI 環境からコマンドライン環境へ．日本乳酸菌学会誌 25: 166-174．
7）	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.
8）	Morita H, Toh H, Oshima K, Murakami M, Taylor TD,
et al. (2009) Complete genome sequence of the probiotic
Lactobacillus rhamnosus ATCC 53103. J Bacteriol 191:
7630-7631.
9）	谷澤靖洋，神沼英里，中村保一，遠野雅徳，大崎研，清水謙多郎，門田幸二（2016）次世代シーケンサーデータの解析手
法：第 7 回ロングリードアセンブリ．日本乳酸菌学会誌 27:101-110．
10）	谷 澤靖洋，神沼英里，中村保一，遠野雅徳，寺田朋子，清水謙多郎，門田幸二（2016）次世代シーケンサーデータの解
析手法：第 8 回アセンブリ後の解析．日本乳酸菌学会誌 27:187-195．
11）	Maizel JV Jr, Lenk RP. (1981) Enhanced graphic matrix
analysis of nucleic acid and protein sequences. Proc Natl

Acad Sci U S A 78: 7665-7669.
12）	Sonnhammer EL, Durbin R. (1995) A dot-matrix programwith dynamic threshold control suited for genomic DNA
and protein sequence analysis. Gene 167: GC1-10.
13）	K rumsiek J, Arnold R, Rattei T. (2007) Gepard: a rapid
and sensitive tool for creating dotplots on genome scale.
Bioinformatics 23: 1026-1028.
14）	Charif D, Thioulouse J, Lobry JR, Perrière G. (2005) Online
synonymous codon usage analyses with the ade4 andseqinR packages. Bioinformatics 21: 545-547.
15）	Cunningham F, Achuthan P, Akanni W, Allen J, AmodeMR, et al. (2019) Ensembl 2019. Nucleic Acids Res 47:
D745-D751.
16）	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.
17）	Sayers EW, Agarwala R, Bolton EE, Brister JR, Canese K,
et al. (2019) Database resources of the National Center for
Biotechnology Information. Nucleic Acids Res 47: D23-D28.
18）	Clough E, Barrett T. (2016) The Gene Expression Omnibus
Database. Methods Mol Biol 1418: 93-110.
19）	Toribio AL, Alako B, Amid C, Cerdeño-Tarrága A, Clarke L,
et al. (2017) European Nucleotide Archive in 2016. Nucleic
Acids Res 45: D32-D36.
20）	Kodama Y, Mashima J, Kosuge T, Kaminuma E, Ogasawara
O, et al. (2018) DNA Data Bank of Japan: 30th anniversary.
Nucleic Acids Res 46: D30-D35.
21）	Zhu Y, Stephens RM, Meltzer PS, Davis SR. (2013) SRAdb:
query and use public next-generation sequencing datafrom within R. BMC Bioinformatics 14: 19.
22）	A fgan 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.
23）	L angmead B, Trapnell C, Pop M, Salzberg SL. (2009)
Ultrafast and memory-efficient alignment of short DNA
sequences to the human genome. Genome Biol 10: R25.
24）	Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B.

44

Jpn. J. Lactic Acid Bact.

Vol. 30, No. 1(2008) Mapping and quantifying mammalian transcriptomes
by RNA-Seq. Nat Methods 5: 621-628.
25）	Robinson MD, McCarthy DJ, Smyth GK. (2010) edgeR: a
Bioconductor package for differential expression analysis ofdigital gene expression data. Bioinformatics 26: 139-140.
26）	Tang M, Sun J, Shimizu K, Kadota K. (2015) Evaluation of
methods for differential expression analysis on multi-groupRNA-seq count data. BMC Bioinformatics 16: 361.

27）	Zhao S, Sun J, Shimizu K, Kadota K. (2018) Silhouette Scoresfor Arbitrary Defined Groups in Gene Expression Data andInsights into Differential Expression Results. Biol Proced
Online 20: 5.
28）	Nagasaki H, Mochizuki T, Kodama Y, Saruhashi S, Morizaki
S, et al. (2013) DDBJ read annotation pipeline: a cloudcomputing-based pipeline for high-throughput analysis of
next-generation sequencing data. DNA Res 20: 383-390.

Methods for analyzing next-generation sequencing data
XIII. RNA-seq analysis (Part 1)
Tomoko Terada1, Mitsuo Sakamoto2, Kentaro Shimizu1, 3,
and Koji Kadota1, 3
1Graduate School of Agricultural and Life Sciences, The University of Tokyo.
2
Microbe Division/Japan Collection of Microorganisms,
RIKEN BioResource Research Center.
3
Collaborative Research Institute for Innovative Microbiology,
The University of Tokyo.

		Abstract

Lactobacillus rhamnosus is a probiotic lactic acid bacterium frequently isolated from humangastrointestinal mucosa of healthy individuals. We first introduce two reports regarding the complete
genome sequence of L. rhamnosus GG. We next compare the two genome sequences by using
dotplot with Gepard and confirm an inverted region. We obtain an RNA-seq dataset (SRP125628
or GSE107337) which examined the acid stress response of this strain and explain the experimental
design. The content of this manuscript is necessary to perform the data analysis which will be
explained next time onwards. Supplementary materials are available online at: http://www.iu.a.utokyo.ac.jp/~kadota/r_seq2.html#about_book_JSLAB.

45