2017年 NGSハンズオン講習会
8月31日

メタゲノム解析森 宙史 (Hiroshi Mori), Ph.D.
国立遺伝学研究所
生命情報研究センター
hmori@nig.ac.jp1

数%ぐらいの菌しか培養できない

2

細菌群集を解析するための様々な実験手法
細菌群集を解析するために使用されてきた実験手法
培養による細菌コロニーのカウント法 ・・・ 培養困難な細菌は解析出来ない染色による細菌の数のカウント法 ・・・ 細菌の数しかわからない
FISH法による特定の細菌の染色法 ・・・ プローブ配列を設計する必要がある
DGGE法による細菌群集の解析法 ・・・ バンドパターンのみであり、細菌群集の全体像をとらえるのは困難
細菌群集が形成するシステムを詳細に解析するためには、これらの手法では断片的な情報しか得られないため、細菌群集についての理解はあまり進んでいなかった

3

16S ribosomal RNA (16S rRNA)・ リボソームの核となるRNAの一つ
・ 全ての細菌が所持

・ 配列間の結合によって高次構造を形成
・ 系統マーカー遺伝子の代表例

・ 100万本以上の配列がデータベースに登録済み
・ 全長1500 base

16S rRNA遺伝子は広範囲の細菌における系統推定を行う上で最適な遺伝子

4

Woese, C. et al. 1990

16S rRNA gene amplicon sequencing analysis(メタ16S解析)
Pre-analysis (Remove Primer, Chimera etc.)
Sequence clustering with species levelby CD-HIT-EST or UCLUST, etc.

DNA extraction

PCR amplification

Taxonomic assignment andComparison between samples

DNA Sequencing

Togo picture gallery by DBCLS islicensed under a Creative Commons
Attribution 2.1 Japan license (c)Who’s there?

5Metagenomic sequencing analysis(メタゲノム解析)
Metagenomic reads

Assemble

MEGAHIT,MetaPlatanus
IDBA-UD etc.

Contig sets

DNA extraction

Gene finding

MGA, MetaGeneMark

Gene sets

DNA Sequencing

Togo picture gallery by DBCLS islicensed under a Creative Commons
Attribution 2.1 Japan license (c)BLASTP

Sample1
Metadata
Sample2
Metadata

Comparative
metagenomics

Who’s there?

Gene Function abundance
Pathway abundance
Pathway reconstruction
Taxonomic abundance
What are they doing?6

メタ16S解析
利点
・安価かつ少量のDNAから系統組成が得られる
・reference配列に依存しない解析も可能
・マシンパワーは少なくて済み、解析ツールも普及（QIIME・mothur等）

欠点
・PCRバイアスの存在
・種以下は分解能に問題あり
・個々の系統の機能が不明

メタゲノム解析
利点
・系統組成と遺伝子機能組成が得られる
・実験によるバイアスが少ない
・優占系統のドラフトゲノムの構築（条件が良ければ可能）

欠点
・reference配列に依存した解析・目的依存で解析手法が変化し、マシンパワーも必要
7

16S rRNA遺伝子配列を用いた
既知Bacteria Phylumの系統樹

Genome sequence available
No genome sequence available

816S rRNA gene amplicon sequencing analysis(メタ16S解析)
Pre-analysis (Remove Primer, Chimera etc.)
Sequence clustering with species levelby CD-HIT-EST or UCLUST, etc.

DNA extraction

PCR amplification

Taxonomic assignment andComparison between samples

DNA Sequencing

Togo picture gallery by DBCLS islicensed under a Creative Commons
Attribution 2.1 Japan license (c)Who’s there?

9PCRに使うプライマーは大きく分けて2種類• 系統特異的なプライマー• 系統Universalなプライマー

両者は何が違うのか？

10

系統(機能)特異的プライマー
例えば、• 腸管出血性大腸菌とそうではない大腸菌を判別したい• Fusarium oxysporumのレースを判別したいなどの病原性の有無の判定を迅速に行いたい場合に使われたりする。
遺伝子レベルで病原性のメカニズムがわかっている生物の場
合、• ある遺伝子を持っているか否か？
• ある遺伝子に一塩基置換があるか無いか？
が、病原性の有無に重要であるような場合には、非常に有効なプライマーが設計可能な場合が多い。
11

系統特異的プライマーの特徴• 系統判定には増幅産物のシーケンスをしなくても大丈夫（PCR後の電気泳動でバンド出るか否か）• どの遺伝子を使うかはバラエティに富む• PCR条件を厳密に検討する必要がある（非特異的増幅の回避が重要）• degenerate primerが少ない
degenerate primerとは？

12Degenerate primer
例: 5′-GTGCCAGCMGCCGCGGTAA-3′• 曖昧（縮重）塩基を使ったプライマー
曖昧塩基塩基1塩基2塩基3R

A

GYCT

S

GCW

ATKGTM

ACB

CGT

DA

GTH

A

CTV

A

CGN

A

C

G

塩基4T

13

Universalプライマー• 幅広い系統群を増幅できるプライマー
用途• 系統の判別• 群集組成を見る
特徴• 増幅産物の系統判別には基本的にシーケンシングが必要• degenerate primerが多い• ターゲットになりうる遺伝子は少数
14

菌類(糸状菌・酵母など）• 28S rRNA遺伝子• 18S rRNA遺伝子• COX1• COX2• rDNA–ITS1 (internal transcribed spacer), rDNA–ITS2
18S rDNA – ITS1 – 5.8S rDNA – ITS2 – 28S rDNA

15

Virus (RNA virusの場合)• Tospovirus Nタンパク質• Comoviridae RNA-dependent RNA polymerase• Tombusviridae RNA-dependent RNA polymerase• Flexiviridae RNA-dependent RNA polymerase、外被タンパク質（Coat protein)• タバコモザイクウイルスなどの棒状ウイルス 外被タンパク質16

細菌• 16S rRNA遺伝子• rDNA–ITS
真核生物: 18S rDNA – ITS1 – 5.8S rDNA – ITS2 – 28S rDNA
原核生物: 16S rDNA – ITS – 23S rDNA – 5S rDNA

17

16S rRNA gene amplicon sequencing analysis(メタ16S解析)
Pre-analysis (Remove Primer, Chimera etc.)
Sequence clustering with species levelby CD-HIT-EST or UCLUST, etc.

DNA extraction

PCR amplification

Taxonomic assignment andComparison between samples

DNA Sequencing

Togo picture gallery by DBCLS islicensed under a Creative Commons
Attribution 2.1 Japan license (c)Who’s there?

18

メタ16S解析パイプラインの例: Ion PGMリードの前処理
Ion PGMリード

fastqファイル

Nを含むリードを除外

長さフィルタ (増幅断片長の中央値±約50 base以内のリードのみ採用)
FASTQからFASTAとQUALITYファイルへの変換
クオリティフィルタ (Average QV >= 25)

TagCleanerによる5’ Primer除去

TagCleanerによる
3’ Primer除去

前処理済みIon PGM
リード19

キメラ除去前処理済みIon PGM
リード

UCLUSTを用いて全サンプルのリードをクラスタリングしOTU化 (Identity 97%,coverage 80%)
OTU代表配列

Reference 16S rRNA
gene database

UCHIME Reference mode
でキメラを検出

UCHIME De novo mode
でキメラを検出

Broad Instituteが提供している
Type Strains + complete or draft genomesの5181本の高精度16Sデータ

両modeでキメラとされた
OTUをキメラと判定、そのOTUを構成する全リードを除去

キメラ除去済みOTU各OTUのサンプルごとの構成本数をまとめたOTU
Table20

系統アサインメント
キメラ除去済みOTURDP Classifierを用いてbootstrap cutoff >= 50%
で各OTUの代表配列を系統アサインメント

各OTUのサンプルごとのリード数を整理したOTU table

Domainレベルの系統組成
Phylumレベルの系統組成

Classレベルの系統組成
Orderレベルの系統組成

Familyレベルの系統組成
Genusレベルの系統組成

各OTUの系統情報とサンプルごとのリード数を整理したOTU table21http://qiime.org/

QIIME citations since publication
Web of Science: ~3300
Google Scholar: ~5050

22https://docs.qiime2.org/2017.7/plugins/available/

23https://docs.qiime2.org/2017.7/tutorials/

24https://mothur.org/wiki/Main_Page25http://www.drive5.com/usearch/

http://www.drive5.com/usearch/manual/pipe_examples.html

26

論文のMethodsの記述としてダメな例• “sequence data were analyzed by QIIME”• “sequence data were analyzed by mothur”• “sequence data were analyzed by USEARCH”としか書かないのはダメ
それらはコマンドの集合体であるため、必ず、どのコマンドを使ってどのようなパラメータ設定で、何をしたのかまで書く。
27

Matsuki T.†, Yahagi K.†,
28Mori H.†,et al, Nature Commun., 2016

http://metagenomics.anl.gov/29

30

Metagenomic sequencing analysis(メタゲノム解析)
Metagenomic reads

Assemble

MEGAHIT,
MetaPlatanus
IDBA-UD etc.

Contig sets

DNA extraction

Gene finding

MGA, MetaGeneMark

Gene sets

DNA Sequencing

Togo picture gallery by DBCLS islicensed under a Creative Commons
Attribution 2.1 Japan license (c)BLASTP

Sample1
Metadata
Sample2
Metadata

Comparative
metagenomics

Who’s there?

Gene Function abundance
Pathway abundancePathway reconstruction
Taxonomic abundance
What are they doing?
31

(山本・森・山田・黒川, 2014 「生命のビッグデータ利用の最前線」 シーエムシー出版より一部改変)

細菌群集
DNA抽出

シーケンス
ライブラリ作成

シーケンシング

塩基配列データのクオリティコントロール

メタデータ記録

メタデータリスト

塩基配列データContig /Scaffold配列

メタゲノム
アセンブル

遺伝子予測

アセンブル結果への配列マッピング

遺伝子のアミノ酸配列データ

機能情報付きアミノ酸配列DB
配列相同性検索(例:KEGG,GenBank nr等) 各遺伝子の機能と群集中の存在量

高精度
塩基配列データ
配列相同性検索

各系統の群集中の存在量

統計学の手法を用いた比較メタゲノム解析

その細菌群集の特異性や他の細菌群集との類似性等の知見

16S rRNA
遺伝子配列DB(例:SILVA, RDP等)

既知のメタゲノム
解析データのDB(例:MG-RAST,
MicrobeDB.jp等)
32

Next Generation Sequencer (NGS) (データはちょっと古いです）Sequencer Name
Specific
Read length Read number /property(base)run
ABI 3730xl

Sanger500–1000

384

Ion Proton

Emulsion PCR

200

80,000,000

MiSeq

Bridge PCR

300

30,000,000

454 GS FLX+

Emulsion PCR

700-1,000

1,000,000

HiSeq 2500/4000

Bridge PCR

150 or 250~3,200,000,000

PacBio RS II /SequelSingle molecule>7,000>50,000

MinION

Single molecule>数kb>10,00033

目的に応じたハードウェアのスペックの目安
原核・真核・メタゲノム・Transcriptomeアセンブル
特に重要なハードウェア: メモリ・CPU
・ 原核生物のゲノムアセンブル 数十GBメモリ (50GB あれば十分)、1CPU 数コア
・ Transcriptomeアセンブル 今のHiSeqのリード1億pairsなら、120GBあれば十分(例えばTrinityなら、100万pairsにつき1GBメモリが目安)(https://github.com/trinityrnaseq/trinityrnaseq/wiki/)・ メタゲノムのアセンブル 組成が単純な群集なら原核生物のアセンブルと同様群集が多様な場合、大量のリードが必要であり、数百GB-数TBメモリが必要な場合もある。 十数コア以上
・ 真核生物のゲノムアセンブル ゲノムサイズ、どれくらいheteroか、等に依存する
数十GB-数TBメモリ、十数コア以上

少なくとも100GB以上のメモリをのせないと、アセンブルは辛い34

目的に応じたハードウェアのスペックの目安
数千本以上の配列相同性検索
特に重要なハードウェア: CPU・ ゲノムやメタゲノムの遺伝子アノテーション
・ 複数ゲノムの比較解析

入力配列が数千本以上になることが多く、計算を並列化して高速化する必要がある。

現在のnr相手のBLASTXやBLASTPは、並列化しても各プロセスで5GBほどメモリを使用するので、メモリもそれなりに必要。

CPU数、コア数が非常に重要、メモリもある程度必要
35

目的に応じたハードウェアのスペックの目安
マッピングツールでのReferenceゲノムへのマッピング
特に重要なハードウェア: ディスク・ Resequencing, RNA-Seq, ChIP-Seq解析におけるマッピング
リードのゲノムへのマッピングは高速でメモリ使用量も少ない。
結果のSAMファイルやBAMファイルが数十から数百GBになったりする。
samtools等でsortしたりすると、その規模のファイルが何個もできる。

マッピングを頻繁にするのなら、ディスクは少なくとも十数TBは必要

36

目的に応じたハードウェアのスペックの目安
1. 原核・真核・メタゲノム・Transcriptomeアセンブル
特に重要なハードウェア: メモリ・CPU

2. 数千本以上の配列相同性検索
特に重要なハードウェア: CPU

3. マッピングツールでのReferenceゲノムへのマッピング
特に重要なハードウェア: ディスク37

午後からの実習で用いるメタゲノムデータ

Backhed F. et al 2015の4 daysの新生児の1サンプルを100万 pairにダウンサンプリングしたデータ

38

(山本・森・山田・黒川, 2014 「生命のビッグデータ利用の最前線」 シーエムシー出版より一部改変)

ヒト乳児腸内細菌群集
DNA抽出

シーケンス
ライブラリ作成

シーケンシング

塩基配列データ

100 万pairsIllumina HiSeq 2000
メタデータ記録

MEGAHIT
Contig配列

メタデータリスト

メタゲノム
アセンブル

Prodigal 遺伝子予測

アセンブル結果への配列マッピング

遺伝子のアミノ酸配列データ

Bowtie 2

機能情報付きBLASTP
アミノ酸配列DB
配列相同性検索(例:KEGG,UniProt KB等) 各遺伝子の機能と群集中の存在量

塩基配列データのクオリティコントロール
Trimmomatic
高精度
塩基配列データ
配列相同性検索

各系統の群集中の存在量

統計学の手法を用いた比較メタゲノム解析

その細菌群集の特異性や他の細菌群集との類似性等の知見

BLASTN

16S rRNA
遺伝子配列DB(例:SILVA, RDP等)

既知のメタゲノム
解析データのDB(例:MG-RAST,MicrobeDB.jp等)
39

今日の実習• ほぼ全てターミナル上で行います• cd• cd Test• TestがWorkingディレクトリです• Test/には、Data/とProgram/があります• コマンドは全てCommandMemo3.txtに書いてありますので、コピペで良いです• マシンのメモリが8GB以下の場合は、Bio-Linuxは4GBで起動しましょう40

今日の実習の参考資料
https://2017-ucsc-metagenomics.readthedocs.io/en/latest/

41

PhiXについて
https://jp.illumina.com/content/dam/illuminamarketing/apac/japan/documents/pdf/2013_ill
umina_techsupport_session16.pdf

42塩基配列の類似性検索ツールにも様々なものが存在する

Referenceと近ければMapper系 (Bowtie 2等)
Referenceと遠いのなら、BLAST系(Chaisson J. M. et al. 2012)43ゲノムアセンブルの二大戦略• Overlap-Layout-Consensus– k-merの共有やローカルアラインメント等でリード間のoverlapを見つけて、短いContigを作成し、さらにContig間をoverlapをもとに結合（layout）。リードのoverlapを領域ごとに集めてマルチプルアラインメント等をしてconsensusをとることでアセンブルする– 例: Celera Assembler, Newbler, Mira, Canu

• de Bruijn Graph– リードをoverlapありのk-merに分割して、多数のリード間のk-merの共有をde Bruijn graphというグラフ構造で表現して、グラフ上で最短経路を見つける問題を解く
44

メタゲノムアセンブルツールの例• IDBA-UD (Peng et al. 2012)– 短いk-merでアセンブルしてContig作成（Contig間のcoverageの差はある程度許容する）。そのContig群を用いて、もう少し長めのk-merでアセンブルしてContig作成。これを繰り返す。最後に、Contig間をまたがるpairリード（paired-endやmate pair)の情報をもとに、scaffoldingする。
– 短いk-merでシーケンスエラー、長いk-merでリピートの問題に対処• MEGAHIT (Li et al. 2015)– Contig作成の方法はIDBA-UDと類似しているが、de Bruijn graphの表現方法が簡素化されているため (http://alexbowe.com/succinctdebruijn-graphs/) 、高速で省メモリ。また、coverageが小さいk-merの扱いについて色々と工夫している。scaffoldingはしない。

• metaSPAdes (Nurk et al. 2017)– Contig作成の方法はIDBA-UDと類似しているが、リードデータ中のstrainレベルの配列多様性をContig/Scaffoldにおいてもできるだけ保つために、サイトに多型があるとContigを分岐する傾向が強い。
45

今日のProkka(https://github.com/tseemann/prokka/blob/master/README.md)• MetaProdigal: CDS prediction• Aragorn: tRNA
• CDSについては、機能アノテーションをするために、次の2 stepを行う。
– BLASTPで、UniProt KB中のEvidence codeがreal protein ortranscript evidenceとなっている、Prokaryote由来のproteinのアミノ酸配列に対して検索– HMMER3でPfamとTIGRFAMに対して検索• 上記2 stepで閾値以上(BLASTならE-value <1e-6) でHitしなければ、hypothetical proteinとする。
46

メタゲノムデータからの系統組成推定
16S rRNA遺伝子を用いる手法
利点: Reference配列DBが充実（RDP, SILVA, GreenGenesなど)
欠点: ゲノム内コピー数の問題

系統マーカータンパク質の遺伝子(rpoB等)を用いる手法
利点: 1ゲノム1コピー
欠点: Reference配列DBが貧弱

リードやcontigのマッピング or k-mer組成を用いる手法
利点: Virus, 原核, 真核生物を同時解析可能
欠点
・Reference配列DBが貧弱
・Genus以下の精度がかなり落ちる
・水平伝播配列の扱い

47

パスウェイデータベース• 先週行われた、AJACS河内の九州大の山西先生の資料が参考になります。
http://motdb.dbcls.jp/?AJACS66
http://motdb.dbcls.jp/?plugin=attach&pcmd=op
en&file=170824Pathway_DB_yamanishi_submit_r.pdf&refer=AJACS66• メタゲノムでは、KEGGのKEGG Orthologyを遺伝子機能の単位として使うことが多い48

(山本・森・山田・黒川, 2014 「生命のビッグデータ利用の最前線」 シーエムシー出版より一部改変)

ヒト乳児腸内細菌群集
DNA抽出

シーケンス
ライブラリ作成

シーケンシング

塩基配列データ

100 万pairs

Illumina HiSeq 2000
メタデータ記録

MEGAHIT
Contig配列

メタデータリスト

メタゲノム
アセンブル

MetaProdigal 遺伝子予測

アセンブル結果への配列マッピング

遺伝子のアミノ酸配列データ

Bowtie 2

機能情報付きBLASTP
アミノ酸配列DB
配列相同性検索(例:KEGG,
UniProt KB等) 各遺伝子の機能と群集中の存在量

塩基配列データのクオリティコントロール
Trimmomatic
高精度
塩基配列データ
配列相同性検索各系統の群集中の存在量

統計学の手法を用いた比較メタゲノム解析

その細菌群集の特異性や他の細菌群集との類似性等の知見

BLASTN

16S rRNA
遺伝子配列DB(例:SILVA, RDP等)

既知のメタゲノム
解析データのDB(例:MG-RAST,MicrobeDB.jp等)
49代表的なメタゲノムデータベース
運営者

NCBI Taxonomy +SRA

配列
データ

メタゲノム/ 系統 遺伝子
メタ16S区
組成 機能組別成

サンプル数(2017 8月)

ゲノム等との統合化

NCBI, USA ○×××

673,079×

JGI, USA××××

24,922×

JGI, USA○（Reads+Contigs)○○（メタゲノムのみ収録）

○

7,982

△

Chicago U. ○USA○○○47,313(303,594)

××

○○

○

74,342×https://www.ncbi.nlm.nih.gov/Taxonomy/GOLDhttps://gold.jgi.doe.gov/IMG/M
https://img.jgi.doe.gov/cgi-bin/m/main.cgi

MG-RAST
http://metagenomics.anl.gov/index.html

EBI-Metagenomics

EBI, EU

https://www.ebi.ac.uk/
metagenomics/

MicrobeDB.jp

NIG,×

○

○

○

60,55150○

微生物統合DB「MicrobeDB.jp」
http://microbedb.jp/• 微生物に関するデータを系統・遺伝子・環境の３つの軸に沿って整理・統合し、フルRDFのDBを構築• 約90億トリプルから構成• 12種類のオントロジー＆ボキャブラリの開発• 公開済みの約6万サンプルのメタ16S・メタゲノムデータ、約1万7千株のゲノム・ドラフトゲノムデータを収録• 195種類のStanzaの開発• 解析プロトコルの標準化および解析パイプラインの開発• 単細胞の真菌類・藻類のゲノムデー
タも整理・統合• 自動更新技術の開発
51

52

MG-RAST v.3 Pipeline
ftp://ftp.metagenomics.anl.gov/data/manual/mg-rast-manual.pdf

53Preprocessing: SolexaQA (Average QV, Length, N, 3’ Trim)
Metagenome or Amplicon: Calculate Shannon entropy of first 20 sequence in reads
Dereplication and DRISEE: Identify duplicate in which first 50 bp identical reads >20 times
Screening: Bowtie2 against specified organism genome to remove host genome
Gene Calling: FragGeneScan (>75 bp)
AA Clustering: UCLUST (AA Identity 90%, representative sequence is the longest one)
Protein Identification: sBLAT (OpenMP parallelzation) against M5nr (GenBank, SEED, IMG,UniProt, KEGG, eggNOG)
Annotation Mapping: SEED Subsystems, IMG terms, COG, eggNOGs, GO
Abundance Profiles: E-value, Identity, Alignment length can be specified
rRNA pipeline
BLAT search against 90% clustered SILVA. Identified reads are then clustered at 97% identity.
Longest sequence is the representative of the cluster. BLAT searched against the M5rna 54(SILVA, Greengenes, RDP)