2016.02.01版

講義室後ろにあるUSBメモリ中のhogeフォルダをデスクトッ
プにコピーしておいてください。

農学生命情報科学
特論I 第3回
東京大学大学院農学生命科学研究科
アグリバイオインフォマティクス教育研究ユニット
門田幸二
kadota@iu.a.u-tokyo.ac.jp

Jun 25, 2014

1

講義予定


教科書

第1回（2014年6月11日）
 西：NSG概論。現状や展望など。講義のみ

第2回（2014年6月18日）
 門田：データベース、データ取得、ファイル形式および変換、前処理
 教科書の1.3節周辺

第3回（2014年6月25日）
 門田：アセンブル、マッピング、カウント情報取得

 教科書の2.3節周辺

第４回（2014年7月2日）
 門田：クラスタリング、データ正規化、実験デザイン、分布（モデル）、発現変動解析
 教科書の3.3節周辺

授業の目標・概要
次世代シーケンサ(NGS)の普及により、以前は主にゲノム解析系で必要とされていた配列解析のためのスキルがトランスクリプトーム解析においても要求される時代になっています。本科目では、様々な局面で応用可能な配列解
析系のスキルアップを目指し、RNAシークエンス（RNA-Seq）に基づく（非モデル生物の）トランスクリプトーム解析を題材とした実習を含む講義を行います。
Jun 25, 20142

エラーの具体例(2014年6月13日)

2013年11月1日のセミナーで見せた結果

エラーの原因はメモリ不足だそうです by 孫堅強氏(2014年6月19日)

Jun 25, 20143

課題遂行時に何人か遭遇したエラーの解説

何人かの方が、作業ディレクトリの変更も正しく行い、SRR609266.fastq.gzファイルもhogeフォルダ中に存在するにも関わらず、入力ファイル読み込み時にエラーに遭遇しました。
この理由は2つ考えられます。1つめは、USBメモリにコピーする際に正しくコピーできていなかった可能性、そして2つめはUSBメモリ中のSRR609266.fastq.gzファイル段階では正しいものであったが、各自のPCにコピーする際に正しくコピーできなかった可能性です。講義中に述べたMD5チェックサム(MD5check sum)でファイルの同一性を確認するのは重要ですね。
Jun 25, 2014

4

Contents（第3回）


アセンブル(Assembly) 2つのアプローチ(two approaches)

Comparative approach (reference-based assembly; resequencing)：同一生物種または近縁種のゲノム配列を利用 de novo approach：過去に配列決定されたものの中に近縁種がない場合
 アルゴリズム（計算手順）


k-mer解析 ゲノム用、トランスクリプトーム用、雑感

マッピング（QuasRパッケージを利用）
 シミュレーションデータを用いたマッピングの基礎
 リアルデータのマッピング（カイコsmall RNA-seqデータ）
 課題



カウント情報取得

Jun 25, 20145

ゲノムアセンブル
Comparative approach 同一生物種または近縁種のゲノム配列を利用するreference-based assembly。

Resequencingともいう。
 ヒトゲノムresequencingやSNP解析系はこちら
一個人のヒトゲノム(Wheeler et al., Nature, 452: 872-876, 2008) 日本人ゲノム(Fujimoto et al., Nat. Genet., 42: 931-936, 2010)
 ENCODE project (ENCODE Project Consortium et al., Nature, 489: 57-74, 2012)

de novo approach 過去に配列決定された生物種以外が主な対象。
 パンダ(Li et al., Nature, 463: 311-317, 2008)
 サンゴ(Shinzato et al., Nature, 476: 320-323, 2011)（NGS由来の比較的短い）配列決定されたリードのみから、目的生物種のゲノム配列を決めること（組み立てること）

Jun 25, 2014

6

参考

Tips


断片化されたゲノム配列

リード(read) Sequencerで読んだ塩基配列のこと



コンティグ(contig)
 異なる複数のリードがACGTの切れ目なくペアードエンド解析

連結されたもの
 右図ではA-Dの四つのコンティグ
Scaffold (supercontig)
 コンティグ間の位置関係を表したもの 「A-D-B-C」ではなく「A-B-C-D」という関係



アセンブル

N50
 得られた複数のコンティグを最も長いコンテ

ィグから順番に連結していったときに
combined total lengthの50%になったときのコンティグの長さ
Jun 25, 2014

A

B

C

D

7

参考

Tips


Coverage（カバレッジ）

赤字の部分、間違って「小さめ」と書いてしまっていたことに2016年2月1日に気づいたので修正しました。この種のミスはいけませんね。失礼しましたm(_ _)m ゲノム解読したいときなどに、解読するために必要とされる指標となる数値。ゲノムサイズ(X)に対する、sequencerで読んだ塩基配列長の和のこと。一般に、この数値が高いほどよい。
Sequence depthという表現と実質的に同じような指標。


アセンブル時に用いるkの値はいくつがいいの？
 複数のkの値を試すようです。2013年ごろから自動的に決めてくれるものが増えたようです。



アセンブル結果の評価基準は？
 よくわかりません。平均コンティグ長やN50が論文の表でよく記述されます。このあたりの数値を大きくするだけなら、kの値を大きめにすればいいです。
 ただ、ゲノムアセンブリの場合には実質的には長ければ長いほどよいという感じみたいです。


アセンブルプログラムを実行して得られる出力ファイルはどんな感じ？
 （基本的に）multi-FASTA形式のファイルです。
>contig1
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT…
>contig2
ACGATGCAGCCTTAACGA…
>contig3…

Jun 25, 20148

El-Metwally et al., PLoS Comput Biol., 9: e1003345, 2013

ゲノムアセンブルの手順
1. 前処理(pre-processing filtering) クオリティの低いリードやコンタミを除去するステップ。塩基置換(substitution)やイ

ンデル(indels; insertion/deletion)を含むリードの除去や補正(error correction)。
 4つのアプローチ：k-mer, suffix tree/array, multiple sequence alignment, hybrid

2. グラフ構築(graph construction)
 前処理後のリードを用いて、リード間のオーバーラップ(overlap)を頼りにつなげていくステップ。シークエンスエラー(sequencing error)と多型(polymorphism)の違いを見るべく、グラフ構築時にエラー補正を行うものもある。
 4つのアプローチ：OLC, de Bruijn graph (k-mer), greedy, hybrid

3. グラフ簡易化(graph simplification)
 グラフ構築後に、複雑化したグラフをシンプルにしていくステップ。連続したノード(nodes; 頂点) やバブルのマージ作業に相当。

4. 後処理(post-processing) コンティグ(contigs)やスカッフォールド(scaffolds)を得るステップ。ミスアセンブリの同定も含む。
大きく分けて4つの手順からなる
Jun 25, 2014

9

様々な戦略があります。
k-merを利用する方法の基本を紹介しますJun 25, 201410

1. 前処理：k-mer


k-mer頻度解析
 NGSデータを入力として、リード長より短いk連続塩基からなる部分文字列を発生させるのが最初のステップ。発生させたk-merの出現頻度情報をもとに、カバレッ
ジ、ゲノムサイズ推定、コンタミリードの除去などを行う
 例1：20塩基長のNGSリードをk=19で分割すると2個のk-merを発生可能
CACCAGGACATGAAGACGCGCACCAGGACATGAAGACGCACCAGGACATGAAGACGCG 例2：20塩基長のNGSリードをk=17で分割すると4個のk-merを発生可能

CACCAGGACATGAAGACGCGCACCAGGACATGAAGACACCAGGACATGAAGACG
CCAGGACATGAAGACGCCAGGACATGAAGACGCG
L塩基長のNGSリードをk-merに分割すると(L – k + 1)個のk-merを発生可能
Jun 25, 201411

項目名は変更する可能性あり

k=19としてk-mer頻度分布を作成。
入力ファイルについて説明します。
Jun 25, 201412

1. 前処理：k-mer

sample32_ngs.fastaは、50塩基からなるランダム
大きく分けて4つの手順からなる塩基配列をリファレンスとしています。20塩基からなる10個の部分配列をNGSリードとしています。
Jun 25, 201413

1. 前処理：k-mer

各リード中のdescription部分の記大きく分けて4つの手順からなる
述が部分配列の位置情報に相当
Jun 25, 201414

Jun 25, 2014k=19としてk-merの種類ごとに頻度情報を取得したのがkmerオブジェクト

15

1回出現したk-merが12個、2回出現したk-merが3個、という解釈をする

Jun 25, 201416

1回出現したk-merが6個、6回出現したk-merが5個、という解釈をする

Jun 25, 201417

入力ファイル中にAAGACという部分文字列は3つ存在するということリマインドですJun 25, 201418

1. 前処理：k-mer

sample32_ngs.fastaは、50塩基からなるランダム塩基配列をリファレンスとしています。20塩基からなる10個の部分配列をリードとしています。総塩基数は200なのでリファレンス配列の4倍の長さ、つまり4X coverageです。

Jun 25, 201419

sample36_ngs.fastaは、10,000塩基からなるランダム塩基配列をリファレンスとしています。80塩基からなる5,000個の部分配列をリードとしています。総塩基数は400,000なのでリファレンス配列の40倍の長さ、つまり40X coverageです。

1. 前処理：k-mer

Jun 25, 201420

(k-merの種類は問わずに)37回出現したk-mer
が2個あったということ

出現回数のmedianは22可能なk-merの種類数は4^21個。実際の種類
数は9,978個でゲノムサ
イズ(=10,000)と酷似

40Xの分布。シークエンスエラーがない場合に相当
Jun 25, 201421

(k-merの種類は問わずに)31回出現したk-mer
が97個あったということ

出現回数のmedianは22可能なk-merの種類数は4^21個。実際の種類
数は9,976個でゲノムサ
イズ(=10,000)と酷似

40Xの分布。シークエンスエラーがない場合に相当。実行ごとに結果が多少異なりますが、どこかで計算をさぼらな
いといけない現実的な理由があるためかも(バグかも)…。

Jun 25, 201422

(k-merの種類は問わずに)59回出現したk-mer
が19個あったということ

出現回数のmedianは48可能なk-merの種類数は4^21個。実際の種類
数は9,978個でゲノムサ
イズ(=10,000)と酷似

100Xの分布。シークエンスエラーがない場合に相当。coverageの増加(40X  100X)に伴い、出現頻度分布x軸方向で右側にシフトしていることがわかる
Jun
25, 201423

(k-merの種類は問わずに)26回出現したk-mer
が5個あったということ

出現回数のmedianは42可能なk-merの種類数は4^31個。実際の種類
数は9,968個でゲノムサ
イズ(=10,000)と酷似

100Xの分布。シークエンスエラーがない場合に相当。k-merの値が変わっても(k=21 31)、k-merの種類数はゲノムサイズとほぼ同じ
Jun
25, 201424

(k-merの種類は問わずに)67回出現したk-mer
が15個あったということ

出現回数のmedianは54可能なk-merの種類数は4^11個。実際の種類
数は9,965個でゲノムサ
イズ(=10,000)と酷似

100Xの分布。シークエンスエラーがない場合に相当。k-merの値が変わっても(k=31  11)、
k-merの種類数はゲノムサイズとほぼ同じだが、実際にはエラーを多く含むので非現実的。
Jun25, 2014

25

Chikhi and Medvedev, Bioinformatics, 30: 31-37, 2014のFig. 1

k-mer解析は様々な場面で利用されます① 大まかなゲノムサイズ推定（とカバレッジ） 現実のNGSデータはシークエンスエラーやリピート配列を含むため正規分布っぽくはならないが主要なピークの位置情報をもとにゲノムサイズを推定可能（らしい）
 リード長とリード数から得られたNGSデータ中の総塩基数情報はわかっているので、得られたゲノムサイズで割ることでカバレッジ(coverage)が分かる

Jun 25, 2014

k = 41, 51, 61

k = 41, 51, 61

k = 41, 51, 61

S. aureus

H. sapiens chr 14

B. impatiens

ゲノムサイズ：2.8 Mb
Coverage： 167X
リード数：5,000,000
リード長：101 bp

ゲノムサイズ：88 Mb
Coverage： 70X
リード数：62,000,000
リード長：101 bp

ゲノムサイズ：250 Mb
Coverage： 247X
リード数：497,000,000
リード長：124 bp

26Chikhi and Medvedev, Bioinformatics, 30: 31-37, 2014のFig. 1

k-mer解析は様々な場面で利用されます② 前処理（フィルタリング） シークエンスエラーやリピート配列由来リードの除去

k = 41, 51, 61

ゲノム
4X
20X

シークエンスエラーに由来する、数回しか出現しないk-mer由来リードは除去対象
Jun 25, 2014

リピート配列に由来する、異常に多く出現するk-mer由来
リードも除去対象

27

Kajitani et al., Genome Res., in pressのFig. 2

k-mer解析は様々な場面で利用されます③ 前処理（エラー補正）
 ヘテロ接合度の高い2倍体ゲノム(highly heterozygous diploid genomes)由来リードの補正
父親由来ゲノム…CACCAGGACATGAAGACGCGTTCA…
…CACCAGGACATCAAGACGCGTTCA…
母親由来ゲノム

シークエンスエラーに由来する、数回しか出現しな
いk-mer由来リードは除去対象Jun 25, 2014

リピート配列に由来する、異常に多く出現する
k-mer由来リードも除去対象CACCAGGACATACCAGGACATG
CCAGGACATGA
父親由来k-mer…
GACATGAAGAC…
TGAAGACGCGTGAAGACGCGTCAAGACGCGTCA
両親間で配列の異なる領域を含むkmerの出現頻度は、配列が同じ領域由来k-merの出現頻度(c)の1/2となる。アセンブルに悪影響を及ぼすため、母親由来ゲノム配列に揃えるなどして補正
28

El-Metwally et al., PLoS Comput Biol., 9: e1003345, 2013

ゲノムアセンブルの手順
1. 前処理(pre-processing filtering) クオリティの低いリードやコンタミを除去するステップ。塩基置換(substitution)やイ

ンデル(indels; insertion/deletion)を含むリードの除去や補正(error correction)。
 4つのアプローチ：k-mer, suffix tree/array, multiple sequence alignment, hybrid

2. グラフ構築(graph construction)
 前処理後のリードを用いて、リード間のオーバーラップ(overlap)を頼りにつなげていくステップ。シークエンスエラー(sequencing error)と多型(polymorphism)の違いを見るべく、グラフ構築時にエラー補正を行うものもある。
 4つのアプローチ：OLC, de Bruijn graph (k-mer), greedy, hybrid

3. グラフ簡易化(graph simplification)
 グラフ構築後に、複雑化したグラフをシンプルにしていくステップ。連続したノード(nodes; 頂点) やバブルのマージ作業に相当。

4. 後処理(post-processing) コンティグ(contigs)やスカッフォールド(scaffolds)を得るステップ。ミスアセンブリの同定も含む。
大きく分けて4つの手順からなる
Jun 25, 201429

様々な戦略があります
が、k-merを利用する方法の基本を紹介します

Jun 25, 201430

El-Metwally et al., PLoS Comput Biol., 9: e1003345, 2013

2.グラフ構築(graph construction)


k-merアプローチ(de Bruijnグラフ)
 リードを全ての可能なk-merに分割し、有向グラフを作成(k=9の例)

リード1: CACCAGGACATGAAGACGCG

リード2: CCAGGACATGAAGACGCGTT

CACCAGGACACCAGGACA
CCAGGACAT…
ATGAAGACGTGAAGACGCGAAGACGCG

CCAGGACATCAGGACATGAGGACATGA…
GAAGACGCGAAGACGCGTAGACGCGTTCACCAGGAC

CCAGGACATGAAGACGCGCAGGACATG…
ACCAGGACA
Jun 25, 2014…

GAAGACGCG

AGACGCGTT

CCAGGACAT○：ノード（node；頂点）→：エッジ（edge；辺）AAGACGCGT
31

El-Metwally et al., PLoS Comput Biol., 9: e1003345, 2013

2.グラフ構築(graph construction)


k-merアプローチ(de Bruijnグラフ)
 同一ノードをマージしてde Bruijnグラフを作成(k=9)

リード2: CCAGGACATGAAGACGCGTT

リード1: CACCAGGACATGAAGACGCGCACCAGGAC

CCAGGACATGAAGACGCGCAGGACATG…

…

GAAGACGCG
AGACGCGTT

CCAGGACATACCAGGACAAAGACGCGTCACCAGGAC

CCAGGACATGAAGACGCG
AGACGCGTT…
ACCAGGACA
Jun 25, 2014CAGGACATGAAGACGCGT
32

El-Metwally et al., PLoS Comput Biol., 9: e1003345, 2013

ゲノムアセンブルの手順
1. 前処理(pre-processing filtering) クオリティの低いリードやコンタミを除去するステップ。塩基置換(substitution)やイ

ンデル(indels; insertion/deletion)を含むリードの除去や補正(error correction)。
 4つのアプローチ：k-mer, suffix tree/array, multiple sequence alignment, hybrid

2. グラフ構築(graph construction)
 前処理後のリードを用いて、リード間のオーバーラップ(overlap)を頼りにつなげていくステップ。シークエンスエラー(sequencing error)と多型(polymorphism)の違いを見るべく、グラフ構築時にエラー補正を行うものもある。
 4つのアプローチ：OLC, k-mer (de Bruijn graph), greedy, hybrid

3. グラフ簡易化(graph simplification)
 グラフ構築後に、複雑化したグラフをシンプルにしていくステップ。連続したノード(nodes; 頂点) やバブルのマージ作業に相当。

4. 後処理(post-processing) コンティグ(contigs)やスカッフォールド(scaffolds)を得るステップ。ミスアセンブリの同定も含む。
大きく分けて4つの手順からなる
Jun 25, 201433

El-Metwally et al., PLoS Comput Biol., 9: e1003345, 2013

3. グラフ簡易化(graph simplification)

連続したノード(nodes; 頂点) やバブルのマージ
リード2: CCAGGACATGAAGACGCGTT

リード1: CACCAGGACATGAAGACGCGCACCAGGAC

CCAGGACATGAAGACGCG
AGACGCGTT…

ACCAGGACACAGGACATGAAGACGCGTCACCAGGACATGAAGACGCGTT

Jun 25, 2014

この2つのリードだけで簡易化した結果

34

Kajitani et al., Genome Res., in pressのSuppl. Fig. 8に相当

3. グラフ簡易化(graph simplification)

連続したノード(nodes; 頂点) やバブルのマージ父親由来ゲノム…CACCAGGACATGAAGACGCGTTCA…
…CACCAGGACATCAAGACGCGTTCA…
母親由来ゲノム

父親由来ゲノム
リード1: CACCAGGACATGAAGACGCG

リード2: CCAGGACATCAAGACGCGTT

CAGGACATGAAGACGCGCACCAGGACATAAGACGCGTT
CAGGACATCAAGACGCG
SNPなど塩基に違いがあればバブル構造になります

Jun 25, 201435

El-Metwally et al., PLoS Comput Biol., 9: e1003345, 2013

ゲノムアセンブルの手順
1. 前処理(pre-processing filtering) クオリティの低いリードやコンタミを除去するステップ。塩基置換(substitution)やイ

ンデル(indels; insertion/deletion)を含むリードの除去や補正(error correction)。
 4つのアプローチ：k-mer, suffix tree/array, multiple sequence alignment, hybrid

2. グラフ構築(graph construction)
 前処理後のリードを用いて、リード間のオーバーラップ(overlap)を頼りにつなげていくステップ。シークエンスエラー(sequencing error)と多型(polymorphism)の違いを見るべく、グラフ構築時にエラー補正を行うものもある。
調べると沢山見つかります 4つのアプローチ：OLC, k-mer (de Bruijn graph), greedy, hybrid

3. グラフ簡易化(graph simplification)
 グラフ構築後に、複雑化したグラフをシンプルにしていくステップ。連続したノード(nodes; 頂点) やバブルのマージ作業に相当。

4. 後処理(post-processing) コンティグ(contigs)やスカッフォールド(scaffolds)を得るステップ。ミスアセンブリの同定も含む。

Jun 25, 201436

ゲノムアセンブル（Linux）
比較的ロングリードの454データ用

非モデル生物やヘテロ接合度の高い生物種用Jun 25, 2014

微生物など小～中規
模ゲノム配列決定用37

教科書p18-22

ゲノムアセンブル（Linux以外）
ゲノムアセンブリ以外にもマッピングなど一通りの解析が可能らしい実験系のヒトはわりと使っている人が多いらしい。アセン
ブルはVelvetだけか?!

RパッケージはありませんJun 25, 201438

教科書p22

アセンブルの評価関連
 様々な試み
 Assemblathon 2 (http://assemblathon.org/) GAGE (http://gage.cbcb.umd.edu/) 自分のゲノムプロジェクトでどの程度のcoverageが必要か？
 アセンブリ結果がどんな感じになるかの見通し どのソフトウェア（とパラメータ）を使うべきか

Jun 25, 2014

最近のアセンブラは大抵GAGEやAssemblathon 2を用いた性能評価結果を示しています39

トランスクリプトームアセンブル

一番よく使われているのはTrinityのようです

Jun 25, 201440

Martin and Wang, Nature Reviews Genet., 12: 671-682, 2011

ゲノム用とトランスクリプトーム用の違い
Sequencing depth (coverage)情報の利用法
 ゲノムの場合（例えば）配列長の10倍読んだデータなら、平均的にゲノムのどの領域も10回程度読まれていると仮定される（10X coverage）
 k-mer出現頻度分布に基づくエラー補正が可能
 多くのアセンブラはcoverage情報をリピート配列の認識に利用


ゲノム
10X

 トランスクリプトーム（RNA-seq）の場合

転写物ごとに大きく異なる：低発現転写物はlow coverage, 高発
現転写物はhigh coverage
PacBioが普及すれば
 アセンブル前の段階でどのk-merがどの転写物由来かはわからないので、k-mer出現頻度の外れ値としてartifactsを除去する トランスクリプトーム用
戦略は（低発現転写物がターゲットの場合には）不可能。ただし はもはや必要なし?!
、low coverageなものはたとえ除去していなくてもアセンブルさ
れにくい。
転写物1
転写物2
転写物3
…


Jun 25, 2014

トランスクリプトーム

41

アセンブルの直観的な理解旧世代シーケンサー（ABI3730など）：～1,000塩基
800塩基程度

一致領域(overlap)大→ 信頼性高いNGS (short-read; Illumina)：～数百塩基
100塩基程度
一致領域(overlap)小
→ 信頼性低いNGS (long-read; PacBio)：～数千塩基
×
×

×

×
×
×

Jun 25, 2014

×
×
×
×

×

×
×
×
×
×

×
×
×
×

×
×

×
×
×
×
×
×
×
×
×

エラーは多いが転
写物配列レベルではアセンブルはほぼ不要なレベル42

（今この瞬間を含む）未来予想図


ゲノム配列決定 小～中規模：Illumina MiSeq、PacBio
 大規模：Illumina HiSeq、PacBio

トランスクリプトーム配列決定
 PacBio (+ Illumina)

発現解析
 Illumina HiSeq

PacBioを用いたトランスクリプトー
ム配列決定論文は既に存在する

Jun 25, 201443

Contents（第3回）


アセンブル(Assembly) 2つのアプローチ(two approaches)

Comparative approach (reference-based assembly; resequencing)：同一生物種または近縁種のゲノム配列を利用 de novo approach：過去に配列決定されたものの中に近縁種がない場合
 アルゴリズム（計算手順）


k-mer解析 ゲノム用、トランスクリプトーム用、雑感

マッピング（QuasRパッケージを利用）
 シミュレーションデータを用いたマッピングの基礎
 リアルデータのマッピング（カイコsmall RNA-seqデータ）
 課題



カウント情報取得

Jun 25, 201444

マッピングの基本的なイメージ

基本的なマッピングプログラム（basic aligner; bowtieなど）を用いた場合
リファレンス配列：ゲノム

あるサンプルのRNA-Seqデータ
mapping

遺伝子1

遺伝子2

count

遺伝子3

遺伝子4

リファレンス配列：トランスクリプトーム

遺伝子1

遺伝子2

遺伝子3

count

遺伝子4

マップされたリードをカウントしたデータ（カウントデータ）がその後の数値解析の基礎情報
Jun 25, 201445

マッピング ＝ 大量高速文字列検索
 マップされる側のリファレンス配列：hoge4.fa

 マップする側のRNA-seqデータ（リードと呼ばれる）：”AGG”

出力ファイル

マッピングプログラムの出力：（どのリードが）リファレンス配列
上のどの位置から転写されたものかという座標情報
Jun 25, 201446

教科書p81

マッピング（準備）

マップされる側のリファレンス配列：ref_genome.fa

コピペで作成Jun 25, 201447

教科書p81

マッピング（準備）

マップされる側のリファレンス配列：ref_genome.fa

chr3とchr5の違いは、2番目と7番目の塩基のみ。
マッピングプログラム
bowtie利用時に、”-m”オプションの違いの把握が可能。
Jun 25, 201448

教科書p83-84

マッピング（準備）

マップする側のRNA-seqデータ：sample_RNAseq1.fa

コピペで作成Jun 25, 201449

教科書p83-84

マッピング（準備）

マップする側のRNA-seqデータ：sample_RNAseq1.fa

許容するミスマッチ数による違いや、マップされるべき場所が完全に把握できるように、リードのdescription行に記述されている
Jun 25, 201450

Langmead et al., Genome Biol., 10: R25, 2009

QuasRパッケージを用いてマッピングBasic alignerの1つであるbowtie (Langmead et al., 2009)を利用
 マッピング時に多くのオプションを指定可能

 “-v”：許容するミスマッチ数を指定するオプション。”-v 0”は、リードがリファレンスに完全一致す

るもののみレポート。”-v 2”は、2塩基ミスマッチまで許容してマップされうる場所を探索。
 “-m”：出力するリード条件を指定するオプション。”-m 1”は、複数個所にマップされるリードを除外して、1か所にのみマップされたリードをレポート。”-m 3”は、合計3か所にマップされるリードまでをレポート。
 “--best --strata”：最も少ないミスマッチ数でマップされるもののみ出力する、という意思表示。
これをつけずに”-v 2 -m 1”などと指定すると、たとえ完全一致(ミスマッチ数0)で1か所にのみマップされるリードがあったとしても、どこか別の場所で1塩基ミスマッチでマップされる個所があれば、マップされうる場所が2か所ということを意味し、そのリードは出力されなくなる。それを防ぐのが主な目的
 ...
デフォルトである程度よきに計らってくれるが...実際の挙動を完全に把握できる状況で様々なオプションを試したい

Jun 25, 201451

教科書p86-89

複数のRNA-seqサンプルを実行できるようにリストファイルとして与える
許容するミスマッチ数は0個（”-v 0”）、1か所にマップされるリードのみ出力（”-m 1”）

Jun 25, 201452

教科書p86-89
入力ファイル中の8リードのうち、マップされたのが5リード、マップされなかったのが3リード。R console画面でなく、QCレ
ポートPDFファイル中にも記述あり。

Jun 25, 201453

QuasRパッケージを用いてマッピング

教科書p86-89

実行後

出力ファイルとして実際に取り扱うのはBAM形式ファイルですJun 25, 201454

マッピング結果の出力ファイル形式


教科書p86-89

ゲノム上のどの位置にどのリードがマッピングされたか（トランスクリプ
トームの場合どの転写物配列上のどの位置にどのリードがマッピング
されたか）を表すファイル形式は複数あります。
 SAM (Sequence Alignment/Map) format


SAMtools (Li et al., Bioinformatics, 25: 2078-2079, 2009)

 BAM (Binary Alignment/Map) format


SAMtools (Li et al., Bioinformatics, 25: 2078-2079, 2009)

 BED (Browser Extensible Data) format
BEDtools (Quinlan et al., Bioinformatics, 26: 841-842, 2010)...
実用上はBAM形式、視覚上はBED形式

Jun 25, 201455

マッピング結果の出力ファイル形式

教科書p86-89

BAM形式ファイル

BED形式ファイル

BEDの最小限の情報は、リードIDを含まない
Jun 25, 201456

マッピングオプションと結果の解釈


教科書p86-89“-m 1 --best --strata -v 0”：0ミスマッチで1か所にのみマップされるリードを出力

マップされなかったのは、
計8リード中3リード
Jun 25, 201457

マッピングオプションと結果の解釈


教科書p86-89“-m 1 --best --strata -v 0”：0ミスマッチで1か所にのみマップされるリードを出力

完全一致でも複数個所にマップされるために落とされた2リード
Jun 25, 201458

マッピングオプションと結果の解釈


教科書p86-89“-m 1 --best --strata -v 0”：0ミスマッチで1か所にのみマップされるリードを出力

1塩基ミスマッチのため落とされたリード
Jun 25, 201459

Nie et al., BMC Genomics, 14: 661, 2013

実データのマッピングを行う


教科書p89-90

カイコゲノムにsmall RNA-seqデータをマッピング

目的：カイコゲノム配列にsmall RNA-seqリードをマップ。アダプター配列除去前後でのマッ
プ率の違いを考察（←これが課題）。hoge –SRP016842フォルダ中に2つともあります。

Jun 25, 201460

教科書p89-90

実データのマッピングを行う


カイコゲノムにsmall RNA-seqデータをマッピング

複数のRNA-seqサンプルを実行できるようにリストファイルとして与える

許容するミスマッチ数は2個（”-v 2”）、1か所にマップされるリードのみ出力（”-m 1”）

Jun 25, 201461

実データのマッピングを行う

教科書p89-90

カイコゲノムファイル

実行後

Jun 25, 2014

ファイルサイズ削減のため、配布したhoge – SRP016842
フォルダ中のファイル群はいくつか除いています62

実データのマッピング結果

教科書p89-90

マッピングに要した時間
は5329.58秒(約90分)

マッピングに用いたプログラムやオ
プション情報入力と出力ファイル情報Jun 25, 201463

実データのマッピング結果

教科書p89-90

アダプター配列除去前後のマッピング結果
QCレポートファイルは実際には1つだけ作成される

Jun 25, 201464

教科書p89-90

実データのマッピング結果
アダプター配列除去前

アダプター配列除去後

おそらくどのマッピングプログラムもこのようなサマリーレポート
ファイルを出力する。上：クオリティ分布、下：塩基組成

塩基組成があたかも同じ種類のものが大量に存在しているように見えるがバグか?!
Jun 25, 201465

アダプター配列除去前

アダプター配列除去後たしかに同じ種類のsmall RNA配列が沢山存在してそう。
念のためsequence logosで確認してみる。
Jun 25, 201466

Sequence logosでも似たような結果。プ
ログラムのバグでないことは確かだろう。

Jun 25, 201467

正しくアダプター配列除去ができていることもわかる

アダプター配列：TGGAATTCTCGGGTGC…

Jun 25, 201468

実データのマッピング結果

アダプター配列除去前

Jun 25, 2014

教科書p89-90

アダプター配列除去後

このsRNA-seqリードは49bp長である。43bp程度以上の比較的長いsRNA
リードの場合、3’側にアダプター配列を含んでいてもその塩基数は短いため、1塩基ミスマッチまで許容するとマップされるということだろう。

69

課題


アダプター配列除去前後のsmall RNA-seqデータを
カイコゲノムにマップし、マップ率を比較する
1.

マッピング前の総リード数を述べよ アダプター配列除去前のSRR609266.fastq.gz：
 アダプター配列除去後のhoge4.fastq.gz：

2.

マッピング後の「マップされたリード数」を述べよ アダプター配列除去前のSRR609266.fastq.gz：
 アダプター配列除去後のhoge4.fastq.gz：

3.

Jun 25, 2014

結果の考察。

70

教科書p90-95

マッピング結果からのカウント情報取得

アノテーション情報を利用する場合
 UCSC Genes, Ensembl Genesなど様々なテーブル名を指定可能

 gene, exon, promoter, junctionなど様々なレベルを指定可能

アノテーション情報がない場合
 マップされたリードの和集合領域を同定したのち、領域ごとのリード数をカウント
 BEDtools (Quinlan et al., 2010)中のmergeBedプログラムを実行して和集合領域

同定後、intersectBedプログラムを実行してリード数をカウントする作業に相当

count

基本的なイメージ
領域1

Jun 25, 2014

2

3

4

71

教科書p90-95

マッピング結果からのカウント情報取得

アノテーション情報を利用する場合
 UCSC Genes, Ensembl Genesなど様々なテーブル名を指定可能

 gene, exon, promoter, junctionなど様々なレベルを指定可能

アノテーション情報がない場合
 マップされたリードの和集合領域を同定したのち、領域ごとのリード数をカウント
 BEDtools (Quinlan et al., 2010)中のmergeBedプログラムを実行して和集合領域

同定後、intersectBedプログラムを実行してリード数をカウントする作業に相当

sample1countsample2

複数サンプルの場合には領域が変わりうるJun 25, 201472

教科書p90-95*_range.txtというカウントデータのファイルが作成される
Jun 25, 201473

教科書p90-95“*.bam”という文字列を”*_range.txt”という文字列に変更したものを出力ファイル名として自動的に生成しているJun 25, 201474

マッピング結果からのカウント情報取得*.bed*_range.txt

カウント数はこちら

Jun 25, 201475

マッピング結果からのカウント情報取得

リストファイル中で指定したサンプル
名がカウントデータ行列の列名となる

Jun 25, 201476

昔よく見かけたカウントデータ取得手段
basic alignerの1つであるBowtieを利用
 最大2塩基ミスマッチまで許容してリファレンス配列の1か所とのみ一致するリード（uniquely mapped reads or unique mapper）数をカウント
 Marioni et al., Genome Res., 18:1509-1517, 2008

 Bullard et al., BMC Bioinformatics, 11:94, 2010
 Risso et al., BMC Bioinformatics, 12:480, 2011
 ReCount (Frazee et al., BMC Bioinformatics, 12:449, 2011)
 …

SpliceMap (Au et al., 2010)などのsplice-aware alignerだと相当時間がかかるという現実的な問題もあるのだろう。講義や講習会では到底無理。
→ ユーザの記憶に残らない → 実際に使われない...

上記情報はshort-readの頃の情報なので既に古いかも…。今はlong-readになっているのでsplice-aware alignerの一種のTophatなどから得られたカウント情報だろう
Jun 25, 201477

定量化：遺伝子レベル ⇔ isoformレベル


全体的な流れとしては遺伝子レベル → isoformレベル


例：新規splice variantの発見（Twine et al., PLoS One, 6: e16266, 2011）

遺伝子セット解析（Gene Ontology解析やパスウェイ解析など）のための基本情報は遺伝子レベルの解像度
 複数エクソン → 遺伝子レベルの要約統計量 exon union method (Mortazavi et al., Nat. Methods, 5: 621-628, 2008)

全てのisoforms間で用いられているexonの情報（union：和集合）を利用

 exon intersection method (Bullard et al., BMC Bioinformatics, 11: 94, 2010)

複数isoforms間で共通して用いられているexonの情報のみ（intersection：積集合）を利用

count情報を得る際に、どのexonの情報を用いるか?

Jun 25, 201478

Garber et al., Nat. Methods, 8: 469-477, 2011のFig. 3c遺伝子のカウント数の定義
算出された生リードカウント結果
exon union method（和集合）の場合：20 reads
Exon intersection method（積集合）の場合：11 reads

様々な思想があり、当然その後の解析結果に影響を及ぼしますJun 25, 201479

教科書p90-95

Unionがデフォルトらしいが、exon-union methodに相当する記述ではないようだ。他にもpairedendやstrand情報など奥が深いのでご注意。

Jun 25, 201480