生態学のためのメタバーコーディングとDNAバーコーディング

田辺晶史
2017年7月26日

はじめに

本書は、筆者が開発した塩基配列アセンブラー「Assams」と、配列のクラスタリングとホスト生物同定のためのパッケージ「Claident」の利用を普及促進するために執筆を開始しました。 これらのプログラムの利用方法のみならず、データの取得方法からDNAバーコーディングによる生物の同定方法まで解説しています。

細菌などの研究においては既に広く利用されているメタバーコーディングですが、細菌以外の分類群でも幅広い応用が考えられます。 例えば、土壌や動植物体内の真菌、水中の動植物プランクトンなどにも適用可能だと思います。 水中の溶け出している大型動植物の環境DNAでも同様の解析が可能でしょう。 本書では、環境DNAやメタゲノムから、DNAバーコーディングに利用される領域を増幅し、比較的長い解読長の第2世代以降のシーケンサーで配列を読み取って、DNAバーコーディングによって配列のホスト生物を特定することで生物の在・不在を調査する方法を解説します。 また、メタバーコーディングに限らず、様々な分野で応用が進みつつあるDNAバーコーディングを用いて生物を同定する方法についても解説しています。

本書はクリエイティブ・コモンズの表示-継承 2.1 日本ライセンスの下で配布します。 このライセンスの下では、原著作者の明示を行う限り、利用者は自由に本書を複製・頒布・展示することができます。 また、原著作者の明示と本ライセンスまたは互換性のあるライセンスの適用を行う限り、本書を改変した二次著作物の作成・配布も自由に行うことができます。 詳しい使用許諾条件を見るには
http://creativecommons.org/licenses/by-sa/2.1/jp/
をチェックするか、クリエイティブ・コモンズに郵便にてお問い合わせください。 住所は:171 Second Street, Suite 300, San Francisco, California 94105, USA です。

本書が皆さんの役に立つことができましたら幸いです。 この機会を与えて下さった京都大学大学院人間・環境学研究科の東樹宏和博士、水産研究・教育機構中央水産研究所の長井敏博士と、本書をお読みの皆さんに感謝します。

凡例

本書ではコンピュータに入力するコマンドやその結果を表記する際に以下のように記述しています。

# コメント > command argument1 \ argument2 \ argument3↓ output of command > command argument1 argument2 argument3↓ output of command

上記の例ではcommand argument1 argument2 argument3という全く同じコマンドを2回実行しており、コマンド実行後にoutput of commandがコマンドにより表示されています。 ここで、#から改行まではコメントを表しており、入力の必要はありません。 行頭の>とそれに続くスペースはコマンドの入力の開始を表しており、↓までがコマンドとオプションの入力内容になります。 >とそれに続くスペースはあくまで入力の開始を示すためのものですので、入力しないで下さい。 ↓は入力の終端を表し、ここでEnterキーを押すことを指示する記号です。 ↓を入力しないようにして下さい。 なお、コマンドとオプションを見やすくするためにコマンドやオプションの途中に改行を意図的に入れることがありますが、そのような改行の直前には\を記してあります。 したがって、\が直前にある改行はコマンドの終端や改行入力の指示を意味しません。 また、表示環境によってはワードラップ機能により筆者の意図しない改行が入ってしまうことがありますが、これもコマンドの終端や改行入力の指示を意味しませんので注意して下さい。

また、本書では様々なファイルを使用しますが、その内容は以下のように記述しています。

| 1行目の内容 | 2行目の内容

この例では、行頭の|とそれに続くスペースはファイル内の行頭を表しており、ファイル作成の際は入力しないように注意して下さい。 これは、ワードラップ機能による筆者の意図しない改行とファイルに入力すべき改行を区別できるようにするためのものです。

第0章 ソフトウェアのインストールと環境設定

本書では、Debian GNU/Linux jessie (以下Debian)、Ubuntu Linux 14.04 LTS (以下Ubuntu)を利用環境として想定しています。 Windows環境の方は、Linuxのどちらかをインストールして環境を整えて下さい。 動作が遅くても構わないなら、CygwinやWindows 10 build 14393以降で提供されているWindows Subsystem for LinuxなどのWindows上で動作するUNIX互換環境にインストールすることも可能です。 Linuxのインストールには、インストール用CD・DVD・USBメモリが使えます。 HDDやSSDが1台しかない場合、Windows用のEaseUS Partition Masterなどのパーティションリサイズが可能なソフトを使ったり、インストール用CD・DVD・USBメモリ内に用意されている機能を使ってWindows用の領域を縮小して容量を空ける必要があります。 HDDまたはSSDを追加してスペースを確保していただいても良いでしょう。 USBメモリやUSB接続HDDにインストールすることも可能です。 Ubuntuは見た目・操作性の異なる何種類かがありますが、標準のUbuntuよりもXubuntuというのがおすすめです。

Macでも、これらのOSがインストールできます。 空きディスクがない場合はディスクユーティリティを用いてLinuxインストール用のスペースを空ける必要があります。 ディスクユーティリティでHDD・SSDをクリックして、MacOS Xが使用している領域を縮小して下さい。 そしてrEFItやrEFIndというソフトをインストールして起動時に起動デバイス・起動OS選択メニューが出るようにします(インストール後に2回の再起動が必要です)。 そうすると、Linuxインストール用CD・DVD・USBメモリから起動できるようになりますので、そこから起動してインストールして下さい。 くれぐれも既存のOS用領域を誤って削除しないようにご注意下さい。 空いているディスク領域がある場合はディスクユーティリティでの操作は必要ありません。 rEFItやrEFIndを導入してインストールCD・DVD・USBメモリから起動し、空いている領域にLinuxをインストールしていただければ結構です。 またMacでも、インストール先をUSBメモリやUSB接続HDDにすることができます。

なお、本書は64bit対応のIntel・AMD製CPU搭載機しか想定していません。 それ以外の環境でも動作するでしょうが、自力で解決していただく必要があります。 導入するLinuxは64bit版にして下さい。 32bit版では大容量メモリを活かすことができません。

0.1 LinuxへのClaidentおよび同定用データベースと必要なプログラムのインストール

以下のコマンドをターミナルかコンソールで実行して下さい。 必要なものが全てインストールされます。 sudoを利用可能なユーザーである必要があります。 途中、sudoの実行時に何回かパスワードを質問されますので、入力して下さい。

> mkdir -p ~/workingdirectory↓ > cd ~/workingdirectory↓ > wget https://www.claident.org/installClaident_Debian.sh↓ > sh installClaident_Debian.sh↓ > wget https://www.claident.org/installOptions_Debian.sh↓ > sh installOptions_Debian.sh↓ > wget https://www.claident.org/installDB_Debian.sh↓ > sh installDB_Debian.sh↓ > wget https://www.claident.org/installUCHIMEDB_Debian.sh↓ > sh installUCHIMEDB_Debian.sh↓ > cd ..↓ > rm -r workingdirectory↓

標準のインストール先は/usr/local以下になります。 また、Permission deniedと言われた直後にパスワードを尋ねられたりしますが、パスワードに答えることで進行する場合は問題ありません。 これは、最初はsudoなしにユーザー権限でインストールを試み、うまく行かなかった場合(ここでエラーになる)に初めてsudoを使って管理者権限でインストールしようとする(ここでパスワードを尋ねられる)ためです。

もしも外部ネットワークへのアクセスにプロキシを設定する必要がある場合は、上記のコマンド実行の前に以下のコマンドを実行して環境変数を設定しておいて下さい。 これにより、外部へはプロキシを経由してアクセスが行われるようになります。

> export http_proxy=http://server.address:portnumber↓ > export https_proxy=$http_proxy↓ > export ftp_proxy=$http_proxy↓ > export all_proxy=$http_proxy↓

なお、ユーザー名とパスワードを用いた認証が必要なプロキシでは、以下のようにして下さい。

> export http_proxy=http://username:password@server.address:portnumber↓ > export https_proxy=$http_proxy↓ > export ftp_proxy=$http_proxy↓ > export all_proxy=$http_proxy↓

0.1.1 バージョンアップの場合には

全てのプログラムとデータベースを更新する場合は、インストールと同様の手順でコマンドを実行して下さい。 この手順では、Assams、Claident、PEAR、VSEARCH、Metaxa、ITSxが/usr/local以下へインストールされ、/usr/local/share/claident以下へNCBI BLAST+と同定用BLASTデータベース・Taxonomyデータベース、その他の依存プログラム等がインストールされます。 NCBI BLAST+やBLASTデータベースが別途システムにインストールされている場合でも、Claidentが利用するBLAST+やデータベースとは共存可能です。

なお、以下のコマンドを利用することで、一部の更新を無効化して他のものだけ更新することができます。

> mkdir -p ~/workingdirectory↓ > cd ~/workingdirectory↓ # Assamsの更新を無効化 > touch .assams↓ # Claidentの更新を無効化 > touch .claident↓ # PEARの更新を無効化 > touch .pear↓ # VSEARCHの更新を無効化 > touch .vsearch↓ # NCBI BLAST+の更新を無効化 > touch .blast↓ # プログラムを更新 > wget https://www.claident.org/installClaident_Debian.sh↓ > sh installClaident_Debian.sh↓ # sff_extractの更新を無効化 > touch .sffextract↓ # HMMerの更新を無効化 > touch .hmmer↓ # MAFFTの更新を無効化 > touch .mafft↓ # Metaxaの更新を無効化 > touch .metaxa↓ # ITSxの更新を無効化 > touch .itsx↓ # オプションプログラムを更新 > wget https://www.claident.org/installOptions_Debian.sh↓ > sh installOptions_Debian.sh↓ # 同定用データベース「overall」の更新を無効化 > touch .overall↓ # 同定用データベースを更新 > wget https://www.claident.org/installDB_Debian.sh↓ > sh installDB_Debian.sh↓ # キメラ検出用データベース「Claident Databases for UCHIME」の更新を無効化 > touch .cdu↓ # キメラ検出用データベース「rdp」の更新を無効化 > touch .rdp↓ # キメラ検出用データベース「silva」の更新を無効化 > touch .silva↓ # キメラ検出用データベース「unite」の更新を無効化 > touch .unite↓ # キメラ検出用データベースを更新 > wget https://www.claident.org/installUCHIMEDB_Debian.sh↓ > sh installUCHIMEDB_Debian.sh↓ > cd ..↓ > rm -r workingdirectory↓

0.1.2 標準以外の場所にインストールする

前述の方法でインストールすると、標準では/usr/local以下へインストールされます。 コマンドは/usr/local/binにインストールされます。 他のプログラム(特に旧バージョン)との共存を図りたい場合は、インストール場所を変更して、使用する場合に環境変数PATHを変更するのがいいでしょう。 以下のコマンドで標準とは異なる場所へインストールすることができます。

> mkdir -p ~/workingdirectory↓ > cd ~/workingdirectory↓ > export PREFIX=インストール先にしたい場所↓ > wget https://www.claident.org/installClaident_Debian.sh↓ > sh installClaident_Debian.sh↓ > wget https://www.claident.org/installOptions_Debian.sh↓ > sh installOptions_Debian.sh↓ > wget https://www.claident.org/installDB_Debian.sh↓ > sh installDB_Debian.sh↓ > wget https://www.claident.org/installUCHIMEDB_Debian.sh↓ > sh installUCHIMEDB_Debian.sh↓ > cd ..↓ > rm -r workingdirectory↓

なお、使用する際には、以下のコマンドを実行する必要があります。

> export PATH=インストールした場所/bin:$PATH↓

毎回上記コマンドを実行するのが面倒な場合は、~/.bash_profile~/.bashrcの末尾に上記コマンドを記述しておけばいいでしょう。

0.1.3 複数のバージョンを共存させる

Claidentは標準のインストール先にインストールすると以前のバージョンがあっても上書きされますが、前述の通り標準以外の場所にインストールすれば複数のバージョンを共存させることができます。 ただし、ユーザーのホームディレクトリ(/home/ユーザー名)の直下もしくは/etc/claidentに作成される設定ファイル.claidentは共存できませんので、本設定ファイルは使用するバージョンに合わせて切り替える必要があります。 ユーザーのホームディレクトリの設定ファイルが優先的に読み込まれますので、複数バージョンを共存させる場合はそれぞれ別のユーザーを用意してそのユーザーのホームディレクトリ以下にインストールしていただければ、ユーザーを切り替えることでバージョンを切り替えることができます。

第1章 次世代シーケンサーによる大量サンプルの塩基配列決定

ここでは、Roche GSシリーズシーケンサーや、Ion PGM、Illumina MiSeqによるタグ付き塩基配列決定法について簡単に説明します。 これらのシーケンサーは、いずれも400bp以上の長さの塩基配列解読能力があり、メタバーコーディング、DNAバーコーディングに適しています。 ただし、MiSeqは両側から解読したもの(ペアエンドリード)を連結する必要があります。 したがって、両側の配列間でオーバーラップができるよう、できれば400bp程度、長くても500bp程度のアンプリコンができるようにPCRしなくてはなりません。 連結せずに片側ずつ解析することもできますが、reverse側配列は質がかなり低いため、あまり有用ではないと思います。

次世代シーケンサーは非常に多くの塩基配列を解読できますが、かつてはランニングコストが高く、1度に多数のサンプル由来配列を由来サンプルがわかるように解読することができませんでした。 しかし、ランニングコストはかなり下がってきた上、塩基配列に由来サンプル識別用の数塩基の「MID (multiplex identifier)タグ」を予め付加して解読することで、多数のサンプルの塩基配列を混合(multiplex)して1度に解読できるようになりました。 これにより、1サンプル当たりのコストも劇的に低下しています。 MIDタグをバーコードと呼ぶことも多いですが、DNAバーコーディングでは利用する塩基配列をバーコードと呼び、紛らわしいのでここでは単にタグに統一します。 なお、インデックス配列と呼ぶ場合もあります。 他の文書を読まれる場合は注意して下さい。

なお、後の解析では、PCRの際に生成されるキメラ配列や読み間違いの多い配列が多様性の把握の際に問題となりますが、PCRのレプリケートを確保して塩基配列決定を行うことで、全レプリケートで共通して見られる配列のみをキメラでなく読み間違いも少ない配列と考えることができます。 これは、キメラのできる配列の組み合わせとキメラ配列の「継ぎ目」は多数あるが、キメラでない配列は非常に限られていてしかもPCR前からあって観測される数がキメラより多くなりやすいため、真の配列は全てのレプリケートで検出されやすいが、キメラはそうではないこと、同様に読み間違いのパターンも無数にあるが真の配列はただ一通りで真の配列が全レプリケートで検出されやすいことを利用するものです。 ソフトウェアによる処理だけではキメラや読み間違いを十分に除去できませんが、この方法を併用することで検出効率を改善できると期待できます。 キメラと読み間違いの検出後、レプリケートは足し合わせて解析することができます。

1.1 タグ・アダプター配列付きプライマーを用いたPCR

タグを配列に付加するには、タグの付いたプライマーを用いてPCRを行うのが最も簡単で確実です。 PCR方式では初期投資としてタグ付きプライマーを購入する必要があり、タグの種類を多くしようとするとこれに最大数十万円を要します。 次世代シーケンサー用のサンプル調製キットでは、予めサンプル内のDNA末端に特定のアダプター配列のDNAが付加されていることを前提としています。 そこで、以下のような配列のプライマーを用いてPCRを行います。 なお、鋳型DNAはメタゲノムでも単一のゲノムでもどちらでも構いません。

5' ― [アダプター配列] ― [タグ配列] ― [本来のプライマー配列] ― 3'

両側にタグを付けた場合は、増幅産物は以下のような形になります。

5' ― [アダプター配列F] ― [タグ配列F] ― [本来のプライマー配列F] ― [解読したい配列] ― [本来のプライマー配列R逆相補] ― [タグ配列R逆相補] ― [アダプター配列R逆相補] ― 3'

片側から解読する場合、読み始めがタグ配列Fで、本来のプライマー配列F、解読したい配列の順で読むことになります。

どのようなタグを使うかは、Hamady et al. (2008)のSupplementなど、既存の文献を見て決めて下さい。 ただし、タグが必要なのは解読開始側だけですので、片側からしか読まない場合には反対側はタグのないプライマーを使います。 両側から読む場合でも、サンプル数が少ない場合は片側だけでいいのですが、解読ミスのためフォワード配列とリバース配列の対応がおかしくなることがあり、片側しかタグがないとそれが検出できませんので両側に付けることを推奨します。

このようなプライマーを用いてPCRを行うと、プライマーは鋳型DNAにY字型にアニールし、アダプター配列とタグ配列が端に付加された増幅産物ができます。 この増幅産物をできるだけ等濃度で混合してシーケンスサンプルを調製し、シーケンサーで塩基配列を解読します。 シーケンスサンプル調整ステップ以降はメーカーから提供されている純正キットとそのプロトコルに従って下さい。 メーカー提供のプロトコルにはいくつか改善すべき点があり、改善によってデータの品質が向上することもあるようです。 詳しくはその筋の方にお問い合わせ下さい。 濃度の測定は、Nanodropなどの分光光度計では十分な精度で行うことができません。 そのため、最低でもThermoFisher (旧Invitrogen)のQubitなどで測定して下さい。 シーケンスサンプルのアダプター配列をターゲットとして定量PCRによって測定するタイプのキットが最も高精度ですが、これはランニングコストが非常に高いためサンプルが多い場合は適用できないと思います。

なお、プライマー配列もサンプルを区別する用途に使うことができますので、例えば植物のrbcLmatKなどの複数領域配列を同時に各サンプルから得ることもできます。 もちろんこれらは塩基配列データを見ることでも区別できます。 PCRでは、伸長時間を長めに取り、サイクル数を少なくするようにして下さい。 シーケンスサンプル調製に必要なDNA量は多くありませんので、サイクル数をそれほど多くする必要はないのです。 サイクル数が多く伸長時間が短いほど、「途中まで合成された産物」が生成され、それが次のサイクルで別の鋳型に基いて伸長する、ということが多くなります。 このようにして生成されたDNAはキメラDNAと呼ばれ、「存在しない新奇生物の発見」や「多様性の過大評価」に繋がります。 Finnzymes社のPhusionや東洋紡のKODなどの高正確性DNAポリメラーゼを用いるのも対策として有効です。 キメラDNAの生成を抑制するためには、熱変性後にアニーリング温度へと下げるときの速度を下げることがよいという報告があります(Stevens et al., 2013)ので、もし温度の変化速度を設定可能なサーマルサイクラーであれば温度変化速度を下げるようにして下さい。 後述するプログラムによる処理によって除去することもできますが、完全ではありませんので生成を抑えるのに越したことはありません。

増幅の難しい鋳型の場合、夾雑物があっても増えやすくする島津製作所のAmpdirectなどの試薬を用いてPCRを行ったり、そもそもDNAを抽出する段階でホモジナイザーやビーズ破砕機を用いて細胞を破砕した上でDNAを抽出することで収量を増やすなどの対策が必要になります。 破砕の前にディープフリーザーで凍結して、凍結した状態で破砕するのも効果的です。 ポリフェノール類や多糖類を除去するキットを用いる必要がある場合もあるかもしれません。 この辺りの課題は従来の方法と違いはありません。 タグ・アダプターの付いていないプライマーでは増幅できるがタグ・アダプター付きのプライマーだとうまくいかないことがありますが、その場合はタグ・アダプターなしのプライマーでPCRした増幅産物を用意し、プライマー除去・精製して(ExoSAP-ITなどでいいでしょう)から、タグ・アダプターの付いたプライマーで数サイクルのPCRを行い、増幅産物にタグ・アダプターを付加するというやり方でうまくいくことがあります。

1.1.1 中間アダプターを用いたコスト抑制方法

タグ・アダプター配列付きプライマーは非常に長いため、高価になりがちです。 しかも、増幅するプライマーセットごとに用意しなくてはならず本数も多くなってしまいます。 そこで、PCRを2回行うことで、1組のタグ・アダプター付きプライマーを使い回す方法を紹介します。 まず、1回目は以下のプライマーセットを用いてPCRを行います。

5' ― [中間アダプター配列] ― [本来のプライマー配列] ― 3'

これで末端に中間アダプター配列の付いた増幅産物ができます。 中間アダプター配列はシーケンスに必要なアダプター配列と異なるもので、非特異的増幅を起こさないものであれば何でも構いません。 そういう配列は各種知られていますので、既に使用されているものを探してきて下さい。 そして、この増幅産物を鋳型として以下のプライマーセットでさらにPCRを行います。

5' ― [本来のアダプター配列] ― [タグ配列] ― [中間アダプター配列] ― 3'

このような2段階PCRによって、1組のタグ・アダプター付きプライマーを使い回すことができるようになりますので、コスト削減になります。 ただし、PCRを繰り返すため、PCRによる人工的な置換やPCRバイアスが増加しますし、解読できる長さも減少します。 その点には注意が必要です。 なお、増幅産物は以下のような形になります。

5' ― [本来のアダプター配列F] ― [タグ配列F] ― [中間アダプター配列F] ― [本来のプライマー配列F] ― [解読したい配列] ― [本来のプライマー配列R逆相補] ― [中間アダプター配列R逆相補] ― [タグ配列R逆相補] ― [本来のアダプター配列R逆相補] ― 3'

IlluminaのNextera XT Index Kitを用いたmultiplex法(Illumina corporation, 2013)はこれをキットで実現しているものです。 この方法で2つのタグ配列と両側からの解読を行う場合、最初に中間アダプター配列Fの後ろから、つまり本来のプライマー配列Fから読み始め、解読したい配列方向へ解読します。 次に、中間アダプター配列R逆相補の後ろからタグ配列R逆相補を解読します。 3番目にアダプター配列Fの後ろからタグ配列Fを解読し、方向を反転させて、最後に中間アダプター配列Rの後ろから本来のプライマー配列R、および解読したい配列(逆ストランド)を解読します。 それぞれがR1、R2、R3、R4を含むファイル名で保存されます。 R1、R2、R3は同じstrandで、R4だけ逆strandになります。 R2はタグ配列R逆相補で、R3がタグ配列Fです。 ややこしいので注意が必要です。 なお、上記のキットではR1とR4のシーケンスプライマーは中間アダプター配列をターゲットとするものになっているため、本来のプライマー配列がR1とR4のデータ配列に入ってしまい、解読したい領域が500bp程度ある場合には長さが不足する可能性がありますが、本来のプライマーをシーケンスプライマーとして利用することでデータ配列にプライマー部分が入らないようすることも可能です。 ただし、この方法では後述するNの挿入によるシーケンス品質の改善ができません。

1.1.2 Nの挿入によるシーケンス品質の改善

Illuminaでは、フローセル上の蛍光を光学センサーで読み取っていますが、メタゲノムを鋳型としたユニバーサルプライマーでのPCR増幅産物は配列が似ているため、フローセル上でDNAが近接していると区別が困難になります。 また、読み始め(特に最初の12塩基)に「ほとんどの配列がAであるため発光点が少なく真っ暗になる」と、異常と判断して解読を停止してしまうことがあるとのことです。 ゲノムショットガンやRNA-seqでは近接しているDNAが非常によく似ている可能性がずっと低いために区別できる上、真っ暗になることもないのですが、同じ遺伝子座の増幅産物(しかも読み始めがほとんど変異のないプライマー領域)だとまずいわけです。 そこで、本来のプライマー配列と中間アダプター配列の間にわざとNNNNNNを入れてしまうことにします。 そうすると、この部分の解読時に近接したDNAも区別することができ、真っ暗になることもないため、シーケンス品質が改善するとのことです(Nelson et al., 2014)。 MiSeqやHiSeqをお使いの方はお試し下さい。 さらに、この長さをサンプルごとに変えてしまうと、解読するゲノム上の位置がずれるため、全く同じ配列が隣接していても区別できるようになるためにシーケンス品質が改善するそうです(Fadrosh et al., 2014)。 これらの対処法により、PhiX標準配列を混ぜる量を減らすことができるので、一度に得られるデータ量も増加するはずです。

第2章 塩基配列データの前処理

GSシリーズシーケンサーやIon PGMでは、最終出力として*.sffというファイルが得られます。 Illumina MiSeqの場合は*.fastqとなります。 ここでは、これらのファイルに含まれる配列をサンプルごとに分別(demultiplex)して別々のファイルに分けた上で、低品質の配列の除去を行う方法を説明します。 どの機種でもメーカーのプログラムの機能によるdemultiplexは品質の値を見ていなかったりするため、Claidentのコマンドclsplitseqを利用してdemultiplexを行います。 以下ではコマンドは全てターミナルかコンソールで実行して下さい。 基本的なターミナルの使い方の知識は持っているという前提で話を進めます。 ターミナルの操作に慣れていない方は、最低でも付録Cの内容を理解できるようになっておいて下さい。

2.1 SRA/DRA/ERAの登録データやdemultiplex済FASTQの変換

Claidentでは、配列の名前として配列ID__ランID__タグID__プライマーIDを、また拡張子を除いた入力ファイル名としてランID__タグID__プライマーIDを仮定しています。 そのため、SRA/DRA/ERAの登録データやdemultiplex済FASTQをそのまま用いることができません。 そこで、climportfastqというコマンドを用いることで、上記の仮定を満たす入力ファイルに変換することができます。 ペアエンドの場合は変換する前に第2.3.3節の内容に従ってフィルタリングと連結を先に行っておきます。 変換するには、以下のような入力ファイルを用意します。

| 配列ファイル名1 ランID__タグID__プライマーID | 配列ファイル名2 ランID__タグID__プライマーID | 配列ファイル名3 ランID__タグID__プライマーID

ランIDは適当にでっち上げていただいて構いません。 プライマーIDは同じプライマーを用いたファイル間で同じになるようにすれば何でも構いません。 タグIDはファイルごとに異なるものにします。 配列ファイル名と同じでも構いません。

上記の入力ファイルが用意できたら、以下のようにclimportfastqを実行することで、出力フォルダに変換後のファイルが作成されます。

> climportfastq \ --numthreads=CPU数 \ 入力ファイル \ 出力フォルダ↓

変換後、シングルエンドデータであれば第2.2.3節に従って低品質配列の除去を行って下さい。

2.2 GSシリーズシーケンサーやIon PGMの場合

2.2.1 SFFからFASTQの生成

まずはSFFファイルからFASTQ形式のファイルを作成します。

> sff_extract -c SFFファイル名↓

-cオプションを付けることで、読み始めのTCAGが除去されてファイルが作成されます。 もし後述するタグ配列の冒頭にTCAGを付けている場合には、このオプションは使用しないで下さい。 SFFファイルがHOGEHOGE.sffだとすると、HOGEHOGE.fastqというFASTQ形式のファイルとして塩基配列は保存されます。 HOGEHOGE.xmlというファイルもできますが、ここでは使いませんので削除して構いません。 出力される塩基配列はタグ配列で始まり、解読開始側プライマー配列、解読対象の塩基配列データと続き、反対側プライマー配列(の逆相補配列)で終わります。 ただし、全域が読めるとは限りませんので、途中で読み終わっている配列も含んでいます。 なお、ここではsff_extractコマンドを用いていますが、読み始めのTCAGを除いた塩基配列のFASTQ形式ファイルが得られれば、シーケンサーメーカーから提供される純正のコマンドを使っていただいても構いません。 逆に後述するclsplitseqに与えるタグ配列の冒頭にTCAGを追加していただいても構いませんが、TCAG部分の品質値まで厳格にチェックされます。

2.2.2 塩基配列をサンプルごとに分ける

次に、タグ配列とプライマー配列に基づいてサンプルごとに配列を別々のファイルに分離(demultiplex)します。 事前に以下のようなタグ配列を記したFASTAファイル、解読開始側プライマー配列を記したFASTAファイルがそれぞれ必要です。

| >タグID | タグ配列 | >examplesample1 | ACGTACGT

| >プライマーID | プライマー配列 | >exampleprimer1 | ACGTACGTACGTACGTACGT

タグ配列には縮重コード表記は使えませんが、プライマー配列には縮重コード表記を利用可能です。 タグ配列、プライマー配列ともに複数の配列を記述しておくことも可能です。 なお、1.1.1節で説明したような中間アダプターを用いたりした場合は、プライマー配列ファイルは以下のような内容で作成して下さい。

| >プライマーID | 中間アダプター配列プライマー配列 | >exampleprimer1 | TGATACTCGATACGTACGTACGTACGTACGTACGT

つまり、プライマー配列ファイルには、タグ配列と解読対象の配列の間にある配列を記述して下さい。

これらのファイルが用意できたら、以下のコマンドで塩基配列をサンプルごとに別ファイルへ分割します。

> clsplitseq \ --runname=ランID \ --tagfile=タグ配列ファイル \ --primerfile=プライマー配列ファイル \ --minqualtag=27 \ --numthreads=使用するCPU数 \ 入力ファイル \ 出力フォルダ↓

ランIDはランごとに異なるように適当に決めて指定して下さい。 通常はシーケンサーのメーカーが適当に付けているはずなので、それを記述します。 ファイル名やファイル中の配列名などに含まれていますが、メーカーによってルールが異なりプログラムから自動的に取得することが面倒なため、ユーザー側で指定する必要があります。 --minqualtagはタグ配列部分の品質値の下限を指定するオプションです。 この値より品質値が低い塩基がタグ配列に含まれる場合は、その配列は出力されません。 品質値27という値は、Roche GSシリーズシーケンサーの読み取りミスをクラスタリング・コンセンサス配列作成によって除去するのに適した値です(Kunin et al., 2010)。 他機種では別の値の方が良い可能性があります。 品質値30がかなり広く使用されていますので、30にしてもいいでしょう。

タグを利用したマルチプレックスシーケンスを行っていない場合、--tagfileオプションが不要となりますが、単に省略しただけだとこの後のClaidentのコマンドで処理可能な配列ファイルが作成されません。 そのため、ダミーでも構いませんので--indexname=タグIDというオプションを付加して実行して下さい。

出力される配列からは、タグとプライマーの配列部分は除去されています。 タグ配列部分の認識は厳密一致です。 1塩基の読み間違いを許すなどのオプションはありません。 プライマー配列の認識は、Needleman-Wunschアルゴリズムを用いたアライメントを行い、14%まで不一致を許容します(閾値は変更可能)。 出力フォルダには、ランID__タグID__プライマーID.fastq.gzという名前の塩基配列ファイルが保存されます。 clsplitseqでは並列処理によって多数のCPUを活かした高速化が可能です。 CPUが4個の場合、--numthreadsオプションに4を指定すれば最も高速に処理ができるはずです。 ただし、OSやディスクへの書き込み速度がボトルネックになって思うように速度が出ないこともあります。 なお、出力ファイルはGZIPで圧縮されていますので、GZIP圧縮FASTQに対応していないプログラムで読み込むにはgzip -dなどで圧縮を解除しておく必要があります。 これ以降のClaidentのコマンドは対応しているためそのまま使えます。

DDBJ Sequence Read Archive (DRA)など、次世代シーケンサーの生データを登録するサイトには、ここで生成されたGZIP圧縮FASTQを登録します。

数ラン分に分けて大量サンプルをシーケンスした場合

clsplitseqによる処理を複数回行います。 ただし、clsplitseqは出力先フォルダが既に存在しているとエラーになり処理できません。 そこで、2回目以降は--appendオプションを付けて実行します。 例えば以下のようにして下さい。

> clsplitseq \ --runname=ランID \ --tagfile=タグ配列ファイル \ --primerfile=プライマー配列ファイル \ --minqualtag=27 \ --numthreads=使用するCPU数 \ 入力ファイル1 \ 出力フォルダ↓ > clsplitseq \ --runname=ランID \ --tagfile=タグ配列ファイル \ --primerfile=プライマー配列ファイル \ --minqualtag=27 \ --numthreads=使用するCPU数 \ --append \ 入力ファイル2 \ 出力フォルダ↓
異なる長さのタグが混在している場合について

clsplitseqでは、タグの長さは全て同一であることを仮定しています。 その仮定を使うことで処理を高速化しています。 そのため、タグの長さが不統一の場合は、長さの異なるタグごとに別々のFASTAファイルを準備し、以下のようにして複数回clcplitseqによる処理を行って下さい。

> clsplitseq \ --runname=ランID \ --tagfile=タグ配列ファイル1 \ --primerfile=プライマー配列ファイル \ --minqualtag=27 \ --numthreads=使用するCPU数 \ 入力ファイル \ 出力フォルダ↓ > clsplitseq \ --runname=ランID \ --tagfile=タグ配列ファイル2 \ --primerfile=プライマー配列ファイル \ --minqualtag=27 \ --numthreads=使用するCPU数 \ --append \ 入力ファイル \ 出力フォルダ↓
リバースプライマー配列の認識と除去

これまでの手順では、リバースプライマー配列がアニールする部分の配列は除去されていません。 これもできれば除去した方が良いでしょう。 それには、まず以下のような内容のリバースプライマー配列のファイルを用意します。

| >プライマーID | プライマー配列 | >exampleprimer1 | TCAGTCAGTCAGTCAGTCAG

配列は複数記述できますが、フォワードプライマーと1対1対応していることを仮定していますので注意して下さい。 フォワードプライマーと配列数が異なるとエラーになります。 配列順序も対応していることを仮定していますので、リバースプライマーは同一でもフォワードプライマーが異なるサンプルがある場合、それぞれ異なるプライマーIDを割り当てて両方記述する必要があります。

ファイルが用意できたら、clsplitseqを以下のように実行して下さい。

> clsplitseq \ --runname=ランID \ --tagfile=タグ配列ファイル \ --primerfile=プライマー配列ファイル \ --reverseprimerfile=リバースプライマー配列ファイル \ --reversecomplement \ --minqualtag=27 \ --numthreads=使用するCPU数 \ 入力ファイル \ 出力フォルダ↓

このコマンドでは前述のコマンドの処理に加えて、さらにリバースプライマー逆相補配列をNeedleman-Wunschアルゴリズムを用いたアライメントによって15%まで不一致を許して探します。 該当する配列が見つかったら、該当箇所の先頭以降を元の配列から除去します。 該当する配列が見つからなかった場合は、他の条件を満たしていればそのまま出力ファイルに保存されます。 ただし、読み間違いなどによってリバースプライマーの付く部位以降が除去しきれないこともあります。 なお、--needreverseprimerというオプションを追加することで、リバースプライマー逆相補配列に一致する部位が見つからなかった配列は出力されなくなります。 最後まで解読しきれた配列だけにしたい場合はこのオプションを付けて実行してみて下さい。

2.2.3 低品質配列の除去

FASTQ形式の塩基配列には、解読時のデータから推定された塩基の品質値がありますので、その情報に基づいて質の低い配列をこの時点である程度除去してしまうことができます。 以下のようにclfilterseqを実行して下さい。

> clfilterseq \ --minqual=27 \ --minquallen=3 \ --minlen=350 \ --maxlen=400 \ --maxplowqual=0.1 \ --numthreads=使用するCPU数 \ 入力ファイル \ 出力ファイル↓

ここで、--minqualは品質値の下限値で、終末端側からこの値以上の品質値が--minquallenの値以上に連続して出現するまで配列を削っていきます。 その上で、--minlenの長さ未満の配列は除去し、--maxlenより長い場合はこの値になるまで終末端を削ります。 残った配列において、--minqual未満の配列が--maxplowqualの割合より多い場合はその配列は除去します。 通常は出力はファイルになりますが、--output=folderというオプションを与えた場合はフォルダを作成して入力ファイルと同じ名前のファイルに出力します。 既存のフォルダに出力する場合は--appendを追加して下さい。

clsplitseqの出力フォルダ内の全配列ファイルに対して同じ処理を行う場合は、以下のようなコマンドを実行して下さい。

> for f in clsplitseqの出力フォルダ/*.fastq.gz↓ do clfilterseq \ --output=folder \ --append \ --minqual=27 \ --minquallen=3 \ --minlen=350 \ --maxlen=400 \ --maxplowqual=0.1 \ --numthreads=使用するCPU数 \ $f \ 出力フォルダ↓ done↓

2.3 Illumina社シーケンサーの場合

2.3.1 BCLからFASTQの生成

Illumina社のシーケンサーでは、標準でdemultiplexする機能がありますが、タグ配列部分の品質値を一切見ていないため、タグ配列の品質が低い配列がファイルに含まれています。 このため、生データであるBCLファイルから、Illuminaが配布しているbcl2fastq (1.x系と2.x系はどちらでも構いません)を用いて、demultiplexしていないFASTQ配列を生成し、clsplitseqでdemultiplexを行います。 bcl2fastqのインストールは付録に書きましたのでそちらを参照願います。 また、Illumina社のシーケンサーでは300bpのシーケンスでは301bp (150bpのシーケンスでは151bp)のデータが出力されるため、末尾の1塩基は除去する必要があります(何故最初から1塩基除いたデータを出力してくれないのかは知りません)。 下記の例ではbcl2fastqでFASTQ配列生成の際に除去するようになっています。

まず、シーケンサー付属マシンからBaseCallsフォルダを含むランデータフォルダをまるごとbcl2fastqをインストールしたPCへコピーします。 なお、BaseCallsフォルダ内またはランデータフォルダ内にSampleSheet.csvが存在する場合は、これをリネームしておく必要があります。

次に、bcl2fastqの1.x系では以下のようにコマンドを実行します。

> cd BaseCallsフォルダの直上フォルダ↓ > configureBclToFastq.pl \ --fastq-cluster-count 0 \ --use-bases-mask Y300n,Y8,Y8,Y300n \ --input-dir BaseCallsフォルダ \ --output-dir 出力フォルダ↓ > make -j4↓

--fastq-cluster-count 0でファイルの自動分割機能を無効にし、--use-bases-mask Y300n,Y8,Y8,Y300nでフォワード側配列を300bp (末端1bpは切り捨て)、インデックス1 (タグ配列R逆相補)を8bp、インデックス2 (タグ配列F)を8bp、リバース側配列を300bp (末端1bpは切り捨て)とし、それぞれ*_R1_001.fastq.gz*_R2_001.fastq.gz*_R3_001.fastq.gz*_R4_001.fastq.gzとして出力するよう設定しています。 これはDual Index (インデックスは各8bp)で両側から各300bpずつ解読する(300PE)設定でシーケンスした場合ですので、異なる設定でシーケンスした場合は値を適宜変更して下さい。 Single Index (長さ6bp)で片側250bpの解読を行った場合は--use-bases-mask Y250n,Y6となり、Dual Index (長さ各8bp)で片側300bpの場合は--use-bases-mask Y300n,Y8,Y8となります。 make -j4を実行すると、4つのCPUを使用してBCLからFASTQの生成が行われます。 生成されたFASTQはGZIPというプログラムで圧縮されたものになります。 ファイル名の末尾の拡張子.gzはそれを示すものです。 ClaidentはGZIPで圧縮されたFASTQをそのまま扱えるため、圧縮したままで構いません。

bcl2fastqの2.x系を使用している場合は以下のようにコマンドを実行します。

> bcl2fastq \ --processing-threads 4 \ --use-bases-mask Y300n,Y8,Y8,Y300n \ --runfolder-dir ランデータフォルダ \ --output-dir 出力フォルダ↓

--processing-threadsは処理に使用するCPU数を指定するオプションです。 CPUが16個あるPCで専有できるなら16を指定して下さい。 --use-bases-maskオプションの意味は1.x系と同じです。 --runfolder-dirにはランデータフォルダを指定します。 1.xでは--input-dirにBaseCallsフォルダを指定していましたが、2.xでは指定方法が異なりますので注意して下さい。

2.3.2 塩基配列をサンプルごとに分ける

FASTQファイルができたら、タグ配列とプライマー配列に基づいてサンプルごとに配列を別々のファイルに分離(demultiplex)します。 事前に以下のようなタグ配列を記したFASTAファイル、解読開始側プライマー配列を記したFASTAファイルがそれぞれ必要です。 ペアエンドの場合は、当然リバース側の配列を記したファイルも必要です。

| >タグID | タグ配列 | >examplesample1 | ACGTACGT

| >プライマーID | プライマー配列 | >exampleprimer1 | ACGTACGTACGTACGTACGT

タグ配列には縮重コード表記は使えませんが、プライマー配列には縮重コード表記を利用可能です。 タグ配列、プライマー配列ともに複数の配列を記述しておくことも可能ですが、フォワードプライマーとリバースプライマーは配列数も配列順序も対応していることを仮定しています。 リバースプライマーは同一でもフォワードプライマーが異なるサンプルがある場合、それぞれ異なるプライマーIDを割り当てて両方記述する必要があります。 タグに関しても同様です。 なお、プライマーの前にNNNNNNを付加した場合は、それも含めてプライマーを記述して下さい。 サンプルごとにNの長さを変化させている場合は最も長い場合だけを記述して下さい。

これらのファイルが用意できたら、以下のコマンドで塩基配列をサンプルごとに別ファイルへ分割します。

> clsplitseq \ --runname=ランID \ --index1file=インデックス1配列(タグ配列R逆相補)ファイル \ --index2file=インデックス2配列(タグ配列F)ファイル \ --primerfile=プライマー配列ファイル \ --reverseprimerfile=リバースプライマー配列ファイル \ --minqualtag=30 \ --numthreads=使用するCPU数 \ 入力ファイル1 \ 入力ファイル2 \ 入力ファイル3 \ 入力ファイル4 \ 出力フォルダ↓

入力ファイルは*_R1_001.fastq.gz*_R2_001.fastq.gz*_R3_001.fastq.gz*_R4_001.fastq.gzの順で記述して下さい。 --index1fileはIlluminaが言うところのインデックス1、つまりタグ配列R逆相補を指定し、--index2fileにはIlluminaが言うところのインデックス2、つまりタグ配列Fを指定して下さい。 逆にしないように注意して下さい。 フォワードプライマーは14%まで、リバースプライマーは15%まで不一致を許容して探します。 なお、プライマーの前にNNNNNNを付加した場合は、--truncateN=enableを付加して下さい。 このオプションにより、プライマー配列先頭のNとデータ配列の対応する部位は削った上で不一致度を算出するようになります。 Nの長さを変化させていても、その部位を除いて不一致度が算出されるようになるためプライマー配列ファイルにはNが最長の場合だけを記述しておけば全てヒットするようになります。 処理が終わったら、各ファイルの配列数を確認し、Illumina純正のプログラムでdemultiplexしたファイルの配列数と比較して下さい。 正しく処理できていれば、全てIllumina純正プログラムで処理した場合より少なくなっているはずです。

本来のプライマー配列をシーケンスプライマーとして用いた場合、データ配列にはプライマー配列に当たる部位が含まれていません。 この場合は--primerfile--reverseprimerfileオプションが必要ないわけですが、配列の名前をこの後のClaidentのコマンドで処理可能な形式とするためにオプションとして--primername=プライマーIDを付加して実行して下さい。 ダミーでも構いません。

タグを利用したマルチプレックスシーケンスを行っていない場合、--index1file--index2fileオプションが不要となりますが、単に省略しただけだとこの後のClaidentのコマンドで処理可能な配列ファイルが作成されません。 そのため、ダミーでも構いませんので--indexname=タグIDというオプションを付加して実行して下さい。

出力されるファイルはランID__タグID__プライマーID.forward.fastq.gzランID__タグID__プライマーID.reverse.fastq.gzになります。 DDBJ Sequence Read Archive (DRA)など、次世代シーケンサーの生データを登録するサイトには、ここで生成されたGZIP圧縮FASTQを登録します。 なお、DRAへの登録では配列の長さが揃っているかどうかを指定する箇所がありますが、長さが可変のプライマー部位の配列を除去していますから揃っていません(プライマーセットが1つであってもアライメント時にギャップが入る可能性があるため揃わない)ので誤って揃っていると指定しないようにして下さい。

2.3.3 フォワード配列とリバース配列の連結

オーバーラップがある場合

ペアエンドでシーケンスを行った場合には、clconcatpairを用いてフォワード配列とリバース配列を連結します。 フォワード配列とリバース配列にオーバーラップがある場合には、以下のように実行することでVSEARCHを呼び出して連結します。

> clconcatpair \ --mode=OVL \ --numthreads=使用するCPU数 \ 入力フォルダ \ 出力フォルダ↓

以上のコマンドでは、入力フォルダから*.forward.fastq*.reverse.fastqを探しだして自動的に連結します。 末尾に.gz.bz2が付く圧縮ファイルでも構いません。 出力ファイルは指定したフォルダ内に*.fastq.gzとして作成されます。

入力ファイル名が*.forward.fastq*.reverse.fastqになっていない場合は、以下のように実行することで1組のファイルの配列を連結可能です。

> clconcatpair \ --mode=OVL \ --numthreads=使用するCPU数 \ 入力ファイル1 \ 入力ファイル2 \ 出力ファイル↓

入力ファイル1にフォワード配列のFASTQファイル、入力ファイル2にリバース配列のFASTQファイルを指定して下さい。出力ファイルは末尾に.gz.bz2を付けないと自動的に圧縮はされません。

オーバーラップがない場合

フォワード配列とリバース配列間にオーバーラップがない場合には、先に以下のようにclfilterseqを用いて低品質部位のトリミングや低品質配列の除去を行います。

> clfilterseq \ --minqual=30 \ --minquallen=3 \ --minlen=100 \ --maxplowqual=0.1 \ --numthreads=使用するCPU数 \ 入力ファイル1 \ 入力ファイル2 \ 出力フォルダ↓

ここで、--minqualは品質値の下限値で、終末端側からこの値以上の品質値が--minquallenの値以上に連続して出現するまで配列を削っていきます。 その上で、--minlenの長さ未満の配列は除去します。 残った配列において、--minqual未満の配列が--maxplowqualの割合より多い場合はその配列は除去します。 この処理では、ペアエンドのどちらか一方でも除去される場合、ペアとなるもう一方も除去します。 出力はフォルダを作成して入力ファイルと同じ名前のファイルになります。 既存のフォルダに出力する場合は--appendを追加して下さい。 フォルダ内の全ての*.forward.fastq*.reverse.fastqに適用する場合は、以下のように実行して下さい。

> for f in `ls *.forward.fastq.gz | grep -P -o '^[^\.]+'`↓ do clfilterseq \ --minqual=30 \ --minquallen=3 \ --minlen=100 \ --maxplowqual=0.1 \ --numthreads=使用するCPU数 \ $f.forward.fastq.gz \ $f.reverse.fastq.gz \ 出力フォルダ↓ done↓

このように低品質配列を除去した後、以下のようにclconcatpairを用いて連結を行います。

> clconcatpair \ --mode=NON \ --numthreads=使用するCPU数 \ 入力フォルダ \ 出力フォルダ↓

この方法では、以下のようなフォワード配列とリバース配列を入力として仮定します。

5' ― フォワード配列 ― 3' 5' ― リバース配列 ― 3'

これらの配列を、clconcatpair --mode=NONでは以下のように連結します。

5' ― リバース配列(逆相補) ― ACGTACGTACGTACGT ― フォワード配列 ― 3'

ノイジー配列とキメラ配列の除去はそのまま通常の連結したペアエンドデータと全く同様に行い、配列のクラスタリングに用いるclclassseqvと代表配列への生配列のマッピングを行うclrecoverseqvの実行時に--paddinglen=16を加えて実行します。 上記の連結では、連結点を後から識別できるように人工配列であるACGTACGTACGTACGTを加えているため、配列の類似度が過大になってしまいます。 そこで、完全一致するはずのACGTACGTACGTACGTの長さ分を配列の類似度計算から除外することで正しい配列類似度に基づくクラスタリングを行います。

ホスト生物の同定をする際には、ACGTACGTACGTACGTを探して配列を分割し、フォワード配列とリバース配列で別々に同定し、第7.3節の内容に従って2つの同定結果を統合します。 一般にフォワード配列の方が品質が高いはずですので、フォワード配列とリバース配列の分子同定能力や変異性について何も事前情報がない場合はフォワード配列に基づく同定結果を優先する設定で統合するのが良いでしょう。 配列の分割はcldivseqを以下のように実行することで行えます。

>cldivseq \ --query=ACGTACGTACGTACGT \ 入力ファイル \ 出力ファイル1 \ 出力ファイル2

出力ファイル1がリバース配列の逆相補配列、出力ファイル2がフォワード配列になります。

PEARを用いて連結する場合

フォワード配列とリバース配列にオーバーラップがある場合には、PEAR (Zhang et al., 2014)を用いて連結することもできます。 PEARの場合はdemultiplexした配列ファイルのあるフォルダで、以下のコマンドを実行して下さい。

> pear \ -p 0.0001 \ -u 0 \ -j 使用するCPU数 \ -f ランID__タグID__プライマーID.forward.fastq.gz \ -r ランID__タグID__プライマーID.reverse.fastq.gz \ -o ランID__タグID__プライマーID↓

正常に処理が終わると、以下のファイルができます。

ランID__タグID__プライマーID.assembled.fastq

連結された配列

ランID__タグID__プライマーID.unassembled.forward.fastq

連結されなかったフォワード配列

ランID__タグID__プライマーID.unassembled.reverse.fastq

連結されなかったリバース配列

ランID__タグID__プライマーID.discarded.fastq

検定で一致していると判定されなかった配列

この後で使用するのはランID__タグID__プライマーID.assembled.fastqだけです。 また、出力は圧縮されていないため、ディスク容量の節約のためにも圧縮しておいた方がいいかもしれません。

なお、フォルダ内の全てのファイルに対してPEARによる連結を連続して行うには、以下のようなコマンドをお使い下さい。

> for f in `ls *.forward.fastq.gz | grep -P -o '^[^\.]+'`↓ do pear \ -p 0.0001 \ -u 0 \ -j 使用するCPU数 \ -f $f.forward.fastq.gz \ -r $f.reverse.fastq.gz \ -o $f↓ done↓
VSEARCHを用いて連結する場合

フォワード配列とリバース配列にオーバーラップがあるなら、直接VSEARCHを用いて連結することも可能です。 以下のようにコマンドを実行して下さい。

> vsearch \ --threads 使用するCPU数 \ --fastq_mergepairs ランID__タグID__プライマーID.forward.fastq.gz \ --reverse ランID__タグID__プライマーID.reverse.fastq.gz \ --fastq_allowmergestagger \ --fastqout ランID__タグID__プライマーID.assembled.fastq↓

増幅産物長が解読長よりも短い場合、フォワード配列の読み終わりがリバース配列の読み始めを超えてしまったりリバース配列の読み終わりがフォワード配列の読み始めを超えてしまいますが、デフォルトではそのような配列が連結されません。 --fastq_allowmergestaggerはこれを許して連結するオプションです。 反対側の読み始めを超えてしまったオーバーハング部分は人工配列なので切り捨てられます。 上記のようなケースの可能性がないライブラリを解読した場合は指定する必要はありません。 PEARでは--fastq_allowmergestaggerを指定した場合と同じ処理がなされます。 PEARと同様に連結できなかった配列を得たい場合は--fastqout_notmerged_fwdオプションと--fastqout_notmerged_revオプションに保存先のファイル名を指定して下さい。 そのほか、最小オーバーラップ長を--fastq_minovlenで、連結後の最小・最大配列長を--fastq_minmergelen--fastq_maxmergelenで、オーバーラップ中の最大許容不一致数を--fastq_maxdiffsで、連結後の最大許容期待エラー数を--fastq_maxeeで指定できます。 clconcatpairからVSEARCHを呼び出す場合のデフォルト設定は--fastq_minovlen 20 --fastq_maxdiffs 20としていますが、VSEARCHのデフォルトでは--fastq_minovlen 10 --fastq_maxdiffs 5となっています。

フォルダ内の全てのファイルに対してVSEARCHによる連結を連続して行うには、以下のようなコマンドをお使い下さい。

> for f in `ls *.forward.fastq.gz | grep -P -o '^[^\.]+'`↓ do vsearch \ --threads 使用するCPU数 \ --fastq_mergepairs $f.forward.fastq.gz \ --reverse $f.reverse.fastq.gz \ --fastq_allowmergestagger \ --fastqout $f.assembled.fastq↓ done↓

2.3.4 低品質配列の除去

FASTQ形式の塩基配列には、解読時のデータから推定された塩基の品質値がありますので、その情報に基づいて質の低い配列をこの時点である程度除去してしまうことができます。 以下のようにclfilterseqを実行して下さい。

> clfilterseq \ --minqual=30 \ --maxplowqual=0.1 \ --numthreads=使用するCPU数 \ 入力ファイル \ 出力ファイル↓

このコマンドでは、入力配列において--minqual未満の配列が--maxplowqualの割合より多い場合はその配列を除去します。 通常は出力はファイルになりますが、--output=folderというオプションを与えた場合はフォルダを作成して入力ファイルと同じ名前のファイルに出力します。 既存のフォルダに出力する場合は--appendを追加して下さい。 Illuminaでペアエンドを連結した場合は、両端は品質値が高く、オーバーラップ部分も一致していれば品質値が高くなりますので、末端から品質値の低い配列を削るのはあまり意味がありません。 連結したペアエンドではこのように品質値が低い配列が一定量を超えた配列を除去するのがよいでしょう。 もちろん、既存のFASTQの品質確認を行うプログラムを利用されてもいいと思います。 そのようなプログラムとしては、FastQC (Andrews, 2010)やPRINSEQ (Schmieder and Edwards, 2011)があります。

PEARで連結した全配列ファイルに対して同じ処理を行う場合は、以下のようなコマンドを実行して下さい。

> for f in *.assembled.fastq↓ do clfilterseq \ --output=folder \ --append \ --minqual=30 \ --maxplowqual=0.1 \ --numthreads=使用するCPU数 \ $f \ 出力フォルダ↓ done↓

VSEARCHの機能を用いて、品質値に基づいた配列全体での期待エラー数に上限を設定したフィルタリングも可能です。 これは以下のようなコマンドで適用できます。

> vsearch \ --threads 使用するCPU数 \ --fastq_filter 入力ファイル \ --fastq_maxee 1.0 \ --fastqout 出力ファイル↓

片側だけのシーケンスデータや、連結していないデータの場合は、Roche GSシリーズシーケンサやIon PGMの場合と同様の処理を行えばいいと思います。 第2.2.3節をご参照下さい。 ただし、品質値や長さの閾値は適宜変更する必要があると思います。

2.4 同じ鋳型を複数回PCRまたは複数回シーケンスした場合について

キメラ配列はPCRによって生成されますので、同一PCR由来配列なのか、異なるPCR由来配列なのかをプログラムに正しく認識させなくてはなりません。 そのため、同じ鋳型を複数回PCRまたは複数回シーケンスした場合には、プログラムがデータを正しく扱えるようにいくつかの操作が必要になります。 なお、ここでは、Nested PCRや、増幅のためのPCRとタグ配列付加のためのPCRはまとめて1回と数えることに注意して下さい。 つまり、ここで「複数回PCR」とは、「増幅産物を鋳型にしてさらにPCRで増幅する」ことではなく、「同じ鋳型を使用したPCR」を複数回行うことです。

2.4.1 同じ鋳型を複数回PCRして同じタグ配列を付加して同じランでシーケンスした場合

実験デザインが間違っており、複数回PCRを活かしたキメラ検出・除去は行えません。 通常の1回のPCRと同様にしか扱えません。

2.4.2 同じ鋳型を複数回PCRして別のタグ配列を付加して同じランでシーケンスした場合

ClaidentではサンプルID=ランID__タグID__プライマーIDですから、同じ鋳型由来のサンプルが複数存在するデータとなります。 わざと別サンプルのままにしておいて、1つ以上のサンプルで出現しないOTUはキメラや読み間違い(PCRの合成ミスも含む)のある配列と判断するというのも良い考えだと思います。 ただ、キメラの生成や読み間違いがランダムであればこれはベストな方法でしょうが、できやすいキメラや起きやすい読み間違いというものがもしあれば(あるはずです)、同一のキメラ配列・読み間違いのある配列が全サンプルに出現する可能性があります。 この方法がどの程度正確かという情報はまだ十分ではありませんので、プログラムによるキメラ検出も併用した方が良いかもしれません。 この方法の効果を調べた研究(Lange et al., 2015)も参照して下さい。 この方法はClaidentでサポートしているため実施方法は後述します。 サンプルごとの各OTUの配列数の集計表をこの後で作成しますが、集計表の加工コマンドclfiltersumにより複数のサンプルを統合して1サンプルにすることが可能です(各OTUの配列数は和になります)。

2.4.3 同じ鋳型を1回PCRして同じタグ配列を付加して複数のランでシーケンスした場合

1回のランでは、たまたま得られる配列が少ないことがあります。 そのため同一サンプルを複数のランでシーケンスし、全ランのデータを使いたいことがあります。 ClaidentではサンプルID=ランID__タグID__プライマーIDですから、同じ鋳型由来のサンプルが複数存在するデータとなります。 そのような場合はそのまま別サンプルとして解析を進めることもできますし、1サンプルにまとめてから解析を進めることもできます。 別サンプルのまま処理を進める方を推奨します。 OTUピッキングのためのクラスタリング時に、複数のランの配列をまとめてクラスタリングを行い、1つのまとまった結果を得て下さい。 そして、clsumclassでのサンプルごとの各OTUの配列数の集計表を作成の際に--runnameオプションを用いてRunIDを置換して1サンプルにまとめるか、集計表作成後に集計表の加工コマンドclfiltersum--runnameオプションを指定して複数のサンプルを統合して1サンプルにまとめて下さい(各OTUの配列数は和になります)。

同じ鋳型由来の複数サンプルを1つのサンプルに統合してからこの後の解析に進む場合は、以下のように処理して下さい。

> clsplitseq \ オプション省略 \ --runname=HOGEHOGE \ 入力ファイル1 \ 出力フォルダ↓ > clsplitseq \ オプション省略 \ --runname=HOGEHOGE \ --append \ 入力ファイル2 \ 出力フォルダ↓

このように--runnameオプションでランIDを置換し、2つの入力ファイルのラン名を同一にしてやることで、2つのファイルからの同一サンプルからの配列がHOGEHOGE__タグID__プライマーID.fastq.gzという1つのファイルに書き出されます。

ClaidentではランIDとタグIDとプライマーIDが同じなら、同一PCR産物=同一サンプルの配列として扱います。 PCRが1回でタグ配列ファイルとプライマー配列ファイルが同じなら、当然タグIDとプライマーIDは同じになっているはずですので、ランIDも同じにしてしまうことで、同一PCR産物の配列として認識させることができるようになります。

2.4.4 同じ鋳型を複数回PCRして同じタグ配列を付加して別々のランでシーケンスした場合

複数ラン分の配列ファイルがあるはずですので、clsplitseqによるサンプルごとの配列の分別とclfilterseqによる低品質配列の除去を別々に実行し、別々のフォルダに出力します。 この後のノイジー配列とキメラ配列の除去も別々に行い、クラスタリングの際に複数ラン分のデータをまとめて入力ファイルとして与えて下さい。 そして、サンプルごとの各OTUの配列数の集計表を作成する際に複数のサンプルを統合して1サンプルにまとめて下さい(各OTUの配列数は和になります)。

2.4.5 同じ鋳型を複数回PCRして異なるタグ配列を付加して別々のランでシーケンスした場合

1回のランでは、たまたま得られる配列が少ないことがあります。 そのため同一サンプルを複数のランでシーケンスし、全ランのデータを使いたいことがあります。 ClaidentではサンプルID=ランID__タグID__プライマーIDですから、同じ鋳型由来のサンプルが複数存在するデータとなります。 そのような場合はそのまま別サンプルとして解析を進めることもできますし、1サンプルにまとめてから解析を進めることもできます。 別サンプルのまま処理を進め、ノイジー配列とキメラ配列の除去の際に「同じ鋳型由来の複数のPCR産物に共通に出現しない塩基配列はキメラもしくはノイジー配列とみなす」方法を適用することを推奨します。 そして、サンプルごとの各OTUの配列数の集計表を作成する際に複数のサンプルを統合して1サンプルにまとめて下さい(各OTUの配列数は和になります)。

clsplitseqは、複数の配列ファイルがあるはずですので以下のように実行して下さい。 ここではファイルが2つの場合を示します。

> clsplitseq \ オプション省略 \ --tagfile=タグ配列ファイル1 \ --primerfile=プライマー配列ファイル \ 入力ファイル1 \ 出力フォルダ↓ > clsplitseq \ オプション省略 \ --tagfile=タグ配列ファイル2 \ --primerfile=プライマー配列ファイル \ --append \ 入力ファイル2 \ 出力フォルダ↓

プライマー配列ファイル、出力フォルダは同じものを与えて下さい。 --runnameオプションでは異なるランIDを指定します。 また、タグ配列ファイルは個別に用意して与えて下さい。 ただし、同じ鋳型由来の配列に付加したタグIDは同じにして下さい。 配列が同じタグであっても、異なる鋳型由来の配列に付加したタグのIDを同一にしてはいけません。 例えば以下の2ファイルを用意したとします。

| >sample1 | ACGTACGT | >sample2 | ATGCATGC

| >sample1 | ATGCATGC | >sample2 | ACGTACGT

この場合、sample1の鋳型には1ラン目にはACGTACGTを付加し、2ラン目ではATGCATGCを付加したことになります。 sample2の鋳型では逆になります。 この後の解析に関しては後述します。

第3章 重複配列カウントによるノイジー配列とキメラ配列の除去

Claidentでは、配列のクラスタリング結果から読み間違いの多いノイジーな配列を推定して除去することができます。 方法はCD-HIT-OTU (Li et al., 2012)とほとんど同じです。 また、旧パイプラインではVSEARCHに実装されたUCHIME (Edgar et al., 2011)アルゴリズムを呼び出して結果を受け取ることで、キメラ配列の検出・除去も可能です。 新パイプラインでは、クラスタリング処理が終わってからキメラ検出・除去処理を適用します。

以下のコマンドでノイジー配列の検出・除去が行えます。 なお、多数のファイルを一度に与えることができますが、できるだけ1ラン分のファイルを与えて下さい。 これは、ランごとに塩基配列決定の質にばらつきがあるためです。 1ラン分の配列が非常に多く(1000万超)時間がかかってしまう場合は、1サンプルごとにこの処理を行って下さい。

> clcleanseqv \ --derepmode=PREFIX \ --primarymaxnmismatch=0 \ --secondarymaxnmismatch=1 \ --pnoisycluster=0.5 \ --numthreads=使用するCPU数 \ 入力ファイル1 \ 中略 \ 入力ファイルN \ 出力フォルダ↓

オプションのうち、--derepmodeは前処理において全長一致(FULLLENGTH)で配列をまとめるか、前方一致(PREFIX)で配列をまとめるかを指定します。 オーバーラップのあるペアエンド配列を連結して得られた配列を処理する場合はFULLLENGTHを、シングルエンド配列を処理する場合はPREFIXを指定して下さい。 --primarymaxnmismatchはプライマリクラスタリングで許容するミスマッチ数で通常は0にします。 --secondarymaxnmismatchはセカンダリクラスタリングで許容するミスマッチ数で、通常は1にします。 ただし、「先頭から末尾までに読み間違いが全くない配列」がほとんどないと思われるようなあまりにもノイズの多いデータでは、--primarymaxnmismatchを1に、--secondarymaxnmismatchを3にしてみて下さい。 それでもダメなら--primarymaxnmismatchを2に、--secondarymaxnmismatchを5にして試します。 --secondarymaxnmismatch--primarymaxnmismatchの2倍より1大きい値に設定して下さい。 --pnoisyclusterはノイジー配列と判定する感度に関するオプションです。 0以上1以下で指定して下さい。 大きいほど感度が上がります。 デフォルト値・推奨値は0.5です。 クラスタリングを最終的に97%以下の閾値で行うのであればこれでいいと思いますが、99%前後でクラスタリングする場合、--pnoisyclusterは0.9くらいにするといいかもしれません。 ただし、レアOTUを捨ててしまう可能性が高くなるので、1サンプル当たりの配列数を多めに読んでおく必要があります。

出力フォルダには、以下のファイル群が保存されます。

parameter.txt

プライマリクラスタをノイジーと判定するのに使用した所属配列数の下限値

primarycluster.denoised.fasta.gz

ノイジーと判定されなかったプライマリクラスタの代表配列

primarycluster.fasta.gz

プライマリクラスタの代表配列

secondarycluster.fasta.gz

セカンダリクラスタの代表配列

ランID__タグID__プライマーID.noisyreads.txt.gz

ノイジー配列と判定された配列のリスト

ランID__タグID__プライマーID.singletons.txt.gz

プライマリクラスタリングでシングルトンになった配列のリスト

他にも多くのファイルが出力されますが、後の解析に必要なものもありますので削除しないようにして下さい。

3.1 PCRレプリケート法によるノイジー配列・キメラ配列の除去

PCRレプリケートを用意してノイジー配列・キメラ配列を検出・除去するには、まず同一鋳型DNA由来のPCRレプリケートはどのサンプルとどのサンプルなのかを記したファイルを用意する必要があります。 以下のような内容で用意して下さい。

| sample1  sample2  sample3 | sample4  sample5 | sample6  sample7

同一の行内にタブ区切りで書かれたサンプルが、同一鋳型由来のPCRレプリケートであることを表しています。 異なる鋳型由来のサンプルは別の行に記述して下さい。 PCRレプリケートは3つ以上あっても構いません。 また、レプリケート数が鋳型ごとに異なっても構いません。

上記ファイルが用意できたら、以下のようにclcleanseqvを実行することでレプリケート間で共通しないプライマリクラスタをキメラもしくはノイジーであると判断して除去します。

> clcleanseqv \ --replicatelist=PCRレプリケート指示ファイル \ --derepmode=PREFIX \ --primarymaxnmismatch=0 \ --secondarymaxnmismatch=1 \ --pnoisycluster=0.5 \ --numthreads=使用するCPU数 \ 入力ファイル1 \ 中略 \ 入力ファイルN \ 出力フォルダ↓

レプリケートが3つ以上であっても、全てのレプリケートで検出されたプライマリクラスタのみがノイジーでもキメラでもないと判断されます。 また、1つのプライマリクラスタが複数の鋳型から検出されていた場合、どれか1つの鋳型でノイジーまたはキメラと判定されていれば、他の鋳型でもノイジーまたはキメラとみなして除去します。 これらの設定は以下のオプションで変更することができます。

--minnreplicate

整数値で指定。 この値以上のレプリケートで観測されたプライマリクラスタはそのサンプルではノイジーでもキメラでもないと見なす。 デフォルト値は2

--minpreplicate

小数値で指定。 プライマリクラスタが観測されたレプリケート数/サンプルの総レプリケート数がこの値以上であれば、そのサンプルではノイジーでもキメラでもないと見なす。 デフォルト値は1

--minnpositive

整数値で指定。 この値以上の配列(サンプルではない)がノイジーまたはキメラと判定されていれば、全サンプルでノイジーまたはキメラと見なす。 デフォルト値は1

--minppositive

小数値で指定。 ノイジーまたはキメラと判定された配列数/プライマリクラスタの総配列数がこの値以上であれば、全サンプルでノイジーまたはキメラと見なす。 デフォルト値は0

--runname

ランIDを指定。 サンプル名に含まれるランIDが全てこれに置換される。 同一サンプル名になったサンプルは統合される。

--minnreplicate--minpreplicateの値に基づく1サンプル内での判定と、--minnpositive--minppositiveの値に基づく全サンプルにまたがる判定が2段階で行われていることに注意して下さい。 前者では2つの条件を両方満たすとノイジーでもキメラでもないと判定され、後者では2つの条件を両方満たすとノイジーまたはキメラであると判定されます。 また、同一のプライマリクラスタが、あるサンプルではノイジーまたはキメラだが、別のサンプルではそうでない、といった判定は行えません。 最終的な判定は全サンプル共通です。 また、PCRレプリケート指示ファイルに記されていないサンプルは判定に使われません。

上記のようにclcleanseqvを実行すると、以下のようなファイルも出力フォルダに保存されます。

primarycluster.chimeraremoved.fasta.gz

キメラと判定されなかったプライマリクラスタの代表配列

primarycluster.cleaned.fasta.gz

ノイジーでもキメラでもないと判定されたプライマリクラスタの代表配列

ランID__タグID__プライマーID.chimericreads.txt.gz

キメラ配列と判定された配列のリスト

第4章 塩基配列クラスタリングに基づくOTUピッキング

4.1 サンプル間での配列クラスタリング

塩基配列のクラスタリングを行ってOTUを作成するには、以下のコマンドを実行して下さい。

> clclassseqv \ --minident=0.97 \ --numthreads=使用するCPU数 \ 入力ファイル1 \ 中略 \ 入力ファイルN \ 出力フォルダ↓

入力ファイルには、clcleanseqvの実行時にPCRレプリケート法によるノイジー配列・キメラ配列の除去を併用した場合はprimarycluster.cleaned.fasta.gzを、併用しなかった場合にはprimarycluster.denoised.fasta.gzを与えて下さい。 オーバーラップのないペアエンドデータをclconcatpairを用いて連結した配列では、--paddinglen=16をオプションとして加えて実行して下さい。

出力フォルダには、以下のファイルが保存されています。

clustered.otu.gz

どの生配列がどのOTUに割り当てられたかを記録しているファイル

clustered.fasta

OTUの代表配列

4.2 OTU代表配列への生配列のマッピング

ノイジー配列・キメラ配列と判定されてしまった配列でも、OTUの代表配列に対して指定した類似度の閾値で貼り付けることが可能なのであれば、復活させることができます。 これにより、データの減少を抑えることができます。

以下のコマンドを実行して下さい。

> clrecoverseqv \ --minident=0.97 \ --centroid=OTU代表配列のファイル \ --numthreads=使用するCPU数 \ 入力ファイル1 \ 中略 \ 入力ファイルN \ 出力フォルダ↓

OTU代表配列のファイルにはclclassseqvの出力したclustered.fastaを用います。 入力ファイルとしては、clcleanseqvの出力したprimarycluster.fasta.gzを使用して下さい。 出力フォルダには、clclassseqvと同様のファイル群が作成されます。 なお、オーバーラップのないペアエンドデータをclconcatpairを用いて連結した配列では、--paddinglen=16をオプションとして加えて実行して下さい。

第5章 OTUピッキング結果の要約と後処理

5.1 集計表の作成

以下のコマンドで、クラスタリング結果から各サンプルにおける各OTUに割り当てられた生配列数を表にすることができます。

> clsumclass \ --output=Matrix \ 入力ファイル \ 出力ファイル↓

ここでの入力ファイルは、クラスタリング結果の出力フォルダに保存されている、clustered.otu.gzです。 出力されるファイルは表5.1のようなタブ区切りのテキストファイルで、Excelなどの表計算ソフトやRで読み込むことができます。 ただし、表計算ソフトは巨大な行列の全てを読み込むことができない可能性がありますので注意して下さい。 このファイルに基づいて、群集生態学的な解析を行うことができます。

samplename OTU1 OTU2 OTU3 OTU4 OTU5 OTU6
sampleA 2371 0 0 12 3 0
sampleB 0 1518 0 25 0 1
sampleC 1398 0 0 8 77 6
sampleD 0 1436 0 10 0 0
sampleE 0 0 1360 0 15 3
sampleF 0 0 977 55 6 8
表 5.1: 集計表の例 — 数値は各サンプルにおける各OTUの観測配列数です。

5.2 集計表からの特定のOTU・サンプルの除去

clsumclassでは、1回だけ出現した配列も含めて全て出力されますが、以下のコマンドを用いてサンプルやOTUを選別することができます。 入力ファイルはclsumclassの出力ファイルです。

> clfiltersum \ オプション \ 入力ファイル \ 出力ファイル↓

使用できるオプションは以下の通りです。

--minnseqotu

整数値で指定。 この値に基づいて「割り当てられた生配列数がどのサンプルでも指定値未満のOTU」を除去します。

--minpseqotu

小数値で指定。 この値に基づいて「割り当てられた生配列数/サンプルの生配列総数がどのサンプルでも指定値未満のOTU」を除去します。

--minntotalseqotu

整数値で指定。 この値に基づいて「割り当てられた生配列総数が指定値未満のOTU」を除去します。

--minnseqsample

整数値で指定。 この値に基づいて「どのOTUの生配列数も指定値未満のサンプル」を除去します。

--minpseqsample

小数値で指定。 この値に基づいて「そのサンプルのそのOTUに割り当てられた生配列数/OTUの生配列総数がどのOTUでも指定値未満のサンプル」を除去します。

--minntotalseqsample

整数値で指定。 この値に基づいて「生配列総数が指定値未満のサンプル」を除去します。

--otu

残したいOTU名をカンマで区切って指定します。

--negativeotu

除去したいOTU名をカンマで区切って指定します。

--otulist

ファイル名を指定。 1行に1つのOTU名を記したテキストファイルとして残したいOTU名を指定します。

--negativeotulist

ファイル名を指定。 1行に1つのOTU名を記したテキストファイルとして除去したいOTU名を指定します。

--otuseq

ファイル名を指定。 FASTA形式の配列ファイルとして残したいOTU名を指定します。

--negativeotuseq

ファイル名を指定。 FASTA形式の配列ファイルとして除去したいOTU名を指定します。

--sample

残したいサンプル名をカンマで区切って指定します。

--negativesample

除去したいサンプル名をカンマで区切って指定します。

--samplelist

ファイル名を指定。 1行に1つのサンプル名を記したテキストファイルとして残したいサンプル名を指定します。

--negativesamplelist

ファイル名を指定。 1行に1つのサンプル名を記したテキストファイルとして除去したいサンプル名を指定します。

--replicatelist

ファイル名を指定。 PCRレプリケート関係にあるサンプルをタブ区切りで同一行に記述する。 出力で1つのサンプルに統合される。

--runname

ランIDを指定。 集計表中のサンプル名に含まれるランIDが全てこれに置換される。 同一サンプル名になったサンプルは統合される。

後述するキメラ配列の除去や他領域配列の除去を行った上で、残った配列を用いて--otuseqなどのオプションを用いたOTUの選別を行う場合は、他の選別より先に行って下さい。

5.3 UCHIMEによるキメラ配列除去

クラスタリング後の配列に対してリファレンスなしでのキメラ配列除去とリファレンスありでのキメラ配列除去を追加で行うことができます。 両方行う場合は、リファレンスなしでのキメラ配列除去を先に行ってから、リファレンスありでのキメラ配列除去を行います。 利用するコマンドはclrunuchimeです。 このコマンドでは、UCHIMEアルゴリズムを実装しているVSEARCHを呼び出して、キメラ配列除去を行います。

5.3.1 リファレンスなしでのキメラ配列除去

clrunuchimeを以下のように実行して下さい。

> clrunuchime \ --otufile=*.otu.gz \ 入力ファイル \ 出力フォルダ↓

出力フォルダの内容については後述します。

5.3.2 リファレンスありでのキメラ配列除去

clrunuchimeを以下のように実行して下さい。 リファレンスなしでのキメラ配列除去を行った場合は、nonchimeras.fastaというファイルがあるはずですので、それを入力ファイルにすればいいでしょう。

> clrunuchime \ --referencedb=リファレンスデータベース名 \ 入力ファイル \ 出力フォルダ↓

UCHIME用のリファレンスデータベース名と内容は以下の通りです。

cdu12s

動物12S用のClaidentデータベース(20170513版)

cducox1

動物COX1(COI)用のClaidentデータベース(20170513版)

cducytb

動物Cyt-b用のClaidentデータベース(20170513版)

cdumatk

植物matK用のClaidentデータベース(20170513版)

cdurbcl

植物rbcL用のClaidentデータベース(20170513版)

cdutrnhpsba

植物trnH-psbA用のClaidentデータベース(20170513版)

rdpgoldv9

細菌16S用のRDP Goldデータベース(v9版)

silva128LSUref

LSU rRNA用のSILVAデータベース(Release 128版)

silva128SSUrefnr99

SSU rRNA用のSILVAデータベース(Nr99 Release 128版)

unite20161201

真菌ITS用のUNITEデータベース(20161201版)

unite20161201untrim

真菌ITS用のUNITEデータベース(20161201トリミングなし版)

unite20161201its1

真菌ITS用のUNITEデータベース(20161201ITS1版)

unite20161201its2

真菌ITS用のUNITEデータベース(20161201ITS2版)

5.3.3 出力フォルダの内容について

出力フォルダに保存されるファイルは以下の4つです。

chimeras.fasta

キメラと判定された配列

nonchimeras.fasta

キメラでないと判定された配列

uchimealns.txt

判定時のアライメント

uchimeout.txt

判定時の親配列やスコアなど

uchimeout.txtの各項目の意味は
http://drive5.com/usearch/manual/uchimeout.html
をご覧下さい。

5.4 OTU配列からの低頻度出現OTU配列の除去

クラスタリングの出力フォルダには、clustered.fastaというファイルとして、最終的に得られた代表配列が保存されています。 これを次章のDNAバーコーディングに用いますが、全サンプルを通して少量しか出現しないものまで含めると数が多すぎることがあります。 そこで、以下のようにして5配列以上出現するものだけを抽出することができます。

> clfilterseq \ --otufile=*.otu.gz \ --minnseq=5 \ 入力ファイル \ 出力ファイル↓

5.5 保存領域配列認識による領域分割

ここまでの配列が複数領域にまたがっている場合(たとえばITS1・5.8S rRNA・ITS2といった風に)、領域ごとに分割した方が良いかもしれません。 領域間かその付近に保存的な領域があれば、それを境界の目印として領域を分割できます。 特にユニバーサルプライマー配列などを使うといいでしょう。 これは以下のコマンドで行うことができます。

> cldivseq \ --query=保存領域の配列 \ --border=start \ 入力ファイル \ 出力ファイル1 \ 出力ファイル2↓

このコマンドでは、指定された配列をNeedleman-Wunschアルゴリズムを用いたアライメントによって15%まで不一致を許して探します。 該当する部位が見つかったら、その左端を境界として配列を2つに分け、それぞれ出力ファイル1と出力ファイル2に保存します。 該当する部位が見つからなかった場合は、そのまま出力ファイル1だけに保存されます。 なお、指定する配列は対象配列と同一のストランドにしておいて下さい。 ストランドが異なる場合は、--reversecomplementオプションを付加して実行することで、--queryに指定された配列の逆相補配列をアライメントに用いるようになります。

--borderオプションをstartではなくendにすることで、検索に一致した部位の右端を境界として分割することができます。 --border=bothにした場合は、検索に一致した部位はどちらのファイルにも出力されません(これがデフォルト設定です)。 検索した配列が見つからなかった場合は出力ファイル1のみに配列が保存され、出力ファイル2には配列が保存されませんが、これでは都合が悪い場合は--makedummyオプションを付加することでダミーの塩基配列が出力ファイル2に保存されます。 分割した各領域の配列で別個にホスト生物を同定した上で同定結果を統合する場合など、配列間で対応がとれていないと困る場合にご利用下さい。 ここでは2分割の例を挙げましたが、分割後にさらに分割を繰り返すことで3分割や4分割も可能です。

5.6 ITSxやMetaxaによるITS・SSU rRNA・LSU rRNA配列の抽出

ITSxは核ITS1・ITS2領域を認識してその部位だけを取り出すことができるプログラムです(Bengtsson-Palme et al., 2013)。 多くの分類群でITSは非常に変異に富んでいますが、ずっと保守的なSSU・LSU rRNAと隣接しているため、これらの配列が残っていると、その配列の影響でITSでは明確に区別できる遠縁な生物のSSU・LSU配列がBLAST検索時にヒットしてしまい、うまく同定できなくなってしまうことがあります。 ITSxでITS領域のみを切り出して同定することでこのような問題に対処できます。

MetaxaはITSxと同様にSSU (12S/16S/18S) rRNAとLSU (26S/28S) rRNAを区別してそれぞれの部位だけを取り出すことができるプログラムです(Bengtsson et al., 2011)。 SSUは核由来かミトコンドリア由来か葉緑体由来か細菌由来かも区別できます。 SSU rRNAは真核生物のメタバーコーディング・DNAバーコーディングに広く利用されている領域ですが、核だけでなくミトコンドリアや葉緑体、混入した細菌にも含まれているため、核のSSUを標的とするPCRを行った場合にも多少はミトコンドリアや細菌のSSU rRNAが混入してしまいます。 そのような配列は群集生態学的解析の障害となるため、このプログラムで前もって除去しておくと後の解析が楽になります。

これらのプログラムで作成したFASTA配列ファイルを第5.2節のclfiltersum--otuseqオプションに指定して集計表を加工することで、この配列ファイルに含まれているOTUだけの集計表を作成することができます。

5.7 非ターゲット領域の探索と削除

ITSxやMetaxaはITSとSSU rRNA領域にしか使用できませんが、もっと汎用的な非ターゲット領域の探索方法があります。 一つはBLASTなどで検索して配列がどの遺伝子のものかを判定してくれるプログラムにかけることで、もう一つは多重整列プログラムを利用することです。 ClustalW2・ClustalX2 (Larkin et al., 2007)やMAFFT (Katoh and Standley, 2013)には、多重整列した結果を系統的に近いものが近くになるように並び替えて出力することができます。 この機能を使って並び替えた上で、多重整列の閲覧ソフトで表示して目で確認すればどれが非ターゲット領域かすぐわかります。 非ターゲット領域のOTUを除去したFASTA配列ファイルを第5.2節のclfiltersum--otuseqオプションに指定して集計表を加工することで、この配列ファイルに含まれているOTUだけの集計表を作成することができます。

第6章 DRAへのデータ登録

クラスタリング後の配列は従来通りDDBJなどに登録すればいいですが、新型シーケンサーの生データはDDBJ Sequence Read Archive (DRA)というところへ登録することになっています(サンガー法シーケンサーの生データはTrace Archive)。 本書にあるように複数サンプルからのDNAにタグを付けてシーケンスした場合、サンプルごとのデータファイルに分割して登録することを求められます。 これはclsplitseqによって分割したFASTQを使えば問題ありません。

また、サンプルについての情報(メタデータ)を記述したXMLファイルを作成しなくてはなりません。 データ登録を受け付けている側でも、XML作成を補助するツールが用意されていますが大量サンプルを扱っている場合は用意するのは大変です。 そこで、DRA登録用XML作成を補助するclmaketsvおよびclmakexmlというコマンドを用意しています。 以下でこれらの使用法を説明します。

DRAへの登録には、DDBJが管理するD-wayでのユーザーアカウント作成と公開鍵の登録が必要です。 詳しくはDRA Handbook
http://trace.ddbj.nig.ac.jp/dra/submission.shtml
をご覧下さい。 登録の際、Submission、Study、Experiment、Sample、Runの概念と対応関係を把握している必要がありますので、その部分を特に注意して読んで理解しておいて下さい。

DRAでは、StudyはBioProjectという研究プロジェクトデータベースに登録して、それを参照することになります。 BioProjectへの登録はBioProject Handbook
http://trace.ddbj.nig.ac.jp/bioproject/submission.html
を参照して予め済ませておいて下さい。

さらに、SampleはBioSampleという研究サンプルデータベースに登録して、それを参照します。 BioSampleへの登録はBioSample Handbook
http://trace.ddbj.nig.ac.jp/biosample/submission.html
を参照して予め済ませておいて下さい。 なお、メタゲノムを鋳型としてユニバーサルプライマーで増幅を行ったシーケンスサンプルの場合、MIMarks-SurveyをMIxSのタイプとして指定します。 サンプルの採集された高度、水深、温度、湿度、pHなど、様々な情報を付加しておくことができます。 後世のためにも、面倒がらずにできるだけ多くの情報を付加しておいて下さい。 なお、各項目の意味は、Genomic Standards Consortiumから提供されているチェックリスト
http://wiki.gensc.org/index.php?title=MIMARKS
を確認して下さい。 これを見てもわからない場合はDDBJに問い合わせて下さい。 特に、地下や海底下のサンプルでは、平均海水面からサンプリング地点までの高さ・深さ、地面からサンプリング地点までの高さ・深さ、平均海水面から地面までの高さ・深さを入れる項目があってややこしいのでご注意下さい。

6.1 XML作成用タブ区切りテキストの作成

DRAに登録するXMLには、登録するFASTQファイルごとにそのファイルに含まれている配列の情報を記述します。 これは非常に複雑で大変なため、一旦単純なタブ区切りテキストに書き出した上で、Excelなどを用いてそれを編集し、編集後のタブ区切りテキストからXMLを作成します。 タブ区切りテキストを生成するコマンドはclmaketsvです。 以下のように使って下さい。

> clmaketsv \ 入力ファイル1 \ 中略 \ 入力ファイルN \ 出力ファイル↓

入力ファイルは登録するFASTQを指定して下さい。 ワイルドカードも使えます。 タブ区切りテキストができたら、表計算ソフトなどに読み込ませて、各セルを埋めていって下さい。 なお、Excelでは特定の文字列を勝手に日時などと認識して変換してしまう機能(無効にする方法はありません)がありますので、十分注意して下さい。 各セル内では、[Hogehoge,Fugafuga]とあった場合はHogehogeまたはFugafugaのどちらかを選んで括弧と選ばなかったものを消して下さい。 <Fill in this cell>とあった場合は、<>内の指示に従ってセルを埋めて下さい(括弧は残さない)。 その他のセルは適当に自分で考えて何とかして下さい。 なお、BioProject IDやBioSample IDは、正式なアクセッション番号がまだ発行されていない場合、DDBJに申請したときのSubmission ID (BioProjectではPSUBから始まり、BioSampleではSSUBから始まるもの)を入力しておけば問題ありません。

6.2 タブ区切りテキストからのDRA登録用XMLファイルの生成

タブ区切りテキストの編集が終わったら、タブ区切りテキストとして保存した上で以下のようにclmakexmlを実行して下さい。

> clmakexml \ タブ区切りテキストファイル \ DRAから割り振られたsubmission-ID↓

登録に必要なXMLファイル群が生成されます。 なお、タブ区切りテキストは一度に複数指定することも可能です。 ただし、各ファイルの最初の3行に記述する内容は、最初のファイルのものしか利用されません。 2ファイル目以降の最初の3行は無視されます。

DRAでは、「Create new submission(s)」によってsubmission-IDが割り振られます。 ユーザーID-0001などとなっているはずです。 DRAから割り振られたsubmission-IDにはこれを指定します。 すると、submission-ID.*.xmlという名前のファイルが3つ生成されます。 これをDRAのXML Upload機能を利用して送信します。 その他もろもろの処理が終わると、登録してあるメールアドレスに割り振られたアクセッション番号が送られてきます。 ただし、論文にはBioProject IDを書くことが多いと思います。

第7章 DNAバーコーディングによる配列のホスト生物同定

DNAバーコーディングは、近年非常に多くの分野で応用が進められている、生物の同定方法です。 しかし、既知配列データベースが不十分な上、それを前提とした同定アルゴリズムが欠けていました。 そこで、筆者は「問い合わせ配列と最近隣配列との間の変異量<分類群内の最大変異量」という規準を考案し、これを実現するアルゴリズムQCauto法(Tanabe and Toju, 2013)をClaidentに実装しました。 ただし、その配列のホストの可能性がある全生物が記載済みで、しかもそのバーコード配列もデータベースに登録済みである場合には、そのような難しいことを考えずに最近隣配列と同種とみなせばよいでしょう。 そちらの方法もClaidentには実装してあります。

7.1 BLAST検索による近隣既知配列群の取得

以下のコマンドで、「問い合わせ配列と最近隣配列との間の変異量<分類群内の最大変異量」という条件を満たすために検討すべき近隣既知配列群(のGenBank ID)を取得できます。

> clidentseq \ --blastdb=overall_genus \ --numthreads=使用するCPU数 \ 入力ファイル \ 出力ファイル↓

入力ファイルにはFASTA形式の塩基配列ファイルを指定します。 --blastdbオプションでは、BLAST検索に用いるデータベース名を指定します。 筆者が用意しているデータベースは以下の通りです。

animals_COX1_genus

動物(Metazoa)のmtDNA COX1配列で属以下の情報があるもの

animals_COX1_species

同上だが種以下の情報があるもの

animals_mt_genus

動物(Metazoa)のmtDNA配列で属以下の情報があるもの

animals_mt_species

同上だが種以下の情報があるもの

eukaryota_LSU_genus

真核生物のLSU (28S) rRNA配列で属以下の情報があるもの

eukaryota_LSU_species

同上だが種以下の情報があるもの

eukaryota_SSU_genus

真核生物のSSU (18S) rRNA配列で属以下の情報があるもの

eukaryota_SSU_species

同上だが種以下の情報があるもの

fungi_ITS_genus

真菌のITS配列で属以下の情報があるもの

fungi_ITS_species

同上だが種以下の情報があるもの

overall_class

NCBI ntの中で綱以下の情報があるもの

overall_order

同上だが目以下の情報があるもの

overall_family

同上だが科以下の情報があるもの

overall_genus

同上だが属以下の情報があるもの

overall_species

同上だが種以下の情報があるもの

plants_matK_genus

緑色植物のcpDNA matK配列で属以下の情報があるもの

plants_matK_species

同上だが種以下の情報があるもの

plants_rbcL_genus

緑色植物のcpDNA rbcL配列で属以下の情報があるもの

plants_rbcL_species

同上だが種以下の情報があるもの

plants_trnH-psbA_genus

緑色植物のcpDNA trnHpsbA配列で属以下の情報があるもの

plants_trnH-psbA_species

同上だが種以下の情報があるもの

prokaryota_16S_genus

原核生物の16S rRNA配列で属以下の情報があるもの

prokaryota_16S_species

同上だが種以下の情報があるもの

prokaryota_all_genus

原核生物の配列で属以下の情報があるもの

prokaryota_all_species

同上だが種以下の情報があるもの

semiall_class

overall_classから脊椎動物とCaenorhabditisDrosophilaを除いたもの

semiall_order

overall_orderから脊椎動物とCaenorhabditisDrosophilaを除いたもの

semiall_family

overall_familyから脊椎動物とCaenorhabditisDrosophilaを除いたもの

semiall_genus

overall_genusから脊椎動物とCaenorhabditisDrosophilaを除いたもの

semiall_species

overall_speciesから脊椎動物とCaenorhabditisDrosophilaを除いたもの

overall_*は非常に巨大なため、多くのメモリを必要としますが、これらは万能なのでどんな配列にも使えますし、想定外の分類群や想定外の遺伝子座の配列の混入があってもとんでもない誤同定を避けることができます。 overall_genusでは属以下の情報がある、つまりそれなりに分類群の情報が信頼できそうな配列を選別していますが、あまりにマイナーな分類群だとそれでは近隣配列が見つからないことがあります。 その場合は、綱以下まで同定された配列のデータベースoverall_classを使ってみて下さい。 その他のデータベースは、overall_*よりもメモリ占有量がずっと少ないため、メモリの少ないコンピュータでも動作します。

7.1.1 キャッシュデータベースの構築による高速化

clidentseqはBLAST検索を何度も繰り返しますが、それなりの時間と大量のメモリを必要とします。 そこで、予め問い合わせ配列と類似した配列の上位1万本(設定で変更可)のヒットした部位のみをBLASTデータベースから取得してキャッシュとなるBLASTデータベースを作成しておき、clidentseqではこれを対象として検索することで、劇的な高速化が可能です。 これを行うのがclmakecachedbで、以下のように用います。

> clmakecachedb \ --blastdb=overall_genus \ --numthreads=使用するCPU数 \ 入力ファイル \ 出力フォルダ↓

clidentseqの実行時には、--blastdbオプションに上記コマンドの出力フォルダを指定することで、キャッシュデータベースが使われます。 ただし、clmakecachedbclidentseqの入力ファイルは同一である必要があります。 メモリの使用量も大幅に減少しますので、overall_*を使用してもメモリ16GB程度で済みます。

7.1.2 参照配列データベースが全種を網羅している場合

問い合わせ配列のホストの可能性がある全生物が記載済みで、しかもそのバーコード配列も全て参照配列データベースに登録済みである場合には、以下のようにして類似度99%以上の上位1位(1位タイを含む)を近隣配列とすればいいでしょう。

> clidentseq \ blastn -task megablast -word_size 16 end \ --method=1,99% \ --blastdb=overall_genus \ --numthreads=使用するCPU数 \ 入力ファイル \ 出力ファイル↓

類似度99%以上の上位1位の配列を探し出す場合、BLAST検索オプションは-task megablast -word_size 16としても見逃すことはまずないでしょうから、ここでは高速化のためにこのように指定しています。

7.2 近隣既知配列群に基づく同定

以下のコマンドで、近隣既知配列群の所属分類群が同一になるまで分類階層を上げていくことで配列を同定します。 これはlowest common ancestor (LCA) algorithmと呼ばれています(Huson et al., 2007)

> classigntax \ --taxdb=overall_genus \ 入力ファイル \ 出力ファイル↓

入力ファイルには前節で作成したclidentseqの出力ファイルを指定して下さい。 --taxdbは既知配列の所属分類群データベース指定オプションです。 BLASTデータベースと同じ名前のものがインストールされていますので、それを指定して下さい。

classigntaxのデフォルト設定では、少なくとも2本以上の近隣既知配列がないと同定することができません。 第7.1.2節のようにclidentseqの実行時に--method=1,99%を指定していると、近隣既知配列が1本しかないために一つも同定できないことになります。 この場合は以下のようにして必要な近隣既知配列数を1にします。

> classigntax \ --taxdb=overall_genus \ --minnsupporter=1 \ 入力ファイル \ 出力ファイル↓

出力ファイルは表7.1のような形のタブ区切りのテキストになっています。

query phylum genus species
seqA Ascomycota Chloridium Chloridium virescens
seqB Ascomycota Chloridium Chloridium virescens
seqC Ascomycota Chloridium Chloridium virescens
seqD Basidiomycota Amanita Amanita fuliginea
seqE Basidiomycota Coltriciella Coltriciella dependens
seqF Basidiomycota Filobasidium Filobasidium uniguttulatum
seqG Basidiomycota Laccaria Laccaria bicolor
seqH Basidiomycota Lactarius Lactarius quietus
seqI Basidiomycota Russula Russula densifolia
seqJ Basidiomycota Russula Russula densifolia
seqK Basidiomycota Russula Russula densifolia
seqL Basidiomycota Russula Russula vesca
seqM Basidiomycota Agaricus
seqN Basidiomycota Amanita
seqO Basidiomycota Amanita
seqP Ascomycota Bisporella
seqQ Ascomycota Capronia
seqR Ascomycota Capronia
seqS Ascomycota Cenococcum
表 7.1: 同定結果の例 — 空欄はunidentified、つまり不明ということです。 横幅を抑えるためいくつかの分類階層は省いてあります。

classigntaxのデフォルト設定では、全ての近隣既知配列の所属分類群が同一になるまで分類階層を上げていきますので、誤同定された配列が混ざっていたりした場合には下位の分類階層はほとんど不明になってしまいます。 これは「ほぼ正しいと思われる」結果を得る=「間違えない」ことをデフォルト設定では優先しているためですが、近隣既知配列が多い場合、多少不一致を許容して誤同定の可能性を増やしてでもできるだけ下位の分類階層まで同定したい場合があります。 そのような場合には以下のように--maxpopposer--minsoratioを指定します。

> classigntax \ --taxdb=overall_genus \ --maxpopposer=0.05 \ --minsoratio=19 \ 入力ファイル \ 出力ファイル↓

classigntaxは結果として採用する分類群に該当する配列をsupporter配列とし、該当しない配列をopposer配列とします。 --maxpopposerにはopposer配列の存在を許容する割合を百分率で指定します。 --minsoratioはsupporter配列数とopposer配列数の比の下限値を指定します。 2つのオプションが必要なのは、その分類階層の情報がない配列が存在し得るためです。 そのような配列はsupporterでもopposerでもないので、2つのオプションが必要になります。 上記の例では、近隣既知配列中、5%までのopposer配列を許容し、opposer配列の19倍以上のsupporter配列を必要とします。

7.3 複数の同定結果の統合

例えば植物のrbcLmatKなど、複数のバーコード領域を併用しなければ同定できないこともありますし、overall_genusoverall_classでの同定結果のいいとこ取りをしたい場合があります。 また、厳密なLCAによる同定結果と5%までopposerを許容するLCAの同定結果の組み合わせでいいとこ取りしたい場合もあるでしょう。 そのようなことが、複数の同定結果を統合することで可能になります。

植物のrbcLmatKを併用して同定する場合は、より深くまで同定できている方の結果を採用するのがいいでしょう。 これは以下のコマンドでできます。

> clmergeassign \ --priority=equal \ --preferlower \ rbcLでの同定結果 \ matKでの同定結果 \ 出力ファイル↓

より保守的に考えるなら、以下のようにして、両方の同定結果が同一か、一方では未同定の場合だけ採用することもできます。 ただし、上の階層で不一致だった場合には、同定結果が一方で未同定でも採用されません。

> clmergeassign \ --priority=equal \ rbcLでの同定結果 \ matKでの同定結果 \ 出力ファイル↓

trnHpsbAも併用して、最も深くまで同定できているものを採用する場合は以下のようにします。

> clmergeassign \ --priority=equal \ --preferlower \ rbcLでの同定結果 \ matKでの同定結果 \ trnH-psbAでの同定結果 \ 出力ファイル↓

overall_genusの結果を優先的に採用し、overall_genusで同定できている分類階層までは一致しているがoverall_classの方がより下位まで同定できている場合だけ採用するには、以下のようにします。

> clmergeassign \ --priority=descend \ overall_genusでの同定結果 \ overall_classでの同定結果 \ 出力ファイル↓

同様に、厳密なLCAによる結果を優先し、厳密なLCAで同定できている分類階層まで一致しているが制約を緩めたLCAによる結果の方がより下位まで同定できている場合にのみ採用するには、以下のようにします。

> clmergeassign \ --priority=descend \ 厳密なLCAでの同定結果 \ 制約を緩めたLCAでの同定結果 \ 出力ファイル↓

引用文献

  • Andrews, S., 2010, “Software distributed by the author at http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.”.
  • Bengtsson-Palme, J., Ryberg, M., Hartmann, M., Branco, S., Wang, Z., Godhe, A., De Wit, P., Sánchez-García, M., Ebersberger, I., de Sousa, F., Amend, A., Jumpponen, A., Unterseher, M., Kristiansson, E., Abarenkov, K., Bertrand, Y. J. K., Sanli, K., Eriksson, K. M., Vik, U., Veldre, V., and Nilsson, R. H., 2013, “Improved software detection and extraction of ITS1 and ITS2 from ribosomal ITS sequences of fungi and other eukaryotes for analysis of environmental sequencing data”, Methods in Ecology and Evolution, 4, No. 10, 914–919.
  • Bengtsson, J., Eriksson, K. M., Hartmann, M., Wang, Z., Shenoy, B. D., Grelet, G.-A., Abarenkov, K., Petri, A., Rosenblad, M. A., and Nilsson, R. H., 2011, “Metaxa: a software tool for automated detection and discrimination among ribosomal small subunit (12S/16S/18S) sequences of archaea, bacteria, eukaryotes, mitochondria, and chloroplasts in metagenomes and environmental sequencing datasets.”, Antonie Van Leeuwenhoek, 100, No. 3, 471–475, Oct.
  • Edgar, R. C., Haas, B. J., Clemente, J. C., Quince, C., and Knight, R., 2011, “UCHIME improves sensitivity and speed of chimera detection.”, Bioinformatics, 27, No. 16, 2194–2200, Aug.
  • Fadrosh, D. W., Ma, B., Gajer, P., Sengamalay, N., Ott, S., Brotman, R. M., and Ravel, J., 2014, “An improved dual-indexing approach for multiplexed 16S rRNA gene sequencing on the Illumina MiSeq platform.”, Microbiome, 2, No. 1,  6.
  • Hamady, M., Walker, J. J., Harris, J. K., Gold, N. J., and Knight, R., 2008, “Error-correcting barcoded primers for pyrosequencing hundreds of samples in multiplex.”, Nature Methods, 5, No. 3, 235–237, Mar.
  • Huson, D. H., Auch, A. F., Qi, J., and Schuster, S. C., 2007, “MEGAN analysis of metagenomic data.”, Genome Research, 17, No. 3, 377–386, Mar.
  • Illumina corporation, 2013, “16S metagenomic sequencing library preparation.”.
  • Katoh, K. and Standley, D. M., 2013, “MAFFT multiple sequence alignment software version 7: improvements in performance and usability.”, Molecular Biology and Evolution, 30, No. 4, 772–780, Apr.
  • Kunin, V., Engelbrektson, A., Ochman, H., and Hugenholtz, P., 2010, “Wrinkles in the rare biosphere: pyrosequencing errors can lead to artificial inflation of diversity estimates.”, Environmental Microbiology, 12, No. 1, 118–123, Jan.
  • Lange, A., Jost, S., Heider, D., Bock, C., Budeus, B., Schilling, E., Strittmatter, A., Boenigk, J., and Hoffmann, D., 2015, “AmpliconDuo: A Split-Sample Filtering Protocol for High-Throughput Amplicon Sequencing of Microbial Communities.”, PLoS One, 10, No. 11, e0141590.
  • Larkin, M. A., Blackshields, G., Brown, N. P., Chenna, R., McGettigan, P. A., McWilliam, H., Valentin, F., Wallace, I. M., Wilm, A., Lopez, R., Thompson, J. D., Gibson, T. J., and Higgins, D. G., 2007, “Clustal W and Clustal X version 2.0.”, Bioinformatics, 23, No. 21, 2947–2948, Nov.
  • Li, W., Fu, L., Niu, B., Wu, S., and Wooley, J., 2012, “Ultrafast clustering algorithms for metagenomic sequence analysis.”, Briefings in Bioinformatics, 13, No. 6, 656–668, Nov.
  • Nelson, M. C., Morrison, H. G., Benjamino, J., Grim, S. L., and Graf, J., 2014, “Analysis, optimization and verification of Illumina-generated 16S rRNA gene amplicon surveys.”, PLoS One, 9, No. 4, e94249.
  • Schmieder, R. and Edwards, R., 2011, “Quality control and preprocessing of metagenomic datasets.”, Bioinformatics, 27, No. 6, 863–864, Mar.
  • Stevens, J. L., Jackson, R. L., and Olson, J. B., 2013, “Slowing PCR ramp speed reduces chimera formation from environmental samples.”, Journal of Microbiological Methods, 93, No. 3, 203–205, Jun.
  • Tanabe, A. S. and Toju, H., 2013, “Two new computational methods for universal DNA barcoding: a benchmark using barcode sequences of bacteria, archaea, animals, fungi, and land plants.”, PLoS One, 8, No. 10, e76910.
  • Zhang, J., Kobert, K., Flouri, T., and Stamatakis, A., 2014, “PEAR: a fast and accurate Illumina Paired-End reAd mergeR.”, Bioinformatics, 30, No. 5, 614–620, Mar.

付録 A Assamsを用いた旧パイプラインでの解析

A.1 重複配列カウントによるノイジー配列とキメラ配列の除去

以下のコマンドでノイジー配列とキメラ配列の検出・除去が行えます。 なお、多数のファイルを一度に与えることができますが、できるだけ1ラン分のファイルを与えて下さい。 これは、ランごとに塩基配列決定の質にばらつきがあるためです1ラン分の配列が非常に多く(100万超)時間がかかってしまう場合は、1サンプルごとにこの処理を行って下さい。

> clcleanseq \ uchime --minh 0.1 --mindiv 0.8 end \ --detectmode=N+C \ --pnoisycluster=0.5 \ --numthreads=使用するCPU数 \ 入力ファイル1 \ 中略 \ 入力ファイルN \ 出力フォルダ↓

オプションのうち、uchimeからendの間はUCHIMEの実行オプションです。 --minh--mindivはキメラ配列と判定する感度に関するオプションです。 --minhは値が小さいほど感度が高くなります。 UCHIMEの作者は0.1以上5以下の範囲にするよう推奨しています。 デフォルト値は0.1です。 --mindivの指定値が0.8の場合、0.8%以下の違いしかない配列間はキメラが形成されていても無視されます。 最終的に97%一致規準で配列をクラスタリングするとしたら、2%程度の違いしかない配列間のキメラはクラスタリング・コンセンサス配列作成によって潰されるため、この値を大きくすることで計算が速くなります。 また、--pnoisyclusterはノイジー配列と判定する感度に関するオプションです。 0以上1以下で指定して下さい。 大きいほど感度が上がります。 デフォルト値・推奨値は0.5です。 クラスタリングを最終的に97%以下の閾値で行うのであればこれでいいと思いますが、99%前後でクラスタリングする場合、--pnoisyclusterは0.9くらいにするといいかもしれません。 ただし、レアOTUを捨ててしまう可能性が高くなるので、1サンプル当たりの配列数を多めに読んでおく必要があります。

各入力ファイルは同一の鋳型を同一のチューブ内でPCRした配列群を含んでいる必要があります。 Claidentがキメラ検出に利用しているUCHIMEでは、元配列数が多いコンセンサス配列ほどPCR当初から存在した=キメラではないと仮定し、元配列数が少ないコンセンサス配列が、より元配列数の多いコンセンサス配列の部分配列の組み合わせで説明できるかどうかを検討します。 したがって、UCHIMEに「同一の鋳型を同一のチューブ内でPCRした配列群内でのコンセンサス配列の元配列数」を正確に伝える必要があるため、1ファイルは同一の鋳型を同一のチューブ内でPCRした配列群でなくてはならないのです。 これは1ラン分のファイルを一度に与えることよりもはるかに優先度が高いので、同一の鋳型を同一のチューブ内でPCRしたサンプルを複数回のランで塩基配列決定した場合などは複数ラン分のデータを入力しても構いません。

UCHIMEによるキメラ判定はそれぞれのサンプルごとに行われ、複数のサンプルで観測されている完全一致配列の場合、どれか1つのサンプルでキメラと判定されていれば、他のサンプルでもキメラとみなして除去します。 ただし、下記のオプションでこの動作を変更することができます。 両オプションの両方を満たせばキメラと判定されます。

--minnpositive

整数値で指定。 この値以上の配列(サンプルではない)がキメラと判定されていれば、全サンプルでキメラと見なす。 デフォルト値は1

--minppositive

小数値で指定。 キメラと判定された配列数/完全一致した全配列数がこの値以上であれば、全サンプルでキメラと見なす。 デフォルト値は0

出力フォルダには、以下のファイル群が保存されます。

ランID__タグID__プライマーID.chimeraremoved.dereplicated.fastq.gz

キメラ配列除去し完全一致配列をまとめた配列

ランID__タグID__プライマーID.chimeraremoved.fastq.gz

キメラ配列除去した配列

ランID__タグID__プライマーID.denoised.dereplicated.fastq.gz

ノイジー配列除去し完全一致配列をまとめた配列

ランID__タグID__プライマーID.denoised.fastq.gz

ノイジー配列除去した配列

ランID__タグID__プライマーID.cleaned.dereplicated.fastq.gz

キメラ配列とノイジー配列を除去し完全一致配列をまとめた配列

ランID__タグID__プライマーID.cleaned.fastq.gz

キメラ配列とノイジー配列を除去した配列

ランID__タグID__プライマーID.chimericreads.txt.gz

キメラ配列と判定された配列のリスト

ランID__タグID__プライマーID.noisyreads.txt.gz

ノイジー配列と判定された配列のリスト

ランID__タグID__プライマーID.singletons.txt.gz

完全一致配列をまとめたときにシングルトンになった配列のリスト

ランID__タグID__プライマーID.uchime.txt.gz

UCHIMEによるキメラ判定結果のログ

ランID__タグID__プライマーID.parameter.txt

まとまった完全一致配列をノイジーと判定するのに使用した所属配列数の下限値

他にも多くのファイルが出力されますが、後の解析に必要なものもありますので削除しないようにして下さい。 キメラ配列除去を無効化してノイジー配列除去のみを行う場合は--detectmode=noise、逆にノイジー配列除去を無効化してキメラ配列除去のみを行う場合は--detectmode=chimeraオプションを付けて実行して下さい。 ただし、--detectmode=noiseで処理した結果を--detectmode=chimeraで処理しても、一度にノイジー配列・キメラ配列除去を行った場合(--detectmode=N+Cを指定した場合またはデフォルト設定)と結果は一致しませんのでご注意下さい。 処理する場合は同時に行う必要があります。 これは、キメラ配列もノイジー配列と判定するための情報を持っており、同時処理の場合にはそれを利用するためです。

A.1.1 PCRレプリケート法によるノイジー配列・キメラ配列の除去

PCRレプリケートを用意してノイジー配列・キメラ配列を検出・除去するには、まず同一鋳型DNA由来のPCRレプリケートはどのサンプルとどのサンプルなのかを記したファイルを用意する必要があります。 以下のような内容で用意して下さい。

| sample1  sample2  sample3 | sample4  sample5 | sample6  sample7

同一の行内にタブ区切りで書かれたサンプルが、同一鋳型由来のPCRレプリケートであることを表しています。 異なる鋳型由来のサンプルは別の行に記述して下さい。 PCRレプリケートは3つ以上あっても構いません。 また、レプリケート数が鋳型ごとに異なっても構いません。

上記ファイルが用意できたら、以下のようにclcleanseqを実行することでレプリケート間で共通しないプライマリクラスタをキメラもしくはノイジーであると判断して除去します。

> clcleanseq \ --replicatelist=PCRレプリケート指示ファイル \ uchime --minh 0.1 --mindiv 0.8 end \ --detectmode=N+C \ --pnoisycluster=0.5 \ --numthreads=使用するCPU数 \ 入力ファイル1 \ 中略 \ 入力ファイルN \ 出力フォルダ↓

レプリケートが3つ以上であっても、全てのレプリケートで検出されたプライマリクラスタのみがノイジーでもキメラでもないと判断されます。 また、1つのプライマリクラスタが複数の鋳型から検出されていた場合、どれか1つの鋳型でノイジーまたはキメラと判定されていれば、他の鋳型でもノイジーまたはキメラとみなして除去します。 これらの設定は以下のオプションで変更することができます。

--minnreplicate

整数値で指定。 この値以上のレプリケートで観測されたプライマリクラスタはそのサンプルではノイジーでもキメラでもないと見なす。 デフォルト値は2

--minpreplicate

小数値で指定。 プライマリクラスタが観測されたレプリケート数/サンプルの総レプリケート数がこの値以上であれば、そのサンプルではノイジーでもキメラでもないと見なす。 デフォルト値は1

--minnpositive

整数値で指定。 この値以上の配列(サンプルではない)がノイジーまたはキメラと判定されていれば、全サンプルでノイジーまたはキメラと見なす。 デフォルト値は1

--minppositive

小数値で指定。 ノイジーまたはキメラと判定された配列数/プライマリクラスタの総配列数がこの値以上であれば、全サンプルでノイジーまたはキメラと見なす。 デフォルト値は0

--runname

ランIDを指定。 サンプル名に含まれるランIDが全てこれに置換される。 同一サンプル名になったサンプルは統合される。

--minnreplicate--minpreplicateの値に基づく1サンプル内での判定と、--minnpositive--minppositiveの値に基づく全サンプルにまたがる判定が2段階で行われていることに注意して下さい。 前者では2つの条件を両方満たすとノイジーでもキメラでもないと判定され、後者では2つの条件を両方満たすとノイジーまたはキメラであると判定されます。 また、同一のプライマリクラスタが、あるサンプルではノイジーまたはキメラだが、別のサンプルではそうでない、といった判定は行えません。 最終的な判定は全サンプル共通です。 また、PCRレプリケート指示ファイルに記されていないサンプルは判定に使われません。

A.2 塩基配列クラスタリングに基づくOTUピッキング

A.2.1 サンプル内での配列クラスタリング

以下のコマンドで、同一の鋳型を同一のチューブ内でPCRした配列群のクラスタリングを行います。 前節の結果得られる、完全一致配列をまとめた結果(ランID__タグID__プライマーID.cleaned.dereplicated.fastq.gz)を入力します。

> clclassseq \ --minident=0 \ --strand=plus \ --numthreads=使用するCPU数 \ 入力ファイル1 \ 中略 \ 入力ファイルN \ 出力フォルダ↓

前節の出力ファイルに対してこの処理を適用する場合、以下のようにコマンドを打てばいいでしょう。

> clclassseq \ --minident=0 \ --strand=plus \ --numthreads=使用するCPU数 \ "前節の出力フォルダ/*.cleaned.dereplicated.fastq.gz" \ 出力フォルダ↓

オプションのうち、--minident=0は不一致が1塩基以下の配列をまとめる設定です(実際には、類似度が(最短配列の長さ-11/最短配列の長さ-10)以上の配列をまとめます。 最長配列の長さが400bpで最短配列の長さが350bpなら、0.997となり、400bpの配列でも1塩基しか不一致は許容されませんので意図通りになりますが、最短配列と最長配列の差が大きすぎると意図した通りになりません。 そのため、最短配列は最長配列の半分より11塩基以上長くなくてはなりません)。 最終的に97%で配列をまとめる場合もこの時点では0にしておくのがおすすめです。 最終的に99%でまとめる場合にもここでは0にしておきます。 --strand=plusは同一ストランドの比較だけでクラスタリングを行うオプションです。 片側からしか読んでいないので、両ストランドの比較をする必要はなく、比較するとかえって異常なクラスタリングをしかねません。 両側から読んでいる場合も、連結した後はストランドは揃っているのでこれでいいはずです。 出力フォルダ内の*.assembled.fastq.gzがクラスタリング後の配列です。

A.2.2 サンプル間での配列クラスタリング

サンプル内クラスタリングデータを用いてサンプル間でのクラスタリングを行い、全体のクラスタリング結果を得るには以下のコマンドを実行します。 1サンプルしかない場合も、上記のサンプル内配列クラスタリングを行ってからこの処理を行って下さい。

> clclassclass \ --minident=0 \ --strand=plus \ --numthreads=使用するCPU数 \ 入力ファイル1 \ 中略 \ 入力ファイルN \ 出力フォルダ1↓ > clreclassclass \ --minident=0.99 \ --strand=plus \ --numthreads=使用するCPU数 \ 出力フォルダ1 \ 出力フォルダ2↓ > clreclassclass \ --minident=0.97 \ --strand=plus \ --numthreads=使用するCPU数 \ 出力フォルダ2 \ 出力フォルダ3↓

ここでの入力ファイルは前節の出力ファイル*.assembled.fastq.gzです。 したがって、実際には以下のように入力することになるでしょう。

> clclassclass \ --minident=0 \ --strand=plus \ --numthreads=使用するCPU数 \ "前節の出力フォルダ/*.assembled.fastq.gz" \ 出力フォルダ1↓ > clreclassclass \ --minident=0.99 \ --strand=plus \ --numthreads=使用するCPU数 \ 出力フォルダ1 \ 出力フォルダ2↓ > clreclassclass \ --minident=0.97 \ --strand=plus \ --numthreads=使用するCPU数 \ 出力フォルダ2 \ 出力フォルダ3↓

1つ目のclclassclassでサンプル内配列クラスタリングでまとめられた配列において、不一致が1塩基以下の配列をまとめます。 2つ目のclreclassclassで、99%以上一致する配列をまとめ直します。 さらに3つ目のclreclassclassで、97%以上一致する配列をまとめ直します。 こうすることで、より高速に、より正確にクラスタリングが行われます。 clreclassclassを使用せずにclclassclassに97%でまとめるよう指示してしまうと、アルゴリズム上「まとめすぎ」が発生しやすくなります。

出力フォルダには、以下のファイルが保存されています。

assembled.contigmembers.gz

どの生配列がどのOTUに割り当てられたかを記録しているファイル

assembled.fastq.gz

OTUのコンセンサス配列

assembled.fasta

OTUのコンセンサス配列

A.2.3 集計表の作成

新パイプラインではclustered.otu.gzを使用しますが、旧パイプラインでは生成されません。代わりにassembled.contigmembers.gzを使用して下さい。

付録 B その他のプログラムのインストール

B.1 bcl2fastqのインストール

Illuminaから提供されている、BCL形式のベースコールデータからFASTQを生成するプログラムbcl2fastqはRTA version 1.18.54より前の機種用のv1.8.4およびRTA version 1.18.54以降の機種用のv2.18.0.12が
http://support.illumina.com/downloads/bcl2fastq_conversion_software.html
からダウンロードできます。 v1.8.4は以下のようにインストールできます。

# alienのインストール $ sudo apt-get install alien # bcl2fastq v1.8.4のダウンロード $ wget \ ftp://webdata:webdata@ussd-ftp.illumina.com/Downloads/Software/bcl2fastq/bcl2fastq-1.8.4-Linux-x86_64.rpm # .rpmから.debへの変換とインストール $ sudo alien -i \ bcl2fastq-1.8.4-Linux-x86_64.rpm # bcl2fastqに必要な追加パッケージのインストール $ sudo apt-get install libxml-simple-perl xsltproc

Ubuntu・Xubuntu・Lubuntuの14.04LTSでは、Perlのバージョンが新しすぎ、より古いバージョンを前提として書かれたbcl2fastqが動作しません。 そこで以下のようにして問題の部分を書き換えます。

$ sudo perl -i -npe \ 's/qw\(ELAND_FASTQ_FILES_PER_PROCESS\)/\("ELAND_FASTQ_FILES_PER_PROCESS"\)/' \ /usr/local/lib/bcl2fastq-1.8.4/perl/Casava/Alignment/Config.pm $ sudo perl -i -npe \ 's/qw\(ELAND_GENOME\)/\("ELAND_GENOME"\)/' \ /usr/local/lib/bcl2fastq-1.8.4/perl/Casava/Alignment/Config.pm

テキストエディタで
/usr/local/lib/bcl2fastq-1.8.4/perl/Casava/Alignment/Config.pm
の747行目と751行目を手動で書き換えても構いません。 qw(HOGEHOGE)("HOGEHOGE")に書き換えて下さい。

v2.18.0.12のインストールをする場合は以下のようにコマンドを実行して下さい。

# rpm2cpioとcpioのインストール > sudo apt-get install rpm2cpio cpio↓ # bcl2fastq2 v2.18.0.12のダウンロード > wget -O bcl2fastq2-v2-18-0-12-linux-x86-64.zip http://bit.ly/2drPH3W↓ # 圧縮ファイルの展開 > unzip -qq bcl2fastq2-v2-18-0-12-linux-x86-64.zip↓ # .rpmを展開してコマンドを抽出 > rpm2cpio \ bcl2fastq2-v2.18.0.12-Linux-x86_64.rpm | cpio -id↓ # コマンドをインストール > sudo mv usr/local/bin/bcl2fastq /usr/local/bin/↓ # その他のファイルをインストール > sudo cp -R usr/local/share/css /usr/local/share/↓ > sudo cp -R usr/local/share/xsl /usr/local/share/↓

付録 C ターミナルコマンド集

C.1 配列を数え上げる

# FASTQの場合 > grep -P -c '^\+\r?\n?$' inputfile↓ # GZIP圧縮FASTQの場合 > gzip -dc inputfile | grep -P -c '^\+\r?\n?$'↓ # BZIP2圧縮FASTQの場合 > bzip2 -dc inputfile | grep -P -c '^\+\r?\n?$'↓ # XZ圧縮FASTQの場合 > xz -dc inputfile | grep -P -c '^\+\r?\n?$'↓ # FASTAの場合 > grep -P -c '^>' inputfile↓ # GZIP圧縮FASTAの場合 > gzip -dc inputfile | grep -P -c '^>'↓ # BZIP2圧縮FASTAの場合 > bzip2 -dc inputfile | grep -P -c '^>'↓ # XZ圧縮FASTAの場合 > xz -dc inputfile | grep -P -c '^>'↓

C.2 配列を閲覧する

# 無圧縮ファイル内容を画面に出力 > cat inputfile↓ # 無圧縮ファイル内容をスクロールしながら見る > less inputfile↓ # GZIP圧縮ファイル内容を画面に出力 > gzip -dc inputfile↓ # GZIP圧縮ファイル内容をスクロールしながら見る > gzip -dc inputfile | less↓ # FASTQの最初の1配列を画面に出力 > head -n 4 inputfile↓ # GZIP圧縮FASTQの最初の1配列を画面に出力 > gzip -dc inputfile | head -n 4↓ # FASTQの特定の配列を検索して画面に出力 > grep -P -A 3 '^\@配列名の正規表現' inputfile↓ # GZIP圧縮FASTQの特定の配列を検索して画面に出力 > gzip -dc inputfile | grep -P -A 3 '^\@配列名の正規表現'↓

C.3 圧縮・展開

# フォルダ内の全.fastq.gzを展開 > for f in *.fastq.gz↓ do gzip -d $f↓ done↓ # フォルダ内の全.fastq.gzを展開(サブフォルダ含む) > for f in `find . -name *.fastq.gz`↓ do gzip -d $f↓ done↓ # 4つのCPUを使用してフォルダ内の全.fastq.gzを展開 > ls *.fastq.gz | xargs -L 1 -P 4 gzip -d↓ # 4つのCPUを使用してフォルダ内の全.fastq.gzを展開(サブフォルダ含む) > find . -name *.fastq.gz | xargs -L 1 -P 4 gzip -d↓ # フォルダ内の全.fastqを圧縮 > for f in *.fastq↓ do gzip $f↓ done↓ # フォルダ内の全.fastqを圧縮(サブフォルダ含む) > for f in `find . -name *.fastq`↓ do gzip $f↓ done↓ # 4つのCPUを使用してフォルダ内の全.fastqを圧縮 > ls *.fastq | xargs -L 1 -P 4 gzip↓ # 4つのCPUを使用してフォルダ内の全.fastqを圧縮(サブフォルダ含む) > find . -name *.fastq | xargs -L 1 -P 4 gzip↓

C.4 抽出・保存

# FASTQから最初の1万配列を抽出 > head -n 40000 inputfile > outputfile↓ # FASTQから最後の1万配列を抽出 > tail -n 40000 inputfile > outputfile↓ # GZIP圧縮FASTQから最初の1万配列を抽出してGZIP圧縮FASTQへ保存 > gzip -dc inputfile | head -n 40000 | gzip -c > outputfile↓ # GZIP圧縮FASTQから最後の1万配列を抽出してGZIP圧縮FASTQへ保存 > gzip -dc inputfile | tail -n 40000 | gzip -c > outputfile↓ # GZIP圧縮FASTQの特定の配列を検索してGZIP圧縮FASTQへ保存 > gzip -dc inputfile | grep -P -A 3 '^\@配列名の正規表現' | gzip -c > outputfile↓