メタゲノム解析(2017NGSハンズオン講習会-2017年8月31日)

Size: px

Start display at page:

Download "メタゲノム解析(2017NGSハンズオン講習会-2017年8月31日)"

いとはあきくぼ
5 years ago
Views:

1 2017 年 NGS ハンズオン講習会 8 月 31 日メタゲノム解析森宙史 (Hiroshi Mori), Ph.D. 国立遺伝学研究所生命情報研究センター hmori@nig.ac.jp 1

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

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

4 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.) DNA extraction Sequence clustering with species level by CD-HIT-EST or UCLUST, etc.

5 16S rrna gene amplicon sequencing analysis ( メタ 16S 解析 ) Pre-analysis (Remove Primer, Chimera etc.) DNA extraction Sequence clustering with species level by CD-HIT-EST or UCLUST, etc. PCR amplification DNA Sequencing Taxonomic assignment and Comparison between samples Togo picture gallery by DBCLS is licensed under a Creative Commons Attribution 2.1 Japan license (c) Who s there? 5

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

DNA extraction Contig sets Gene finding MGA, MetaGeneMark

DBCLS is licensed under a Creative Commons Attribution 2.

Comparative metagenomics Gene Function abundance Pathway

6 Metagenomic sequencing analysis ( メタゲノム解析 ) Metagenomic reads Assemble MEGAHIT, MetaPlatanus IDBA-UD etc. DNA extraction Contig sets Gene finding MGA, MetaGeneMark Gene sets DNA Sequencing BLASTP Togo picture gallery by DBCLS is licensed under a Creative Commons Attribution 2.1 Japan license (c) Sample1 Metadata Sample2 Metadata Comparative metagenomics Gene Function abundance Pathway abundance Pathway reconstruction Taxonomic abundance Who s there? What are they doing? 6

7 メタ 16S 解析利点安価かつ少量の DNA から系統組成が得られる reference 配列に依存しない解析も可能マシンパワーは少なくて済み解析ツールも普及 (QIIME mothur 等 ) 欠点 PCR バイアスの存在種以下は分解能に問題あり個々の系統の機能が不明メタゲノム解析利点系統組成と遺伝子機能組成が得られる実験によるバイアスが少ない優占系統のドラフトゲノムの構築 ( 条件が良ければ可能 ) 欠点 reference 配列に依存した解析目的依存で解析手法が変化しマシンパワーも必要 7

8 16S rrna 遺伝子配列を用いた既知 Bacteria Phylum の系統樹 Genome sequence available No genome sequence available 8

9 16S rrna gene amplicon sequencing analysis ( メタ 16S 解析 ) Pre-analysis (Remove Primer, Chimera etc.) DNA extraction Sequence clustering with species level by CD-HIT-EST or UCLUST, etc. PCR amplification DNA Sequencing Taxonomic assignment and Comparison between samples Togo picture gallery by DBCLS is licensed under a Creative Commons Attribution 2.1 Japan license (c) Who s there? 9

10 PCRに使うプライマーは大きく分けて2 種類系統特異的なプライマー系統 Universalなプライマー両者は何が違うのか? 10

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

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

13 Degenerate primer 例 : 5 -GTGCCAGCMGCCGCGGTAA-3 曖昧 ( 縮重 ) 塩基を使ったプライマー曖昧塩基塩基 1 塩基 2 塩基 3 塩基 4 R A G Y C T S G C W A T K G T M A C B C G T D A G T H A C T V A C G N A C G T 13

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

15 菌類 ( 糸状菌酵母など ) 28S rrna 遺伝子 18S rrna 遺伝子 COX1 COX2 rdna ITS1 (internal transcribed spacer), rdna ITS2 18S rdna ITS1 5.8S rdna ITS2 28S rdna 15

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

17 細菌 16S rrna 遺伝子 rdna ITS 真核生物 : 18S rdna ITS1 5.8S rdna ITS2 28S rdna 原核生物 : 16S rdna ITS 23S rdna 5S rdna 17

18 16S rrna gene amplicon sequencing analysis ( メタ 16S 解析 ) Pre-analysis (Remove Primer, Chimera etc.) DNA extraction Sequence clustering with species level by CD-HIT-EST or UCLUST, etc. PCR amplification DNA Sequencing Taxonomic assignment and Comparison between samples Togo picture gallery by DBCLS is licensed under a Creative Commons Attribution 2.1 Japan license (c) Who s there? 18

19 メタ 16S 解析パイプラインの例 : Ion PGM リードの前処理 Ion PGM リード fastq ファイル N を含むリードを除外長さフィルタ ( 増幅断片長の中央値 ± 約 50 base 以内のリードのみ採用 ) FASTQ から FASTA と QUALITY ファイルへの変換クオリティフィルタ (Average QV >= 25) TagCleaner による 5 Primer 除去 TagCleaner による 3 Primer 除去前処理済み Ion PGM リード 19

20 キメラ除去前処理済み 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 Table 20

系統アサインメントキメラ除去済み OTU RDP Classifier を用いて

21 系統アサインメントキメラ除去済み OTU RDP Classifier を用いて bootstrap cutoff >= 50% で各 OTU の代表配列を系統アサインメント各 OTU のサンプルごとのリード数を整理した OTU table Domainレベルの系統組成 Phylumレベルの系統組成 Classレベルの系統組成 Orderレベルの系統組成 Familyレベルの系統組成 Genusレベルの系統組成各 OTU の系統情報とサンプルごとのリード数を整理した OTU table 21

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

23 23

24 24

25 25

http://www.drive5.com/usearch/ http://www.

26 26

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

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

29 29

30 30

31 Metagenomic sequencing analysis ( メタゲノム解析 ) Metagenomic reads Assemble MEGAHIT, MetaPlatanus IDBA-UD etc. DNA extraction Contig sets Gene finding MGA, MetaGeneMark Gene sets DNA Sequencing BLASTP Togo picture gallery by DBCLS is licensed under a Creative Commons Attribution 2.1 Japan license (c) Sample1 Metadata Sample2 Metadata Comparative metagenomics Gene Function abundance Pathway abundance Pathway reconstruction Taxonomic abundance Who s there? What are they doing? 31

32 細菌群集 ( 山本森山田黒川, 2014 生命のビッグデータ利用の最前線シーエムシー出版より一部改変 ) DNA 抽出シーケンスライブラリ作成シーケンシング塩基配列データメタデータ記録塩基配列データのクオリティコントロールメタデータリスト Contig / Scaffold 配列メタゲノムアセンブル高精度塩基配列データ機能情報付きアミノ酸配列 DB ( 例 :KEGG, GenBank nr 等 ) 遺伝子予測遺伝子のアミノ酸配列データ配列相同性検索各遺伝子の機能と群集中の存在量アセンブル結果への配列マッピング各系統の群集中の存在量統計学の手法を用いた比較メタゲノム解析その細菌群集の特異性や他の細菌群集との類似性等の知見配列相同性検索 16S rrna 遺伝子配列 DB ( 例 :SILVA, RDP 等 ) 既知のメタゲノム解析データのDB ( 例 :MG-RAST, MicrobeDB.jp 等 ) 32

33 Next Generation Sequencer (NGS) ( データはちょっと古いです ) Sequencer Name Specific property Read length (base) ABI 3730xl Sanger Read number / run Ion Proton Emulsion PCR ,000,000 MiSeq Bridge PCR ,000, 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 / Sequel Single molecule >7,000 >50,000 MinION Single molecule > 数 kb >10,000 33

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

35 目的に応じたハードウェアのスペックの目安数千本以上の配列相同性検索特に重要なハードウェア : CPU ゲノムやメタゲノムの遺伝子アノテーション複数ゲノムの比較解析入力配列が数千本以上になることが多く計算を並列化して高速化する必要がある現在の nr 相手の BLASTX や BLASTP は並列化しても各プロセスで 5GB ほどメモリを使用するのでメモリもそれなりに必要 CPU 数コア数が非常に重要メモリもある程度必要 35

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

37 目的に応じたハードウェアのスペックの目安 1. 原核真核メタゲノム Transcriptome アセンブル特に重要なハードウェア : メモリ CPU 2. 数千本以上の配列相同性検索特に重要なハードウェア : CPU 3. マッピングツールでの Reference ゲノムへのマッピング特に重要なハードウェア : ディスク 37

38 午後からの実習で用いるメタゲノムデータ Backhed F. et al 2015 の 4 days の新生児の 1 サンプルを 100 万 pair にダウンサンプリングしたデータ 38

39 ヒト乳児腸内細菌群集 ( 山本森山田黒川, 2014 生命のビッグデータ利用の最前線シーエムシー出版より一部改変 ) DNA 抽出シーケンスライブラリ作成シーケンシング塩基配列データメタデータ記録メタデータリスト Contig 配列 Illumina HiSeq 2000 MEGAHIT メタゲノムアセンブル塩基配列データのクオリティコントロール Trimmomatic 高精度塩基配列データ 100 万 pairs 機能情報付きアミノ酸配列 DB ( 例 :KEGG, UniProt KB 等 ) Prodigal 遺伝子予測遺伝子のアミノ酸配列データ配列相同性検索各遺伝子の機能と群集中の存在量 BLASTP アセンブル結果への配列マッピング Bowtie 2 各系統の群集中の存在量統計学の手法を用いた比較メタゲノム解析その細菌群集の特異性や他の細菌群集との類似性等の知見配列相同性検索 BLASTN 16S rrna 遺伝子配列 DB ( 例 :SILVA, RDP 等 ) 既知のメタゲノム解析データのDB ( 例 :MG-RAST, MicrobeDB.jp 等 ) 39

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

41 今日の実習の参考資料 41

42 PhiX について umina_techsupport_session16.pdf 42

43 塩基配列の類似性検索ツールにも様々なものが存在する Reference と近ければ Mapper 系 (Bowtie 2 等 ) Reference と遠いのなら BLAST 系 (Chaisson J. M. et al. 2012) 43

44 ゲノムアセンブルの二大戦略 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

45 メタゲノムアセンブルツールの例 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 の表現方法が簡素化されているため ( 高速で省メモリまた coverage が小さい k-mer の扱いについて色々と工夫している scaffolding はしない metaspades (Nurk et al. 2017) Contig 作成の方法は IDBA-UD と類似しているがリードデータ中の strain レベルの配列多様性を Contig/Scaffold においてもできるだけ保つためにサイトに多型があると Contig を分岐する傾向が強い 45

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

47 メタゲノムデータからの系統組成推定 16S rrna 遺伝子を用いる手法利点 : Reference 配列 DB が充実 (RDP, SILVA, GreenGenes など ) 欠点 : ゲノム内コピー数の問題系統マーカータンパク質の遺伝子 (rpob 等 ) を用いる手法利点 : 1 ゲノム 1 コピー欠点 : Reference 配列 DB が貧弱リードや contig のマッピング or k-mer 組成を用いる手法利点 : Virus, 原核, 真核生物を同時解析可能欠点 Reference 配列 DB が貧弱 Genus 以下の精度がかなり落ちる水平伝播配列の扱い 47

48 パスウェイデータベース先週行われた AJACS 河内の九州大の山西先生の資料が参考になります en&file=170824pathway_db_yamanishi_submit _r.pdf&refer=ajacs66 メタゲノムでは KEGG の KEGG Orthology を遺伝子機能の単位として使うことが多い 48

49 ヒト乳児腸内細菌群集 ( 山本森山田黒川, 2014 生命のビッグデータ利用の最前線シーエムシー出版より一部改変 ) DNA 抽出シーケンスライブラリ作成シーケンシング塩基配列データメタデータ記録メタデータリスト Contig 配列 Illumina HiSeq 2000 MEGAHIT メタゲノムアセンブル塩基配列データのクオリティコントロール Trimmomatic 高精度塩基配列データ 100 万 pairs 機能情報付きアミノ酸配列 DB ( 例 :KEGG, UniProt KB 等 ) MetaProdigal 遺伝子予測遺伝子のアミノ酸配列データ配列相同性検索各遺伝子の機能と群集中の存在量 BLASTP アセンブル結果への配列マッピング Bowtie 2 各系統の群集中の存在量統計学の手法を用いた比較メタゲノム解析その細菌群集の特異性や他の細菌群集との類似性等の知見配列相同性検索 BLASTN 16S rrna 遺伝子配列 DB ( 例 :SILVA, RDP 等 ) 既知のメタゲノム解析データのDB ( 例 :MG-RAST, MicrobeDB.jp 等 ) 49

50 代表的なメタゲノムデータベース NCBI Taxonomy + SRA h.gov/taxonomy/ GOLD 運営者配列データメタゲノム / メタ 16S 区別系統組成遺伝子機能組成サンプル数 ( 月 ) NCBI, USA 673,079 JGI, USA 24,922 ゲノム等との統合化 IMG/M gi-bin/m/main.cgi MG-RAST l.gov/index.html EBI-Metagenomics metagenomics/ JGI, USA Chicago U. USA (Reads + Contigs) ( メタゲノムのみ収録 ) 7,982 47,313 (303,594) EBI, EU 74,342 MicrobeDB.jp NIG, 60,551 50

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

52 52

53 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

54 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)

GWB

GWB NGS データ解析入門 Web セミナー : De Novo シークエンス解析編 1 NGS 新規ゲノム配列解析の手順シークエンス遺伝子領域の検出アセンブルデータベース検索 2 解析ワークフローと使用ソフトウェアシークエンスデータのインポートクオリティチェック前処理コンティグ配列の作成 CLC Genomics Workbench 遺伝子領域の検出 Blast2GO PRO データベース検索