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

解説

次世代シーケンサーデータの解析手法
第 7 回 ロングリードアセンブリ
谷澤 靖洋 1, 2、神沼 英里 2*、中村 保一 2、遠野 雅徳 3、
大崎 研 4、清水 謙多郎 5、門田 幸二 5*1

東京大学大学院新領域創成科学研究科
国立遺伝学研究所生命情報研究センター
3
農業・食品産業技術総合研究機構 畜産研究部門
4
トミーデジタルバイオロジー株式会社
5
東京大学大学院農学生命科学研究科
2

完全なゲノム配列の再現を目指す上で最も有効な手段は、ロングリードデータの利用である。第 7 回は、ロングリードの代表格である PacBio データを用いた、乳酸菌ゲノム配列決定について述べる。PacBio デー
タの概観、PacBio データ用のアセンブラである HGAP（Hierarchical Genome Assembly Process）の実行、得られたコンティグの概観・検証（クオリティスコア分布、ドットプロット、および BLAST）、重複領域の除去（環状化）など、一連の手順を解説する。ウェブサイト（R で）塩基配列解析（URL：http://
www.iu.a.u-tokyo.ac.jp/~kadota/r_seq.html） 中 に 本 連 載 を ま と め た 項 目（URL：http://www.iu.a.utokyo.ac.jp/~kadota/r_seq.html#about_book_JSLAB）が存在する。ウェブ資料（以下、W）や関連ウェ
ブサイト、および DDBJ Pipeline 実行結果ファイルなどを効率的に活用してほしい。
Key words：genome assembly, long read, DDBJ Pipeline

PacBio のファイル形式とデータ解析の概要のはこのためである。尚、原著論文中に記載されている
PacBio データに相当する DRR024500（第 6 回の W2-2）は、

第 7 回は、乳酸菌（Lactobacillus hokkaidonensis LOOC260T）
1）

第 6 回原稿執筆時に問題（後述）があったことが判明し、の PacBio データを解析する。デー

2016 年 1 月末に DRR054113-054116 に差し替えになって

タ取得に用いられた PacBio RS II は、SMRT cells と呼ばいる 3）。2016 年 3 月現在、公共 DB 上に DRR024500 は存れるセル単位でシーケンスを行う 2）［W1-1］。リード数を在しない。しかし、第 6 回で用いた Illumina MiSeq デー

増やしたい場合はセル数を増やせばよく、乳酸菌論文 1）

タ（DRR024501）を足がかりにすれば、DRP002401［W2-3］

では 4 セル分のデータを得ている。公共データベース（以や DRA002643［W2-4］のような上位階層の ID を経由し下、DB） 上 で は、 セ ル ご と に 異 な る Run ID が 割 り 振て目的の PacBio データ（DRR054113-054116）に辿り着られる。4 セル分の乳酸菌 PacBio データに対して、4 つくことができる［W2-5］。このデータの場合は、セルあたの DRR ID（DRR054113-054116） が 割 り 振 ら れ て い る

りのリード数が全て 163,482 個になっており、4 セル分をゲノム配列決定論文

合わせると 163,482×4＝653,928 リードになる。
＊

To whom correspondence should be addressed.
Phone ： +81-3-5841-2395
Fax ： +81-3-5841-1136
E-mail ： ekaminum@nig.ac.jp (for DDBJ Pipeline)
E-mail ： kadota@bi.a.u-tokyo.ac.jp (for the others)

PacBio データの解析は、全体として Illumina に代表される他の NGS データ解析とは異なるアプローチが必要である。PacBio の生データは、bax.h5 という拡張子のついた特殊な形式のファイルで出力される。正確には、1 セル分につき、3 つの bax.h5 ファイルと 1 つの bas.h5 ファイ— 101 —

Vol. 27, No. 2		

Japanese Journal of Lactic Acid Bacteria

ルが生成される。そして代表的な PacBio 用 de novo アセらも想像できるように、リード数が非常に少ない（915 個）。

ンブラ（正確にはアセンブリのためのパイプライン）であ

配列長は一定ではなくばらついている（923-8,076 bp）。

る HGAP 4） は、PacBio が提供する SMRT Analysis とい

特に読み始め（最初の 1-3 番目まで）のクオリティスコ

うソフトウェアの一部として提供されており、bax.h5 ファ

アの上下動が激しい。また、配列長分布を眺めてみると、イルのみを入力として受け付ける［W2-6］。その一方で、全リードの約 90% は配列長が 4,000 bp 以下であり、5,000

連載第 3 回5）の「ファイルの解凍、圧縮、概観」でも述

べた通り、公共 DB からエンドユーザが取得可能な形式は、

bp 以上の長いリードは十数個程度しかないことが読み取れる［W3-4］。

sra ファイルと FASTQ ファイルの 2 種類のみである（2016

実際の数値情報を調べたい場合には、例えば（R で）塩

年 4 月 23 日現在）。つまり現状では、HGAP プログラム

基配列解析の項目「前処理 | クオリティチェック | 配列長

実行に必要な bax.h5 ファイルを公共 DB から取得するこ分布を調べる」を参考にすればよい。bzip2 圧縮ファイルとはできない。

から gzip 圧縮ファイルを作成しているのは、FASTQ ファ

本稿では、1 セル分（DRR054113）の bax.h5 ファイル

イルを読み込む際に用いる R パッケージ ShortRead 9） のを公開しているが、ファイルサイズは約 2.4GB に達する

readFastq 関数が bzip2 圧縮ファイルに対応していないた［W2-7］。4 セル分全てだと 10GB 弱になるため、実質的めである［W3-5］。ここでは、配列長ごとの出現頻度を計にノート PC レベルでは解析不可能である。また、bax.h5

算する R スクリプトを実行し［W3-6］、5,000 bp 以上のリーファイルは、テキストではなくバイナリ形式である。このド数が実際には 16 個だったことを確認した［W3-7］。

ため、Linux コマンドの組み合わせで指定したリード数か

次に、DRA で表示されているリード数（163,482 個）とらなるサブセットを抽出し（第 6 回の W3）、再び bax.h5

FastQC 実行結果として得られたリード数（915 個）の極ファイルとして保存するような作業はおそらく誰も行っ

端な違いについて解説する。これは、大まかには DRA で6）でも利用可能であ提供されている sra ファイルと FASTQ ファイルのリードる。結論を先に述べると、前処理などを一切行わない生数の違いに起因する。異なるリード数に関する議論自体の bax.h5 ファイルを DDBJ Pipeline にアップロードし、は、以前に Illumina HiSeq データでも行ったが、ここま

DDBJ Pipeline 上で HGAP を実行するのが最も簡便であでの極端な違いではなかったことを思い出してほしい（第る。理由は、① PacBio データの de novo アセンブリには

3 回の W24）。つまり、sra ファイルが 135,073,834 リード、

256GB 程度のメモリを搭載した Linux マシンが必要であFASTQ ファイルが 134,755,996 リードであり、1% 以下のること、② SMRT Analysis のインストールの難易度が高違いであった。DRA 提供の PacBio の sra ファイルは、生

いこと、③解析プロトコルのステップを理解した上で自分データファイル（3 つの bax.h5 ファイルと 1 つの bas.h5

でプログラムを動作させるのが大変であることなどが挙げファイル）を基に作成している。そして、（PacBio データ

られる。

に限らず）DRA 提供の FASTQ ファイルは、sra ファイていない。HGAP は、DDBJ Pipeline

PacBio データの解析自体は、上述の推奨手順（DDBJ

ルを入力として作成している。全体として、PacBio の場

Pipeline に bax.h5 ファイルをアップロードして実行）に

合は「bax.h5  sra  FASTQ」という不可逆的な流れで、従う限り、極めて簡単である。世間一般の評価同様、少な

sra および FASTQ ファイルが作成される。

くとも乳酸菌を含むバクテリアのゲノム配列決定目的でDRR054113 の sra ファイルは約 1.4GB であり、FASTQ

は、我々の経験上も申し分のない結果が得られる点を強調

ファイルの約 2.3MB とは大きく異なる［W3-8］。これはしておきたい。FastQC 7）によるクオリティチェックも基

bzip2 圧縮ファイルのほうが sra ファイルに比べてサイズ

本的に不必要ではあるが、PacBio の特徴や HGAP アルゴ

が小さい

リズム（計算手順）の理解にも役立つため、周辺情報ととによる影響のほうが圧倒的に大きい。DRA で提供されてもに以下で述べる。

いる FASTQ ファイルのリード数（915 個）が少ないのは、5）という一般的な特徴よりも、リード数の違い

FASTQ ファイル作成時に用いる fastq-dump プログラム
PacBio データの概観 1

実行時に、フィルタリングおよびトリミングのオプション
をつけているためである［W3-9］。前述の DRR024500 で
8）

公共 DB の 1 つである DDBJ SRA（以下、DRA） よ

登録されていたデータの問題点も、ファイル変換時のオプ

り、DRR054113 の FASTQ フ ァ イ ル を ダ ウ ン ロ ー ド しション指定に起因するものであった。フィルタリング・ト［W3-1］、FastQC（ver. 0.11.4）を実行する［W3-2］。リー

リミングされていない状態の PacBio データの特徴を眺め

ド数、配列長分布、クオリティスコア分布のいずれにおいるべく、① fastq-dump を含む SRA Toolkit のインストーても、これまで解説してきた Illumina HiSeq 2000 のデー

ルおよび実行、②オリジナル（163,482 個）に近いリード

タ（SRR616268; 連載第 4 回の W8）や Illumina MiSeq の情報を含む FASTQ ファイルを入力として FastQC を実行

データ（DRR024501; 第 6 回の W4）と全般的に異なるこする。

とがわかる［W3-3］。約 2.4MB というファイルサイズか— 102 —

Jpn. J. Lactic Acid Bact.		

SRA Toolkit（ver. 2.5.7）のインストールと利用

Vol. 27, No. 2

生データの特徴をより正確に把握するためである。このFASTQ フ ァ イ ル（163,380 リ ー ド ） の ク オ リ テ ィ ス コ

sra ファイル から FASTQ ファイルを 作 成 する fastq-

ア は 10 弱 で あ る［W6-3］。PacBio RS II の 塩 基 配 列 決

dump プログラムを利用すべく、NCBI が提供する SRA

定精度は約 87%、つまりエラー率は 13%（p err＝0.13）でToolkit のインストールを行う。本家サイトでは、OS ごと

ある 10）。クオリティスコアは -10×log10（p err）として計の tar.gz ファイルを wget でダウンロードし、tar コマン

算されるため 11）、エラー率 13% のときのスコアは、-10×

ドで解凍するやり方が示されているが［W4-1］、ここではlog10（0.13）＝8.86 となる［W6-4］。DRA から直接ダウンロー

第 4 回で紹介した apt-get コマンドを用いるやり方を示す。

ドした FASTQ ファイル（915 リード）の FastQC 実行結「sudo apt-get install ソフトウェア名 」が基本的なやり方

果は、PacBio の公式スペック以上のクオリティスコア（>

ではあるものの、まずはここで指定するソフトウェア名を20）となっており、一言で言えばきれいすぎる［W3-3］。

把握しなければならない。基本戦略は apt-cache コマンド

著者らは、公称値に近い 10 弱というスコアを見るほうがの利用である。
「apt-cache search キーワード 」として、

むしろ心穏やかである。

ソフトウェア名の一部を構成するであろうと思われる部

配列長は 116 ～ 28,874 bp の範囲となっており［W6-3］、

分文字列をキーワードとして指定すればよい。ここでは、

全リード（163,380 個）の半数以上が 1,000 bp 未満、そし「apt-cache -n search SRA 」として、SRA というキーワーてごく少数のリードが 10,000 bp 以上である［W6-5］。具ドを含むソフトウェア名のみ（-n）リストアップさせ、目体的な数値情報は R で確認することができ、85,134 リー的のソフトウェア名が sra-toolkit であることを見出していド（52.1%）が 1,000 bp 未満であり、5,985 リード（3.66%）

でインストーる［W4-2］。
「sudo apt-get install sra-toolkit 」

が 10,000 bp より長いという結果が得られる［W6-6］。しルを行うメリットは、パスを通すところまで行ってくれる

かしこれは本当の意味での生リード情報であり、アダプ

点にある［W4-4］。実用上は、本家サイトのプログラム名

ター配列を含んでいる。このこと自体は他の NGS テクノを参考にしながら apt-cache で検索し、可能な限り apt-

ロジーと同じではあるが、リードの末端部分にアダプター

get を利用してインストールを行うのが効率的であろう。

を含む Illumina とは異なり、PacBio は SMRTbell と呼ば

DRA で表示されているリード数（163,482 個）に近い情れる鉄アレイのような形のライブラリをテンプレートとし報を含む FASTQ ファイルを得るためには、fastq-dumpてシークエンスを行うため、リード内部にもアダプターをコマンドをデフォルトで実行すればよい［W5-6］。他の含む［W6-7］。このアダプター配列をトリムした後のリー

有用なオプションとしては、--gzip や --bzip2 が挙げらド（adapter-trimmed reads） は サ ブ リ ー ド（subreads）

れる［W5-7］。デフォルトの出力は非圧縮 FASTQ ファと呼ばれる。

イルであるが、これらのオプションをつけることで gzip

PacBio RS II は、セル（SMRT cell）あたり 150,000 個

圧縮や bzip2 圧縮ファイルとして出力させることができのウェルを搭載している。ウェルとは、1 分子（1 つのる。FastQC など、多くのプログラムは圧縮ファイルを

DNA 断片）が入る程度の小さい穴のようなものだと解釈

そのまま入力として与えることができるため、ファイルすればよく、Zero-Mode Waveguides（ZMWs）とも呼ば

サイズ削減の観点からも圧縮ファイルとして出力させてれる。「150,000 ZMWs in PacBio RS II」などと表現されおくとよい。残念ながら、今回の fastq-dump 実行結果

るのはこのためである。各ウェル上に 1 分子のテンプレー

は 163,380 リードであり、DRA で表示されている 163,482

トが入ってシークエンスされるため、理論上ウェル数の分

リードと同じ結果は得られなかった［W5-6; W2-5］。まだけしかリード数は得られない

た、DRA の FAQ に記載されている FASTQ 作成時のオド数の理論値の上限は「150,000 リード（正確には 150,292

プション通りに実行しても、DRA から直接ダウンロード

リード）/ セル」である。DRA で表示されているリード

した FASTQ ファイルのリード数（915 個）にはならなかっ数（163,482 個；W2-5）や、fastq-dump 実行結果としてた［W5-2］。これは用いたプログラムのバージョンの違い得られたリード数（163,380 個；W6-3）は、明らかに理論に起因するのかもしれないが、あくまでも sra  FASTQ

値を超えている。この理由は、PacBio RS II のセルには、ファイル作成時の話であり、PacBio データの推奨解析手

解析サンプルのリード情報を得るためのウェル以外に、コ順（DDBJ Pipeline に bax.h5 ファイルをアップロードしントロール用のウェルや、レーザーが正しく ZMW に照て実行）とは無関係である。

射されているかを確認するためのウェルが搭載されている12）。つまり、得られるリーためである。マイクロアレイにも、ポジティブコントロー
PacBio データの概観 2

ルやネガティブコントロール用のスポット（プローブ）が搭載されている。また、Illumina でシークエンスを行う際

DRA で表示されているリード数（163,482 個）に近いにも、PhiX というコントロール用サンプルが含められる。

163,380 リードからなる FASTQ ファイルを得た主な理由

一般には、これらのコントロールデータの詳細情報までは

は、このファイルを入力として FastQC を実行し、PacBio

記載されないと考えれば納得できるだろう。

— 103 —

Vol. 27, No. 2		

Japanese Journal of Lactic Acid Bacteria

乳酸菌ゲノム配列決定論文 1） では、クオリティスコア

は 6,000 bp 以上の長さのサブリードをシードとして用い（Read Quality＝80）、および配列長（500 bp）でリードのるという意味である。このデータの場合は、4 セル分合わフィルタリングを行った後、そのサブリードを入力としてせて平均 4,000 bp のサブリードが 163,376 個ある 1）。トー

HGAP による de novo アセンブリが行われた。4 セル分のタルの塩基数を乳酸菌の平均ゲノムサイズで割ると、4,000

PacBio データの場合、フィルタリング前の（コントロー

bp×163,376 個 / 2,500,000 bp＝261.4 X となり、25X を大ルを除く）リード数は 150,292×4＝601,168 個となる。上きく超える。本稿では 1 セル分しか用いなかったが、それ記前処理（フィルタリングおよびアダプター除去）後に得でも 261.4 / 4＝65.4X の coverage となる。リード全体でられたサブリード数は 163,376 個であり、配列数としては

65.4X あれば、seed となりうる十分に長いリードでも 25X

約 1/4 に激減している［W6-8］。これは乳酸菌データに限っ分確保でき、結果として 4 セルの場合と遜色のないコンた話ではなく、概ね一般的な減少度合いである。

ティグ数（4 contigs）が得られたということなのであろう。
HGAP 法の実行には、PacBio データの特徴、アルゴリズ

DDBJ Pipeline（HGAP の実行）

ムの大まかな理解に加えて、解析サンプルの推定ゲノムサ
イズおよび coverage を把握しておいたほうが結果の解釈

1 セ ル 分（DRR054113） の 3 つ の bax.h5 フ ァ イ ル をも容易になるだろう。

DDBJ Pipeline にアップロードして、HGAP を実行する

尚、乳酸菌原著論文当時は、P4-C2 試薬を用いて平均

一連の手順を示す。ファイルのアップロード自体は、第4,000 bp のサブリードが得られたが、現在の一般的な手順

6 回の W5 で作成した paired-end Illumina MiSeq データ（P6-C4 試薬、20,000 bp ライブラリ、6 時間 movie）で行（DRR024501）の FTP 経由での手順（第 6 回の W14）とうと、8,000－10,000 bp 程度の長さのサブリードが得られ

同 じ で あ る［W7-2］。 登 録（registration） 作 業 は、bax.

る。HGAP の推奨 coverage（リード全体の coverage）は、

h5 ファイルごとに独立して行う必要がある（2016 年 3 月

データにも依存するが概ね 60-100X である 13-14）。
HGAP の計算時間は、DDBJ Pipeline の実体である遺伝

28 日現在；W7-9）。非常に面倒ではあるが、近い将来一研スパコンの混み具合にもよるが、概ね半日～ 1 日であ括登録が可能になるものと期待される。
HGAP 法 4） のアセンブリ手順は、大まかには次のようる［W9-1］。4 セル分で実行した原著論文 1） の解析時は、な流れで行われる：①最も長いサブリード群をシード（リHGAP（Protocol2; ver. 2.0.0）で行い、7 contigs を得た。

ファレンス配列）として用い、②シードに短いサブリード

今回 1 セル分のみで実行した HGAP（Protocol3; ver. 2.2.0）

をマップし、③エラー補正を行った preassembled reads

実行結果は、4 contigs であった［W9-2］。この違いがプを作成し、④それらを用いてアセンブリの本番を行う。

ログラムのバージョンによるものかセル数によるものかま

HGAP 実行時に指定するパラメータは 2 つある。1 つは推では調査していないが、Platanus 15）実行結果の 52 配列（300定ゲノムサイズ（GenomeSize）、そしてもう 1 つはシー

bp 未満の配列を除く；第 6 回の W20-7）に比べて、少なド（seeds）として用いるリードの最短配列長（Minimum

くとも配列数の点ではよい結果が得られている。尚、このSeed Length） で あ る［W8-2］。 本 稿 で は、PacBio デ ー

データの正解は、全て環状で 3 contigs（1 chromosome＋

タのみのが手元にある場合を想定して乳酸菌の平均的なゲ

2 plasmids）、2,400,586 bp（約 2.4MB）である 1）。

ノムサイズ（2.5MB; 2,500,000 bp）を与えたが、Illumina
MiSeq データから得られた推定ゲノムサイズ（約 2.4MB;

HGAP 実行結果の概観と前処理

2,400,000 bp） を 与 え て み て も い い だ ろ う 3）［W8-3］。
PacBio データは、ランダムな位置でシークエンスエラー

DDBJ Pipeline で HGAP（Protocol3; ver. 2.2.0） を 実が起こる。そのため、②で短いリードをマップして、ポジ行した結果（result.zip）には、計 4 つのファイルが含まションごとに出現する塩基の多数決ルールを適用することれる［W9-4］。エンドユーザが欲しい最終結果ファイル

でエラー補正ができる。X に相当する推定ゲノムサイズをは、polished_assembly.fasta である。Linux の基本コマン

与えることで、coverage の計算を行うことができる。

ド［W9-5］や R［W9-6］を駆使して、このファイルの2 つ め の パ ラ メ ー タ で あ る Minimum Seed Length

全体像を大まかに把握し、ある程度予想を立てる。具体は、基本的にデフォルトの 6,000 bp、および Automatic

的には、最も長い 2,289,497 bp のコンティグは、乳酸菌のEstimation で よ い。PacBio の 真 骨 頂 は、 リ ピ ー ト 配 列

平均ゲノムサイズ（約 2.5MB）に近いことから、染色体を超えうる長さのリードを得られる点にある。それゆ（chromosome）だろうと予想した。それ以外の 3 つのコ

え、①で選択するシードは、安定的にアセンブルできる

ンティグ（86,892 bp, 45,853 bp, and 11,372 bp）は、染色25X を超える範囲で、できるだけ長いサブリードに限定

体の一部、プラスミド、ミスアセンブルのいずれかであろするほうがよい。それを自動的に算出するのが Automaticう。もしコンティグが環状であれば、プラスミドである可Estimation である。6,000 bp という値は、25X という条件

能性が高い。そこでまずは、コンティグごとに環状かどうを満たす最短サブリード長が 6,000 bp 未満だった場合にかのチェックを行う。

— 104 —

Jpn. J. Lactic Acid Bact.		

図 1．45,853 bp からなるトリム前の sequence3 の特徴（a）クオリティスコア分布。横軸：塩基ポジション、縦軸：PHRED スコア（高いほどよい）。最初の 1,223 bp までと最後の 1,335 bp がスコア 0。（b）ドットプロッ
ト。最初と最後の約 5,000 bp が重複している（類似している）ことが分かる。
— 105 —

Vol. 27, No. 2

Vol. 27, No. 2		

Japanese Journal of Lactic Acid Bacteria

コンティグ数が少なくコンティグごとに作業を行う場合
は、multi-FASTA ファイルを分割し、コンティグ数分だけ single-FASTA ファイルを作成しておいたほうが効率的な場合もある。ここでは、ファイル分割手段として R を用いるやり方、および自作プログラム（fastaLengthFilter.
py；第 6 回の W12）と Linux コマンドを組み合わせたやり方を示した［W10］。どちらが正解ということはなく、自分の感性に合う手段を用いればよい。
一 般 に、HGAP ア セ ン ブ リ 結 果 と し て 得 ら れ る コ ン
ティグの末端部分のクオリティは、中央部分に比べて低い。HGAP 実行結果には、FASTA ファイル（polished_assembly.fasta）だけでなく FASTQ ファイル（polished_assembly.fastq）も含まれる。ここでは、FASTQ ファイ
ルを入力として、コンティグごとのクオリティスコア分
布を眺めておく。例えば 2 番目に短いコンティグ（45,853
bp; sequence3.fq）のスコア分布の場合（図 1a; W11-9）、
最初の 1,223 bp までと最後の 1,335 bp が連続してスコア0 になっており、中央の 1,224 bp から 44,518 bp までの計
43,295 bp（＝44,518－1,224＋1）がスコア 1 以上になっていると判断できる［W11-12］。これは、コンティグが環状であることを確認したあとに、どの部分をトリミングするかについての合理的な指針を与えるものでもある。
配列のドットプロット
ドットプロット（dot plot）は、比較したい 2 つの配列の類似度を視覚的に評価するために古くから用いられている描画手段である

16）。基本的には、比較する 2 つの配列

を x 軸 y 軸にそれぞれ並べて、同一塩基部分をハイライトさせるだけである。ここでは、単純な塩基配列を用いた

図 2．
（a）仮想環状コンティグ配列 ACTCGTCAGAACTC のドットプロット。灰色で示した最後の 4 塩基が最初の 4塩基と同じで、重複領域に相当する。（b）環状なので重複除去手段は計 5 通り。＊のついた前後 2 塩基づつトリムするやり方が推奨。

ドットプロットの解釈の基礎を述べ、実際の環状コンティ
グ（ ま た は ゲ ノ ム ） の 実 例 を sequence3 で 示 す。Bio-

グの選択肢として、4 塩基重複の場合は計 5 通り存在する

Linux には dotter 17） というドットプロット用プログラム

が、アスタリスク（＊）のついた両末端から概ね同数の塩がプレインストールされているが［W12-1］、これは環状基をトリミングするのが推奨である（図 2b）。その理由は、チェックの本番で用いる。まず seqinr 18）という R パッケー

図 1a で示すようなコンティグ末端部分のクオリティスコ

ジ中の dotPlot 関数を用いて、ドットプロットの基本形を

ア分布を眺めれば納得できるであろう。この単純な例の場示す。環状かどうかのチェックの場合は、同一の配列を並合は、結論としてはどの選択肢を採用しても同じ結果になべて比較する。同一配列間のドットプロットの主な特徴

るものの、数千塩基におよぶ実際のプラスミドコンティグ

は、必ず対角線上に位置する塩基が同じになるという点での重複部分が完全一致となる可能性はほぼゼロという現実

ある［W12-6］。環状の場合は、コンティグの両末端の配的な問題に当てはめるとよい。

列がほぼ同じになる。これは、ドットプロット上で対角線

数万～数百万塩基におよぶ配列のドットプロットの作成と平行の位置に、重複塩基数分だけの長さの一致部分がハは、seqinr パッケージの dotPlot 関数では現実的に難しい

イライトされることで判別できる［W12-8］。

ため［W14-1］、Bio-Linux にプレインストールされてい例えば、ランダムな塩基配列 “ACTCGTCAGA” が真のる dotter プログラムを用いる［W12-1］。HGAP アセンブ

環状ゲノムだと仮定すると、de novo アセンブリによってリ結果として得られた計 4 コンティグのうち、2 番目に短得られるコンティグは “ACTCGTCAGAACTC” のようないコンティグ（45,853 bp; sequence3.fa）が環状かどうか感じとなる。この場合、灰色の 4 塩基分の重複がドットプを確認したい場合は、「dotter sequence3.fa sequence3.fa」

ロット上でハイライトされ（図 2a）、重複除去（塩基のトと打てばよい［W14-2］。ドットプロット全体を眺めると、

リミング）が環状化作業に相当する。末端塩基のトリミン

対角線と平行の直線が配列末端部分に現れていることがわ— 106 —

Jpn. J. Lactic Acid Bact.		

Vol. 27, No. 2

かる（図 1b）。大まかには、最初と最後の約 5,000 bp が重と、45,853 bp か ら な る sequence3.fa の［1, 4884 bp］ と複しており、この重複除去がアセンブル後の後処理に相当［40967, 45853 bp］の範囲が 99% 一致となっていることがする［W14-4］。もちろんより正確な一致領域を調べる必わかる［W15-9］。重複除去の選択肢は、以下に示すよう要はあり、著者らはこの目的のために BLAST 19）を用いる。

に概ね 4,884 通り存在する：

配列相同性検索（BLAST）

・
［4885, 45853 bp］を残す（最初の 4884 bp、および最後の 0 bp 分をトリム）

BLAST（Basic local alignment search tool） は、 局 所・
［4884, 45852 bp］を残す的なアラインメント（local alignment）を行うプログラム（最初の 4883 bp、および最後の 1 bp 分をトリム）である。問い合わせ（query）配列と呼ばれる手元の塩基・
［4883, 45851 bp］を残す

配列またはアミノ酸配列の特徴や性質を把握する目的で、（最初の 4882 bp、および最後の 2 bp 分をトリム）

GenBank などの公共塩基配列 DB に対して検索を行った…

ことがある人は多いだろう。プログラム内部では、query・
［3, 40968 bp］を残す

配列と似たものがあるかどうかを DB 配列に対して検索（最初の 2 bp、および最後の 4885 bp 分をトリム）

し、指定した閾値を満たす query 側と DB 側の部分配列

・
［2, 40967 bp］を残すのアラインメント結果が出力される。ここでは、query 側（最初の 1 bp、および最後の 4886 bp 分をトリム）

および DB 側の配列を sequence3.fa として、（Bio-Linux

・
［1, 40966 bp］を残すにプレインストールされている）BLAST を実行する 20）。

（最初の 0 bp、および最後の 4887 bp 分をトリム）

同 一 配 列 間 の 比 較 な の で、 ト ッ プ ヒ ッ ト（top hit） はsequence3 の全長配列間で 100% 一致のアラインメントと概ねとした理由は、［1, 4884 bp］と［40967, 45853 bp］

なる。今詳細に調べたいアラインメント結果は、セカンドの範囲のアラインメント結果には Gap が含まれており、ヒットの「最初と最後の約 5,000 bp の重複配列」である。

この Gap の取り扱いに関する不確定要素があるためでこれらの予想は、図 1b のドットプロットを事前に眺めてある。範囲と塩基数の関係や計算法を含めて混乱しがち

おけば立てられる。ドットプロットは、コンティグの全体なところではあるが、始端側の［1, 4884 bp］の塩基数像の理解に役立つだけでなく、BLAST 実行結果の理解のは（4884－1＋1）＝4884 bp と 計 算 し、 終 端 側 の［40967,

助けにもなる。補足情報的な位置づけではあるが、なるべ

45853 bp］の塩基数は（45853－40967＋1）＝4887 bp と計く併用するといいだろう。

算する。アラインメント結果の始端側と終端側の塩基数がBLAST の 実 行 は、 ① DB 側 配 列 の BLAST 用 DB へ

異なるのは、結論としては Gap 数の違いに起因するための変換、および② query 配列の相同性検索の 2 ステップ

であり、気にしなくてよい。実際に我々が行う重複除去で 完 了 す る［W15］。 具 体 的 に は、 ① で は makeblastdb

は、概ね両側から同数程度の塩基のトリムである。そのコ マ ン ド で DB 側 配 列 で あ る sequence3.fa を 入 力 と し理由は、アセンブリ結果として得られるコンティグの末て、BLAST 用 DB（インデックスファイル）を作成する

端部分のクオリティは、中央部分に比べて低いためであ［W15-1］。②では、query 側と DB 側の配列の種類（塩基る（図 1a; W11-9）。重複塩基数が 4900 bp 程度であるこ

配列またはアミノ酸配列）や目的に応じて、以下に示す 5とを踏まえ、BLAST アラインメント結果のセカンドまたつのプログラムを使い分ける：はサードヒットの 2400-2500 番目付近を眺め、Gap やミスマッチのない領域でトリミング領域を決定する［W16-1］。

blastn：query 側、DB 側がともに塩基配列。

こ こ で は、tail と cut コ マ ン ド を 組 み 合 わ せ て［2450,

blastp：query 側、DB 側がともにアミノ酸配列。

43422 bp］ の 範 囲 を 抽 出 し、（43422－2450＋1）＝40,973

blastx：query 側は塩基配列、DB 側はアミノ酸配列。

bp の長さの環状コンティグ（sequence3_trimmed.fa）とquery 配列をアミノ酸配列に翻訳して検索。

して出力した［W16-3］。

tblastn：query 側はアミノ酸配列、DB 側は塩基配列。
確

DB 配列をアミノ酸配列に翻訳して検索。

認tblastx：query 側、DB 側がともに塩基配列。両方をア
ミノ酸配列に翻訳して検索。

動作確認（この場合、トリミングが正しく行えたか
どうかの確認）は、複数の手段で行うことを推奨する。

今回は query 側と DB 側がともに塩基配列のため、最BLAST のアラインメント結果からトリム後の塩基配列も一般的な blastn を利用する。計算自体はほぼ一瞬で終の最初と最後の部分がわかるので［W16-2］、トリム後の了し、指定した名前の BLAST 結果ファイル（sequence3_

FASTA フ ァ イ ル（sequence3_trimmed.fa） の 最 初 と 最

blast.txt）が作成される［W15-2］。BLAST 結果を眺める

後を眺めておくのは基本であろう［W16-4］。トリム後の— 107 —

Vol. 27, No. 2		

Japanese Journal of Lactic Acid Bacteria

図 3．40,973 bp からなるトリム後の sequence3 の特徴（a）クオリティスコア分布。図 1a で見られた両端の低クオリティ領域がなくなっていることがわかる。
（b）ドットプロット。図 1b で見られた右上と左下の対角線（配列類似）領域がなくなっていることがわかる。
— 108 —

Jpn. J. Lactic Acid Bact.		

Vol. 27, No. 2

FASTQ ファイル（sequence3_trimmed.fq）を作成して、ことで、より詳細なアセンブリ結果の評価を行うこともでトリム後のクオリティスコア分布（図 3a; W17-2）をトリ

きる。

ム前のもの（図 1a; W11-9）と比較しておくのも重要でおわりにあろう。FASTQ 形式ファイル（sequence3.fq）について
は、塩基配列情報からなる 2 行目とクオリティ情報からな

第 7 回は、ロングリードの代表格である PacBio データ

る 4 行目についてのみ「cut -v 2450-43422」を実行すれの特徴を述べ、DDBJ Pipeline 上での HGAP アセンブリ

ばよい［W17-1］。これは、FASTQ 形式の中身を理解しプログラムの実行、および後処理の一部を示した。現状でていれば、FASTQ ファイルについても任意の領域を手持

は、公共 DB で PacBio の生データ（bax/bas.h5 形式ファちのスキルのみで抽出可能という例でもある。45,853 bp

イル）が提供されていない。また、多くの PacBio 用の解のコンティグから、最初の 2,449 bp および最後の（45,853

析プログラムは、FASTQ ファイルではなく生データファ－43,423＋1）＝2,431 bp をトリムしているので、両側ともイルを入力とするため、一般の研究者が公共 PacBio デーにほぼ同数の塩基をトリムしていることになる。この視点

タにアクセスして解析しづらい状況にあることは間違いなでトリム前後の分布を比較することで、クオリティの低いい。将来的には、公共 DB からの PacBio 生データの提供や、両末端をうまくトリムできていると判断することができ

FASTQ ファイルを入力とした PacBio 用プログラムの提る。トリム後の配列でドットプロットを作成し、重複領域

供が本格化するかもしれない。本稿のデータ取得に用いたが消えていることまで確認しておけば間違いないだろう

NGS 機器は、セルあたり 150,000 ZMWs の PacBio RS II［W17-3］。

System で あ っ た が、2015 年 10 月 に 1,000,000 ZMWs の本 稿 で は ド ッ ト プ ロ ッ ト の 作 成 や BLAST の 実 行 をSequel System がリリースされている。Sequel System がBio-Linux 上で行ったが、例えば NCBI の BLAST ウェブ

出力する生データは、BAM 形式ファイルである。SMRT

サービスを利用してもよい［W18］。今回重複除去して得

Analysis 3.1 では、BAM 形式ファイルを取り扱うとともられた 40,973 bp の長さの環状コンティグ（sequence3_に、bax/bas.h5  BAM コンバータも提供されている。

trimmed.fa; W16-3）は、原著論文中では 40,971 bp のプ

アセンブラについては、今回は DDBJ Pipeline 上で実行

ラスミド（accession number: AP014682）に相当する。2可能な HGAP 法を利用したが、MHAP 法

bp の違いは Gap 由来だと思われる。86,892 bp からなる

ど新しいプログラムも続々と開発されている。環状化（重sequence2 についても、sequence3 と同様の枠組みで重複複除去）については、今回は概念の説明を含めて手作業で除去を行えるので、ぜひチャレンジしてみてほしい。この行ったが、Circlator 22） など自動的に行うプログラムも提sequence2 について重複除去した結果は、原著論文中では案されている。

21）

や Falcon な

81,630 bp の プ ラ ス ミ ド（accession number: AP014681）に相当する。最も長い染色体候補の sequence1（2,289,497

謝辞bp）と最も短い sequence4（11,372 bp）の関係性につい

本連載の一部は、科学技術振興機構 バイオサイエンスては、アノテーション（遺伝子領域および機能予測）結データベースセンター（NBDC）、および情報・システ

果を含めた解釈が必要になるため、乳酸菌に特化したアム研究機構 国立遺伝学研究所（遺伝研）との共同研究

ノテーションツール DFAST（http://dfast.nig.ac.jp）を含（2012-2088, 2013-2070） の 成 果 に よ る も の で す。 ま た、

め次回以降述べる予定である。尚、Illumina MiSeq デー

JSPS 科研費 JP25712032, JP15K06919 の助成を受けたもの

タ（DRR024501）も利用可能なこの乳酸菌データの場合です。DDBJ Pipeline における計算処理は、遺伝研が有すは、マッピング（マップする側が MiSeq データ、マップ

るスーパーコンピュータシステムを利用して行われていされる側が重複除去後の HGAP アセンブリ結果）を行うます。

		参 考 文 献1）	Tanizawa Y, Tohno M, Kaminuma E, Nakamura Y, Arita
M. (2015) Complete genome sequence and analysis of
Lactobacillus hokkaidonensis LOOC260T, a psychrotrophic
lactic acid bacterium isolated from silage. BMC Genomics
16: 240.
2）	Flusberg BA1, Webster DR, Lee JH, Travers KJ, Olivares
EC, et al. (2010) Direct detection of DNA methylation during
single-molecule, real-time sequencing. Nat Methods 7: 461465.
3）	谷澤 靖洋，神沼 英里，中村 保一，清水 謙多郎，門田 幸二（2016）次世代シーケンサーデータの解析手法：第 6 回ゲノムアセンブリ．日本乳酸菌学会誌 27: 41-52．

4）	Chin CS, Alexander DH, Marks P, Klammer AA, Drake J, etal. (2013) Nonhybrid, finished microbial genome assemblies
from long-read SMRT sequencing data. Nat Methods 10:563-569.
5）	孫 建強，三浦 文，清水 謙多郎，門田 幸二（2015）次世代
シーケンサーデータの解析手法：第 3 回 Linux 環境構築から
NGS データ取得まで．日本乳酸菌学会誌 26: 32-41．
6）	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.
7）	Andrews S. (2015) FastQC a quality control tool for highthroughput sequence data, http://www.bioinformatics.

— 109 —

Vol. 27, No. 2		

babraham.ac.uk/projects/fastqc/8）	Mashima J, Kodama Y, Kosuge T, Fujisawa T, Katayama
T, et al. (2016) DNA data bank of Japan (DDBJ) progress
report. Nucleic Acids Res 44: D51-57.
9）	Morgan M, Anders S, Lawrence M, Aboyoun P, Pagès H,
et al. (2009) ShortRead: a bioconductor package for input,
quality assessment and exploration of high-throughput
sequence data. 25: 2607-2608.
10）	R hoads A, Au KF. (2015) PacBio Sequencing and Its
Applications. Genomics Proteomics Bioinformatics 13: 278289.
11）	孫 建強，湯 敏，清水 謙多郎，門田 幸二（2015）次世代シー
ケンサーデータの解析手法：第 4 回クオリティコントロー
ルとプログラムのインストール．日本乳酸菌学会誌 26: 124132．
12）	パ ッ ク マ ン の 挑 戦 ― PacBio シ ー ク エ ン サ ー ― （http://pacbiobrothers.blogspot.jp/2013/05/hgap.html）
13）	U tturkar SM, Klingeman DM, Land ML, Schadt CW,
Doktycz MJ, et al. (2014) Evaluation and validation of denovo and hybrid assembly techniques to derive highquality genome sequences. Bioinformatics 30: 2709-2716.
14）	L iao YC, Lin SH, Lin HH. (2015) Completing bacterial
genome assemblies: strategy and performance comparisons.
Sci Rep 5: 8747.
15）	Kajitani R, Toshimoto K, Noguchi H, Toyoda A, Ogura Y, etJapanese Journal of Lactic Acid Bacteriaal. (2014) Efficient de novo assembly of highly heterozygous
genomes from whole-genome shotgun short reads. Genome
Res 24: 1384-1395.
16）	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.
17）	Sonnhammer EL, Durbin R. (1995) A dot-matrix programwith dynamic threshold control suited for genomic DNA
and protein sequence analysis. Gene 167: GC1-10.
18）	Charif D, Thioulouse J, Lobry JR, Perrière G. (2005) Online
synonymous codon usage analyses with the ade4 andseqinR packages. Bioinformatics 21: 545-547.
19）	Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. (1990)
Basic local alignment search tool. J Mol Biol 215: 403-410.
20）	Zhang Z, Schwartz S, Wagner L, Miller W. (2000) A greedy
algorithm for aligning DNA sequences. J Comput Biol. 7:203-214.
21）	Berlin K, Koren S, Chin CS, Drake JP, Landolin JM, et al.
(2015) Assembling large genomes with single-molecule
sequencing and locality-sensitive hashing. Nat Biotechnol.
33: 623-630.
22）	Hunt M, Silva ND, Otto TD, Parkhill J, Keane JA, et all. (2015)
Circlator: automated circularization of genome assemblies
using long sequencing reads. Genome Biol. 16: 294.

Methods for analyzing next-generation sequencing data
VII. long-read assembly
Yasuhiro Tanizawa1, 2, Eli Kaminuma2, Yasukazu Nakamura2,
Masanori Tohno3, Ken Osaki4, Kentaro Shimizu5,
and Koji Kadota5
Graduate School of Frontier Sciences, The University of Tokyo.
2Center for Information Biology, National Institute of Genetics.
3
Institute of Livestock and Grassland Science,National Agriculture
and Food Research Organization.
4
TOMY Digital Biology, Co., Ltd.
5
Graduate School of Agricultural and Life Sciences, The University of Tokyo.
1		

AbstractLong-read sequencing represented by Pacific Biosciences’ single-molecule real-time (SMRT)
technology has been widely used for microbial genomes. We overview an analysis procedure of

Lactobacillus hokkaidonensis LOOC260T genome using the so-called “PacBio“ data. We describe (i) the
characteristics of PacBio data, (ii) genome assembly using the HGAP program provided in the webbased NGS analysis tool called “DDBJ Pipeline”, (iii) the survey of assembled contigs such as quality
scores and circularization. We also describe a typical strategy of finishing the circular contig (plasmid).

— 110 —