農学生命情報科学特論I

Transcription

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

2 講義予定 第 1 回 (2015 年 6 月 16 日 ) NGS の普及により 以前は主にゲノム解析系で必要とされていた配列解析のためのスキルがトランスクリプトーム解析においても要求される時代になっています 本科目では 様々な局面で応用可能な配列解析系のスキルアップを目指し RNA-Seq に基づくトランスクリプトーム解析を題材とした講義を行います データベース データ取得 ファイル形式 Quality Control 教科書の 1.3 節周辺 第 2 回 (2015 年 6 月 23 日 ) Quality Control k-mer 解析 トリミング ( アダプター配列除去 ) 第 3 回 (2015 年 6 月 30 日 ) フィルタリング アセンブル マッピング カウント情報取得 教科書の 2.3 節周辺 第 4 回 (2015 年 7 月 7 日 ) クラスタリング 実験デザイン 分布 ( モデル ) 発現変動解析 教科書の 3.3 節周辺 教科書 2

3 Contents 先週の課題やエラーについて 課題 1 の解説 アダプター配列除去 (QuasR) 実行時のエラーは門田のミス Paired-end データの取り扱い フィルタリング (filtering) 基本形 : リード数が同じ paired-end の場合 (QuasR パッケージ ) 発展形 : リード数が異なる場合 (ShortRead パッケージ ) アセンブル (assemble) ゲノム用とトランスクリプトーム用 Rockhopper2( バクテリアのトランスクリプトーム用 ) マッピング (mapping) 基礎 オプション 仮想データ説明 マッピング本番 出力ファイル形式 オプションと結果の関係 カウント情報取得 実データ解析 QuasR でマッピング トリミング 高精度なアセンブルとマッピングへ 3

4 k-mer 解析 上級 1 のコードをテンプレートにして Position 7 で特異的に出現している AAGAGCA を R で確認してみよう Head 関数出力結果は最初の 6 個 ( つまり position 1-6) までしか表示させていないので ここでは n=8 として position 1-8 まで表示させるようにしている たしかに position 7 で期待値 (Exp = 90.09) の 倍の出現回数になっていることが分かる 4

5 課題 1 上級 1 のコードをテンプレートにして 1 AGAGCAC 2GAGCACA 3 長さを含めて任意 の k-mer 解析を行い 得られた結果を示せ 1 5

6 課題 1 の解説 1AGAGCAC の k-mer 解析時にエラーが出たと思います この原因は AGAGCAC が 1 つも存在しないポジションがあった場合に そこが NA となるからです mean 関数実行時に要素中に NA が 1 つでもあると計算できないという罠があったのです 6

7 課題 1 の解説 Obs と打ち込んで 何が起こっているのかを悟る 画面に見えているだけでも AGAGCAC が存在しない position が少なくとも 2 ヶ所あることがわかる 7

8 課題 1 の解説 最もシンプルな対処法は mean 関数中の na.rm オプションを TRUE ( デフォルトは FALSE) にするやり方 NA の要素を読み飛ばすオプションなので 平均値の計算が若干不正確になるが AGAGCAC が position 8 で相対的に高頻度に出現しているということを調べるという点では全く問題ない 茂木朋貴氏他数名提供情報 8

9 課題 1 の解説 平均値 ( この場合 ) さえ計算できれば Obs/Exp を計算可能 つまり プロファイルを描画することができる NA となっている場所が 2 ヶ所あることがわかる 9

10 課題 1 の解説 もっとよりよい方法は 1NA となった要素に 0 を入れること is.na 関数は 要素が NA のときに TRUE そうでないときに FALSE を返す 野間口達洋氏他数名提供情報 10

11 課題 1 の解説 0 代入前後で出現頻度の平均値が から に少し下がっているのがわかります 11

12 課題 1 の解説 こんな感じになります 2015 年 6 月 23 日の課題を該当部分のみ差し替えたいヒトは 回収箱に入れてかまいません 課題用紙は後ろにあります 先週提出分と突き合わせて評価しますので 考察部分など差し替える必要がない部分は空欄で構いません 1 12

13 バージョンアップ ひっそりと修正しています レポート中で問題点の指摘や解決法を示して頂いた方々に感謝 m( )m

14 Contents 先週の課題やエラーについて 課題 1 の解説 アダプター配列除去 (QuasR) 実行時のエラーは門田のミス Paired-end データの取り扱い フィルタリング (filtering) 基本形 : リード数が同じ paired-end の場合 (QuasR パッケージ ) 発展形 : リード数が異なる場合 (ShortRead パッケージ ) アセンブル (assemble) ゲノム用とトランスクリプトーム用 Rockhopper2( バクテリアのトランスクリプトーム用 ) マッピング (mapping) 基礎 オプション 仮想データ説明 マッピング本番 出力ファイル形式 オプションと結果の関係 カウント情報取得 実データ解析 QuasR でマッピング トリミング 高精度なアセンブルとマッピングへ 14

15 アダプター配列除去 1 param_nrec のところをいくら変えてもエラーが出続ける問題に多くの受講生が遭遇しました 理由はほぼ間違いなく 単純に私が preprocessreads 関数中に明示的に param_nrec を与えていなかったからです 結果として それが適切に反映されずにデフォルトの 1,000,000 がいつまでも使われていたということでしょう 2 中村浩正氏提供情報 15

16 アダプター配列除去 休み時間にでも動作確認してみてください きっとうまくいくことでしょう

17 Paired-end の取扱い 2015 年 6 月 23 日の講義資料

18 Paired-end の取扱い QuasR (ver ) では まだ pairedend データのアダプター配列除去には対応できていないようだといいましたが その後すぐに開発者が回避策 (workaround) を教えてくれました Thanks to Dr. Stadler 18

19 回避策 基本戦略は single-end として別々にアダプター配列除去を行うが 出力ファイルのリード数が変わらないようにする です 乳酸菌 paired-end RNAseq データ (SRR626268) のアダプター配列除去を行う一連の手順の基本形は 例題 3-5 です 例題 3 では Forward 側のリードファイルを入力として FastQC でレポートされた TruSeq Adapter, Index 3 の除去を行っています 約 1 分 19

20 回避策 例題 4 では 例題 3 の出力ファイルを入力として もう 1 つレポートされていた TruSeq Adapter, Index 2 の除去を行っています FastQC 実行環境にない場合でも 用いた実験プロトコルが分かれば 候補となる adapter や primer 配列の一部で文字列検索すればある程度わかります もちろん この場合は Index 配列部分を含めたりしたほうがいいとは思います 約 1 分 20

21 回避策 例題 5 では Reverse 側のリードファイルを入力として FastQC でレポートされた Illumina Single End PCR Primer 1 の除去を行っています Primer 配列の場合は 逆相補鎖 ( reverse complement) を作成してから与えないとうまく除去できなかったことから このようにしています FastQC を実行できる環境にあれば overrepresented sequences 項目の possible source のところが No Hit になったかどうかでうまくいったかどうかの判定ができます 約 1 分 21

22 Contents 先週の課題やエラーについて 課題 1 の解説 アダプター配列除去 (QuasR) 実行時のエラーは門田のミス Paired-end データの取り扱い フィルタリング (filtering) 基本形 : リード数が同じ paired-end の場合 (QuasR パッケージ ) 発展形 : リード数が異なる場合 (ShortRead パッケージ ) アセンブル (assemble) ゲノム用とトランスクリプトーム用 Rockhopper2( バクテリアのトランスクリプトーム用 ) マッピング (mapping) 基礎 オプション 仮想データ説明 マッピング本番 出力ファイル形式 オプションと結果の関係 カウント情報取得 実データ解析 QuasR でマッピング トリミング 高精度なアセンブルとマッピングへ 22

23 フィルタリング : 基本形 リード数の揃った paired-end データを入力として 許容する N の数と最短配列長でフィルタリングするやり方を示します

24 フィルタリング : 基本形 この例題で用いているファイルは hoge4.fastq.gz と hoge5.fastq.gz でも構いません 中身は同じです 24

25 フィルタリング : 基本形 リード長が 18 塩基未満を捨て 1 塩基だけ N を許容するという閾値でフィルタリングをした結果 998,208 リードが残ったということ 約 2 分 25

26 フィルタリング : 発展形 基本戦略は paired-end を singleend として独立に取り扱い 最後に共通リードを抽出 例えば 1~3 あたりの任意のフィルタリングを singleend として独立に実行 最後に 4 でリード数の異なる paired-end のファイルを入力として 共通リードを出力 ここでは 1 の例題 6-8 を実行して得られたリード数の異なる paired-end ファイルを入力として 4 を実行するやり方を示します 個人見解としては 4 の入力ファイルの段階では FASTQ ではなく FASTA 形式でいいと思います 一連のフィルタリングでそこそこのクオリティが保証されているので わざわざクオリティ情報を保持しておく必要はないだろうという思想 26

27 フィルタリング : 発展形 1 参考 Small RNA-seq データのときはリードの右側にあるアダプター配列除去でしたが 乳酸菌 RNA-seq データのときは左側にアダプターがついているので Rpattern ではなく Lpattern になっている点に注意 2 27

28 フィルタリング : 発展形 参考 出力予定ファイルが存在する状態で実行するとエラーが出るので注意 28

29 フィルタリング : 発展形 参考 1 こんな感じになります 2 ここで見えている出力ファイルのサイズは 既に存在していた hoge5.fastq.gz のファイルサイズであり ここでのコピペ実行結果で得られたものではない点に注意

30 フィルタリング : 発展形 参考 hoge5.fastq.gz を hoge フォルダから削除して再度実行 30

31 フィルタリング : 発展形 参考 得られた hoge5.fastq.gz を入力として 例題 7 を実行 31

32 フィルタリング : 発展形 参考 例題 8 の paired-end reverse 側も実行 出力ファイルのサイズは 63,442,904 bytes 999,136 リード 32

33 乳酸菌 RNA-seq データ 100 万リードのオリジナル?! ファイル SRR616268sub_1.fastq.gz:107 bp, 74,906,576 bytes SRR616268sub_2.fastq.gz:93 bp, 67,158,462 bytes 100 万リードのアダプター除去後のファイル SRR616268sub_trim2_1.fastq.gz(hoge4.fastq.gz): 71,343,605 bytes SRR616268sub_trim2_2.fastq.gz : 63,524,791 bytes 1 リード数の異なる paired-end のファイルを入力として 共通リードを出力するやり方を示します 107 bp 93 bp hoge2_1.fastq.gz: 998,208 リード 71,194,586 bytes hoge2_2.fastq.gz: 998,208 リード 63,375,473 bytes アダプター除去およびフィルタリング後のファイル SRR616268sub_trim_1.fastq.gz(hoge7.fastq.gz): 998,658 リード 71,227,695 bytes SRR616268sub_trim_2.fastq.gz(hoge8.fastq.gz): 999,136 リード 63,442,904 bytes 1 33

34 Tips:file.info 関数 100 万リードのオリジナル?! ファイル SRR616268sub_1.fastq.gz:107 bp, 74,906,576 bytes SRR616268sub_2.fastq.gz:93 bp, 67,158,462 bytes 100 万リードのアダプター除去後のファイル SRR616268sub_trim2_1.fastq.gz(hoge4.fastq.gz): 71,343,605 bytes SRR616268sub_trim2_2.fastq.gz : 63,524,791 bytes file.info 関数でもファイル名 サイズなど様々な情報をみることができます 107 bp 93 bp hoge2_1.fastq.gz: 998,208 リード 71,194,586 bytes hoge2_2.fastq.gz: 998,208 リード 63,375,473 bytes アダプター除去およびフィルタリング後のファイル SRR616268sub_trim_1.fastq.gz(hoge7.fastq.gz): 998,658 リード 71,227,695 bytes SRR616268sub_trim_2.fastq.gz(hoge8.fastq.gz): 999,136 リード 63,442,904 bytes 1 34

35 フィルタリング : 発展形 1 入力はリード数の異なる pairedend ファイル 出力はリード数が同じ paired-end ファイルです description 行の記述形式次第ではエラーが出ると思いますが 適宜変更してください 35

36 フィルタリング : 発展形 最初の 4 リード分の description 部分を表示 ここで見えているもので判断すると SRR は 2 つの入力ファイル中に存在する そして 少なくとも forward 側のファイル (*_1.fastq.gz) 中の 3 リード (SRR ) は reverse 側のファイル (*_2.fastq.gz) 中には存在しないことがわかる 36

37 フィルタリング : 発展形 共通リード抽出用の基本情報として リード ID の積集合 (intersection) を得るのが基本戦略 しかし description 部分の情報をそのまま使うと length=107 と length=93 という部分文字列の違いにより 1 同一リードを同一文字列として認識できないという問題がある 1 37

38 フィルタリング : 発展形 1 2 description 部分を眺め 部分文字列で同一リードを認識できる箇所を探す ここでは 1 description 部分を区切り文字 スペース ( ) で 分断し 21 番目の要素を取り出せば 赤下線部分のみの部分文字列で評価できるという戦略のもとでプログラムを組んでいる 38

39 フィルタリング : 発展形 Forward 側 うまくリード ID (SRR ) のみにできていることがわかります 39

40 フィルタリング : 発展形 Reverse 側 うまくリード ID (SRR ) のみにできていることがわかります 40

41 参考 Tips: コピペで使いまわす わざわざ別のオブジェクト名に置き換えているので 一見ややこしいですが 黒枠内のコードを使いまわせるので 間違いを減らせるというメリットがあります 41

42 フィルタリング : 発展形 得られた共通リード数 (998,428) は forward 側 (998,658) および reverse 側 (999,136) より少ないので妥当 42

43 フィルタリング : 発展形 fastq1 オブジェクトの 1, 13, 17, 22 番目の要素が共通 43

44 フィルタリング : 発展形 is.element 関数は 右側の集合の中の要素が存在する位置を TRUE それ以外を FALSE として返す 共通である fastq1 オブジェクトの 1, 13, 17, 22 番目の要素が TRUE になっていることがわかる 44

45 フィルタリング : 発展形 共通リードのみなので リード数が揃っていることがわかる 45

46 乳酸菌 RNA-seq データ 100 万リードのオリジナル?! ファイル SRR616268sub_1.fastq.gz:107 bp, 74,906,576 bytes SRR616268sub_2.fastq.gz:93 bp, 67,158,462 bytes 100 万リードのアダプター除去後のファイル SRR616268sub_trim2_1.fastq.gz(hoge4.fastq.gz): 71,343,605 bytes SRR616268sub_trim2_2.fastq.gz : 63,524,791 bytes hoge2_1.fastq.gz: 998,208リード 71,194,586 bytes 2 hoge2_2.fastq.gz: 998,208リード 63,375,473 bytes アダプター除去およびフィルタリング後のファイル SRR616268sub_trim_1.fastq.gz(hoge7.fastq.gz): 998,658 リード 71,227,695 bytes SRR616268sub_trim_2.fastq.gz(hoge8.fastq.gz): 999,136 リード 63,442,904 bytes hoge1_1.fastq.gz: 998,428 リード 63,446,240 bytes hoge1_2.fastq.gz: 998,428 リード 56,019,410 bytes 1はデフォルトの 許容するN 数が2 最短配列長が14 のときの結果 2は若干厳しめの 許容するN 数が1 最短配列長が18 のときの結果なので 2のほうがリード数が若干少ない 107 bp 1 93 bp 46

47 Contents 先週の課題やエラーについて 課題 1 の解説 アダプター配列除去 (QuasR) 実行時のエラーは門田のミス Paired-end データの取り扱い フィルタリング (filtering) 基本形 : リード数が同じ paired-end の場合 (QuasR パッケージ ) 発展形 : リード数が異なる場合 (ShortRead パッケージ ) アセンブル (assemble) ゲノム用とトランスクリプトーム用 Rockhopper2( バクテリアのトランスクリプトーム用 ) マッピング (mapping) 基礎 オプション 仮想データ説明 マッピング本番 出力ファイル形式 オプションと結果の関係 カウント情報取得 実データ解析 QuasR でマッピング トリミング 高精度なアセンブルとマッピングへ 47

48 アセンブル : ゲノム用 1 ゲノム用 2 非モデル生物やヘテロ接合度の高い生物種用 3 微生物など小 ~ 中規模ゲノム配列決定用 アセンブルのアルゴリズム周辺は で web ページ内を検索

49 アセンブル : 転写物用 転写物 ( トランスクリプトーム ) 用

50 Rockhopper 1Windows のヒトは.exe ファイル Macintosh のヒトは.dmg ファイルで Rockhopper を起動 hoge - Rockhopper フォルダ中に Rockhopper_Results というフォルダがないヒトは ダブルクリックで起動したときに自動的に作成されます 3 最初 An error occurred while opening the file 的なポップアップが出現しますが 30 秒ほど待つと 4 Rockhopper の GUI が起動します

51 Rockhopper 1( 見づらいが )DE NOVO と赤字で書いている部分をクリック 2 入力ファイルを指定すべく BROWSE ボタンを押す

52 Rockhopper 1 ファイルの場所を指定 2Desktop 3hoge

53 乳酸菌 RNA-seq データ 100 万リードのオリジナル?! ファイル SRR616268sub_1.fastq.gz:107 bp, 74,906,576 bytes SRR616268sub_2.fastq.gz:93 bp, 67,158,462 bytes 100 万リードのアダプター除去後のファイル オリジナルの 100 万リードからなる生データファイルで実行 ( の発展課題と同じ ) 107 bp SRR616268sub_trim2_1.fastq.gz(hoge4.fastq.gz): 71,343,605 bytes SRR616268sub_trim2_2.fastq.gz : 63,524,791 bytes 93 bp hoge2_1.fastq.gz: 998,208 リード 71,194,586 bytes hoge2_2.fastq.gz: 998,208 リード 63,375,473 bytes アダプター除去およびフィルタリング後のファイル SRR616268sub_trim_1.fastq.gz(hoge7.fastq.gz): 998,658 リード 71,227,695 bytes SRR616268sub_trim_2.fastq.gz(hoge8.fastq.gz): 999,136 リード 63,442,904 bytes hoge1_1.fastq.gz: 998,428 リード 63,446,240 bytes hoge1_2.fastq.gz: 998,428 リード 56,019,410 bytes 53

54 Rockhopper:SE まずは single-end (SE) データとして実行 1 入力ファイルを選び 2 開く 3SUBMIT 約 1 分

55 Rockhopper:SE 1 入力ファイルは 1,000,000 リードだったので 一定の基準で内部的にフィルタリングしているのだろう 2 アセンブル後のリファレンス配列に対してマッピングした結果 1,964 リードがマップされたということだろう 3 アセンブルによって得られた転写物配列数は 8 個しかなかった 4 その平均配列長 (121 bp) と 5 中央値 (106 bp) 6 アセンブルされたトータルの塩基数は 974 これは 3 4 の結果 (8 121=968) とほぼ同じなので まあよしとしよう

56 乳酸菌 RNA-seq データ 100 万リードのオリジナル?! ファイル SRR616268sub_1.fastq.gz:107 bp, 74,906,576 bytes SRR616268sub_2.fastq.gz:93 bp, 67,158,462 bytes 100 万リードのアダプター除去後のファイル SRR616268sub_trim2_1.fastq.gz(hoge4.fastq.gz): 71,343,605 bytes SRR616268sub_trim2_2.fastq.gz : 63,524,791 bytes 次に reverse 側を single-end (SE) データとして de novo assemble 107 bp 93 bp hoge2_1.fastq.gz: 998,208 リード 71,194,586 bytes hoge2_2.fastq.gz: 998,208 リード 63,375,473 bytes アダプター除去およびフィルタリング後のファイル SRR616268sub_trim_1.fastq.gz(hoge7.fastq.gz): 998,658 リード 71,227,695 bytes SRR616268sub_trim_2.fastq.gz(hoge8.fastq.gz): 999,136 リード 63,442,904 bytes hoge1_1.fastq.gz: 998,428 リード 63,446,240 bytes hoge1_2.fastq.gz: 998,428 リード 56,019,410 bytes 56

57 Rockhopper:SE まずは single-end (SE) データとして実行 1 入力ファイルを選び 2 開く 3SUBMIT 約 1 分

58 Rockhopper:SE 1 入力ファイルは 1,000,000 リード 内部的にフィルタリングした結果 983,854 リードになったのだろう 2 アセンブル後のリファレンス配列に対してマッピングした結果 710,393 リードがマップされたということだろう 3 アセンブルによって得られた転写物配列数は 424 個 4 その平均配列長 (436 bp) と 5 中央値 (228 bp) 6 アセンブルされたトータルの塩基数は 185,233 これは 3 4 の結果 ( =184,864) とほぼ同じなので まあよしとしよう

59 Rockhopper:PE Paired-end データとしてアセンブル 1 念のため一旦消す 2 Options Clear cache というのがあったので こちらも念のためやってみる 数分フリーズ状態になりました 3Rockhopper アプリごと再起動しても キャッシュに残っているようなので Clear cache をやるのが無難かも

60 Rockhopper:PE しばらく経つと左のような画面になるので 1 paired-end の 1 つ目の forward 側のファイルを指定 2 双方向矢印 のところにカーソルをもっていくと Using paired-end reads; と書いてあるので ここをクリック 3OK を押す

61 Rockhopper:PE 1reverse 側のファイルを選択して 2 開く Paired-end の 2 つのファイルが見られるわけではないが 3SUBMIT

62 Rockhopper:PE 約 2 分で終わるが 残念ながらうまくアセンブルできていないようだ 後にこの原因は forward ファイル中の 3 側の adapter/primer 配列が主原因と判明 62

63 Rockhopper:PE それゆえ 1 Options - Parameter settings で 2Strand specific や Reverse complement reads のあたりをいじっても状況に変化はない

64 アセンブル結果まとめ いろいろやった結果のまとめです 64

65 Contents 先週の課題やエラーについて 課題 1 の解説 アダプター配列除去 (QuasR) 実行時のエラーは門田のミス Paired-end データの取り扱い フィルタリング (filtering) 基本形 : リード数が同じ paired-end の場合 (QuasR パッケージ ) 発展形 : リード数が異なる場合 (ShortRead パッケージ ) アセンブル (assemble) ゲノム用とトランスクリプトーム用 Rockhopper2( バクテリアのトランスクリプトーム用 ) マッピング (mapping) 基礎 オプション 仮想データ説明 マッピング本番 出力ファイル形式 オプションと結果の関係 カウント情報取得 実データ解析 QuasR でマッピング トリミング 高精度なアセンブルとマッピングへ 65

66 マッピング基礎 教科書 p81-89 基本イメージ マップされたリード数 = 発現量 ではないが マップされたリード数のカウント情報は 発現量推定の基本情報 基本的なマッピングプログラム (bowtie など ) を用いた場合 リファレンス配列 : ゲノム count T1 サンプルの RNA-Seq データ mapping 遺伝子 1 遺伝子 2 遺伝子 3 遺伝子 4 リファレンス配列 : トランスクリプトーム count 遺伝子 1 遺伝子 2 遺伝子 3 遺伝子 4 66

67 マッピング基礎 教科書 p81-89 マップされる側のリファレンス配列 :hoge4.fa マップする側のRNA-seqデータ ( リードと呼ばれる ): AGG マッピング = 大量高速文字列検索 という捉え方でよい マッピングプログラムの出力は ( どのリードが ) リファレンス配列上のどの位置から転写されたものかという座標情報 出力ファイル 67

68 マッピング 1 教科書 p81-89 R パッケージとして提供されているものは圧倒的に少ないですが QuasR パッケージ (Gaidatzis et al., 2015) のおかげで Windows OS 上でもマッピングを気軽に利用できるようになりました これは内部的に Bowtie (Langmead et al., 2009) プログラムを利用しています 68

69 オプション QuasR パッケージは内部的に Bowtie を利用 マッピング時に多くのオプションを指定可能 -v : 許容するミスマッチ (mismatch) 数を指定するオプション -v 0 は リードがリファレンスに完全一致するもののみレポート -v 2 は 2 塩基ミスマッチまで許容してマップされうる場所を探索 -m : 出力するリード条件を指定するオプション -m 1 は 複数個所にマップされるリードを除外して 1 か所にのみマップされたリードをレポート -m 3 は 合計 3 か所にマップされるリードまでをレポート --best --strata : 最も少ないミスマッチ数でマップされるもののみ出力する という意思表示 これをつけずに -v 2 -m 1 などと指定すると たとえ完全一致 ( ミスマッチ数 0) で 1 か所にのみマップされるリードがあったとしても どこか別の場所で 1 塩基ミスマッチでマップされる個所があれば マップされうる場所が 2 か所ということを意味し そのリードは出力されなくなる それを防ぐのが主な目的... 教科書 p81-89 オプションは デフォルトである程度よきに計らってくれるが... 実際の挙動を完全に把握できる状況で様々なオプションを試したい 69

70 教科書 p81-89 仮想リファレンス配列 マップされる側のリファレンス配列 :ref_genome.fa chr3 と chr5 の違いは 2 番目と 7 番目の塩基のみ 主に -m オプション (-m 1:1 ヶ所にマップされるリードのみ出力 ) の違いの把握が可能 70

71 教科書 p81-89 仮想リファレンス配列 マップされる側のリファレンス配列 :ref_genome.fa コピペで作成しています 71

72 教科書 p81-89 仮想 RNA-seq リード マップする側の RNA-seq データ :sample_rnaseq1.fa 許容するミスマッチ数による違いや マップされるべき場所が完全に把握できるように リードの description 行に情報を記載 72

73 教科書 p81-89 仮想 RNA-seq リード マップする側の RNA-seq データ :sample_rnaseq1.fa コピペで作成しています 73

74 Contents 先週の課題やエラーについて 課題 1 の解説 アダプター配列除去 (QuasR) 実行時のエラーは門田のミス Paired-end データの取り扱い フィルタリング (filtering) 基本形 : リード数が同じ paired-end の場合 (QuasR パッケージ ) 発展形 : リード数が異なる場合 (ShortRead パッケージ ) アセンブル (assemble) ゲノム用とトランスクリプトーム用 Rockhopper2( バクテリアのトランスクリプトーム用 ) マッピング (mapping) 基礎 オプション 仮想データ説明 マッピング本番 出力ファイル形式 オプションと結果の関係 カウント情報取得 実データ解析 QuasR でマッピング トリミング 高精度なアセンブルとマッピングへ 74

75 マッピング本番 やってみましょう

76 マッピング本番 やってみましょう 1 2 複数の RNA-seq サンプルを実行できるようにリストファイルとして与える 許容するミスマッチ数は 0 個 ( -v 0 ) 1 か所にマップされるリードのみ出力 ( -m 1 ) 76

77 マッピング本番 hoge mapping_kiso1 に最低限必要な 3 つの入力ファイルが揃っていることを確認してコピペ ファイヤーフォールか何かでアクセス許可関連のポップアップが出たら OK を押す 77

78 出力ファイル形式 出力ファイルとして実際に取り扱うのは 1BAM 形式ファイルです 2 赤枠部分はランダムに発生させる文字列部分なので ヒトによって異なります

79 出力ファイル形式 実用上は BAM 形式 視覚上は BED 形式 ゲノム上のどの位置にどのリードがマッピングされたか ( トランスクリプトームの場合どの転写物配列上のどの位置にどのリードがマッピングされたか ) を表す出力ファイル形式は複数あります 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) 79

80 出力ファイル形式 BED の最小限の情報は リード ID を含まない BAM 形式ファイル BED 形式ファイル 80

81 オプションと結果の関係 マップされなかったのは 赤枠の 8 リード中 3 リード -m 1 --best --strata -v 0 :0 ミスマッチで 1 か所にのみマップされるリードを出力 81

82 オプションと結果の関係 完全一致でも複数個所にマップされるために落とされたのは 2 リード -m 1: 1 ヶ所にマップされるリードのみ出力 -m 1 --best --strata -v 0 :0 ミスマッチで 1 か所にのみマップされるリードを出力 82

83 オプションと結果の関係 1 か所にのみマップされるがミスマッチのため落とされたのは 1 リード -v 0: 許容するミスマッチ数は 0 -m 1 --best --strata -v 0 :0 ミスマッチで 1 か所にのみマップされるリードを出力 83

84 Contents 先週の課題やエラーについて 課題 1 の解説 アダプター配列除去 (QuasR) 実行時のエラーは門田のミス Paired-end データの取り扱い フィルタリング (filtering) 基本形 : リード数が同じ paired-end の場合 (QuasR パッケージ ) 発展形 : リード数が異なる場合 (ShortRead パッケージ ) アセンブル (assemble) ゲノム用とトランスクリプトーム用 Rockhopper2( バクテリアのトランスクリプトーム用 ) マッピング (mapping) 基礎 オプション 仮想データ説明 マッピング本番 出力ファイル形式 オプションと結果の関係 カウント情報取得 実データ解析 QuasR でマッピング トリミング 高精度なアセンブルとマッピングへ 84

85 カウント情報取得 アノテーション情報を利用する場合 UCSC known Genes, Ensembl Genes など様々なテーブル名を指定可能 gene, exon, promoter, junction など様々なレベルを指定可能 アノテーション情報がない場合 R でのアノテーション情報利用は TranscriptDb が基本 マップされたリードの和集合領域を同定したのち 領域ごとのリード数をカウント BEDtools (Quinlan et al., 2010) 中の mergebed プログラムを実行して和集合領域同定後 intersectbed プログラムを実行してリード数をカウントする作業に相当 count 領域

86 カウント情報取得 アノテーション情報を利用する場合 UCSC known Genes, Ensembl Genes など様々なテーブル名を指定可能 gene, exon, promoter, junction など様々なレベルを指定可能 アノテーション情報がない場合 アノテーション情報がない場合の戦略は 複数サンプルの場合には領域が変わりうる Cufflinks を知っているヒトは cuffmerge と同じイメージだと思えばよい マップされたリードの和集合領域を同定したのち 領域ごとのリード数をカウント BEDtools (Quinlan et al., 2010) 中の mergebed プログラムを実行して和集合領域同定後 intersectbed プログラムを実行してリード数をカウントする作業に相当 sample1 count sample2 86

87 カウント情報取得 1 アノテーション情報がない場合のやり方です

88 カウント情報取得 1 hoge mapping_kiso1 に入力ファイルは揃っています マッピング結果ファイルが既に存在しますが log ファイルを見て同じオプションで実行した結果があるかどうかを自動判定します 再度マッピングを実行する必要がない場合は 以前に得られていた bam ファイルを入力として 速やかに 1 のところの作業に移行します 1 88

89 カウント情報取得 1 出力ファイルは何も指定していませんが *_range.txt という名前のファイルが作成されます これは.bam という名前のファイルを内部的に入力として読み込み その文字列中の.bam を _range.txt に置換したものを出力ファイル名として自動作成しているからそうなります 89

90 カウント情報取得 1.bed ファイルと *_range.txt ファイルを見比べると理解が深まるでしょう *_range.txt ファイルの一番右側の列がカウント情報です *.bed *_range.txt 90

91 カウント情報取得 2 複数の RNA-seq リードファイルのマッピングからカウントデータ取得までを一気に行う例です

92 カウント情報取得 2 hoge mapping_kiso2 に入力ファイルは揃っています 92

93 カウント情報取得 2 コピペ 出力ファイル群のうち 主に取り扱うのは カウントデータを含むファイル (hoge4_count.txt) 93

94 カウント情報取得 2 リストファイル中で指定したサンプル名 (sample1 と sample2) がカウントデータ行列の列名となります 94

95 Contents 先週の課題やエラーについて 課題 1 の解説 アダプター配列除去 (QuasR) 実行時のエラーは門田のミス Paired-end データの取り扱い フィルタリング (filtering) 基本形 : リード数が同じ paired-end の場合 (QuasR パッケージ ) 発展形 : リード数が異なる場合 (ShortRead パッケージ ) アセンブル (assemble) ゲノム用とトランスクリプトーム用 Rockhopper2( バクテリアのトランスクリプトーム用 ) マッピング (mapping) 基礎 オプション 仮想データ説明 マッピング本番 出力ファイル形式 オプションと結果の関係 カウント情報取得 実データ解析 QuasR でマッピング トリミング 高精度なアセンブルとマッピングへ 95

96 実データ解析 アノテーション情報がない場合のやり方です 1 96

97 実データ解析 hoge mapping_paired1 に解析したい fastq.gz ファイルを移動させて実行 約 8 分 97

98 実データ解析 1 このあたりで約 1 分 リファレンス配列のインデックス作成部分 約 3MB の乳酸菌ゲノムなので短時間で終わる ヒトゲノムだと数十分から数時間といったところか 2 のところで約 7 分 オリジナルの総リード数は 1.4 億 1/100 以下の 100 万リードだから 10 分足らずで終わる

99 実データ解析 1 マッピングの実体である qalign 関数実行前後の時間を計測し /60 = 7.05 分という結果を得ている 1 99

100 実データ解析 1 マッピングに用いられたオプション ( パラメータ ) 情報 2 これがメインのマッピング結果ファイルの名前 ( いわゆる BAM ファイル ) 3 マッピング結果の概観 forward 側と reverse 側合わせて 200 万リード中 8,204 リードしかマップされてない!

101 実データ解析 1list.files を用いて ここまでの作業で作成されたファイルを概観 元は黒枠の 4 つのファイルのみだったが マッピングが終わった時点で数多くのファイルが作成されていることが分かる 2 マッピング結果を取扱うときの基本は BAM ファイル 3 マッピングに用いられたオプション ( パラメータ ) 情報は この log ファイルからも閲覧可能

102 QC レポート作成 1sub 関数を用いて黒下線で見られるファイル名中の拡張子部分.bam を _QC.pdf に置換したファイル名を作成して out_f に格納している 拡張子部分周辺の文字列だけを変更する程度の出力ファイル名であれば このように自動的に作成することができる 数十秒程度 1 102

103 BED ファイル作成 1 これは BAM ファイルがバイナリ形式でそんなものだと慣れてしまえば不必要 単純に視覚的に見やすいように bam ファイルを読み込んで BED 形式ファイルを作成しているところ ( つまりファイル形式の変換 ) 2 得られた BED 形式ファイルの中身 マップされたリード数が 8,204 個であったことから BED ファイルの行数が 8,204 行なのは妥当

104 QC レポート概観 1QC レポートを眺める 最初のページは quality score 分布 2forward 側と 3 reverse 側で左右に表示されている カイコ small RNA-seq データ (single-end) で左側のみにしか表示されていなかったときは違和感を覚えたが paired-end データで眺めると至極妥当な配置であったことが分かる

105 QC レポート概観 2 3 1QC レポートの 2page 目 2Forward のリードの左側は確かにガタガタしていたので妥当だが 3 右側もこんなになってたっけ?! と思いつつ ( 後にこれが犯人と判明 ) 4 4page 目 リード全体の mapped/unmapped の割合を表示 確かに 200 万リード中 8,204 リードしかマップされてなかったので妥当

106 QC レポート概観 15page 目 200 万リード中 8,204 リードマップされていたが paired-end なので ここでは分母を 8,204/2 = 4, e+03 と表現しているようだ 1 ヶ所にのみマップされたのは全体の 91.7% 26page 目 これはおそらくマップされたリードの中でどのポジションにミスマッチがあったかを示す部分 3 と 4 を眺めて 確かに FastQC で眺めていたのは最初の 50 bp 分で adapters/primers の除去は 5 側のみだった!

107 アセンブル結果再考 このぶざまな状況を打破する糸口が完璧に見つかった! 107

108 Contents 先週の課題やエラーについて 課題 1 の解説 アダプター配列除去 (QuasR) 実行時のエラーは門田のミス Paired-end データの取り扱い フィルタリング (filtering) 基本形 : リード数が同じ paired-end の場合 (QuasR パッケージ ) 発展形 : リード数が異なる場合 (ShortRead パッケージ ) アセンブル (assemble) ゲノム用とトランスクリプトーム用 Rockhopper2( バクテリアのトランスクリプトーム用 ) マッピング (mapping) 基礎 オプション 仮想データ説明 マッピング本番 出力ファイル形式 オプションと結果の関係 カウント情報取得 実データ解析 QuasR でマッピング トリミング 高精度なアセンブルとマッピングへ 108

109 課題 1, 2, 3 課題 1: Rockhopper のアセンブル結果 を QuasR マッピング結果の QC レポート と絡めて説明せよ 課題 2: どうすればマッピング結果を改善できるのか ( マップ率を上げることができるのか ) 考えられる戦略を述べよ 課題 3: 今回の結果を踏まえ リファレンスゲノムがない ( マッピング結果がない ) 場合の高精度なアセンブル戦略について述べよ 109

110 課題 4 課題 4:1 は FastQC の Per base sequence content 2 は QuasR の QC レポートファイル中の同様な結果である 赤点線枠部分に違いが見られるが この理由について主に描画方法の観点から考察せよ よく分からないヒトは別の観点でもよい もちろん色の違いは本質的ではない

111 Rockhopper リトライ これだけ結果が変わります かなりイケてます 111

112 QuasR リトライ これだけ結果が変わります かなりイケてます 112