特論I - PDF 無料ダウンロード

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

講義予定第 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) に基づく ( 非モデル生物の ) トランスクリプトーム解析を題材とした実習を含む講義を行います教科書 2

Contents( 第 2 回 ) イントロダクション (Introduction) NGS データ概観 (PacBio と Illumina) NGS データベース (DB) データ形式 (FASTQ 形式 ) SRAdb パッケージを用いたデータ取得エラーへの対処前処理 (Pre-processing) qrqc パッケージを用いた Quality Control (QC) アダプター配列除去基本戦略 (girafe パッケージ ) 昔は正常に動作していたのにという例 (QuasR パッケージ ) アダプター除去を含む様々なフィルタリングの組合せ (ShortRead パッケージ ) 課題 3

トランスクリプトームとはある特定の状態の組織や細胞中に存在する全 RNA( 転写物 transcripts) の総体様々なトランスクリプトーム解析技術マイクロアレイ DNA マイクロアレイ Affymetrix GeneChip タイリングアレイなど配列決定に基づく方法 EST SAGE CAGE 次世代シーケンサ (RNA-seq) など電気泳動に基づく方法 Differential Display AFLP HiCEP など調べたい組織でどの遺伝子がどの程度発現しているのかを一度に観察 4

トランスクリプトームとはある状態のあるサンプル ( 例 : 目 ) のあるゲノムの領域ヒト遺伝子 1 遺伝子 2 遺伝子 3 遺伝子 4 遺伝子全体 ( ゲノム ) どの染色体上のどの領域にどの遺伝子があるかは調べる個体 ( 例 : ヒト ) が同じなら不変 ( 目だろうが心臓だろうが ) AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA 転写物全体 ( トランスクリプトーム ) 遺伝子 1 は沢山転写されている ( 発現している ) 遺伝子 4 はごくわずかしか転写されてない 5

トランスクリプトームとはある状態のあるサンプル ( 例 : 目 ) のあるゲノムの領域光刺激ヒト遺伝子 1 遺伝子 2 遺伝子 3 遺伝子 4 遺伝子全体 ( ゲノム ) どの染色体上のどの領域にどの遺伝子があるかは調べる個体 ( 例 : ヒト ) が同じなら不変 ( 目だろうが心臓だろうが ) AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA 転写物全体 ( トランスクリプトーム ) 遺伝子 2 は光刺激に応答して発現亢進遺伝子 4 も光刺激に応答して発現亢進 6

トランスクリプトーム情報を得る手段教科書 p9 光刺激前 (T1) の目のトランスクリプトーム遺伝子 1 遺伝子 2 遺伝子 3 遺伝子 4 これがいわゆる遺伝子発現行列光刺激後 (T2) の目のトランスクリプトーム遺伝子 1 遺伝子 2 遺伝子 3 遺伝子 4 教科書 p9 の図 1-8 に示してあるように実際には遺伝子 = 転写物ではない点に注意! マイクロアレイ RNA-seq 7

トランスクリプトーム取得 (NGS) 次世代シーケンサー (Illumina 社の場合 ) 光刺激前 (T1) の目のトランスクリプトーム配列決定教科書 p9 ペアードエンド法断片配列の両末端が数百塩基以内の対の 2 種類の配列が得られる数百塩基程度に断片化シングルエンド法約 50-250 塩基アダプター配列を両末端に付加シングルエンド法の場合アダプター 1 アダプター 2 数百塩基程度 8

様々な NGS プラットフォーム Roche:GS FLX+ System Illumina:HiSeq System Life Technologies:SOLiD System Pacific Biosciences:PacBio RS II System 9

実際のデータ PacBio のロングリードデータも出始めています 10

実際のデータ比較的新しい論文のリードごとの塩基配列情報は見られるものの FASTQ ファイルがまだ生成されてなくてダウンロードができないこともある 12

実際のデータ DDBJ SRA (DRA) がだめな場合は NCBI SRA にトライ 13

実際のデータリードごとの塩基配列情報を見る場合はここ 14

実際のデータ最初のリードの塩基配列が表示されます 15

実際のデータ 2 番目のリードを表示リードによって長さが異なる 16

実際のデータチェックを入れるとクオリティスコアも表示される 17

実際のデータ Illumina 社の Genome Analyzer SRA061145:Marioni et al., Genome Res., 18: 1509-1517, 2008 PacBio データの長さがよく分かりますサンプルデータ 1, 36 bp Applied Biosystems 社の SOLiD4 System SRA000306:Cloonan et al., Nat. Methods, 5: 613-619, 2008 サンプルデータ 5, 25-35 bp Illumina 社の HiSeq 2000 SRA062939:Chan et al., Hum. Mol. Genet., 22: 2662-2675, 2013 サンプルデータ 25, 100 bp 18

教科書 p12-16 Contents( 第 2 回 ) イントロダクション (Introduction) NGS データ概観 (PacBio と Illumina) NGS データベース (DB) データ形式 (FASTQ 形式 ) SRAdb パッケージを用いたデータ取得エラーへの対処前処理 (Pre-processing) qrqc パッケージを用いた Quality Control (QC) アダプター配列除去基本戦略 (girafe パッケージ ) 昔は正常に動作していたのにという例 (QuasR パッケージ ) アダプター除去を含む様々なフィルタリングの組合せ (ShortRead パッケージ ) 課題 19

NGS 用データベース SRA, DRA, ENA が公式?! な NGS data repository 2 次 DB の DBCLS SRA も面白いサイトです 20

DBCLS SRA Nakazato et al., PLoS One, 2013 全体をざっくりと知りたい場合によく利用しています 21

教科書 p14 NGS データ取得は全体像の理解が大事 Chan et al., Hum. Mol. Genet., 22: 2662-2675, 2013 論文中に NGS データは Gene Expression Omnibus (GEO) 中に GSE42960 で登録したという記載あり NCBI SRA 内で GSE42960 で検索 22

NGS データ取得は全体像の理解が大事原著論文を読むことで Illumina HiSeq 2000 を使っていること 2 群間比較用データであることは既知ニコチンアミド処理群 (Nico) 対未処理群 (Untreated; UT) 生データをダウンロードすると 9GB 程度になる 23

NGS データ取得は全体像の理解が大事 GSM1054024, SRX210739, SRR649760, ID だらけです 24

NGS データ取得は全体像の理解が大事ある論文のデータ全体を指し示す ID の大元が GSE42960 でそれに付随する実験 ID の SRX210739 ラン ID の SRR649760 などの全体をまとめたデータをメタデータという Illumina の場合はフローセルと呼ばれるスライドグラス程度の大きさのものを用いて一度に 8 サンプル分の sequencing が可能この実験をラン (RUN) という 25

NGS データ取得は全体像の理解が大事リードごとの塩基配列情報を見る場合はここ 26

NGS データ取得は全体像の理解が大事ダウンロードはここだが... FASTQ 形式ではなく SRA 形式ファイルなので非推奨 27

http://rgm22.nig.ac.jp/mediawiki-ogareport/index.php/raw_data_archiving/sharing_at_ddbj 様々なファイル形式情報量 :SRA-full > SRA-lite > FASTQ > FASTA SRA-full: 塩基配列クオリティ情報 Intensity 情報など画像以外の全て SRA-lite:SRA-fullからIntensity 情報を除いて軽量化したもの FASTQ: 塩基配列とクオリティ情報のみからなるもの FASTA: 塩基配列のみからなるものファイルサイズ (SRA-full : SRA-lite : FASTQ : FASTA) 6 : 3 : 2 : 1 例 :SRA-fullはFASTQの約 3 倍 FASTQ 形式ファイルの利用が基本 28

FASTA 形式と FASTQ 形式 FASTA 形式 1 行目 : > ではじまる一行の description 行 2 行目 : 配列情報教科書 p11-13 >SEQ_ID GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT FASTQ 形式 1 行目 : @ ではじまる 1 行の description 行 2 行目 : 配列情報 3 行目 : + からはじまる 1 行 ( の description 行 ) 4 行目 : クオリティ情報 @SEQ_ID GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT +!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65 http://en.wikipedia.org/wiki/fastq_format 29

公共 DB からデータを取得する場合 ENA Sequence Read Archive (ERA; 欧 ) FASTQ 形式でダウンロード可能 NCBI Sequence Read Archive (SRA; 米 ) SRA 形式でダウンロード可能 DDBJ Sequence Read Archive (DRA; 日 ) FASTQ 形式と SRA-Lite 形式でダウンロード可能 ENA を概観しながらサンプル数と得られるファイル数の違いなどを認識論文中の情報を頼りに SRA または SRP ID を入手する手続きメタデータ ( 全体像 ) 情報を把握 30

NGS データ取得は全体像の理解が大事 Chan et al., Hum. Mol. Genet., 22: 2662-2675, 2013 論文中に NGS データは Gene Expression Omnibus (GEO) 中に GSE42960 で登録したという記載あり ENA で GSE42960 で検索 31

NGS データ取得は全体像の理解が大事 Chan et al., Hum. Mol. Genet., 22: 2662-2675, 2013 論文中に NGS データは Gene Expression Omnibus (GEO) 中に GSE42960 で登録したという記載あり実質的にどちらでもよい 32

NGS データ取得は全体像の理解が大事 Chan et al., Hum. Mol. Genet., 22: 2662-2675, 2013 ENA は全体像をつかみやすい計 6 個の FASTQ 形式ファイルになるようだ 33

NGS データ取得は全体像の理解が大事 Chan et al., Hum. Mol. Genet., 22: 2662-2675, 2013 SRA, DRA, ENA どれでもよいが論文から得られる GSE ID を頼りに SRA ID または SRP ID 情報を入手 (SRA062939 または SRP017580) するとともにメタデータ情報を把握すべしダウンロードして得られるのは Run ID から始まる SRR という名前のファイル 34

R で NGS データ取得教科書 p71-73 SRAdb (Zhu et al., BMC Bioinformatics, 14: 19, 2013) SRAdb パッケージを利用して R 経由で NCBI SRA から FASTQ ファイル群をダウンロード可能数時間かかるのでやらないで! 36

R で NGS データ取得 SRA または SRP ID を与えることでコピペで FASTQ ファイルをダウンロード可能です 37

ここまでのコマンドで NCBI SRA 中の全てのメタデータ情報を含んだ SRAmetadb.sqlite の gzip 圧縮ファイルがダウンロードされる東大有線 LAN で 20 分程度 38

ダウンロード後は自動で解凍解凍後のファイルは 6GB 程度 39

SRAmetadb.sqlite ファイルのダウンロードは結構大変だが一度ダウンロードしておけば格納されている SRA メタデータ情報は利用可能 40

ENA でみられるメタデータ情報が R Console 画面上でも見られます 41

ENA でみられるメタデータ情報が R Console 画面上でも見られます 42

1 ライブラリ名 2 総リード数 3 ファイル名 4 ファイルサイズ 1 2 3 4 1 2 3 4 43

k オブジェクトは 6 行 19 列からなる 1 ライブラリ名 2 総リード数 3 ファイル名 4 ファイルサイズ 44

計 6 個の gzip 圧縮 FASTQ ファイルが得られます 45

教科書 p11 FASTQ ファイル実例 gzip 圧縮ファイル解凍後の SRR633902_1.fastq をテキストエディタで表示 4 行で 1 つのリード情報を表すのが FASTQ 形式です 46

教科書 p11 FASTQ ファイル実例 gzip 圧縮ファイル解凍後の SRR633902_1.fastq をテキストエディタで表示解凍時にファイルが壊れていますなどというメッセージが出たらダウンロードに失敗していると解釈すべし 47

教科書 p11 FASTQ ファイル実例 gzip 圧縮ファイル解凍後のSRR633902_1.fastq をテキストエディタで表示 MD5 チェックサムが王道ですが私はペアエンドリードの場合はもう片方のファイルサイズと比較して判断します基本思考停止して再ダウンロード 48

Quality Control FASTQ 形式ファイルを入力として全体像を眺める作業 FastQC が有名だが R パッケージもいくつかある 50

qrqc パッケージで全体像を眺める qrqc はシンプルだが若干動作が不安定な印象使われた実績はあり 51

hoge フォルダ中にあり通常利用時は (Windows のヒトは ) 右クリックでファイルの保存 52

R colsole 画面上でペーストするとエラーは出るが SRR633902_1_sub-report というフォルダは作成される 53

縦軸 :Phred quality score 横軸 : 塩基のポジション 50 bp あたりから 100 bp にかけてクオリティ値が下がってますがこれが一般的な傾向です 54

SRR633902_1_sub.fastq 入力の FASTQ 形式ファイル中の ASCII コードで 1 文字表記にされた Phred スコア ( 赤矢印部分 ) を変換して表示している教科書 p10-13 に対応表などの記載あり 55

Quality Control 作業内容フィルタリング (filtering) クオリティ値の低い塩基やリードの除去 rrna や trna の除去トリミング (trimming) 最初の 35 塩基のみ利用など重複除去 (de-duplication) コンタミ (contamination) バーコード配列 (barcoding) アダプター配列除去 (adapter removal) 実験デザインや使用する機器にもよるが様々な前処理が行われます 56

Kraken Davis et al., Methods, 63: 41-49, 2013 Linux と Mac 用のみ Table 1 Table 1 を見るといろいろできるように見えるが 57

Review 論文だと Paszkiewicz et al., Front Genet., 5: 157, 2014 Table 1 評価項目によって印象は変わりますね 58

Quality Control 作業内容フィルタリング (filtering) クオリティ値の低い塩基やリードの除去 rrna や trna の除去トリミング (trimming) 最初の 35 塩基のみ利用など重複除去 (de-duplication) コンタミ (contamination) バーコード配列 (barcoding) アダプター配列除去 (adapter removal) Davis et al., Methods, 63: 41-49, 2013 特にアダプター配列除去 (adapter removal) は small RNA sequencing (srna-seq) の場合にマップ率に多大な影響を及ぼします Fig. 2 59

アダプター配列除去 girafe パッケージのデフォルト設定はイマイチですが感覚をつかむ上では便利なのでそれを利用して説明します 60

アダプター配列除去のイメージアダプター配列 :CATCGATCCTGCAGGCTAGAGACAGAT FASTQ 形式ファイル :SRR037439.fastq 3 側の 2 塩基が除去される 3 側の 5 塩基が除去される 3 側の 0 塩基が除去される 3 側の 4 塩基が除去されるアラインメント時に指定するパラメータ ( 一致に何点不一致に ) 次第で結果が変わる 3 側の 2 塩基が除去される 61

アダプター配列除去のイメージアダプター配列 :CATCGATCCTGCAGGCTAGAGACAGAT FASTQ 形式ファイル :SRR037439.fastq 3 側の 2 塩基が除去される 3 側の 5 塩基が除去される 3 側の 0 塩基が除去される 3 側の 4 塩基が除去されるこの結果は一致に +1 点不一致に -1 点を与えた場合です具体的にどういう計算をしているのだろう? 3 側の 2 塩基が除去される 62

アダプター配列除去のイメージ一塩基づつずらしたアラインメントの overlap の範囲で一致 (+1), 不一致 (-1) の総和を計算し最も得点の高かったものを採用している NNNNNNNNNNNNNNNAGACAGTTGATTTAGCATAG score(case1): -1 CATCGATCCTGCAGGCTAGAGACAGAT score(case2): -2 CATCGATCCTGCAGGCTAGAGACAGAT score(case3): -1 CATCGATCCTGCAGGCTAGAGACAGAT score(case4): -4 CATCGATCCTGCAGGCTAGAGACAGAT score(case5): +3 CATCGATCCTGCAGGCTAGAGACAGAT score(case6): -6 CATCGATCCTGCAGGCTAGAGACAGAT score(case7): -5 CATCGATCCTGCAGGCTAGAGACAGAT score(case8): -4 CATCGATCCTGCAGGCTAGAGACAGAT score(case9): -5 CATCGATCCTGCAGGCTAGAGACAGAT Case5(3 側の 5 塩基を除去 ) のスコアが最大 63

アダプター配列除去 girafe パッケージのデフォルト設定はイマイチですが感覚をつかむ上では便利なのでそれを利用して説明します 64

FASTQ ファイル読み込み後のリード塩基配列情報は sread 関数で抽出可能 table 関数を用いて配列長分布を調べているこの場合 35 bp のものが 500 個あったということ 65

入力ファイルをちゃんと読み込めていることがわかります 66

アダプター配列除去後のリード塩基配列情報は sread 関数で抽出可能 table 関数を用いて配列長分布を調べているこの場合 19 bp のものが 3 個など 67

table 関数を用いて配列長分布を調べているトリム前は 500 リードの配列長はすべて 35bp だったがトリム後に 19bp 長になっているものが 3 つ存在するそれを調べる 68

特定の条件を満たすリードを調べる配列長が 19 bp のものの位置情報を取得しその数を確認 obj が TRUE となる要素のみに対して塩基配列と description 情報を表示最も多くアダプター配列を含むリード ID を特定できたアダプター配列除去アルゴリズムの詳細を知ることで girafe のデフォルトパラメータがイマイチであることを知る 69

アダプター配列除去のイメージ一塩基づつずらしたアラインメントのoverlapの範囲で一致 (+1), 不一致 (-1) の総和を計算し最も得点の高かったものを採用している CNNNNNNNNNNNNNTGTGTCCTTGCCGTTGCAGGT score(case1): -1 CATCGATCCTGCAGGCTAGAGACAGAT score(case2): -2 CATCGATCCTGCAGGCTAGAGACAGAT score(case3): -1 CATCGATCCTGCAGGCTAGAGACAGAT score(case4): -4 CATCGATCCTGCAGGCTAGAGACAGAT score(case5): -1 CATCGATCCTGCAGGCTAGAGACAGAT score(case6): -4 CATCGATCCTGCAGGCTAGAGACAGAT score(case7): -3 CATCGATCCTGCAGGCTAGAGACAGAT score(case8): -6 CATCGATCCTGCAGGCTAGAGACAGAT score(case9): -5 CATCGATCCTGCAGGCTAGAGACAGAT score(case10): -2 score(case11): -9 score(case12): -8 score(case13): -3 score(case14): -10 CATCGATCCTGCAGGCTAGAGACAGAT CATCGATCCTGCAGGCTAGAGACAGAT CATCGATCCTGCAGGCTAGAGACAGAT CATCGATCCTGCAGGCTAGAGACAGAT CATCGATCCTGCAGGCTAGAGACAGAT CATCGATCCTGCAGGCTAGAGACAGAT score(case15): -9 score(case16): +2 CATCGATCCTGCAGGCTAGAGACAGAT ミスマッチ数が7 個! ( 黒矢印の数 ) 70

? 関数名で詳細な使用法を学ぶ一致に +1 不一致に -1 を与えているところ一致に +1 不一致に -1 を与えて 2 塩基以上のオーバーラップがないとトリムしない設定にしているのがデフォルト 71

2013 年 11 月ごろはうまくいっていましたが 2014 年 6 月に試すとエラーが出ていましたという例 73

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

アダプター配列除去 ( 推奨のやり方 ) hoge フォルダ中のファイルを解凍すれば実行できますが見るだけにして 75

全部で 11,928,428 リード配列長は 49 bp 76

総リード数は不変だがアダプター配列除去によって配列長にバリエーションができたことがわかる最短の 18 bp のものが 55,342 リード最長の 47 bp のものが 63,998 リード 77

N を含むリードがちゃんと消えていることがわかる最短の 18 bp のものが 55,264 リード最長の 47 bp のものが 63,763 リード 78

配列長の範囲を 20-30 bp に限定すると 2,619,892 リードに減ることがわかる最短の 20 bp のものが 517,002 リード最長の 30 bp のものが 58,713 リード 79

出力は FASTA 形式にしているアダプター配列や各種前処理後はクオリティスコア情報はいらないだろうという思想主なメリットはファイルサイズ 80

アダプター配列除去 ( 推奨のやり方 ) readdnastringset 関数は gzip 圧縮ファイルも読み込み可能 gzip 圧縮ファイルとして保存することも可能 81

Tips rcode_adapter.txt N を一つでも含むリードの除去を行うステップを省く場合 (# を左端に入れればそのコマンドは実行されない ) 82

課題 1 と 2 原著論文中のアダプター配列除去後の配列長分布は右表のとおりであった rcode_adapter.txt( の一部 ) Nie et al., BMC Genomics, 14: 661, 2013 1. 右表と同じように 18-44 塩基の範囲内にある small RNA リードのみを抽出するためにはどこをどう変更すればよいか示せ 2. 指定した範囲に含まれる総リード数を示せ 83

課題 3 右表に示されているように原著論文中のアダプター配列除去を含むフィルタリング後の総リード数は 11,691,441 個であった以下に様々な条件で得られた総リード数を示す条件 1-1( 許容するミスマッチ数 =1; N を含んでもよい ):11,619,415 個 Nie et al., BMC Genomics, 14: 661, 2013 条件 1-2( 許容するミスマッチ数 =1; N を全く含まない ):11,599,894 個条件 2-1( 許容するミスマッチ数 =0; N を含んでもよい ):11,357,039 個条件 2-2( 許容するミスマッチ数 =0; N を全く含まない ):11,338,479 個自分でもいくつか試し結果を簡単に考察せよ原著論文も明確に条件を記述しているわけではないので細かな違いは気にしなくてよい 84

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