3rd-jikken-ngs

(1)

生命情報実験A

次世代シークエンサーのデータを用いた

ゲノム解析

慶應義塾大学理工学部生命情報学科

榊原康文、佐藤健吾

(2)

ねらい

• これからの生命科学において要となるツールである次世代シークエンサー（NGS）が産生するデータを用いたゲノム解析を体験する。 • 次世代シークエンサー - 長所：高速かつ低コスト - 短所：得られる一本一本の配列が短い（＝ショートリード） • 例：ヒトゲノムの解読 2003年キャピラリーシークエンサー 13年、3000億円以上 2012年次世代シークエンサー 10日間、70万円 Illumina GAIIx

(3)

解析対象 – 納豆菌

-• 納豆菌(Bacillus subtilis natto)は有用物質を作る

- γ-PGA、ナットウキナーゼ、エラスターゼ、ポリアミン

• 枯草菌(Bacillus subtilis)の近縁種

• BEST195株ゲノムの決定 [Nishito et al., 2010]

• 遺伝的多様性を持つ菌種ミツカン(愛知) オシキリ食品(北海道) 萬歳食品(宮城) 黒石納豆(青森) あづま食品(栃木) タカノフーズ(茨城) 水戸納豆(茨城) やぐちフーズ(埼玉) 鎌倉山納豆(神奈川) ヤマダフーズ(秋田) 高畠納豆(山形) 奥野食品(三重) 旭松食品(大阪) 丸美屋(熊本) ホンコン(中国) ソウル(韓国) チェンマイ(タイ) キネマ(ネパール) http://cache.cart-imgs.fc2.com/user_img/freeformat/2500_1_35.jpg

(4)

解析対象 – 納豆菌

-• アジア各国の納豆に似た食品から採取した５株 - 韓国株（KorC1） - ラオス株（LaoA1） - ミャンマー株（MyaA2） - ネパール株（NepD5） - タイ株（ThaB）

Mapping and variant calling

Filtered-corrected short reads sequenced from the eight B. subtilis strains were mapped to the BEST195 genome sequence using BWA. An average of 86.2% of reads were mapped to the BEST195 genome with an average of 193.6-fold coverage across the entire genome. Based on the mapping result, SNPs and INDELs were detected for all strains using GATK, and effect impacts of each variant on a genome were estimated by SnpEff. The statistics of the mapping and variant calls for each strain are summarized in Table C in S1 File.

To score and vectorize detected variants (see Material and Methods), we performed PCA and hierarchical clustering analysis, and the results are shown in Fig 3. The principal compo-nents obtained in PCA were transformations of variant score vectors by a linear combination that was chosen to maximize the variance of the score vectors of all eight strains. As shown in Fig 3-A, the first principal component (contributing rate: 51%) indicates a feature of the non-Japanese strains, and the second principal component (contributing rate: 19%) can be regarded as a feature to distinguish strain LaoA1 from the other non-Japanese strains. The Japanese

strains converged and formed a small cluster. For the first principal component, principal scores are high for genes BSNT_09336, BSNT_09102, and BSNT_09338, which mean that these genes contribute to the first principal component. Although they all are annotated as

Fig 3. The results of PCA and hierarchical clustering based on variant vectors. (A) Biplot of principal component analysis based on variant vectors. The dots show the eightB. subtilis strains, and the upper left image is an enlarged image focused on the three Japanese strains located near (0, 0). The fist principal component features the non-Japanese strains, and the second principal component can be regarded as a feature to distinguish strain LaoA1 and the other non-Japanese strains. (B) Hierarchical clustering of the eight B. subtilis strains based on the Euclidean distance between variant scores of each strain using the furthest neighbor method. The different cluster indicates that strains have different variant score patterns. (C)

Geographical location of each country. doi:10.1371/journal.pone.0141369.g003

Comparative Genome Analysis of B. subtilis from Fermented Foods

PLOS ONE | DOI:10.1371/journal.pone.0141369 October 27, 2015 10 / 21

(5)

解析の流れ

– サンプル調製とシークエンシング

-• ペアエンドリード - 断片化されたショートリードの両端を読む平均インサート長 500bp ペアエンドリード Illumina GAIIx シークエンシング (アジア株５株) 納豆菌ゲノムゲノム抽出サンプル調製断片化 500 bp ペアエンドリード

(6)

解析の流れ de novo アセンブリ

-• 参照するゲノムがない場合の解析方法推定ゲノム（スキャフォルド） _{アノテーション}遺伝子枯草菌納豆菌 BEST195 アジア株 • オーソログ遺伝子 • リピート領域の解析ペアエンドリード（アジア株５株）アセンブリ SPAdes [Bankevich et al, 2012] 遺伝子予測 glimmer [Salzberg et al, 1998] 比較ゲノム Murasaki [Popendorf et al, 2010]

(7)

解析の流れ – 変異解析

-• 参照するゲノムがある場合の解析方法マッピング Bowtie2 [Langmead et al, 2012] 多型の同定 VarScan2 [Koboldt et al, 2012] 変異影響度アノテーション snpEff [Cingolani et al, 2012] BEST195ゲノム _{BEST195ゲノム} A A A A A C • 同義置換 • 非同義置換 • 非コード領域の置換 • … • 遺伝子の機能 • 表現型の解析ペアエンドリード（アジア株５株）

(8)

進め方

• アジア５株のうち、どれか１株を解析する。 - 学科整理番号 % 5 = • 0 ⇒ 韓国株（KorC1） • 1 ⇒ ラオス株（LaoA1） • 2 ⇒ ミャンマー株（MyaA2） • 3 ⇒ ネパール株（NepD5） • 4 ⇒ タイ株（ThaB） • 各々の計算は、ウェブツールGalaxyを操作して実行する。 • 計算結果は、Windows上のツールGMV、IGVを用いて視覚化する。

(9)

Galaxy

• ペンシルバニア州立大学で開発されているWebベースのゲノム解析プラットホーム

• コンピュータに詳しくない研究者でもバイオインフォマティクス研究が可能に！

(10)

準備（Galaxyアカウント作成）

• 実験用Galaxyサーバにアクセスする。

- http://galaxy.dna.bio.keio.ac.jp/

• ユーザ登録を行う。

(11)

準備（Galaxyアカウント作成）

• ユーザ登録を行う。

- 中央上部の「User」→「Register」

ニックネーム

(12)

解析の流れ de novo アセンブリ

(13)

準備（de novo アセンブリ）

• まずはログインする。

- 先ほど作成した「ユーザ名」「パスワード」でログインする。

(14)

• 担当する株のリード配列をHistoryにインポートする

- 「Shared Data」→「Data Libraries」→「natto」→「DNA」

- 「to History」をクリック→「de novo assembly」選択→「Import」

リード配列をインポート

(15)

FASTAフォーマット

• 元々はFASTAというプログラムで使われていた配列フォーマットだが、他のプログラムでも広く使われている。 - 「〜.fa」とか「〜.fasta」というファイル名であることが多い。 • １行目： “>”で始まるヘッダ • ２行目以降：実際の配列

>AP011541 Bacillus subtilis subsp. natto BEST195 DNA, complete genome. ATCTTTTTCGGCTTTTTTTAGTATCCACAGAGGTTATCGACAACATTTTCACATTACCAA

CCCCTGTGGACAAGGCTTTTTCAACAGGTTGTCCGCTTTGTGGATAAGATTGTGACAACC ATTGCAAGCTCTCGTTTATTTTGGTATTATATTTGTGTTTTAACTCTTGATTACTAATCC TACCTTTCCTCTTTATCCACAAAGTGTGGATAAGTTGTGGATTGATTTCACACAGCTTGT GTAGCAGGTTGTCCACAAGTTGTGAAATTTGTCGAAAAGCTATTTATCTACTATATTATA

(16)

FASTQフォーマット

• NGSから出力される配列フォーマット - 「〜.fq」とか「〜.fastq」というファイル名であることが多い。 • １行目： “@”で始まるヘッダ • ２行目：塩基配列 • ３行目： “+” • ４行目：塩基配列のクオリティスコア @HWUSI-EAS1730:24:FC:1:1:2719:1156 1:N:0:TAGCTT GAGGTTAACGGCACATTTCGCGCCAACCATTCCTGCGGACACGATTCNCATATGACAA + HHGGHHHHHHHHDHHGHHHHHHHHDHHHGHGHHBHHHHFHHGHHG@B#B@AEAAEHDB (より詳しくは： https://en.wikipedia.org/wiki/FASTQ_format)

(17)

解析の流れ de novo アセンブリ

(18)

SPAdes

• アセンブリを実行する

この値をいろいろ変えるとアセンブリ結果のパフォーマンス（N50等）が変わる

(19)

SPAdes

• アセンブリを実行する

株の名前_R1.fastq 株の名前_R2.fastq

(20)

ゲノムアセンブリ

(21)

SPAdes

• アセンブリ結果を表示する

- 各配列の長さとカバレッジ（１塩基あたりマッピングされているリード数）

(22)

アセンブリの評価

• 推定ゲノム長 - アセンブリにより得られたコンティグ（スキャフォルド）の長さの和 • 最長コンティグ（スキャフォルド）長 - 最も長いコンティグ（スキャフォルド）長 • N50長、N90長 - コンティグ（スキャフォルド）を長さ順につなげていった時、長さの累積が推定ゲノム長の50(90)%を越えた時のコンティグ（スキャフォルド）の長さ - この長さ以上のコンティグ（スキャフォルド）の長さを足すと、推定ゲノム長の50(90)%を越える推定ゲノム長推定ゲノム長の50% 最長コンティグ長 N50長

(23)

SPAdes stats

(24)

SPAdes stats

• アセンブリ結果を評価する

(25)

Summary Statistics

- 推定ゲノム長、最大スキャフォルド長

Scaffold statsテーブルの第２カラムという意味

(26)

Summary Statistics

(27)

解析の流れ de novo アセンブリ

(28)

Filter SPAdes output

(29)

Glimmer3

(30)

Convert Glimmer to GFF

• 遺伝子予測結果を他のプログラムで使える形式に変換する

同じ名前のデータが二つあるが、番号が若い方を選ぶ

(31)

GFFフォーマット

• General Feature Format

• ゲノムアノテーション（遺伝子など）を記述するためのフォーマット

natto0 confirmed CDS 539263 540387 . + . ID=BSNT_00864 natto0 confirmed CDS 777821 778975 . + . ID=BSNT_01234 natto0 confirmed CDS 1283805 1285061 . + . ID=BSNT_02038 natto0 confirmed CDS 1474240 1475478 . - . ID=BSNT_02338 natto0 confirmed CDS 1545714 1546838 . + . ID=BSNT_02447

配列名開始座標終了座標ストランドその他の属性

（IDなど）タイプ

(32)

解析の流れ de novo アセンブリ

(33)

枯草菌と納豆菌BEST195のゲノム配列を準備

• 「Shared Data」 → 「Data Libraries」→「natto」からインポート

• 「Shared Data」→「Data Libraries」→「natto」→「Reference」

枯草菌

納豆菌 BEST195

(34)

GenBankフォーマット

• 配列アノテーションを記述するフォーマットの一つ • GenBankで使用されている

LOCUS AP011541 4091591 bp DNA linear HTG 25-MAR-2010 DEFINITION Bacillus subtilis subsp. natto BEST195, *** SEQUENCING IN PROGRESS

***, 12 ordered pieces. ACCESSION AP011541

SOURCE Bacillus subtilis subsp. natto BEST195 ORGANISM Bacillus subtilis subsp. natto BEST195

Bacteria; Firmicutes; Bacillales; Bacillaceae; Bacillus. REFERENCE 1

AUTHORS Nishito,Y., Osana,Y., Hachiya,T., Popendorf,K., Toyoda,A., Fujiyama,A., Itaya,M. and Sakakibara,Y.

TITLE Whole genome assembly of a natto production strain Bacillus subtilis natto from very short read data

JOURNAL BMC Genomics (2010) In press FEATURES Location/Qualifiers

source 1..4091591

/organism="Bacillus subtilis subsp. natto BEST195" /mol_type="genomic DNA" /strain="BEST195" /sub_species="natto" /db_xref="taxon:645657" gene 410..1750 /gene="dnaA" /locus_tag="BSNT_00001" CDS 410..1750 /gene="dnaA" /locus_tag="BSNT_00001" /note="from glimmer orf00001" /codon_start=1

/transl_table=11

/product="chromosomal replication initiation protein" /protein_id="BAI83445.1" /db_xref="GI:291482370" /translation="MENILDLWNQALAQIEKKLSKPSFETWMKSTKAHSLQGDTLTIT APNEFARDWLESRYLHLIADTIYELTGEELSIKFVIPQNQDVEDFMPKPQVKKAVKED TSDFPQNMLNPKYTFDTFVIGSGNRFAHAASLAVAEAPAKAYNPLFIYGGVGLGKTHL MHAIGHYVIDHNPSAKVVYLSSEKFTNEFINSIRDNKAVDFRNRYRNVDVLLIDDIQF

(35)

Murasaki

(36)

結果をダウンロード

• 比較ゲノム解析の結果をGMVで可視化するために、Murasakiと

(37)

GMV (Murasaki Viewer)

• GMVをダウンロードして適当な場所で解凍

http://www.dna.bio.keio.ac.jp/lecture/jikken/data/gmv-win.zip

• Murasakiの実行結果が入っているzipファイルを解凍して実行 • 展開先/output/test.anchors を開く

(38)

GMV (Murasaki Viewer)

• Glimmer3による遺伝子予測結果を貼り付ける [File]-[Load Annotation File]から

(39)

GMV (Murasaki Viewer)

• アンカーを並べ替える

[Edit]-[Sort Sequence Sources by Anchors]から

(40)

GMV (Murasaki Viewer)

�� _{�� }�� _{�� }�� _�� _�� _�� _�� _��_�� _�� _�� _� _��_�� _�� _� _�� _�� _�� _� ��

(41)

解析の流れ – 変異解析

(42)

準備（変異解析）

• 新しいHistoryを作る

ここから「Create New」を選ぶ

(43)

• 担当する株のリード配列をHistoryにインポートする

- 「Shared Data」→「Data Libraries」→「natto」→「DNA」

- 「to History」をクリック→「variant analysis」選択→「Import」

リード配列をインポート

(44)

解析の流れ – 変異解析

(45)

Bowtie2

• ショートリードを参照ゲノムにマッピングする

株の名前_R1.fastq 株の名前_R2.fastq

(46)

Bowtie2

• ショートリードを参照ゲノムにマッピングする

平均インサート長＋200

(47)

(48)

SAM Tools flagstat

(49)

SAM/BAMフォーマット

• 配列をゲノムにマッピングしたときに生成されるアラインメントのフォーマット

- SAMフォーマット：プレーンテキスト形式

- BAMフォーマット： SAMフォーマットをバイナリ形式したもの

HWUSI-EAS1730:24:FC:1:51:11857:4571 99 natto0 1 42 58M = 398 455 ATCTTTTTCGGCTTTTTTTAGTATCCACA GAGGTTATCGACAACATTTTCACATTACC IIHHIIIIHIIIIIIIIIIDE8GGGIIIDFHHIGHIHIEIGIHIIGG:GGGDDGGIGH AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i: 0 NM:i:0 MD:Z:58 YS:i:-1 YT:Z:CP RG:Z:KorC1

HWUSI-EAS1730:24:FC:1:112:3847:12416 163 natto0 1 23 4M8I46M = 400 457 ATGTGGATATCTTTTTCGGCTTTTTTTAG TATCCACAGAGGTTATCGACAACATTTTC @DBDGGGG?GGIIIIIIBIIIFIIIIHHIIG?G?GGEDGIHHHIIGDGGGDII@DGGG AS:i:-34 XN:i:0 XM:i:1 XO:i: 1 XG:i:8 NM:i:9 MD:Z:2C47 YS:i:0 YT:Z:CP RG:Z:KorC1

(50)

マッピングの評価

(51)

解析の流れ – 変異解析

(52)

Mpileup call variants

(53)

VarScan2

(54)

VCFフォーマット

• Variant Call Format

- 配列の多型を記述する

#CHROM POS ID REF ALT QUAL FILTER INFO ..

AP011541 1562 . T A . PASS ADP=132; ..

AP011541 1617 . C A . PASS ADP=143; ..

AP011541 1992 . A C . PASS ADP=129; ..

AP011541 2185 . A T . PASS ADP=116; ..

(55)

解析の流れ – 変異解析

(56)

snpEff

(57)

多型の種類

• 一塩基多型（SNP） - 同義置換（synonymous coding） - 非同義置換（non-synonymous coding） • 挿入・欠失（INDEL） - フレームシフト（frame-shift） • 遺伝子領域以外の多型

ATG AAT TGC AGC ACC ...

M N C S T

SNP

ATG AAT TGC ATT ACC ...

M N C I T

SNP

ATG AAT TGG CAG TAC C....

M N W Q Y

挿入

ATG AAT TGC AGT ACC ... M N C S T 参照ゲノム Moderate Low High Effect Modifier 詳しくは http://snpeff.sourceforge.net/SnpEff_manual.html#eff を参照

(58)

変異解析結果

(59)

変異解析結果

(60)

結果を可視化

• 下のリンクをクリックすると、IGVにマッピング結果と変異解析の結果を表示することができる。

(61)

IGV (Integrative Genomics Viewer)

• IGVをダウンロードして適当な場所で解凍

http://www.dna.bio.keio.ac.jp/lecture/jikken/data/IGV_2.3.71.zip

• IGV.batをクリックして実行

• リファレンスゲノムをダウンロード

(62)

IGV (Integrative Genomics Viewer)

遺伝子アノテーションマッピング結果

(63)

レポート課題について

• ウェブページにある「レポート課題の進め方」をよく読む • 割り当てられた株に関して以下を報告する。 - アセンブリの推定ゲノム長、最長コンティグ長、N50長 - BEST195に対するマッピング率 - 同定された変異の種類と数 • Murasakiで計算した比較ゲノムの結果を使って、トランスポゾン（“transposase”で検索）を探す。遺伝子領域に入り込んで破壊していないかを調べる。 • 比較ゲノムやリシークエンシングの結果を使って、納豆のねばねばに関係する遺伝子を中心に調べる。