農学生命情報科学特論I - PDF 無料ダウンロード

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

講義予定第 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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

乳酸菌 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

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

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

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

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

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

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

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

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

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

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

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

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

乳酸菌 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

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

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

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

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

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

乳酸菌 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 万リードからなる生データファイルで実行 (2015.05.12 の発展課題と同じ ) 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

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

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

乳酸菌 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

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

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

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

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

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

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

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

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

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

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

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

オプション 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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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