特論I

Size: px

Start display at page:

Download "特論I"

れれごみぶち
5 years ago
Views:

1 版講義室後ろにある USB メモリ中の hoge フォルダをデスクトップにコピーしておいてください農学生命情報科学特論 I 第 3 回東京大学大学院農学生命科学研究科アグリバイオインフォマティクス教育研究ユニット門田幸二 kadota@iu.a.u-tokyo.ac.jp Jun 25,

講義予定第 1 回 (2014 年 6 月 11 日 ) 西 :NSG 概論現状や展望など講義のみ第 2 回 (2014 年 6 月 18 日 ) 門田 : データベースデータ取得ファイル形式および変換前処理教科書の 1.3 節周辺第 3 回 (2014 年 6 月 25 日 ) 門田 : アセンブルマッピングカウント情報取得教科書の 2.

2 講義予定第 1 回 (2014 年 6 月 11 日 ) 西 :NSG 概論現状や展望など講義のみ第 2 回 (2014 年 6 月 18 日 ) 門田 : データベースデータ取得ファイル形式および変換前処理教科書の 1.3 節周辺第 3 回 (2014 年 6 月 25 日 ) 門田 : アセンブルマッピングカウント情報取得教科書の 2.3 節周辺第 4 回 (2014 年 7 月 2 日 ) 門田 : クラスタリングデータ正規化実験デザイン分布 ( モデル ) 発現変動解析教科書の 3.3 節周辺授業の目標概要次世代シーケンサ (NGS) の普及により以前は主にゲノム解析系で必要とされていた配列解析のためのスキルがトランスクリプトーム解析においても要求される時代になっています本科目では様々な局面で応用可能な配列解析系のスキルアップを目指し RNA シークエンス (RNA-Seq) に基づく ( 非モデル生物の ) トランスクリプトーム解析を題材とした実習を含む講義を行います教科書 Jun 25,

3 エラーの具体例 (2014 年 6 月 13 日 ) 2013 年 11 月 1 日のセミナーで見せた結果エラーの原因はメモリ不足だそうです by 孫堅強氏 (2014 年 6 月 19 日 ) Jun 25,

課題遂行時に何人か遭遇したエラーの解説何人かの方が作業ディレクトリの変更も正しく行い SRR609266.fastq.

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

5 Contents( 第 3 回 ) アセンブル (Assembly) 2 つのアプローチ (two approaches) Comparative approach (reference-based assembly; resequencing): 同一生物種または近縁種のゲノム配列を利用 de novo approach: 過去に配列決定されたものの中に近縁種がない場合アルゴリズム ( 計算手順 ) k-mer 解析ゲノム用トランスクリプトーム用雑感マッピング (QuasR パッケージを利用 ) シミュレーションデータを用いたマッピングの基礎リアルデータのマッピング ( カイコ small RNA-seq データ ) 課題カウント情報取得 Jun 25,

6 ゲノムアセンブル Comparative approach 同一生物種または近縁種のゲノム配列を利用する reference-based assembly Resequencing ともいうヒトゲノム resequencing や SNP 解析系はこちら一個人のヒトゲノム (Wheeler et al., Nature, 452: , 2008) 日本人ゲノム (Fujimoto et al., Nat. Genet., 42: , 2010) ENCODE project (ENCODE Project Consortium et al., Nature, 489: 57-74, 2012) de novo approach 過去に配列決定された生物種以外が主な対象パンダ (Li et al., Nature, 463: , 2008) サンゴ (Shinzato et al., Nature, 476: , 2011) (NGS 由来の比較的短い ) 配列決定されたリードのみから目的生物種のゲノム配列を決めること ( 組み立てること ) Jun 25,

7 Tips リード (read) Sequencer で読んだ塩基配列のことコンティグ (contig) 異なる複数のリードが ACGT の切れ目なく連結されたもの右図では A-D の四つのコンティグ Scaffold (supercontig) N50 コンティグ間の位置関係を表したもの A-D-B-C ではなく A-B-C-D という関係得られた複数のコンティグを最も長いコンティグから順番に連結していったときに combined total length の 50% になったときのコンティグの長さ断片化されたゲノム配列ペアードエンド解析アセンブル A B C D 参考 Jun 25,

8 Tips Coverage( カバレッジ ) ゲノム解読したいときなどに解読するために必要とされる指標となる数値ゲノムサイズ (X) に対する sequencer で読んだ塩基配列長の和のこと一般にこの数値が高いほどよい Sequence depth という表現と実質的に同じような指標アセンブル時に用いる k の値はいくつがいいの? 複数の k の値を試すようです 2013 年ごろから自動的に決めてくれるものが増えたようですアセンブル結果の評価基準は? 参考よくわかりません平均コンティグ長や N50 が論文の表でよく記述されますこのあたりの数値を大きくするだけなら k の値を大きめにすればいいですただゲノムアセンブリの場合には実質的には長ければ長いほどよいという感じみたいですアセンブルプログラムを実行して得られる出力ファイルはどんな感じ? ( 基本的に )multi-fasta 形式のファイルです赤字の部分間違って小さめと書いてしまっていたことに 2016 年 2 月 1 日に気づいたので修正しましたこの種のミスはいけませんね失礼しました m( )m Jun 25, 2014 >contig1 GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT >contig2 ACGATGCAGCCTTAACGA >contig3 8

9 ゲノムアセンブルの手順 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) を得るステップミスアセンブリの同定も含む El-Metwally et al., PLoS Comput Biol., 9: e , 2013 大きく分けて 4 つの手順からなる Jun 25,

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

11 1. 前処理 :k-mer k-mer 頻度解析 NGSデータを入力としてリード長より短いk 連続塩基からなる部分文字列を発生させるのが最初のステップ発生させたk-merの出現頻度情報をもとにカバレッジゲノムサイズ推定コンタミリードの除去などを行う例 1:20 塩基長のNGSリードをk=19で分割すると2 個のk-merを発生可能 CACCAGGACATGAAGACGCG CACCAGGACATGAAGACGC ACCAGGACATGAAGACGCG 例 2:20 塩基長の NGS リードを k=17 で分割すると 4 個の k-mer を発生可能 CACCAGGACATGAAGACGCG CACCAGGACATGAAGAC ACCAGGACATGAAGACG CCAGGACATGAAGACGC CAGGACATGAAGACGCG L 塩基長の NGS リードを k-mer に分割すると (L k + 1) 個の k-mer を発生可能 Jun 25,

12 項目名は変更する可能性あり k=19 として k-mer 頻度分布を作成入力ファイルについて説明します Jun 25,

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

14 1. 前処理 :k-mer 大きく分けて 4 つの手順からなる各リード中の description 部分の記述が部分配列の位置情報に相当 Jun 25,

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

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

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

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

19 1. 前処理 :k-mer sample32_ngs.fasta は 50 塩基からなるランダム塩基配列をリファレンスとしています 20 塩基からなる 10 個の部分配列をリードとしています総塩基数は 200 なのでリファレンス配列の 4 倍の長さつまり 4X coverage です Jun 25,

20 1. 前処理 :k-mer sample36_ngs.fastaは 10,000 塩基からなるランダム塩基配列をリファレンスとしています 80 塩基からなる5,000 個の部分配列をリードとしています総塩基数は400,000なのでリファレンス配列の40 倍の長さつまり40X coverageです Jun 25,

21 (k-mer の種類は問わずに )37 回出現した k-mer が 2 個あったということ出現回数の median は 22 可能な k-mer の種類数は 4^21 個実際の種類数は 9,978 個でゲノムサイズ (=10,000) と酷似 40X の分布シークエンスエラーがない場合に相当 Jun 25,

22 (k-mer の種類は問わずに )31 回出現した k-mer が 97 個あったということ出現回数の median は 22 可能な k-mer の種類数は 4^21 個実際の種類数は 9,976 個でゲノムサイズ (=10,000) と酷似 40Xの分布シークエンスエラーがない場合に相当実行ごとに結果が多少異なりますがどこかで計算をさぼらないといけない現実的な理由があるためかも ( バグかも ) Jun 25,

23 (k-mer の種類は問わずに )59 回出現した k-mer が 19 個あったということ出現回数の median は 48 可能な k-mer の種類数は 4^21 個実際の種類数は 9,978 個でゲノムサイズ (=10,000) と酷似 100Xの分布シークエンスエラーがない場合に相当 coverageの増加 (40X 100X) に伴い出現頻度分布 x 軸方向で右側にシフトしていることがわかる Jun 25,

24 (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,

(=10,000) と酷似 100Xの分布シークエンスエラーがない場合に相当 k-merの値が変わっても

25 (k-mer の種類は問わずに )67 回出現した k-mer が 15 個あったということ出現回数の median は 54 可能な k-mer の種類数は 4^11 個実際の種類数は 9,965 個でゲノムサイズ (=10,000) と酷似 100Xの分布シークエンスエラーがない場合に相当 k-merの値が変わっても (k=31 11) k-merの種類数はゲノムサイズとほぼ同じだが実際にはエラーを多く含むので非現実的 Jun 25,

26 Chikhi and Medvedev, Bioinformatics, 30: 31-37, 2014 の Fig. 1 k-mer 解析は様々な場面で利用されます 1 大まかなゲノムサイズ推定 ( とカバレッジ ) 現実の NGS データはシークエンスエラーやリピート配列を含むため正規分布っぽくはならないが主要なピークの位置情報をもとにゲノムサイズを推定可能 ( らしい ) リード長とリード数から得られた NGS データ中の総塩基数情報はわかっているので得られたゲノムサイズで割ることでカバレッジ (coverage) が分かる k = 41, 51, 61 k = 41, 51, 61 k = 41, 51, 61 S. aureus ゲノムサイズ :2.8 Mb Coverage: 167X リード数 :5,000,000 リード長 :101 bp H. sapiens chr 14 ゲノムサイズ :88 Mb Coverage: 70X リード数 :62,000,000 リード長 :101 bp B. impatiens ゲノムサイズ :250 Mb Coverage: 247X リード数 :497,000,000 リード長 :124 bp Jun 25,

27 Chikhi and Medvedev, Bioinformatics, 30: 31-37, 2014 の Fig. 1 k-mer 解析は様々な場面で利用されます 2 前処理 ( フィルタリング ) シークエンスエラーやリピート配列由来リードの除去 k = 41, 51, 61 ゲノム 4X 20X シークエンスエラーに由来する数回しか出現しないk-mer 由来リードは除去対象リピート配列に由来する異常に多く出現するk-mer 由来リードも除去対象 Jun 25,

28 k-mer 解析は様々な場面で利用されます 3 前処理 ( エラー補正 ) ヘテロ接合度の高い 2 倍体ゲノム (highly heterozygous diploid genomes) 由来リードの補正 Kajitani et al., Genome Res., in press の Fig. 2 父親由来ゲノム CACCAGGACATGAAGACGCGTTCA CACCAGGACATCAAGACGCGTTCA 母親由来ゲノムシークエンスエラーに由来する数回しか出現しない k-mer 由来リードは除去対象リピート配列に由来する異常に多く出現する k-mer 由来リードも除去対象 CACCAGGACAT ACCAGGACATG CCAGGACATGA 父親由来 k-mer GACATGAAGAC TGAAGACGCGT GAAGACGCGTC AAGACGCGTCA 両親間で配列の異なる領域を含むk- merの出現頻度は配列が同じ領域由来 k-merの出現頻度 (c) の1/2となるアセンブルに悪影響を及ぼすため母親由来ゲノム配列に揃えるなどして補正 Jun 25,

29 ゲノムアセンブルの手順 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) を得るステップミスアセンブリの同定も含む El-Metwally et al., PLoS Comput Biol., 9: e , 2013 大きく分けて 4 つの手順からなる Jun 25,

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

31 El-Metwally et al., PLoS Comput Biol., 9: e , グラフ構築 (graph construction) k-mer アプローチ (de Bruijn グラフ ) リードを全ての可能な k-mer に分割し有向グラフを作成 (k=9 の例 ) リード 1: CACCAGGACATGAAGACGCG リード 2: CCAGGACATGAAGACGCGTT CACCAGGAC ACCAGGACA CCAGGACAT ATGAAGACG TGAAGACGC GAAGACGCG CCAGGACAT CAGGACATG AGGACATGA GAAGACGCG AAGACGCGT AGACGCGTT CACCAGGAC CCAGGACAT GAAGACGCG CAGGACATG GAAGACGCG AGACGCGTT ACCAGGACA CCAGGACAT : ノード (node; 頂点 ) : エッジ (edge; 辺 ) AAGACGCGT Jun 25,

32 El-Metwally et al., PLoS Comput Biol., 9: e , グラフ構築 (graph construction) k-mer アプローチ (de Bruijn グラフ ) 同一ノードをマージして de Bruijn グラフを作成 (k=9) リード 1: CACCAGGACATGAAGACGCG リード2: CCAGGACATGAAGACGCGTT CACCAGGAC CCAGGACAT GAAGACGCG CAGGACATG GAAGACGCG AGACGCGTT ACCAGGACA CCAGGACAT AAGACGCGT CACCAGGAC CCAGGACAT GAAGACGCG AGACGCGTT AAGACGCGT ACCAGGACA CAGGACATG Jun 25,

33 El-Metwally et al., PLoS Comput Biol., 9: e , 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,

34 El-Metwally et al., PLoS Comput Biol., 9: e , グラフ簡易化 (graph simplification) 連続したノード (nodes; 頂点 ) やバブルのマージリード 1: CACCAGGACATGAAGACGCG リード 2: CCAGGACATGAAGACGCGTT CACCAGGAC CCAGGACAT GAAGACGCG AGACGCGTT ACCAGGACA CAGGACATG AAGACGCGT CACCAGGACATGAAGACGCGTT この 2 つのリードだけで簡易化した結果 Jun 25,

35 Kajitani et al., Genome Res., in press の Suppl. Fig. 8 に相当 3. グラフ簡易化 (graph simplification) 連続したノード (nodes; 頂点 ) やバブルのマージ父親由来ゲノム CACCAGGACATGAAGACGCGTTCA CACCAGGACATCAAGACGCGTTCA 母親由来ゲノム父親由来ゲノムリード1: CACCAGGACATGAAGACGCG リード2: CCAGGACATCAAGACGCGTT CAGGACATGAAGACGCG CACCAGGACAT AAGACGCGTT CAGGACATCAAGACGCG SNP など塩基に違いがあればバブル構造になります Jun 25,

El-Metwally et al., PLoS Comput Biol., 9: e1003345, 2013 ゲノムアセンブルの手順 1.

36 El-Metwally et al., PLoS Comput Biol., 9: e , 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,

37 ゲノムアセンブル (Linux) 比較的ロングリードの 454 データ用非モデル生物やヘテロ接合度の高い生物種用微生物など小 ~ 中規模ゲノム配列決定用 Jun 25,

38 教科書 p18-22 ゲノムアセンブル (Linux 以外 ) ゲノムアセンブリ以外にもマッピングなど一通りの解析が可能らしい実験系のヒトはわりと使っている人が多いらしいアセンブルは Velvet だけか?! R パッケージはありません Jun 25,

アセンブルの評価関連様々な試み Assemblathon 2 (http://assemblathon.

39 アセンブルの評価関連様々な試み Assemblathon 2 ( GAGE ( 自分のゲノムプロジェクトでどの程度の coverage が必要か? アセンブリ結果がどんな感じになるかの見通しどのソフトウェア ( とパラメータ ) を使うべきか教科書 p22 最近のアセンブラは大抵 GAGE や Assemblathon 2 を用いた性能評価結果を示しています Jun 25,

40 トランスクリプトームアセンブル一番よく使われているのは Trinity のようです Jun 25,

41 ゲノム用とトランスクリプトーム用の違い Sequencing depth (coverage) 情報の利用法ゲノムの場合 ( 例えば ) 配列長の 10 倍読んだデータなら平均的にゲノムのどの領域も 10 回程度読まれていると仮定される (10X coverage) k-mer 出現頻度分布に基づくエラー補正が可能多くのアセンブラは coverage 情報をリピート配列の認識に利用トランスクリプトーム (RNA-seq) の場合 Martin and Wang, Nature Reviews Genet., 12: , 2011 転写物ごとに大きく異なる : 低発現転写物は low coverage, 高発現転写物は high coverage アセンブル前の段階でどのk-merがどの転写物由来かはわからないので k-mer 出現頻度の外れ値としてartifactsを除去する戦略は ( 低発現転写物がターゲットの場合には ) 不可能ただし low coverageなものはたとえ除去していなくてもアセンブルされにくい転写物 1 転写物 2 転写物 3 ゲノム 10X PacBio が普及すればトランスクリプトーム用はもはや必要なし?! Jun 25, 2014 トランスクリプトーム 41

42 アセンブルの直観的な理解旧世代シーケンサー (ABI3730 など ):~1,000 塩基 800 塩基程度一致領域 (overlap) 大信頼性高い NGS (short-read; Illumina):~ 数百塩基 100 塩基程度一致領域 (overlap) 小信頼性低い NGS (long-read; PacBio):~ 数千塩基 Jun 25, 2014 エラーは多いが転写物配列レベルではアセンブルはほぼ不要なレベル 42

43 ( 今この瞬間を含む ) 未来予想図ゲノム配列決定小 ~ 中規模 :Illumina MiSeq PacBio 大規模 :Illumina HiSeq PacBio トランスクリプトーム配列決定 PacBio (+ Illumina) 発現解析 Illumina HiSeq PacBio を用いたトランスクリプトーム配列決定論文は既に存在する Jun 25,

44 Contents( 第 3 回 ) アセンブル (Assembly) 2 つのアプローチ (two approaches) Comparative approach (reference-based assembly; resequencing): 同一生物種または近縁種のゲノム配列を利用 de novo approach: 過去に配列決定されたものの中に近縁種がない場合アルゴリズム ( 計算手順 ) k-mer 解析ゲノム用トランスクリプトーム用雑感マッピング (QuasR パッケージを利用 ) シミュレーションデータを用いたマッピングの基礎リアルデータのマッピング ( カイコ small RNA-seq データ ) 課題カウント情報取得 Jun 25,

45 マッピングの基本的なイメージ基本的なマッピングプログラム (basic aligner; bowtie など ) を用いた場合リファレンス配列 : ゲノム count あるサンプルの RNA-Seq データ mapping 遺伝子 1 遺伝子 2 遺伝子 3 遺伝子 4 リファレンス配列 : トランスクリプトーム count 遺伝子 1 遺伝子 2 遺伝子 3 遺伝子 4 Jun 25, 2014 マップされたリードをカウントしたデータ ( カウントデータ ) がその後の数値解析の基礎情報 45

46 マッピング = 大量高速文字列検索マップされる側のリファレンス配列 :hoge4.fa マップする側の RNA-seq データ ( リードと呼ばれる ): AGG 出力ファイル Jun 25, 2014 マッピングプログラムの出力 :( どのリードが ) リファレンス配列上のどの位置から転写されたものかという座標情報 46

47 マッピング ( 準備 ) マップされる側のリファレンス配列 :ref_genome.fa 教科書 p81 コピペで作成 Jun 25,

48 マッピング ( 準備 ) マップされる側のリファレンス配列 :ref_genome.fa 教科書 p81 chr3 と chr5 の違いは 2 番目と 7 番目の塩基のみマッピングプログラム bowtie 利用時に -m オプションの違いの把握が可能 Jun 25,

49 マッピング ( 準備 ) マップする側の RNA-seq データ :sample_rnaseq1.fa 教科書 p83-84 コピペで作成 Jun 25,

50 マッピング ( 準備 ) マップする側の RNA-seq データ :sample_rnaseq1.fa 教科書 p83-84 許容するミスマッチ数による違いやマップされるべき場所が完全に把握できるようにリードの description 行に記述されている Jun 25,

51 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 か所ということを意味しそのリードは出力されなくなるそれを防ぐのが主な目的... Langmead et al., Genome Biol., 10: R25, 2009 デフォルトである程度よきに計らってくれるが... 実際の挙動を完全に把握できる状況で様々なオプションを試したい Jun 25,

52 教科書 p86-89 複数の RNA-seq サンプルを実行できるようにリストファイルとして与える許容するミスマッチ数は 0 個 ( -v 0 ) 1 か所にマップされるリードのみ出力 ( -m 1 ) Jun 25,

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

54 QuasR パッケージを用いてマッピング教科書 p86-89 実行後 Jun 25, 2014 出力ファイルとして実際に取り扱うのは BAM 形式ファイルです 54

55 マッピング結果の出力ファイル形式教科書 p86-89 ゲノム上のどの位置にどのリードがマッピングされたか ( トランスクリプトームの場合どの転写物配列上のどの位置にどのリードがマッピングされたか ) を表すファイル形式は複数あります SAM (Sequence Alignment/Map) format SAMtools (Li et al., Bioinformatics, 25: , 2009) BAM (Binary Alignment/Map) format SAMtools (Li et al., Bioinformatics, 25: , 2009) BED (Browser Extensible Data) format... BEDtools (Quinlan et al., Bioinformatics, 26: , 2010) 実用上は BAM 形式視覚上は BED 形式 Jun 25,

56 マッピング結果の出力ファイル形式 BAM 形式ファイル教科書 p86-89 BED 形式ファイル Jun 25, 2014 BED の最小限の情報はリード ID を含まない 56

57 マッピングオプションと結果の解釈 -m 1 --best --strata -v 0 :0 ミスマッチで 1 か所にのみマップされるリードを出力教科書 p86-89 マップされなかったのは計 8 リード中 3 リード Jun 25,

58 マッピングオプションと結果の解釈 -m 1 --best --strata -v 0 :0 ミスマッチで 1 か所にのみマップされるリードを出力教科書 p86-89 完全一致でも複数個所にマップされるために落とされた 2 リード Jun 25,

59 マッピングオプションと結果の解釈 -m 1 --best --strata -v 0 :0 ミスマッチで 1 か所にのみマップされるリードを出力教科書 p 塩基ミスマッチのため落とされたリード Jun 25,

60 Nie et al., BMC Genomics, 14: 661, 2013 実データのマッピングを行うカイコゲノムに small RNA-seq データをマッピング教科書 p89-90 目的 : カイコゲノム配列に small RNA-seq リードをマップアダプター配列除去前後でのマップ率の違いを考察 ( これが課題 ) hoge SRP フォルダ中に 2 つともあります Jun 25,

61 実データのマッピングを行うカイコゲノムに small RNA-seq データをマッピング教科書 p89-90 複数の RNA-seq サンプルを実行できるようにリストファイルとして与える許容するミスマッチ数は 2 個 ( -v 2 ) 1 か所にマップされるリードのみ出力 ( -m 1 ) Jun 25,

62 実データのマッピングを行う教科書 p89-90 カイコゲノムファイル実行後 Jun 25, 2014 ファイルサイズ削減のため配布した hoge SRP フォルダ中のファイル群はいくつか除いています 62

実データのマッピング結果教科書 p89-90 マッピングに要した時間は 5329.

63 実データのマッピング結果教科書 p89-90 マッピングに要した時間は秒 ( 約 90 分 ) マッピングに用いたプログラムやオプション情報入力と出力ファイル情報 Jun 25,

64 実データのマッピング結果教科書 p89-90 アダプター配列除去前後のマッピング結果 QC レポートファイルは実際には 1 つだけ作成される Jun 25,

65 実データのマッピング結果教科書 p89-90 アダプター配列除去前アダプター配列除去後おそらくどのマッピングプログラムもこのようなサマリーレポートファイルを出力する上 : クオリティ分布下 : 塩基組成 Jun 25, 2014 塩基組成があたかも同じ種類のものが大量に存在しているように見えるがバグか?! 65

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

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

68 正しくアダプター配列除去ができていることもわかるアダプター配列 :TGGAATTCTCGGGTGC Jun 25,

69 実データのマッピング結果教科書 p89-90 アダプター配列除去前アダプター配列除去後 Jun 25, 2014 この srna-seq リードは 49bp 長である 43bp 程度以上の比較的長い srna リードの場合 3 側にアダプター配列を含んでいてもその塩基数は短いため 1 塩基ミスマッチまで許容するとマップされるということだろう 69

70 課題アダプター配列除去前後の small RNA-seq データをカイコゲノムにマップしマップ率を比較する 1. マッピング前の総リード数を述べよアダプター配列除去前の SRR fastq.gz: アダプター配列除去後の hoge4.fastq.gz: 2. マッピング後のマップされたリード数を述べよアダプター配列除去前の SRR fastq.gz: アダプター配列除去後の hoge4.fastq.gz: 3. 結果の考察 Jun 25,

マッピング結果からのカウント情報取得アノテーション情報を利用する場合 UCSC Genes, Ensembl Genes など様々なテーブル名を指定可能 gene, exon, promoter, junction など様々なレベルを指定可能アノテーション情報がない場合教科書 p90-95

71 マッピング結果からのカウント情報取得アノテーション情報を利用する場合 UCSC Genes, Ensembl Genes など様々なテーブル名を指定可能 gene, exon, promoter, junction など様々なレベルを指定可能アノテーション情報がない場合教科書 p90-95 マップされたリードの和集合領域を同定したのち領域ごとのリード数をカウント BEDtools (Quinlan et al., 2010) 中の mergebed プログラムを実行して和集合領域同定後 intersectbed プログラムを実行してリード数をカウントする作業に相当領域 count 基本的なイメージ Jun 25,

72 マッピング結果からのカウント情報取得アノテーション情報を利用する場合 UCSC Genes, Ensembl Genes など様々なテーブル名を指定可能 gene, exon, promoter, junction など様々なレベルを指定可能アノテーション情報がない場合教科書 p90-95 マップされたリードの和集合領域を同定したのち領域ごとのリード数をカウント BEDtools (Quinlan et al., 2010) 中の mergebed プログラムを実行して和集合領域同定後 intersectbed プログラムを実行してリード数をカウントする作業に相当 sample1 count sample2 複数サンプルの場合には領域が変わりうる Jun 25,

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

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

75 マッピング結果からのカウント情報取得 *.bed *_range.txt カウント数はこちら Jun 25,

76 マッピング結果からのカウント情報取得リストファイル中で指定したサンプル名がカウントデータ行列の列名となる Jun 25,

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

78 定量化 : 遺伝子レベル isoform レベル全体的な流れとしては遺伝子レベル isoform レベル例 : 新規 splice variant の発見 (Twine et al., PLoS One, 6: e16266, 2011) 遺伝子セット解析 (Gene Ontology 解析やパスウェイ解析など ) のための基本情報は遺伝子レベルの解像度複数エクソン遺伝子レベルの要約統計量 exon union method (Mortazavi et al., Nat. Methods, 5: , 2008) 全ての isoforms 間で用いられている exon の情報 (union: 和集合 ) を利用 exon intersection method (Bullard et al., BMC Bioinformatics, 11: 94, 2010) 複数 isoforms 間で共通して用いられている exon の情報のみ (intersection: 積集合 ) を利用 count 情報を得る際にどの exon の情報を用いるか? Jun 25,

79 遺伝子のカウント数の定義算出された生リードカウント結果 Garber et al., Nat. Methods, 8: , 2011 の Fig. 3c exon union method( 和集合 ) の場合 :20 reads Exon intersection method( 積集合 ) の場合 :11 reads 様々な思想があり当然その後の解析結果に影響を及ぼします Jun 25,

80 教科書 p90-95 Union がデフォルトらしいが exon-union method に相当する記述ではないようだ他にも pairedend や strand 情報など奥が深いのでご注意 Jun 25,

農学生命情報科学特論I

農学生命情報科学特論I 2015.07.01 版 USB メモリ中の hoge フォルダをデスクトップにコピーしておいてください前回 (6/23) の hoge フォルダがデスクトップに残っているかもしれないのでご注意ください農学生命情報科学特論 I 第 3 回大学院農学生命科学研究科アグリバイオインフォマティクス教育研究プログラム門田幸二 ( かどたこうじ ) kadota@iu.a.u-tokyo.ac.jp