Rでゲノム・トランスクリプトーム解析

06.03.05 版実習用 PC のデスクトップ上に hoge フォルダがありますこの中に解析に必要な入力ファイルがありますネットワーク不具合時はローカル環境で html ファイルを起動して各自対応してください R で塩基配列解析 : ゲノム解析からトランスクリプトーム解析まで東京大学大学院農学生命科学研究科アグリバイオインフォマティクス教育研究プログラム門田幸二 ( かどたこうじ ) kadota@iu.a.u-tokyo.ac.jp http://www.iu.a.u-tokyo.ac.jp/~kadota/

自己紹介学歴および職歴少数のスタッフで行っているアグリバイオの活動のみで基本的に手一杯ここ数年でさらに研究 << 教育のヒトに現在研究は片手間以下限界以下のスタッフ数でアグリバイオの本務を行っているため精神状態をなるべく平静に保つべく優先順位の低い活動には関与しません 00 年 3 月東京大学大学院農学生命科学研究科博士課程修了 00 年 4 月産業技術総合研究所 CBRC 003 年月放射線医学総合研究所先端遺伝子発現研究センター 005 年月 ~ 東京大学大学院農学生命科学研究科アグリバイオインフォマティクス人材養成プログラム ( 科学技術振興調整費 : 004/0-009/3) アグリバイオインフォマティクス教育研究プログラム ( 特別教育研究経費 : 009/4~04/3) アグリバイオインフォマティクス教育研究プログラム他大学の学生や社会人も受講できる希少なバイオインフォ教育プログラム科目以上の合格者数

主な活動基本スタンスは優先順位とエフォート基本独裁一匹狼ロビー活動なし門田教への勧誘なし信者になっても ( オールフリー派なので w) メリットゼロ受益者が金と時間をかけずに効率的に学べる教材整備が最優先東大アグリバイオの大学院講義 ( バイオインフォ全般 ) R を中心としたハンズオン講義 ( 平成 6 年度 ~) 受講人数が多い ( 最大 30 名 ) のでクラウド ( ウェブツール ) 系実習は実質的に不可能講義補助員 (TA) が数名のみなので Linux 系実習も困難 NBDC/ 東大アグリバイオ /HPCI の NGS ハンズオン講義 (NGS に特化 ) Linux を中心としたハンズオン講義 ( 平成 6 年度 ~) 受講人数は多い ( 最大 7 名 ; おそらくアグリバイオ本体に次ぐ規模 ) が受講生の意識レベルが高く ( きっちり予習をやるヒトが多数派 ) 環境構築済みノート PC 数 TA 数が充実しているため本格的な Linux 実習が成立しうる日本乳酸菌学会誌の NGS 連載 Linux を中心とした自習用教材 ( 平成 6 年度 ~) バクテリア ( 乳酸菌 ) データを主に Bio-Linux 上で解析するノウハウを提供第 6 回 (06 年 3 月予定 ) 分以降は DDBJ Pipeline( ウェブツール ) の利用法も紹介データ取得インストール実行に時間がかかるものも自習なので時間を気にせずにできるハンズオン講義よりも心穏やかその他研究 ( 発現変動解析精度向上のためのアルゴリズム開発や評価 ) HPCI 講習会バイオインフォマティクス実習コースの講師丸日だが上記の主要 3 項目に比べれば心穏やか 3

Contents イントロダクション (R で ) 塩基配列解析アグリバイオ NGS ハンズオン講習会日本乳酸菌学会の NGS 連載 HPCI 講習会の PC 環境ゲノム解析 NGS データ解析戦略 DDBJ Pipeline と R の関係用語説明 de novo アセンブリ実行および結果を R で解析塩基配列解析基礎 ( 塩基ごとの出現頻度解析 ) 各種テクニックや注意事項 R コードの解説塩基配列解析基礎 ( 基本情報取得 ) 塩基配列解析基礎 3( 配列長でフィルタリング ) アノテーショントランスクリプトーム配列プロモーター配列取得 4

(R で ) 塩基配列解析 03 年秋以降の講義資料や連載原稿の PDF を簡単な解説つきで公開講義資料系は年以上昔のものは参考程度ウェブサイトが見づらいとか見栄えに関する要望は無視 ( 優先順位が閾値以下 ) http://www.iu.a.u-tokyo.ac.jp/~kadota/r_seq.html 5

(R で ) 塩基配列解析 Linux 系の教材日本乳酸菌学会誌の NGS 連載第 4 回第 5 回第 6 回は 06 年 3-4 月ごろ公開予定 6

(R で ) 塩基配列解析 04 年 4 月刊行の R 本トランスクリプトーム解析全般の基礎知識的なところはこの本の第章をご覧ください 7

(R で ) 塩基配列解析講習会講義講演資料の PDF 時系列 ( 新古 ) 順にリストアップ 8

アグリバイオこれら 3 科目の講義資料を順番にみていくとよい科目名 : 農学生命情報科学特論 I 内容 : 公共 DB チェックサム QC 前処理アセンブリマッピング RPKM 発現変動など実施日 :05.06.6 05.06.3 05.06.30 05.07.07 科目名 : 機能ゲノム学内容 : データ取得正規化クラスタリング発現変動解析多重比較問題機能解析など実施日 :05.05. 05.05.9 05.05.6 05.06.09 科目名 : ゲノム情報解析基礎内容 :R の基礎 GC 含量計算や CpG 解析上流配列解析 R のバージョンの違いなど実施日 :05.04.07 05.04.4 05.04. 9

アグリバイオ例えば特論 I の第 4 回講義資料は講義資料をクリックすればよいが科目名 : 農学生命情報科学特論 I 内容 : 公共 DB チェックサム QC 前処理アセンブリマッピング RPKM 発現変動など実施日 :05.06.6 05.06.3 05.06.30 05.07.07 科目名 : 機能ゲノム学内容 : データ取得正規化クラスタリング発現変動解析多重比較問題機能解析など実施日 :05.05. 05.05.9 05.05.6 05.06.09 科目名 : ゲノム情報解析基礎内容 :R の基礎 GC 含量計算や CpG 解析上流配列解析 R のバージョンの違いなど実施日 :05.04.07 05.04.4 05.04. 0

NGS ハンズオン講習会 NGS ハンズオン講習会のほうが R については基本的なところをきっちり抑えているのでやも自習してください

NGS ハンズオン講習会 NGS ハンズオン講習会へのリンクもあり

NGS ハンズオン講習会大元 (NBDC) のサイトへのリンク NBDC のサイトのほうが見やすいのは全日程終了から数か月後に整形して公開するから当然です動画 ( 統合 TV と YouTube) も公開されている 3

NGS ハンズオン講習会 R は 05 年 7 月 9-30 日に開催動画はこちら 4

乳酸菌 NGS 連載 NGS 連載関連主に全体像原稿およびウェブ資料関係赤枠の個別の回のところで原稿中のプログラムへのリンクやコピペ用 Linux コマンドなどを利用可能 5

HPCI 講習会の PC 環境実習用 PC 環境はの手順に従って R および必要なパッケージのインストールが完了している状態です自分の PC で復習したい場合はを参考にして自力で環境構築してください 6

HPCI 講習会の PC 環境具体的な順番は R 本体のインストール各種 R パッケージのインストールです 3 の基本的な利用法の習得は HPCI 講習会の枠組みでは必須ではありません 3 7

HPCI 講習会 HPCI 講習会のバイオインフォマティクス実習コースの中の一部が門田担当 8

HPCI 講習会具体的には R を使った NGS 解析を基礎から学ぶのうちの塩基配列解析 ( 特にゲノム解析とトランスクリプトーム解析部分 ) が門田の担当 9

NGS データ解析戦略解析受託企業に外注 :Linux コマンドを知らなくてもよいクラウド ( ウェブツール ):Linux コマンドを知らなくてもよい DDBJ Pipeline (Nagasaki et al., DNA Res., 0: 383-390, 03) Illumina BaseSpace Galaxy (Goecks et al., Genome Biol., : R86, 00) Linux コマンドを駆使 ( 旧来型 ) なるべく自力で解析自分の置かれている環境予算ポリシーによっても異なるどの選択肢でも基本正解 R は主に統計解析部分で使われている Linux コマンドや NGS 解析用プログラムのインストールなどを練習しスパコンを使いこなす NBDC/ 東大アグリバイオ /HPCI の NGS ハンズオン講習会の方向性

DDBJ Pipeline と R 解析受託企業に外注 :Linux コマンドを知らなくてもよいクラウド ( ウェブツール ):Linux コマンドを知らなくてもよい DDBJ Pipeline (Nagasaki et al., DNA Res., 0: 383-390, 03) Illumina BaseSpace Galaxy (Goecks et al., Genome Biol., : R86, 00) Linux コマンドを駆使 ( 旧来型 ) なるべく自力で解析 DDBJ Pipelineだけで全てのNGS 解析ができるわけではない Rもまた然り特にRでは ( 門田の知る限り )de novoアセンブリは不可能現実を知りうまく使い分けるべし DDBJ Pipeline 上でde novoアセンブリを行った結果の解釈や確認をrで行い塩基配列解析基礎のスキルがあってよかったと思った実例を紹介 Linux コマンドや NGS 解析用プログラムのインストールなどを練習しスパコンを使いこなす NBDC/ 東大アグリバイオ /HPCI の NGS ハンズオン講習会の方向性

DDBJ Pipeline DDBJ Pipeline では主にマッピングや de novo アセンブリができる特に後者ができるのは非常に有難い 3 新規アカウント作成から de novo アセンブリまでの詳細については乳酸菌連載第 6 回ウェブ資料を参照 3 Nagasaki et al., DNA Res., 0: 383-390, 03 3

NGS データ乳酸菌 (Lactobacillus hokkaidonensis LOOC60 T ) ゲノム解読論文 (PMID: 5879859) Illumina MiSeq データ (DRR450) と PacBio データ (DRR04500) を併用することで complete genome を得ることができたという内容尚 DRR04500 は登録内容の誤りが判明し 06 年月末に削除され DRR0543-0546 に差し替えられている Tanizawa et al., BMC Genomics, 6: 40, 05 4

NGS データ Full text リンク先で全文を見られる Availability of supporting data という項目をよく眺めると生データが DDBJ Sequence Read Archive (DDBJ SRA; 略して DRA) に DRR04500 と DRR450 という ID で登録されていることがわかる Tanizawa et al., BMC Genomics, 6: 40, 05 5

NGS データ Genome sequencing and de novo assembly という項目を見ると paired-end で 5,94,60 リードと書いてある一応公共 DB(DRA) 上で確認する Tanizawa et al., BMC Genomics, 6: 40, 05 6

NGS データ Genome sequencing and de novo assembly という項目を見ると paired-end で 5,94,60 リードと書いてある 3DRR0450 という ID のほうは 4,97,30 リードと書いてある 5,94,60 / =,97,30 であるウェブサイト上の数値は single-end としてのリード数と考えれば妥当 3 4 Tanizawa et al., BMC Genomics, 6: 40, 05 7

用語 : リード断片化されたゲノム配列リードとは Sequencer で読んだ塩基配列のこと黒矢印の一本一本がリードに相当する single-end の場合 paired-end の場合 8

用語 :single-end 断片化されたゲノム配列断片化された配列の片側のみを読む場合を single-end というこの場合は右向き矢印のみ single-endの場合 paired-end の場合 9

用語 :paired-end 断片化された配列の両側から読む場合を paired-endという右向き矢印と3 左向き矢印のリードが読まれることになるそれぞれを forward 側リード reverse 側リードなどと呼ぶ断片化されたゲノム配列 single-end の場合 paired-endの場合 forward 側 3 reverse 側 30

用語 :paired-end Illumina MiSeqデータ (DRR450) の場合 forward 側 reverse 側ともに矢印の長さが 50 bp 矢印の本数( リード数 ) が計 5,94,60 個 ( 約 594 万 ; 片側のみで約 97 万 ) に相当断片化されたゲノム配列 single-end の場合 paired-end の場合 3

DDBJ SRA (DRA) DRA で Illumina MiSeq データ (DRR0450) を概観 Paired-end の FASTQ ファイルをダウンロードする場合は forward 側と 3reverse 側のつに分割されますをクリック 3 3

DDBJ SRA (DRA) forward 側 :DRR0450_.fastq.bz reverse 側 :DRR0450_.fastq.bz のような感じ DRA の場合は bzip 圧縮 FASTQ ファイルをダウンロード可能乳酸菌ゲノム配列決定論文ではこのデータを入力として de novo アセンブリが行われた Tanizawa et al., BMC Genomics, 6: 40, 05 33

用語 : コンティグ入力 :paired-end FASTQ ファイル ( 通常は )paired-end のリードファイルを入力として de novo アセンブリプログラムを実行した結果として得られる異なる複数のリードが ( ACGT の切れ目なく ) つなげられたもの contiguous sequence( 連続的な配列 ) という意味通常元のリード長よりも長くなる Assembly( コンティグの作成 ) contig contig contig3 contig4 contig5 34

用語 :scaffold 入力 :paired-end FASTQ ファイル得られたコンティグにリードをマップし Assembly( コンティグの作成 ) contig contig contig3 contig4 contig5 Scaffold contig contig contig3 contig4 contig5 35

用語 :scaffold 入力 :paired-end FASTQ ファイル得られたコンティグにリードをマップしペアの情報を頼りにコンティグ間に N を入れて連結したもの supercontig ともいう scaffold の数は contig の数よりも少なくなる尚 N を入れた部分を gap という Assembly( コンティグの作成 ) contig contig contig3 contig4 contig5 Scaffold NNNNN NNNN NNNNNN scaffold scaffold 36

用語 :gap close 入力 :paired-end FASTQ ファイル得られた scaffolds にリードをマップし Assembly( コンティグの作成 ) contig contig contig3 contig4 contig5 Scaffold NNNNN NNNN NNNNNN scaffold Gap close scaffold NNNNN NNNN NNNNNN scaffold scaffold 37

用語 :gap close 入力 :paired-end FASTQ ファイル得られた scaffolds にリードをマップし gap 周辺にマップされたリードの塩基で N を置換 gap の N がなくなり閉じていく (close) ので gap close という ( おそらく ) Assembly( コンティグの作成 ) contig contig contig3 contig4 contig5 Scaffold NNNNN NNNN NNNNNN scaffold Gap close scaffold CA A T C GG G TA A NNNCA TNNC GGNNTA scaffold scaffold 38

de novo アセンブリ入力 :paired-end FASTQ ファイル最も有名な NGS データ用 de novo ゲノムアセンブリプログラムである Velvet (Zerbino and Birney, Genome Res., 008) は Step までを実行比較的最近開発された Platanus (Kajitani et al., Genome Res., 04) は Step3 までを実行してくれる Step: Assembly contig contig contig3 contig4 contig5 Step: Scaffold Velvet NNNNN NNNN NNNNNN scaffold scaffold Step3: Gap close Platanus CA A T C GG G TA A NNNCA TNNC GGNNTA scaffold scaffold 39

乳酸菌論文は乳酸菌 (Lactobacillus hokkaidonensis LOOC60 T ) ゲノム解読論文では Illumina MiSeq データ (DRR450) の de novo アセンブリプログラムとして Platanus (ver..) を利用している Tanizawa et al., BMC Genomics, 6: 40, 05 40

DDBJ Pipeline DDBJ Pipeline では主にマッピングや de novo アセンブリができる特に後者ができるのは非常に有難い 3 新規アカウント作成から de novo アセンブリまでの詳細については乳酸菌連載第 6 回ウェブ資料を参照ここでは説明は必要最小限にして R のハンズオンへと移行する 3 Nagasaki et al., DNA Res., 0: 383-390, 03 4

DDBJ Pipeline で Platanus DDBJ Pipeline のプログラム選択画面 Velvet や Platanus を選択可能 DDBJ Pipeline: Nagasaki et al., DNA Res., 0: 383-390, 03 Platanus: Kajitani et al., Genome Res., 4: 384-395, 04 43

DDBJ Pipeline で Platanus De novo アセンブリの一般的な手順がわかっていれば赤枠内の Step-3 の説明の意味がなんとなくわかる DDBJ Pipeline は基本的にボタンをポチポチ押していくだけ DDBJ Pipeline: Nagasaki et al., DNA Res., 0: 383-390, 03 Platanus: Kajitani et al., Genome Res., 4: 384-395, 04 44

DDBJ Pipeline で Platanus アセンブリ終了後の画面 Platanus 実行結果ファイル (platanusresult.zip) をダウンロードして解凍したのが Platanus: Kajitani et al., Genome Res., 4: 384-395, 04 45

DDBJ Pipeline で Platanus アセンブリ終了後の画面 Platanus 実行結果ファイル (platanusresult.zip) をダウンロードして解凍したのが hoge フォルダ中の platanusresult Platanus: Kajitani et al., Genome Res., 4: 384-395, 04 46

DDBJ Pipeline で Platanus 入力 :paired-end FASTQ ファイル一般的な de novo アセンブリの手順を知っておけばファイル名から最終的な結果が 3out_gapClosed.fa だと認識できる Step: Assembly contig contig contig3 contig4 contig5 Step: Scaffold 3 NNNNN NNNN NNNNNN scaffold scaffold Step3: Gap close CA A T C GG G TA A NNNCA TNNC GGNNTA scaffold scaffold 47

塩基配列解析基礎入力 :paired-end FASTQ ファイル ( アセンブリ実行結果の )multi-fasta ファイルを読み込んで塩基ごとの出現頻度解析ができる Step 実行後 (out_contig.fa) は N がなく Step 実行後 (out_scaffold.fa) に N ができて 3Step3 実行後 (out_gapclosed.fa) に N が減るのだろうと妄想できるそれを自力で確認することでアルゴリズムの理解を深めることができる Step: Assembly contig contig contig3 contig4 contig5 Step: Scaffold 3 NNNNN NNNN NNNNNN scaffold scaffold Step3: Gap close CA A T C GG G TA A NNNCA TNNC GGNNTA scaffold scaffold 50

塩基ごとの出現頻度解析 ( アセンブリ実行結果の )multi- FASTA ファイルを読み込んで塩基ごとの出現頻度解析を行う項目 5

塩基ごとの出現頻度解析例題 7 が Platanus の Step3 実行後のファイル (out_gapclosed.fa) を入力とするものなのでそのままコピペできて便利これを実行します 5

塩基ごとの出現頻度解析つまり Platanus 実行結果ファイル (platanusresult.zip) をダウンロードし解凍して得られた platanusresult フォルダ中の out_gapclosed.fa を入力として塩基ごとの出現頻度解析を行う 53

R の起動と作業ディレクトリ変更ファイルディレクトリの変更 3 デスクトップ hoge - platanusresult を指定する 3 54

作業ディレクトリ変更ヒトによって若干見栄えは違うだろうが 5-7 が同じになればよい 3 5 4 6 7 55

getwd() 作業ディレクトリ変更の確認です getwd() と打ち込んで確認ののように見えていれば OK 56

list.files() R 上で現在の作業ディレクトリ中のファイルを眺めるのが list.files() GUI 画面上で眺めている platanusresult フォルダ中のものと同じものが見えていることがわかる 57

list.files() 参考ファイルが存在しないフォルダ上で list.files() とやると character(0) という結果になる 58

塩基ごとの出現頻度解析当たり前ですが解析したいディレクトリ ( またはフォルダ ) を正しく指定できていなければエラーに遭遇しますまた解析したいファイルが存在しない状態でもエラーが出ます今は解析したい入力ファイル (out_gapclosed.fa) が R Console 画面上でも 3 見えているのでエラーなく動くはずです 3 59

基本はコピペ一連のコマンド群をコピーして R Console 画面上でペースト Windows のヒトは CTRL と ALT キーを押しながらコードの枠内で左クリックすると全選択できますトリプルクリックでも OK Macintosh はよくわかりません 60

途中経過と終了後コピペ直後と実行後の状態エラーなく実行できたときはこんな感じになります一見何も変化がないように見えますが 6

結果の解説解析結果 ( 塩基ごとの出現頻度情報 ) は hoge7.txt というファイルに保存されている list.files() とやると確かに自分が出力ファイル名として指定した 3hoge7.txt が存在することがわかる 6

結果の解説もちろん出力ファイル (hoge7.txt) は手の届く場所 ( つまり作業ディレクトリ内 ) にある getwd() や 3 現在時刻を表示する date() はただの確認用 4 エクセルで眺めるとこんな感じ 4 3 63

R 上で眺める赤枠程度の情報量ならエクセルなどをわざわざ開くまでもなく R 上で眺めればよい例えばここでは出力ファイル名を out_f というオブジェクト名で取り扱っている 3out_f と打てば対応関係がわかる 3 64

R 上で眺める R コードの最後の部分がファイルに保存するところ out_f に書き込んでいるのは out というオブジェクトの情報 3out の中身を見れば hoge7.txt と同じ情報を得られる 3 65

sum で総塩基数を得る out オブジェクトは数値ベクトル sum は数値ベクトルの総和を計算する関数 out に対して実行した結果 (,356,09) は入力ファイル (out_gapclosed.fa) の総塩基数を調べていることと同義 66

sum で総塩基数を得る DDBJ Pipeline 実行結果画面上の数値と同じ入力ファイル (out_gapclosed.fa) は DDBJ Pipeline 上で Platanus という de novo アセンブリプログラムを実行した結果だったことを思い出そう 67

目的をおさらい入力 :paired-end FASTQ ファイル ( アセンブリ実行結果の )multi-fasta ファイルを読み込んで塩基ごとの出現頻度解析ができる Step 実行後 (out_contig.fa) は N がなく Step 実行後 (out_scaffold.fa) に N ができて 3Step3 実行後 (out_gapclosed.fa) に N が減るのだろうと妄想できるそれを自力で確認することでアルゴリズムの理解を深めることができる Step: Assembly contig contig contig3 contig4 contig5 Step: Scaffold 3 NNNNN NNNN NNNNNN scaffold scaffold Step3: Gap close CA A T C GG G TA A NNNCA TNNC GGNNTA scaffold scaffold 68

目的をおさらい入力 :paired-end FASTQ ファイル ( アセンブリ実行結果の )multi-fasta ファイルを読み込んで塩基ごとの出現頻度解析ができる Step 実行後 (out_contig.fa) は N がなく Step 実行後 (out_scaffold.fa) に N ができて 3Step3 実行後 (out_gapclosed.fa) に N が減るのだろうと妄想できるそれを自力で確認することでアルゴリズムの理解を深めることができるを調べるにはどうすればいいか? Step: Assembly contig contig contig3 contig4 contig5 Step: Scaffold NNNNN NNNN NNNNNN scaffold scaffold Step3: Gap close CA A T C GG G TA A NNNCA TNNC GGNNTA scaffold scaffold 69

入力ファイルを変更テンプレートの out_gapclosed.fa を out_scaffold.fa に変更すればよい 7

入力ファイルを変更適当なテキストエディタ ( ここでは EmEditor) に例題をコピペし必要最小限の変更を施したところ 7

変更後のコードをコピペコピペ 73

ありがちなミスこれはエラーメッセージです w エラーの理由は出力予定ファイル (hoge7.txt) を開くことができないというもの Permission denied( 権限が与えられていない ) はアク禁みたいなものです Tips: ワードパッドやメモ帳で開く分にはエラーは出ないようです 74

ありがちなミスエラーの原因はエクセルで hoge7.txt を開いているから閉じて再実行すればエラーは出なくなる 75

再実行エクセルを閉じて再実行した結果エラーは出ていないことがわかる out オブジェクトの中身を見ると確かに N がある! 76

納得できる結果入力 :paired-end FASTQ ファイル入力ファイル (out_scaffold.fa) のイメージはのような感じなので N が 49 個あったという結果は合理的 Step: Assembly contig contig contig3 contig4 contig5 Step: Scaffold NNNNN NNNN NNNNNN scaffold scaffold Step3: Gap close CA A T C GG G TA A NNNCA TNNC GGNNTA scaffold scaffold 77

ありがちなミス最終行の部分で改行をキチンと含めないとハマる 78

ありがちなミス最終行の部分で改行をキチンと含めないと最後の write.table 関数部分が実行されないつまりファイルが作成されません 79

実際の利用時は hoge フォルダ直下にある rcode.txt のような無駄なコメントを除いてスリムにした一連のスクリプトを作成しておき一気にコピペ 80

一気に結果を得る hoge フォルダ直下にある rcode.txt のような無駄なコメントを除いてスリムにした一連のスクリプトを作成しておき一気にコピペコピペ後に自分が指定した出力ファイルができていることを確認 8

結果のまとめ result_step*.txt の結果をまとめたものが 8

結果の解釈入力 :paired-end FASTQ ファイル Step 実行後は N が 0 Step 実行後に N が 49 個生成されたということはいくつかの contigs がまとめられて scaffolds になったのだろう 3 Step3 で N が 0 個になったのはおそらくたまたまうまくいっただけ 49 個よりも減ったということが重要で gap close がうまく機能したと判断できる Step: Assembly contig contig contig3 contig4 contig5 Step: Scaffold 3 NNNNN NNNN NNNNNN scaffold scaffold Step3: Gap close CA A T C GG G TA A NNNCA TNNC GGNNTA scaffold scaffold 83

コード内部の理解は重要入力 :paired-end FASTQ ファイル R を使うことでアセンブリプログラムの内部挙動の把握や理解ができる他の例は配列数 (contig 数や scaffold 数と書くと説明しづらいので配列数に統一 ) 配列数は Step Step で減り Step Step3 では不変だろうと予想 Step: Assembly contig contig contig3 contig4 contig5 Step: Scaffold NNNNN NNNN NNNNNN scaffold scaffold Step3: Gap close CA A T C GG G TA A NNNCA TNNC GGNNTA scaffold scaffold 84

R コードの解説配列数の把握の仕方の前に Tips list.files() 実行時に pattern オプションをつけて任意の文字列を含むもののみ表示させることが可能ここでは out_ という文字列を含むもの ( ファイル ) のみ表示させている入力ファイルの存在確認 86

R コードの解説配列数の把握の仕方赤枠部分をコピペ 87

R コードの解説 4 Biostrings という R パッケージを library 関数で読み込んで Biostrings パッケージが提供する関数群を利用可能な状態にする (Biostrings が提供する )readdnastringset 関数を用いて 3FASTA 形式の 4 入力ファイルを読み込んだ結果を 5fasta というオブジェクト ( ものという理解でよい ) に格納 4 3 5 88

Tips: 配列数 fasta オブジェクトの中身を表示 ( ここでの目的の ) 配列数は 7 個スカラー値として配列数情報のみ取り出したい場合は 3 ベクトルの要素数を調べる length 関数を利用する 3 89

答え合わせ DDBJ Pipeline 実行結果の数値 (7 個 ) と同じことがわかります最長の配列 (Maximum contig size; 57,78 bp) と最短の配列 (Minimum contig size; 0 bp) も R 上で把握できます 90

Tips: 配列長配列長の情報は (DNAStringSet という形式で保持されている )fasta オブジェクト中の width 列の位置に相当する 9

Tips: 配列長配列長情報は width(fasta) とやることで数値ベクトルとして取り出すことができるこの程度 (7 個 ) の配列数ならパッと見で最長と 3 最短のものを確認できるが 3 9

Tips: 配列長ベクトル演算の基本関数を駆使して全貌を把握する上矢印キーを回押して以前打ち込んだコマンドを出すなど上下左右の矢印キーを駆使して効率的に打ち込むべしやの数値はつ前のスライドには存在しないがこれは 3summary 関数実行結果として表示させる有効数字のデフォルトが 4 桁だから summary(width(fasta), digits=6) とすればの 57700 が正しく 5778 と表示されるようになる 3 93

Tips:description 部分 description 行部分は names という関数を用いることで ( 文字列 ) ベクトルとして取り扱うことができるここでは 3:4 という指定を行って最初の 4 個分のみ表示させている 3 94

Tips: 塩基配列部分但しこのノリは塩基配列部分には通用しない w seq という関数は別の意味を持つこと fasta オブジェクトの主要な中身がこの塩基配列情報であるためと理解すればよい 95

Tips: 塩基配列部分どうしても文字列ベクトルなどで取り出したい場合は as.character 関数を使うが DNAStringSet 形式の fasta オブジェクトをそのまま用いて各種塩基配列解析を行うのが通常のやり方 96

alphabetfrequency 塩基ごとの出現頻度解析の中核となっている関数は alphabetfrequency 実行結果である hoge の中身は数値行列この段階で塩基配列解析から数値解析に切り替わる塩基の種類には多型がある例えば 3M (A or C) 4K (G or T) など門田はシロイヌナズナのゲノム配列で ACGTN 以外のものを見た記憶あり 3 4 98

dim dim 関数で行数と列数を把握 alphabetfrequency は配列ごとに結果を返しているので 7 行からなると解釈 8 列であることから塩基の種類数は 8 個と解釈 3 行列の一部要素の取り出し例 3 99

is.element は出現頻度情報を得たい塩基の種類を指定するところほとんどの場合 ACGTN のみで事足りるのでこのようにしている is.element 関数は条件判定 ( 集合演算 ) を行っている行列 hoge の列名 (column names) からなるベクトルの中から 3param_base で指定された要素が存在する場所を TRUE そうでないところを FALSE と評価するのが 4is.element 関数 3 3 4 00

条件を満たす列のみ行列の subsetting は [ 行, 列 ] で指定する [, 列 ] で列のみの指定 [ 行, ] で行のみの指定となる hoge[, obj] は obj ベクトルの TRUE となっている列の位置のみ取り出すことに相当する hoge[:, ] で hoge 行列の最初の行分のみ表示 3hoge[:, obj] の合わせ技でさらに param_base で指定した塩基のみを出力できるようになる 3 0

colsums alphabetfrequency 実行結果は配列ごとに各塩基の出現頻度を計算しているそのため hoge は 7 行分の要素からなる colsums は行列データを入力として列ごとに総和 (column sum) を計算する関数 colsums を適用することで配列ごとではなくファイル全体の出現頻度を得ることができる ( 今得たい情報はこれ ) 0

colsums は最初の行分のみで列ごとの総和を計算する場合ではエラーとなっている colsums の入力として与えている hoge[, obj] は最初の行分のみからなるつまり入力が次元の行列データではなく次元のベクトルになってしまっているため 3 行頭に # をつけており実際にはこのコードは動作していない 3 03

apply が一般的かもこのコードで実際に動かしているのは apply 関数を用いるほう結果は colsums と同じおそらく行列演算で行ごとや列ごとに何かを行うときには一般に apply 関数を用いるので一応示した 04

apply の説明 apply は入力データに対して列ごと ( 行ごとの場合はここをにする ) に 3 総和を計算する sum 関数を適用するみたいな指定を行う colsums だと sum を計算することしかできないが apply の場合は 3 のところの関数名を mean, median, max などいろいろ自在に変更できる 3 05

as.matrix 実はこの入力ファイルの場合は as.matrix という関数をつけなくても 3 つけたときと同じ結果が得られるつけている理由は apply(as.matrix( ),, sum) などとして行ごとに sum 関数を適用したいときに配列数が複数の場合でも単数の場合でも統一的にエラーなく処理できたという記憶があったから 3 06

as.matrix 挙動の違いは入力データの行列が行しかない ( 配列数がつの ) 場合に出てくる複数行からなる ( 配列数がつ以上の ) 場合と比べればエラーメッセージの違いがわかります 07

思考停止するべからず as.matrix をつけてエラーメッセージが出てないからといってこれが正しいわけではないことに気づこう 08

思考停止するべからずの実行結果である 050 という数値は単純に番目の配列の長さ今調べたいのは塩基ごとの出現頻度情報なので 3 が正解! 3 09

思考停止するべからず少なくともこのサンプルコードは配列数がつしかない場合にはうまく動かないことが既知欲しい結果が ( この場合は ) 数値ベクトルになっていない段階でおかしいと思えるようになりましょう一般論としては得られる結果をイメージし特にイメージと異なる場合に疑いの目で結果を眺めよう 0

R コードの解説赤下線のように沢山のオプションを駆使している sep= t は区切り文字を指定するオプション t はタブ区切りの意味 row.names=t は行の名前 (row names) を TRUE にせよという意味ここが T になると 3 赤枠部分の情報が追加される FALSE にするとこの列は消える 4 col.names=f は col.names=t にしたときに無意味なヘッダー行が含まれるのが嫌だったのでこうしているだけ 3 4

目的をおさらい入力 :paired-end FASTQ ファイル配列数は Step Step で減り Step Step3 では不変だろうと予想 Step: Assembly contig contig contig3 contig4 contig5 Step: Scaffold NNNNN NNNN NNNNNN scaffold scaffold Step3: Gap close CA A T C GG G TA A NNNCA TNNC GGNNTA scaffold scaffold

目的をおさらい 3 配列数は Step Step で減り Step Step3 では不変だろうと予想 3hoge フォルダ直下の rcode.txt は配列数をカウントする必要最小限のコード 349 7 7 で予想通りの結果 3

塩基配列解析基礎の情報を一気に得る項目もあります 5

塩基配列解析基礎 FASTA 形式ファイルを読み込んで各種情報を得る項目 6

塩基配列解析基礎 rcode3.txt 3 例題の入力ファイル名部分を out_gapclosed.fa に変更した 3rcode3.txt をコピペで実行 7

塩基配列解析基礎 rcode3.txt 出力ファイル (hoge.txt) を開かずに write.table 関数で書きだしている tmp の中身を表示 8

塩基配列解析基礎 rcode3.txt DDBJ Pipeline 上の Platanus 実行結果と完全一致 GC 含量情報なども得られる 9

塩基配列解析基礎 rcode3.txt 配列数の算出法 length(fasta) や最短配列長 min(width(fasta)) も前のスライドで解説したものと同じです 0

塩基配列解析基礎 3 このアセンブル結果の最短配列長は 0 bp 通常アセンブル結果ファイルから一定の配列長 ( 例 :300 bp) 未満のものは除去される

塩基配列解析基礎 3 FASTA 形式ファイルを読み込んで指定した配列長以上のもののみ残して FASTA 形式ファイルで出力する項目例題 5 は out_gapclosed.fa を読み込んで 300 bp 以上の配列のみ hoge5.fasta ファイルに保存するスクリプト

塩基配列解析基礎 3 赤枠部分をコピペ入力ファイルを読み込んだ直後の fasta オブジェクトは 37 個の配列からなる赤下線で見えているものが 300 bp 未満なのでフィルタリングされる 3 3

塩基配列解析基礎 3 width(fasta) は配列長情報からなる数値ベクトル 300 bp という閾値情報からなる param_length で条件判定した結果が obj に格納されている 4

塩基配列解析基礎 3 param_length 以上 (>=) という条件を満たすものが TRUE そうでないものが FALSE 5

塩基配列解析基礎 3 オリジナルの 7 配列からなる fasta オブジェクトの中から obj が TRUE となる (300 bp 以上の ) 配列は 5 個 6

塩基配列解析基礎 3 こういう上書きはアリですもちろん fasta みたいな別名にしてもいいがヒトゲノム配列などを取り扱うときにはノート PC レベルではメモリ的に厳しくなります 7

塩基配列解析基礎 3 writexstringset 関数を使えば fastaオブジェクトの中身を指定したファイルに書きだすことができるのXStringSetのXは何でもよいみたいな意味 fastaがdnastringssetという形式で格納されていることアミノ酸配列 (Amino Acids) を格納する形式として AAStringSetという形式が存在することからそれらを同じ関数で統一的に取り扱えるようにするため 8

塩基配列解析基礎 3 出力ファイルは FASTA 形式で保存した hoge5.fasta 行あたりの塩基数を 50 個に指定している 9

アノテーションアノテーション ( 遺伝子注釈付け ) はアセンブル後の配列を入力として与えどこ ( 座標 ) にどんな遺伝子 (gene symbols; gene names; products) がありどんな Gene Ontology ID や KEGG Pathway 上に存在するかなどを得る作業広範囲 KEGG 系 3 バクテリアに特化などいろいろあります 3 3

アノテーションアノテーションファイルの形式は GFF/GTF が有名 3

GFF/GTF 形式ファイルの例 GFF3 形式 ( シロイヌナズナ ; TAIR0_GFF3_genes.gff) 他に refflat 形式など様々なファイル形式が存在します GTF 形式 ( ゼブラフィッシュ ; Danio_rerio.Zv9.75.gtf) 33

GFF の読み込み読み込み段階でコケる読み込みはうまくいったがその後の解析段階でコケるなど Linux 上での解析同様一筋縄ではいきません過去の受講生など多方面からの情報提供のおかげでだいぶ分かってきました 34

GFF の読み込み例題 7 ここで用いている GFF 形式の入力ファイルは 3 から取得しました 3 をクリックしたつもりでよい w 3 35

Ensembl 解説 GFF ファイルはここから取得の gzip 圧縮ファイルをダウンロードして解凍したものが入力ファイル 3 のあたりがバージョン番号概ね月単位でバージョン番号が上がっていく 3 36

Ensembl 解説でこのゲノムの全貌をある程度把握可能原著論文の情報なども合わせることで chromosome and plasmids 環状ゲノムであることも認識可能 3 でゲノム配列も取得できる 3 Tanizawa et al., BMC Genomics, 6: 40, 05 37

Ensembl 解説いろんなものがあって私はよくわかりませんが GFF ファイルと一緒に取り扱いたいときには GFF ファイルと似た名前のを採用します Tanizawa et al., BMC Genomics, 6: 40, 05 38

GFF の読み込み例題 7 が読み込みの基本形 GenomicFeatures というパッケージが提供する 3makeTxDbFromGFF 関数を用いて GFF ファイルを読み込んで TxDb という独特の形式で取り扱えるようにする入力ファイルは 4 hoge L.hokkaidonensis 中にある 4 3 4 39

GFF の読み込み入力ファイルが GFF ver.3 という形式になっていないみたいな警告メッセージが出ているが読み込んだ後の txdb オブジェクトは大丈夫そうだたぶん 40

GFF の読み込み若干自信がないのは GFF ファイル読み込み後ので見えている数値と Ensembl ウェブサイト上で見られる数値が一致していないことに由来,344 や,4 はプラスミドを含むものなのか詳細は不明 4

転写物配列取得 multi-fasta ファイル ( ゲノム配列情報 ) と GFF ファイル ( アノテーション情報 ) を同時に読み込むことで例えばトランスクリプトーム ( 転写物 ) 配列情報を一気に取得することも可能例題 5 43

転写物配列取得は GFF ファイル情報を保持した txdb オブジェクトから transcripts という関数を用いて抽出したい転写物の座標情報を取得した結果を hoge に保存している 44

転写物配列取得 GFF ファイルの見方がよくわかっていなくてもうまく読み込めているらしいことはわかる 45

転写物配列取得 in_f で指定したゲノム配列情報はここで登場ゲノム配列からで指定した座標情報の塩基配列を 3(Biostrings パッケージが提供する )getseq 関数を用いて取得 4(Rsamtools パッケージが提供する )FaFile 関数は getseq 関数利用時に必要なおまじない 4 4 3 46

転写物配列取得 getseq 実行後の fasta オブジェクトが欲しいトランスクリプトーム配列情報ではあるが 47

転写物配列取得の fasta オブジェクトをそのまま FASTA 形式で保存するとで見えているがままの description 情報が書きだされるつまりすべて Chromosome になってしまう 48

転写物配列取得赤枠部分で行っているのは description 部分の記述内容を Chromosome_start_end としてどこの座標由来の塩基配列かがわかるようにしている paste は文字列を sep オプションで指定した文字を間に挟んで連結する関数 3 の例をみれば挙動がわかると期待 3 49

転写物配列取得 description 部分が変わっていることがわかるこれを眺めるだけで出力ファイルをみなくてもうまくいっていると判断できる ( と油断していると時々落とし穴があるので注意 ) 50

プロモーター配列取得基本的にはプロモーター配列取得もトランスクリプトームのときと同じノリ例題 0 5

プロモーター配列取得例題 0 は転写開始点上流 00 bp 下流 0 bp の領域を取得するコードはその領域情報これだけだと 3 確かに 0 bp 分の領域になっていることの確認しかできない 3 53

プロモーター配列取得例題 0 は転写開始点上流 00 bp 下流 0 bp の領域を取得するコード元となっている転写開始点情報を transcripts(txdb) で strand 情報も含めて比較するとよくわかる 54

プロモーター配列取得例題 0 は転写開始点上流 00 bp 下流 0 bp の領域を取得するコード最後まで実行した結果転写開始点上流 00 bp 下流 0 bp の領域を取得するコードなので配列長が全て 0 bp になっており妥当 55

失敗例例題は転写開始点上流 00 bp 下流 0 bp の領域を取得するコード例題 0 との違いは上流の塩基配列数のみ 56

失敗例例題は転写開始点上流 00 bp 下流 0 bp の領域を取得するコード例題 0 との違いは上流の塩基配列数のみ 57

失敗例例題は転写開始点上流 00 bp 下流 0 bp の領域を取得するコード例題 0 との違いは上流の塩基配列数のみ 58

失敗例例題は転写開始点上流 00 bp 下流 0 bp の領域を取得するコード例題 0 との違いは上流の塩基配列数のみ hoge は取得したいプロモーター配列の座標情報 3getSeq を実行するとエラーが出る 3 59

思考停止するべからず例題は転写開始点上流 00 bp 下流 0 bp の領域を取得するコード例題 0 との違いは上流の塩基配列数のみ hoge は取得したいプロモーター配列の座標情報 3getSeq を実行するとエラーが出る 4 fasta と打つと何か出力されるがうまく取れているわけではない! 3 4 60

思考停止するべからず 3 例題は転写開始点上流 00 bp 下流 0 bp の領域を取得するコードこの fasta オブジェクトの中身は 3 この R Console 画面上で以前に行っていた例題 0( 転写開始点上流 00 bp 下流 0 bp の領域を取得するコード ) 実行時に作成されたものが残っているだけであるその証拠に 4 ここが 0 4 6

大事な計算時は Rを再起動し真っ新な状態でコピペするのが一番スッキリ私はいつもコレ普段はいいえを押して作業スペースを保存せずに終了させるがここでは3キャンセルにしておく 3 6

オブジェクトの表示現在利用可能なオブジェクトの表示は ls() 63

オブジェクトの消去 ( この R Console 画面上で利用可能な ) 全オブジェクトの消去は rm(list = ls()) 64

例題再実行 Tips ディレクトリの変更は setwd でも可能お約束のコマンドはのような任意のファイルに書き込んでおいて R 起動直後に無条件でコピペしておく ( のが門田の習慣 ) 66

例題再実行例題再実行結果 getseq 実行部分でコケるところまでは一緒 3 以降は fasta オブジェクトがないのでそれを用いる部分は軒並みエラー祭りになっていることがわかる重要なのはエラーの原因を正確に把握すること 3 67

エラー原因解説エラーを把握すべく欲しい領域 hoge と 3getSeq 実行部分を再度表示 3 68

エラー原因解説エラーの原因はに書かれているように 6 番目のレコード (Chromosome:77844-78053) が切り捨てられている (truncated) というものこれは 3 最後の転写開始点 (77853 番目の塩基 ) の上流 00 bp から下流 0 bp の領域に相当する 3 69

エラー原因解説 seqinfo 関数を使うことで in_f で指定したファイルの配列長情報 (77985 bp) を取得可能 3 で取得しようとしていた領域の一部が存在しないことが原因である 3 70

例題が推奨コード例題は取得予定の座標が存在するかどうかを判定し存在しないものをフィルタリングする部分を追加したコード ( 甲斐政親氏提供 ) 今問題となっているの領域が除去 (filter out) されればよいそれを行ってくれる 7

例題が推奨コード 3 が取得予定の座標が存在するかどうかを判定し存在しないものをフィルタリングする部分 ( 甲斐政親氏提供 ) 3 7

例題が推奨コード赤枠のフィルタリング実行後の状態エラーの原因であった 6 番目の領域がなくなっていることがわかる 73

例題が推奨コード例題を最後まで実行した結果 74

例題は多数つの項目内にも多数の例題がありますうまく動かないままわかってて放置してあるもの気づいていないもの作成当時はうまく動いていたが R 本体のバージョンが上がってからうまく動かなくなっているもの条件判定が不十分なものなど玉石混交です 75

Contents トランスクリプトーム解析イントロダクション : 簡単な原理基本イメージ様々な解析目的解析データ : 乳酸菌 (L. casei A) QuasRでマッピング ( 基礎 ): コード各部の説明と結果の解釈 QuasRでマッピング ( 応用 ): オプションを指定して実行カウント情報取得, サンプル間クラスタリング (TCC) 発現変動解析 (TCC) M-A plot モデル分布統計的手法 3 群間比較 (TCCによるANOVA 的な解析 ) 遺伝子間クラスタリング (MBCluster.Seq) 3 群間比較 (TCCによるANOVA 的な解析 + MBCluster.Seqでのパターン分類 ) 76

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

トランスクリプトーム解析ある状態のあるサンプル ( 例 : 目 ) のあるゲノムの領域遺伝子遺伝子遺伝子 3 遺伝子 4 遺伝子全体 ( ゲノム ) どの染色体上のどの領域にどの遺伝子があるかは調べる個体 ( 例 : ヒト ) が同じなら不変 ( 目だろうが心臓だろうが ) 働いているRNAの種類や量を調べるのが目的光刺激ヒト 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 転写物全体 ( トランスクリプトーム ) 遺伝子は光刺激に応答して発現亢進遺伝子 4 も光刺激に応答して発現亢進 mrna 78

トランスクリプトーム解析光刺激前 (T) の目のトランスクリプトーム遺伝子遺伝子遺伝子 3 遺伝子 4 状態の異なる複数サンプルのデータを取得して解析するのが一般的サンプル間比較光刺激後 (T) の目のトランスクリプトーム遺伝子遺伝子遺伝子 3 遺伝子 4 79

トランスクリプトーム解析光刺激前 (T) の目のトランスクリプトーム遺伝子遺伝子遺伝子 3 遺伝子 4 具体的な目的はやの発現変動遺伝子同定などこれがいわゆる遺伝子発現行列光刺激後 (T) の目のトランスクリプトーム遺伝子遺伝子遺伝子 3 遺伝子 4 80

データ取得光刺激前 (T) の目のトランスクリプトーム遺伝子遺伝子遺伝子 3 遺伝子 4 現在は NGS の利用が主流 NGS を用いた RNA の配列決定 (sequencing) なので RNA-seq これがいわゆる遺伝子発現行列光刺激後 (T) の目のトランスクリプトーム遺伝子遺伝子遺伝子 3 遺伝子 4 Togo picture gallery by DBCLS is Licensed under a Creative Commons 表示. 日本 (c) 8

RNA-seq 概略入力 : サンプルの RNA 出力 : 大量塩基配列データ入力 : 抽出された RNA 断片化出力 : 塩基配列 NGS で配列決定アダプター付加 Togo picture gallery by DBCLS is Licensed under a Creative Commons 表示. 日本 (c) 8

RNA-seq 概略入力 : 抽出された RNA NGS の出力はリードと呼ばれる 00 塩基程度の短い配列が延々と続く巨大なファイル各矢印がつのリードに相当この段階ではまだどのリードがどの転写物由来かは不明 ( なので灰色一色 ) 断片化出力 : 塩基配列 NGS で配列決定アダプター付加 Togo picture gallery by DBCLS is Licensed under a Creative Commons 表示. 日本 (c) 83

遺伝子転写物ある状態のあるサンプル ( 例 : 目 ) のあるゲノムの領域遺伝子遺伝子遺伝子 3 遺伝子 4 遺伝子全体 ( ゲノム ) どの染色体上のどの領域にどの遺伝子があるかは調べる個体 ( 例 : ヒト ) が同じなら不変 ( 目だろうが心臓だろうが ) 赤枠部分の表現は本当は不正確昔は実験機器の解像度が事実上遺伝子レベルだった遺伝子発現解析という表現はその名残りヒト AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA AAAAAAA 転写物全体 ( トランスクリプトーム ) 遺伝子は沢山転写されている ( 発現している ) 遺伝子 4 はごくわずかしか転写されてない mrna 84

遺伝子転写物ある状態のあるサンプル ( 例 : 目 ) のあるゲノムの領域遺伝子遺伝子遺伝子 3 遺伝子 4 ある遺伝子領域から転写 (transcription) されている転写物 (transcript) は種類とは限らないヒト 85

遺伝子転写物ある状態のあるサンプル ( 例 : 目 ) のあるゲノムの領域遺伝子遺伝子遺伝子 3 遺伝子 4 ある遺伝子領域から転写 (transcription) されている転写物 (transcript) は種類とは限らない例えば遺伝子の領域では 3 種類の真の転写物が存在しそのうち種類は既知とするヒト遺伝子領域 exon exon exon3 既知転写物既知転写物未知転写物真の転写物情報 86

遺伝子転写物ある状態のあるサンプル ( 例 : 目 ) のあるゲノムの領域遺伝子遺伝子遺伝子 3 遺伝子 4 実際の細胞内 ( 例 : 目のサンプル ) での発現情報 ( 働いている度合い ) がのような感じだったとするヒト遺伝子領域高発現 exon exon exon3 既知転写物既知転写物低発現未知転写物中発現真の転写物情報真の発現情報 87

遺伝子転写物ある状態のあるサンプル ( 例 : 目 ) のあるゲノムの領域遺伝子遺伝子遺伝子 3 遺伝子 4 NGS 機器を用いて転写されている mrna 配列決定 (RNA-seq) をした結果のイメージヒト遺伝子領域高発現 exon exon exon3 既知転写物既知転写物低発現未知転写物中発現真の転写物情報真の発現情報 RNA-seqで得られるリード情報 ( 色は不明 ) 88

データ解析の出発点トランスクリプトーム (RNA-seq) データ解析の出発点は RNA-seq データファイル RNA-seq データ 89

データ解析の出発点トランスクリプトーム (RNA-seq) データ解析の出発点は RNA-seq データファイルゲノム配列情報 RNA-seq データ 90

データ解析の出発点遺伝子遺伝子遺伝子 3 遺伝子 4 遺伝子領域トランスクリプトーム (RNA-seq) データ解析の出発点は RNA-seq データファイルゲノム配列情報 3 ゲノム上のどこにどんな遺伝子 exon 転写物が存在するかというアノテーション情報 exon exon exon3 既知転写物既知転写物 RNA-seq データ 9

解析結果のイメージ RNA-seq データゲノム配列情報 3 アノテーション情報を利用して 4 未知転写物 ( 新規 isoform) の同定ができる遺伝子遺伝子遺伝子 3 遺伝子 4 遺伝子領域 exon exon exon3 既知転写物既知転写物未知転写物 4 RNA-seq データ 9

解析結果のイメージ遺伝子遺伝子遺伝子 3 遺伝子 4 RNA-seq データゲノム配列情報 3 アノテーション情報を利用して 4 未知転写物 ( 新規 isoform) の同定ができる 5 転写物の発現量 ( 働いている度合い ) 推定も原理的に可能遺伝子領域 5 高発現 exon exon exon3 既知転写物既知転写物低発現未知転写物中発現 RNA-seq データ 93

具体的な戦略は? RNA-seq データゲノム配列情報 3 アノテーション情報を利用して 4 未知転写物 ( 新規 isoform) の同定ができる遺伝子遺伝子遺伝子 3 遺伝子 4 遺伝子領域 exon exon exon3 既知転写物既知転写物未知転写物 4 RNA-seq データ 94

具体的な戦略 RNA-seq データ中の本本のリード ( 横棒 ) がゲノム上のどの領域から転写されたのかを調べる文字列検索と本質的に同じでありこれがマッピングという作業に相当するゲノム RNA-seq データ 95

具体的な戦略 RNA-seq データ中の本本のリード ( 横棒 ) がゲノム上のどの領域から転写されたのかを調べる文字列検索と本質的に同じでありこれがマッピングという作業に相当するゲノム RNA-seq データ 96

具体的な戦略リードの長さが初期は 35 塩基程度だったが現在は 50 塩基程度まで伸びているそのおかげでリードを分割してマッピングすることもできるゲノム RNA-seq データ 97

具体的な戦略分割してマップされたリードは大抵の場合複数のエクソン (exon) をまたぐリードでありジャンクションリード (junction read) と呼ばれるジャンクションリードゲノム exon exon exon3 RNA-seq データ 98

具体的な戦略既知遺伝子 ( 転写物 ) の座標情報と比較することで答え合わせも可能ジャンクションリードゲノム exon exon exon3 既知転写物既知転写物アノテーション情報 ( 既知遺伝子座標情報 ) RNA-seq データ 99

具体的な戦略同様にして他のジャンクションリードも既知転写物と比較することでジャンクションリードゲノム exon exon exon3 既知転写物既知転写物アノテーション情報 ( 既知遺伝子座標情報 ) RNA-seq データ 00

具体的な戦略未知転写物 ( 新規 isoform) の同定も原理的に可能未知転写物?! ジャンクションリードゲノム exon exon exon3 既知転写物既知転写物アノテーション情報 ( 既知遺伝子座標情報 ) RNA-seq データ 0

様々な解析目的トランスクリプトーム配列取得ゲノム配列既知の場合参考新規転写物同定などに相当がメインプログラム多くのメインプログラム内部でや 3 のサブプログラムが動作する例えば Bowtie-Tophat-Cufflinks パイプラインは Cufflinks 内部でジャンクションリードもマップ可能な Tophat( や Tophat) が動作しておりそのさらに内部で基本マッピングプログラムである 3Bowtie( や Bowtie) が動作している 3 03

様々な解析目的トランスクリプトーム配列取得ゲノム配列未知の場合参考トランスクリプトーム配列の de novo アセンブリに相当多くのプログラムは発現量 (FPKM 値 ) も出力してくれます 04

様々な解析目的発現量の正確な推定参考転写物の発現量を正確に推定したい場合は専用のプログラムを使うべし RSEM が有名 3Sailfish も高速なアルゴリズムとして有名 4TIGER は日本語で質問できる上最近の手法比較論文 (Kanitz et al. Genome Biol., 05) でも高評価でおススメ 3 4 05

様々な解析目的発現変動解析 ( 群間比較 ) 参考群間比較で反復あり ( 複製あり ) データの場合は edger 反復なしの場合は DESeq を内部的に用いて頑健な結果を返す TCC がおススメ反復の有無に応じて内部的に用いるパッケージを自動で切り替える 06

様々な解析目的発現変動解析 (3 群間比較 ) 参考 3 群間比較で反復あり ( 複製あり ) データの場合は edger 反復なしの場合は DESeq を内部的に用いて頑健な結果を返す TCC がおススメ (Tang et al., BMC Bioinformatics, 05) 反復の有無に応じて内部的に用いるパッケージを自動で切り替える 07

解析データ ( 乳酸菌 ) マップする側 (paired-end RNA-seq データ ;SRR6668) オリジナルデータ (Illumina HiSeq 000 で取得 ) の情報リード長 :forward 側は 07 bp reverse 側は 93 bp リード数 : ともに 34,755,996 リード ( 約.35 億 ) データ量 :bzip 圧縮状態で計約 5GB 非圧縮 FASTQ で計約 80GB 下記手順実行後のデータ ( 計約 0MB) をマッピングに利用. 最初の00 万リードのみ抽出 ( 計 00 万リード ). forward 側 :3 側 7 bpをトリム後にアダプターを除去 998,649リード 3. reverse 側 :3 側 bp をトリム後にアダプターを除去 999,33 リード 4. 両方で存在するリードのみ抽出 998,5 = 計,997,04 リード forward 側 (SRR6668sub_trim3_.fastq.gz) reverse 側 (SRR6668sub_trim3_.fastq.gz) マップされる側 (Lactobacillus casei A) ゲノムサイズ :,907,89 bp 遺伝子数 :,799 個 Lactobacillus casei A 株のデータ乳酸菌 NGS 連載第 3 回最後のほうでダウンロードしたものと基本的に同じトリムの理由は第 5 回でわかる Lactobacillus_casei_a.GCA_000309565..5.dna.chromosome.Chromosome.fa Lactobacillus_casei_a.GCA_000309565..5.chromosome.Chromosome.gff3 07 bp 93 bp 00 bp 9 bp 08

参考教科書 p-3 FASTA 形式と FASTQ 形式 FASTA 形式がわかるヒトはそれに quality 情報のみが追加されたものという理解でよい FASTA 形式行目 : > ではじまる一行の description 行行目 : 配列情報 >SEQ_ID GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT FASTQ 形式行目 : @ ではじまる行の description 行行目 : 配列情報 3 行目 : + からはじまる行 ( の description 行 ) 4 行目 : クオリティ情報 @SEQ_ID GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT +!''*((((***+))%%%++)(%%%%).***-+*''))**55CCF>>>>>>CCCCCCC65 http://en.wikipedia.org/wiki/fastq_format 09

解析データ ( 乳酸菌 ) 3 DRAで見たマップする側のpaired-end RNA-seqデータ accession 番号は SRR6668 リード数はに示されている大元は3SRAという形式で保存されている DRAでは FASTQ 形式ファイルを提供している SRAからFASTQファイルを作成する際に 4でみえているようなNが多くクオリティが低いリードが除去されるため FASTQファイルで取り扱う際にはで見えているリード数よりも若干少なくなる 4 0

解析データ ( 乳酸菌 ) マップされる側のゲノム配列データとアノテーションデータ (GFF ファイル ) 情報は Ensembl から取得バージョンは年に数回程度以上の頻度で上がっている印象基本は最新だが最も重要なのはゲノム配列とアノテーション情報のバージョンが同じであること

マッピング基礎 QuasR は general なパッケージ精度云々というよりは Windows の R 環境でマッピングを可能にしてくれたという点で感謝 m( )m 前処理マッピング QC レポートカウント情報取得まで可能なのでこのあたりの全貌を大まかに理解するうえで便利 3

マッピング基礎 QuasR は basic なマッピングプログラムである bowtie (Langmead et al., Genome Biol., 009) とジャンクションリードをマッピング可能な SpliceMap (Au et al., Nucleic Acids Res., 00) を選択可能ここではバクテリア ( 乳酸菌 ) ゲノムにマップするので高速な basic aligner を利用例題 4

マッピング基礎マップされる側のリファレンス配列マップする側はリストファイル ( タブ区切りテキストファイル ) として与えるリストファイルの中身は 3paired-end の場合はこんな感じ 3 5

マッピング基礎リストファイルとして与えるメリットは複数サンプルの場合を考えればよいのファイルは行からなるが 3 行目以降に同様な形式で追加するサンプル数分だけ行を増やしていけばよい SampleName 列はカウントデータ ( 後述 ) を得るときにここに記載さいたサンプル名 ( ここでは namae_paired) が列名として使われる 6

解析データ ( 乳酸菌 ) マップする側 (paired-end RNA-seq データ ;SRR6668) オリジナルデータ (Illumina HiSeq 000 で取得 ) の情報リード長 :forward 側は 07 bp reverse 側は 93 bp リード数 : ともに 34,755,996 リード ( 約.35 億 ) データ量 :bzip 圧縮状態で計約 5GB 非圧縮 FASTQ で計約 80GB 下記手順実行後のデータ ( 計約 0MB) をマッピングに利用. 最初の00 万リードのみ抽出 ( 計 00 万リード ). forward 側 :3 側 7 bpをトリム後にアダプターを除去 998,649リード 3. reverse 側 :3 側 bp をトリム後にアダプターを除去 999,33 リード 4. 両方で存在するリードのみ抽出 998,5 = 計,997,04 リード forward 側 (SRR6668sub_trim3_.fastq.gz) reverse 側 (SRR6668sub_trim3_.fastq.gz) マップされる側 (Lactobacillus casei A) ゲノムサイズ :,907,89 bp 遺伝子数 :,799 個マップする側のリストファイルとマップされる側の FASTA 形式ファイル Lactobacillus_casei_a.GCA_000309565..5.dna.chromosome.Chromosome.fa Lactobacillus_casei_a.GCA_000309565..5.chromosome.Chromosome.gff3 7

マッピング基礎デスクトップ hoge L.casei_A_genome フォルダ内に必要なファイルはあります作業ディレクトリをそこに変更して 3list.files() 3 8

マッピング基礎赤枠部分が必要最小限 QuasR パッケージを読み込んで 3 主要関数である qalign を実行するだけで BAM 形式のマッピング結果ファイルを得ることができるの赤枠部分をコピペこのプログラムの機能のいくつかが Windows ファイアウォールでブロックされていますというアラートウィンドウが出ることもあるが 4 その場合はキャンセルボタンを押す ( でないと先に進めない ) 約 4 分 3 4 9

解説 R Console 画面上で見えているのはのあたりマッピング部分の所要時間は 04.96 秒計,997,04 リード中 3693,500 個がマップされ 4,303,54 個がマップされなかったことがわかる 3 4 0

解析データ ( 乳酸菌 ) マップする側 (paired-end RNA-seq データ ;SRR6668) オリジナルデータ (Illumina HiSeq 000 で取得 ) の情報リード長 :forward 側は 07 bp reverse 側は 93 bp リード数 : ともに 34,755,996 リード ( 約.35 億 ) データ量 :bzip 圧縮状態で計約 5GB 非圧縮 FASTQ で計約 80GB 下記手順実行後のデータ ( 計約 0MB) をマッピングに利用. 最初の00 万リードのみ抽出 ( 計 00 万リード ). forward 側 :3 側 7 bpをトリム後にアダプターを除去 998,649リード 3. reverse 側 :3 側 bp をトリム後にアダプターを除去 999,33 リード 4. 両方で存在するリードのみ抽出 998,5 = 計,997,04 リード forward 側 (SRR6668sub_trim3_.fastq.gz) reverse 側 (SRR6668sub_trim3_.fastq.gz) マップされる側 (Lactobacillus casei A) ゲノムサイズ :,907,89 bp 遺伝子数 :,799 個 seqlength の数値はマップされる側のリファレンス配列の総塩基数 Lactobacillus_casei_a.GCA_000309565..5.dna.chromosome.Chromosome.fa Lactobacillus_casei_a.GCA_000309565..5.chromosome.Chromosome.gff3

解説 qalign 関数実行時にマッピングプログラム (bowtie or SpliceMap) の指定を行わなかった理由はデフォルトが basic aligner (unspliced aligner ともいう ) の bowtie であることを知っていたから qalign 関数上では splicedalignment = FALSE として表現される主なマッピング結果である BAM ファイルは 3 拡張子が.bam だがこれはバイナリファイルなので中身を眺めても意味不明 ( 爆 ) 3

解説作業フォルダを眺めると確かに.bam ファイルが作成されているマッピング後の解析は基本的に BAM 形式ファイルを入力として取り扱う 3 赤下線部分の文字列はヒトそれぞれ 3 3

QC レポート PDF 入力データのクオリティスコアやマッピング結果の概要などを PDF ファイルとして出力できます最後の部分が _QC.pdf となる PDF ファイルが作成されます 4

QC レポート PDF 確かに最後の部分が _QC.pdf となる PDF ファイルが作成されていることが分かりますこれを解説します 5

PDF の一部を解説リードのポジションごとのクオリティスコア分布と塩基組成 6

PDF の一部を解説マップされたリードの割合計,997,04 リード (e+06) 中 334.7% がマップされた 4 693,500/,997,04 = 0.3473 なので妥当 3 4 7

PDF の一部を解説 4 3 マップされたリードのうちか所にのみマップされたリード (uniquely mapped reads) が 87.7% 3 複数個所にマップされたリード (non-unique) が.3% 43.47e+05 は片側のみで考えているのかもしれないいずれにせよ 5693500 0.877 = 60800 個程度はゲノム中のか所にのみマップされたと解釈 5 8

BED ファイル作成 BED 形式ファイルはバイナリファイルのため中身を解釈できない BAM ファイルをテキストファイルに変換して可視化したものという理解でよい赤枠部分をコピペ.bed というファイルが作成される 3for というループを回しているのは複数サンプルに対応するため 3 9

BED 概観作成された BED ファイルをエクセルで眺めるマッピング結果はどの配列上の 3 どこから (start) 4 どこまで (end) の場所にリードがマップされたかが最低限わかればよいマップされたリードの総数は 693,500 個だったのでヘッダー行 ( つまり行目に列名情報 ) がないこのファイルの場合 5693,500 行となるのは妥当 3 4 5 30

マッピング応用 QuasR はオプション無指定でもマッピングを行ってくれるがログファイル (QuasR_log_be8cf3864.txt) を眺めても -v で指定する許容するミスマッチ数が何だったのか不明例題 3bowtie (basic aligner のつ ) 利用時に - m -best -strata v 0 オプションをつけて実行する 3 3

マッピング応用以前のマッピング結果が残っている場所で実行してよいが念のためオブジェクトの全消去をやっておこう 33

マッピング応用終了後の状態マップされたリードの総数 (494,9 個 ) はデフォルトの結果 (693,500 個 ) に比べて減っているデフォルトオプション実行時の -v で指定する許容するミスマッチ数は以上だったのだろうと推測可 34

マッピング応用赤枠部分が -m -best -strata v 0 オプションをつけて実行した結果ファイル QC レポート PDF 中の 3 マップされたリードの割合 4 計,997,04 リード (e+06) 中 54.8% がマップされた 494,9/,997,04 = 0.478 なので妥当 3 5 4 35

マッピング応用マップされたリード ( 片側 47,456 リードで計 494,9 リード ) のうちか所にのみマップされたリード (uniquely mapped reads) が 86.6% 3 複数個所にマップされたリード (non-unique) が 3.4% この結果は直感的にオカシイ理由は -m としてか所にのみマップされるリードを出力しているつもりだから 3 36

? 関数名 Bowtie 実行時のオプションを眺めるべく alignmentparameter のところを詳細に調査 qalign 関数の詳細を調べたいときは 3?qAlign と打つ数秒後に html マニュアルが開く 3 37

関数マニュアル qalign 関数の html マニュアルページ上部 Usage ( 基本的な利用法 ) Arguments( オプションの説明 ) 3Details( 詳細情報 ) 4Value( 返り値 ; どんな結果が返ってくるか ) などの情報が見られるこれはかなり難しい例なので最初のうちは sum, mean, alphabetfrequency など挙動を完全に把握できている関数のマニュアルを眺めて慣れておくとよい 38

関数マニュアル調べたい alignmentparameter のところをつぎはぎで表示結論として -m としてか所にのみマップされるリードを出力しているつもりなのになぜ non-unique が 3.4% 含まれるという結果になるのか理解できないもしかしたら複数個所にマップされるリードはランダムにどこかか所の結果が返されるということなのだろうか? 少なくとも門田はマニュアルの説明だけでは挙動を完全にイメージできないプログラムのバグかもしれないし門田の勘違いかもしれないがこれ以上は深追いしない 39

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

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

カウント情報取得アノテーション情報を利用するやり方の例題をやってみましょう 43

違いを説明マッピング基礎の例題と基本的に同じ違いは GFF 形式のアノテーションファイルの指定カウントデータ取得のレベルを指定 3 アノテーションファイルの読み込み 3 44

違いを説明マッピング結果 out とアノテーション情報 txdb を読み込んで param_reportlevel( この場合 gene レベル ) で指定したレベルのカウントデータを取得する qcount 実行部分および結果の保存の部分 45

コピペで実行同じデフォルトオプションで実行したマッピング基礎の例題の結果がある作業ディレクトリ上で念のためオブジェクトの全消去を行ってからコピペ分弱 46

無駄を省く分弱で終わる理由はマッピングを行っていないから QuasR はまずログファイルを調べるそして以前に同じ入力ファイル同じオプションで実行した結果を見つけたら (QuasR_log_be8cf3864.txt) ( あなたが実行したい ) 全てのマッピング結果が見つかったよとなり以前の結果を読み込んでくれる 47

無駄を省くマッピング基礎の例題の結果と同じく 693,500 リードがマップされた out オブジェクトを利用可能なので 3qCount 関数でカウントデータを得ることができる 3 48

qcount qcount 関数実行によって得られたカウントデータ情報を含む count オブジェクトの行数と列数は 77 行列 3count オブジェクトの最初の 6 行分を head 関数で表示 3 49

qcount 列目は配列長列目が目的のカウント情報 3param_reportlevel で指定していたのは 4gene それゆえ 5 行名は遺伝子名の通し番号のようになっているのだろう 4 3 5 50

qcount 遺伝子名の通し番号っぽいものは txdb オブジェクトの元情報である 3in_f3 で読み込んだアノテーションファイル中に書かれている筈 3 5

GFF を眺めて確認 GFF ファイルをエクセルで眺めているの文字列が out_f で指定した出力ファイル (hoge.txt) 中の 3 行名として使われているのだろう 3 5

GFF を眺めて確認基本テクニックを駆使して配列長 350 bpの LCAA_067を確認同じ配列長になっており妥当完全に欠番のない通し番号になっていれば3count 行列の67 行目という指定でもイケるわけではない 4tailで納得 3 4 53

カウント情報取得さきほどまでは例題の解説次は例題 4 55

カウント情報取得例題 4 はマッピング応用の例題と同じくマッピング時のオプションを明示的に指定している許容するミスマッチ数を 0(-v 0) にしたときのマップされたリード数 ( 494,9 個 ) はデフォルトの結果 (693,500 個 ) よりも少なかったカウント行列データの数値は全体的に少なくなると予想できるのでそれを確認するのが主目的 56

カウント情報取得コピペで実行し全体像を見られるところを表示マップされたリード数は確かに 494,9 個だった 3 カウント行列データの数値は確かに全体的に少ない 3 57

比較許容するミスマッチ数を 0(-v 0) にしたときの結果のほうが確かにデフォルトの結果に比べて少ない 58

出力ファイル例題と例題 4 のもうつの違いは配列長 (hoge4_length.txt) と 3 カウント情報 (hoge4_count.txt) を別々に出力している点 3 59

出力ファイルカウント情報ファイル (hoge4_count.txt) の中身はこんな感じ配列情報の列を含まないものが一般的なカウントデータ以降の統計解析の入力ファイルとして利用される 60

カウントデータ解析カウントデータこのデータ (Blekhman et al., 00) は 3 種類の生物種間比較 0,689 genes 8samples ヒト (HS) チンパンジー (PT) アカゲザル (RM) 生物種ごとにメス 3 匹オス 3 匹発現変動遺伝子 (DEG) 同定サンプル間クラスタリング Blekhman et al., Genome Res., 0: 80-89, 00 6

サンプル間クラスタリング 0,689 遺伝子 8 サンプルの biological replicates のみからなるカウントデータ (Blekhman et al., 00) のサンプル間クラスタリングデータの取得や整形については 05.07.9 の講義資料を参考 63

入力ファイルデスクトップ上の hoge フォルダ中に sample_blekhman_8.txt が存在するはずファイルが存在しないヒトは右クリック 3 対象をファイルに保存でデスクトップ上の hoge に保存 3 Blekhman et al., Genome Res., 0: 80-89, 00 64

0,689 genes 入力ファイルこのデータは 3 種類の生物種間比較ヒト (Homo sapiens; HS) チンパンジー (Pan troglodytes; PT) アカゲザル (Rhesus macaque; RM) 生物種ごとにメス 3 匹オス 3 匹雄雌を考慮しなければ biological replicates ( 生物学的な反復 ) は 6 ヒト (Homo sapiens; HS) チンパンジー (Pan troglodytes; PT) アカゲザル (Rhesus macaque; RM) メス (Female) オス (Male) メスオスメスオス Blekhman et al., Genome Res., 0: 80-89, 00 65

コピペで実行一連のコマンド群をコピーして R Console 画面上でペーストブラウザが Internet Explorer の場合は CTRL と ALT キーを押しながらコードの枠内で左クリックすると全選択できます 66

実行結果エラーなく実行できると右下のような画面になっているはずです入力ファイル情報を格納した行列 data の行数が 0,689 列数が 8 となっていることがわかります 67

400 ピクセル出力ファイル出力ファイルのサイズを指定しているのでこんな感じになる hoge8.png 700 ピクセル 68

結果の解釈ヒト (HS) オス3 匹 (M, M, M3) メス3 匹 (F, F, F3) チンパンジー (PT) オス3 匹 (M, M, M3) メス3 匹 (F, F, F3) アカゲザル (RM) オス3 匹 (M, M, M3) メス3 匹 (F, F, F3) hoge8.png 3 生物種間全体で眺めるとヒト (HS) とチンパンジー (PT) はよく似ているヒト (HS) チンパンジー (PT) アカゲザル (RM) 69

結果の解釈ヒト (HS) オス3 匹 (M, M, M3) メス3 匹 (F, F, F3) チンパンジー (PT) オス3 匹 (M, M, M3) メス3 匹 (F, F, F3) アカゲザル (RM) オス3 匹 (M, M, M3) メス3 匹 (F, F, F3) hoge8.png 群間比較 ( 発現変動遺伝子検出 ; DEG 検出 ) を行ったときに HS vs. RM で得られる DEG 数のほうが 3 HS vs. PT で得られる DEG 数よりも多そうヒト (HS) チンパンジー (PT) アカゲザル (RM) 3 70

結果の解釈ヒト (HS) オス3 匹 (M, M, M3) メス3 匹 (F, F, F3) チンパンジー (PT) オス3 匹 (M, M, M3) メス3 匹 (F, F, F3) アカゲザル (RM) オス3 匹 (M, M, M3) メス3 匹 (F, F, F3) hoge8.png 同一生物種でクラスターを形成している RMM は外れサンプルっぽいヒト (HS) チンパンジー (PT) アカゲザル (RM) 7

結果の解釈ヒト (HS) オス3 匹 (M, M, M3) メス3 匹 (F, F, F3) チンパンジー (PT) オス3 匹 (M, M, M3) メス3 匹 (F, F, F3) アカゲザル (RM) オス3 匹 (M, M, M3) メス3 匹 (F, F, F3) hoge8.png ヒト (HS) とアカゲザル (RM) はメスとオスのサンプルが入り混じっているこれらの生物種内でメス群 vs. オス群の群間比較を行っても DEG はほとんど検出されないだろうヒト (HS) チンパンジー (PT) アカゲザル (RM) 7

結果の解釈ヒト (HS) オス3 匹 (M, M, M3) メス3 匹 (F, F, F3) チンパンジー (PT) オス3 匹 (M, M, M3) メス3 匹 (F, F, F3) アカゲザル (RM) オス3 匹 (M, M, M3) メス3 匹 (F, F, F3) hoge8.png チンパンジー (PT) に限っていえばメス 3 匹がクラスターを形成しているのでメス群 vs. オス群の群間比較結果として多少なりとも DEG が検出されるだろうヒト (HS) チンパンジー (PT) アカゲザル (RM) 73

この解釈の仕方は原著論文あり但し TCC パッケージが提供する clustersample 関数を用いた結果以外では責任を持ちませんまた運悪く例外データセットもあるかも Tang et al., BMC Bioinformatics, 6: 36, 05 74

論文レベルの図です手法比較原著論文 (Tang et al., 05) の Additional file 6 に全く同じ図もありつまりコピペで作成したクラスタリング結果がそのまま論文の図に使えるということ 3 書き方は赤下線のような感じでよい 3 Tang et al., BMC Bioinformatics, 6: 36, 05 75

HS vs. RM HS vs. RM の群間比較を TCC パッケージ (Sun et al., 03) で行ってみよう 77

サブセット抽出例題黒枠内をコピペ 78

サブセット抽出ここで取得したいサブセットの列番号やグループ情報を指定発現変動解析に用いるサブセットは 0,689 genes 4 samples のデータ 3 正しくヒト (HS) vs. アカゲザル (RM) になっている 3 79

サブセット抽出入力ファイル (sample_blekhman_8.txt) を眺めるなどして該当サンプルの列の位置を把握していることが前提 80

FDR 赤枠部分までをコピペ q < 0.05 を満たす遺伝子数は,488 個 False discovery rate (FDR) = 0.05 はこの閾値を満たす,488 個を発現変動遺伝子 (Differentially Expressed Genes; DEGs) とみなすと,488*0.05 = 4.4 個は偽物であることを意味する有意水準 (false positive rate; FPR)5% と同じような位置づけであり FDR5% というのは許容する偽物 (non-deg) 混入割合に相当する 8

FDR 最後までコピペ q < 0.30 を満たす遺伝子数は 4,786 個 FDR = 0.30 なので 4,786*0.30 =,435.8 個は偽物で残りの 70% は本物だと判断する 8

DEG 数の見積もり FDR 閾値が比較的緩めのところを眺め 0,689 genes 中 3,300 個程度が本物の DEG と判断する 83

樹形図と一致今比較しているのは HS vs. RM クラスタリング結果からもこれらの発現プロファイルの類似度が低い ( 距離が遠い ) ので妥当 84

M-A plot これが M-A plot 発現変動遺伝子 (DEG) と判定されたものが多数存在することがわかる param_fdr で指定した閾値 (0.05) を満たす遺伝子群がマゼンタ色で表示されている 86

M = log G - log G - - 0 M-A plot DEG が存在しないデータの M-A plot を眺めることで縦軸の閾値のみに相当する倍率変化を用いた DEG 同定の危険性が分かります群間比較用横軸が全体的な発現レベル縦軸が log 比からなるプロット名前の由来はおそらく対数の世界での縦軸が引き算 (Minus) 横軸が平均 (Average) G 群 < G 群 G 群で高発現 G 群 = G 群 G 群 > G 群 G 群で高発現 3 4 5 A = (log G + log G)/ 低発現全体的に高発現 Dudoit et al., Stat. Sinica, : -39, 00 87

DEG 検出結果基本的にはこれが解析結果 88

DEG 検出結果 G(HS) 群 G(RM) 群位は RM 群 (G 群 ) で高発現の DEG p-value とその順位 G 群で高発現 M-A plot の A 値と M 値 q-value FDR 閾値判定結果 q-value < 0.05 を満たす DEG が non-deg が 0 G 群で高発現 89

DEG 検出結果 G(HS) 群 G(RM) 群位も RM 群 (G 群 ) で高発現の DEG p-value とその順位 G 群で高発現 M-A plot の A 値と M 値 q-value FDR 閾値判定結果 q-value < 0.05 を満たす DEG が non-deg が 0 G 群で高発現 90

DEG 検出結果 G(HS) 群 G(RM) 群 3 位は HS 群 (G 群 ) で高発現の DEG p-value とその順位 G 群で高発現 M-A plot の A 値と M 値 q-value FDR 閾値判定結果 q-value < 0.05 を満たす DEG が non-deg が 0 G 群で高発現 9

DEG 検出結果 G(HS) 群 G(RM) 群指定した FDR 閾値 (0.05) をギリギリ満たす,488 位の遺伝子 p-value とその順位 G 群で高発現 M-A plot の A 値と M 値 q-value FDR 閾値判定結果 q-value < 0.05 を満たす DEG が non-deg が 0 G 群で高発現 9

様々な M-A plot 例題 -5 のコピペで作成しました 93

分布やモデル ( 当たり前だが )FDR 閾値を緩めると得られる DEG 数は増える傾向にあることがわかる例題 6 のコピペで作成 HS vs. RM HS vs. RM HS vs. RM 73 DEGs (FDR 0.0%),488 DEGs (FDR 5%) 5,435 DEGs (FDR 40%) 厳しい FDR 閾値緩い少ない DEG 数多い 95

分布やモデル重要 : 黒の分布は non-deg の分布に相当 HS vs. RM HS vs. RM HS vs. RM 73 DEGs (FDR 0.0%),488 DEGs (FDR 5%) 5,435 DEGs (FDR 40%) 厳しい FDR 閾値緩い少ない DEG 数多い 96

統計的手法とは HS vs. RM HS vs. PT 同一群 ( 下段 ) の分布は異なる群 ( 上段 ) の non-deg 分布とよく一致する同一群内のばらつきの分布 (non-deg 分布 ) 以外のものが DEG と判定されるのが統計的手法の結果,488 DEGs 578 DEGs ヒト (HS) チンパンジー (PT) HS vs. HS PT vs. PT RM vs. RM RM vs. RM3 アカゲザル (RM) 7 DEGs 6 DEGs 4 DEGs 0 DEGs 97

統計的手法とは例題 3 同一群内のばらつきの分布 (non- DEG 分布 ) から遠く離れたところに位置するものは 0 に近い p-value 同一群内の遺伝子のばらつきの程度を把握し帰無仮説に従う分布の全体像を把握しておく ( モデル構築 ) non-deg のばらつきの程度を把握しておくことと同義実際に比較したい群の遺伝子のばらつきの程度が non-deg 分布のどのあたりに位置するかを評価 ( 検定 ) HS vs. RM 98

統計的手法とは例題同一群内のばらつきの分布 (non- DEG 分布 ) のど真ん中に位置するものはに近い p-value 同一群内の遺伝子のばらつきの程度を把握し帰無仮説に従う分布の全体像を把握しておく ( モデル構築 ) non-deg のばらつきの程度を把握しておくことと同義実際に比較したい群の遺伝子のばらつきの程度が non-deg 分布のどのあたりに位置するかを評価 ( 検定 ) 99

結果の比較 ( 倍変化 ) 倍率変化 (fold-change; FC) での DEG 検出結果下段の同一群内比較でも多数の偽陽性が検出されている例題 3 をベースに作成 HS vs. RM HS vs. PT 5,966 DEGs 5,077 DEGs ヒト (HS) チンパンジー (PT) HS vs. HS PT vs. PT RM vs. RM RM vs. RM3 アカゲザル (RM),7 DEGs 3,375 DEGs 3,806 DEGs 3,00 DEGs 300

結果の比較 (FDR) HS vs. RM HS vs. PT 統計的手法 (TCC) も多少偽陽性が存在するが倍率変化 (FC) ほど凶悪ではないことがわかるまた高発現側の DEG は FC と比較的よく一致していることがわかる先人が FC のみで比較的信頼性の高い結果を得てきた理由がよくわかる ( 高発現側を信頼するという経験則 ),488 DEGs 578 DEGs ヒト (HS) チンパンジー (PT) HS vs. HS PT vs. PT RM vs. RM RM vs. RM3 アカゲザル (RM) 7 DEGs 6 DEGs 4 DEGs 0 DEGs 30

0,689 genes 3 群間比較このデータは 3 種類の生物種間比較 : ヒト (Homo sapiens; HS) チンパンジー (Pan troglodytes; PT) アカゲザル (Rhesus macaque; RM) どこかの群間で発現変動している遺伝子を検出するやり方を示すヒト (Homo sapiens; HS) チンパンジー (Pan troglodytes; PT) アカゲザル (Rhesus macaque; RM) メス (Female) オス (Male) メスオスメスオス Blekhman et al., Genome Res., 0: 80-89, 00 303

3 群間比較論文 3 群間比較用に特化した手法選択のガイドライン反復ありデータの場合は ( 内部的に edger の関数を用いた )TCC 反復なしの場合は ( 内部的に DESeq を用いた )TCC がおススメ Tang et al., BMC Bioinformatics, 6: 36, 05 304

データ正規化周辺 RPM (Mortazavi et al., Nat. Methods, 5: 6-68, 008) RPKM(Reads per kilobase of exon per million mapped reads) の長さ補正を行わないバージョン Reads Per Million mapped reads の略 TMM 正規化 (Robinson and Oshlack, Genome Biol., : R5, 00) Trimmed Mean of M values の略 edger パッケージに実装されている発現変動遺伝子 (DEG) のデータ正規化時の悪影響を排除すべく M-A plot 上で周縁部にあるデータを使わずに ( トリムして ) 正規化係数を決定する方法 TbT 正規化 (Kadota et al., Algorithms Mol. Biol., 7: 5, 0) TMM 法の改良版で TMM-baySeq-TMM という 3 ステップで正規化を行う方法 st step で得られた TMM 正規化係数を用いて nd step (bayseq) で DEG 同定を行い 3rd step (TMM) では DEG を排除した残りのデータで TMM 正規化 DEG の影響を排除しつつもできるだけ多くの non- DEG データを用いて頑健に正規化係数を決めるという思想 (DEG elimination strategy 提唱論文 ) ideges 正規化 (Sun et al., BMC Bioinformatics, 4: 9, 03) TCC パッケージの原著論文参考発現変動解析 ( サンプル間比較 ) 時に重要となる sequence depth 周辺の正規化法の進展 03.03.07 の HPCI セミナースライドから拝借 DEG elimination strategy (DEGES) を一般化しより高速且つ頑健にしたもの TbT は複製ありのデータのみにしか対応していなかったが複製なしデータにも対応 ideges/edger 正規化法 : 複製ありデータ正規化用 TMM-(edgeR-TMM) n パイプライン ideges/deseq 正規化法 : 複製なしデータ正規化用 DESeq-(DESeq-DESeq) n パイプライン 305

TbT 正規化法参考 TCC パッケージに実装している基本コンセプトの原著論文 TbT 正規化法の説明 04.07. のイルミナウェビナー時のスライドの拝借本来の目的である発現変動遺伝子 (DEG) 自体がデータ正規化時に悪影響を与えるので DEG 候補を除去して正規化を行うほうがよいこと (DEG Elimination Strategy) を提唱した論文既存の正規化法は比較するグループ間で DEG 数に偏りがない (unbiased DE) 場合にはうまく正規化できるが偏りがある場合 (biased DE) にはうまく正規化できないことを示した TbT 法の実体は edger パッケージ中の TMM 正規化法実行 bayseq パッケージ中の DEG 検出法実行および 3DEG 候補を除去した残りの non-deg 候補のみを用いた TMM 正規化法実行の 3 ステップを基本とする TMM-baySeq-TMM パイプライン出力は正規化後の結果 ( 正確には正規化係数 ) なので TbT 正規化後に任意の DEG 検出法を適用することで一連の発現変動解析が終了することになる例えば TbT 正規化法実行後に edger 中の DEG 検出法を適用する一連の手順は TMM-baySeq-TMM-edgeR に相当し原著論文中では edger/tbt と略記している論文中では TbT にした理由を論理的に書いたが本音は ToT に近いものということで TMM と bayseq を採用提案したマルチステップの正規化パイプラインは第および第 3 ステップを繰り返して実行することでより頑健な正規化を実現可能であることも示しているこれが図 3 で説明している iterative TbT approach に相当するものであり TMM-(baySeq-TMM)n とも表現できる例えば iterative TbT 正規化法実行後に edger 中の DEG 検出法を適用する一連の手順は TMM-(baySeq- TMM)n-edgeR に相当する n = 0 の場合は TMM-edgeR となりこれは edger パッケージ中のオリジナルの手順と同じである 306

TCC 参考 TbT 論文の考え方を一般化し R パッケージとしてまとめたという論文 TbT 正規化法の説明 04.07. のイルミナウェビナー時のスライドの拝借 TbT は DEG Elimination Strategy (DEGES; でげす ) に基づく一つの正規化パイプラインにすぎないこと第ステップの bayseq による DEG 同定ステップが律速であり高速化が課題であったことそして各ステップにおいて他の方法が原理的に適用可能であることなどを述べている第ステップの DEG 同定法を edger 中のものに置き換えると TMM-edgeR-TMM という正規化パイプラインになるこれは全て edger パッケージ中の関数のみで成立するため DEGES/edgeR と略記しているまた DEGES 正規化後に edger 中の DEG 同定法を適用する一連の解析手順は DEGES/edgeR-edgeR または TMM-edgeR-TMM-edgeR と表記できるこれは実質的に edger パッケージ中のオリジナルの解析手順を回繰り返して行っていることと同義である ( ただし第 3 ステップの TMM は検出された DEG 候補以外のデータで実行される ) それが実質的に TCC は例えば iterative edger という理解でよいと主張する根拠である TbT 論文中で言及した iterative TbT に相当するものはこの論文中では iterative DEGES ( 略して ideges) と称している例えば ideges/edger-edger は TMM-(edgeR-TMM)n-edgeR に相当する n= は DEGES/edgeR-edgeR に相当する n が以上の場合が ideges に相当するが n の数を増やしてもその分計算コストがかかる一方で実質的に n=3 程度で頭打ちになることを論文中で示しているそれゆえ iterative DEGES のデフォルトは n=3 としている compcoder (Soneson, C., Bioinformatics, 04) 中でもデフォルトはそうなっている 307

0,689 genes 3 群間比較 sample_blekhman_8.txt を入力としてヒト (HS)6 サンプルチンパンジー (PT)6 サンプルアカゲザル (RM)6 サンプルの 3 群間比較を行いますどこかの群間で発現変動している遺伝子を検出するやり方です各群のサンプルは全て別個体です例えばヒトの場合は 6 人分のデータ (6 biological replicates) であり人のサンプルを 6 個に分割したデータ (6 technical replicates) ではありませんヒト (Homo sapiens; HS) チンパンジー (Pan troglodytes; PT) アカゲザル (Rhesus macaque; RM) Blekhman et al., Genome Res., 0: 80-89, 00 308

3 群間比較例題 7 3 入力は sample_blekhman_8.txt 出力は hoge7.txt のみ M-A plot はない 3 Sun et al., BMC Bioinformatics, 4: 9, 03 Tang et al., BMC Bioinformatics, 6: 36, 05 309

3 群間比較ここで各群のサンプル数 ( 列数 ) を指定 FDR 閾値を指定するところだが出力ファイルの読み取り方が分かっていれば気にしなくてもよい Sun et al., BMC Bioinformatics, 4: 9, 03 Tang et al., BMC Bioinformatics, 6: 36, 05 30

3 群間比較コピペで実行した結果の最後の部分を表示約分赤枠部分は 4 種類の FDR 閾値を満たす遺伝子数を表示している例えばは 5% の偽物混入を許容すると 7,47 個が DEG と判定されるということ 3

DEG 数の見積もり FDR 閾値が比較的緩めのところを眺め 0,689 genes 中 8,000 個程度がどこかの群間で発現変動している本物の DEG と判断する 3

出力ファイル解説出力ファイル (hoge7.txt) の中身を解説正規化後のデータ 3 統計解析結果 (p 値 q 値順位情報など ) 4 発現変動順にソートされた結果 4 3 33

コードの解説正規化後のデータを取り出す部分の説明データ正規化 ( 正確には正規化係数を得た ) 後の tcc オブジェクトを入力として getnormalizeddata 関数を用いて正規化後のデータを取得した結果を normalized オブジェクトに格納し 3 出力の一部として組み込んでいる 3 3 34

コードの解説 ( このウェブページではお約束的にそうしているので ) 番左側に行名をしているが実はの result オブジェクトの 3 番左側が同じ情報なのでなくてもよい 3 35

コードの解説 rank( 順位 ) 情報は p 値をベースに計算している p 値が低いほど発現変動の度合いが高いことを意味する実際 3 位の発現パターンはチンパンジー (PT) 群で 4 桁以上のカウント数でそれ以外の群では桁となっており妥当 3 36

ANOVA 的な解析位はアカゲザル (RM) 群で高発現のパターン 3 位はヒト (HS) 群で高発現のパターン ( どの群間で違いがあるかは問わずに ) どこかの群間で発現に差があるものを検出する枠組み (ANOVA 的な解析 ) なのでこのような結果になる 37

コードの解説の段階ではまだ発現変動順にはソートされておらず cbind 関数を用いて列方向で結合した (column bind) 結果を tmp オブジェクトに格納しているだけ実際 result オブジェクトは入力ファイルの遺伝子名の並びになっている 38

コードの解説 ( 赤字のコメント部分をよく見ればわかるが w) 発現変動順にソートしている部分がココ tmp[order(x), ] は x の並びで行をソートするお約束的な書き方発現変動順にソートしたくない場合は行頭に # をつけてコメントアウトすればよい 39

コードの解説の部分が灰色になっているのはあくまでも 0.05, 0.0, 0.0, 0.30 という FDR 閾値を満たす遺伝子数をざっと表示させて概要 (DEG があるかどうかあったとしたらどれくらいか ) を知りたいだけだから 30

コードの解説の部分は出力ファイル中の q.value 列を見ているのと同じ 3

コードの解説例えば 5% FDR threshold というのは q-value (q 値 ) が 0.05 未満 (or 以下 ) という条件判定を行っていることと同義これは有意水準 5% というのが実際の手順として p-value(p 値 ) < 0.05 という条件判定を行っていることを思い出せば納得できる 3

コードの解説番右側の estimateddeg という名前の列が q-value (q 値 ) < 0.05 のところでから 0 に切り替わっているのは param_fdr のところで 0.05 を指定していたからの指定部分を気にしなくてよいと最初に解説したのは q-value 列の取り扱い方法を理解していればそれで充分だからです 33

パターン分類 ( どの群間で違いがあるかは問わずに ) どこかの群間で発現に差があるものを検出する枠組み (ANOVA 的な解析 ) なのでこのような結果になるのはしょうがいないとして G(HS) 群で高発現のものが x 個 G(PT) 群で高発現のものが y 個みたいなものが欲しい! 34

遺伝子クラスタリング ( 個人的には結論として非推奨だが ) 遺伝子間クラスタリングを行っておきどの遺伝子がどのパターンに属するかという情報を MBCluster.Seq などを用いて得ておき特定の FDR 閾値を満たす遺伝子サブセットの分類分けを行えばよいおそらくこれも王道なので一応紹介 36

遺伝子クラスタリング MBCluster.Seq 単体での利用はこちら例題 4 3K-means クラスタリングの一種なので 3 クラスター数を指定する 50 など比較的大きめの値にしても non-redundant にしてくれるので最終的に得られるクラスター数は ( データの性質にもよるが通常は ) 減る 3 37

遺伝子クラスタリング通常のクラスタリングの際は入力データのどの列がどの群かという情報を与えない MBCluster.Seq は同一群内のバラツキを考慮してくれるので列のラベル情報を与えている Model-based clustering というのは ( 発現変動解析時に non-deg の分布を見積もるのと同様に ) 同一群内のバラツキを超えた意味のある発現パターンを見積もって返してくれると理解すればよい 38

コピペ (CTRL + ALT + 左クリックまたはトリプルクリックでコードを全選択して ) コピペ約 8 分 39

結果の説明 800 500 ピクセルの hoge4.png 内部的に乱数を発生させているので見栄えはヒトによって異なる 330

結果の説明下から順に G(HS) G(PT) G3(RM) 順序はのテキストファイル ( 後述 ) と見比べることでわかる 3cluster と 0 は G3 群で高発現 4cluster 6 と 8 は G3 群で低発現パターンのものだみたいに解釈する全体を眺めることで各クラスタを構成するメンバー数 ( 遺伝子数 ) の概要をつかめる例えば 5cluster 9 を構成する遺伝子数が最多とか 3 5 4 33

3 結果の説明テキストファイル (hoge4.txt) の中身番右側の列がクラスター番号情報 3 例題 4 の場合出力ファイルはクラスター番号順にソートされている確かにデンドログラム ( 樹形図 ) でみた通り cluster は G3(RM) 群で高発現パターンになっている 33

コードの解説テキストファイル (hoge4.txt) 作成部分のコード data オブジェクトは入力ファイルを読み込んだ直後のものなので正規化前のデータ 3cls$cluster という数値ベクトルが入力データの遺伝子の並び順にどの遺伝子がどのクラスターに属するかを示した情報に相当する出力がクラスター番号順になっている理由は 4cls$cluster の並びでソートしているから門田亡き後もこのようにコードの中身を自力で解読できるようになっておけば大丈夫 4 3 333

クラスタごとの遺伝子数 hoge4.png を眺めることで cluster 9 を構成する遺伝子数が最多とかが一応わかるが各クラスタを構成するメンバー数 ( 遺伝子数 ) を正確に調べるやり方を伝授 334

クラスタごとの遺伝子数入力データの遺伝子の並び順にどの遺伝子がどのクラスターに属するかを示した cls$cluster という数値ベクトルを入力として table 関数を実行した結果が欲しいものですを眺めると確かに cluster 9 のメンバー数が最多 (7,67 個 ) であることからもこの結果が妥当であることが分かりますこの結果も乱数を発生させているのでヒトによって結果は異なる 3 は理解を助ける補足情報 3 335

TCC + MBCluster.Seq 応用の項目の例題 8 が TCC でどこかの群間で発現変動する遺伝子と MBCluster.Seq を組み合わせたものですコピペ約 0 分 337

TCC + MBCluster.Seq 800 500 ピクセルの hoge8.png 内部的に乱数を発生させているので見栄えはヒトによって異なるこの結果の場合は 3cluster 7 が non-deg パターンであり最多の遺伝子数から構成される 3 338

TCC + MBCluster.Seq テキストファイル (hoge8.txt) の中身 TCC の結果部分は不変一番右側の列の数値はおそらくヒトによって異なる 339

クラスタごとの遺伝子数コードの最後の部分を表示全遺伝子 (0,689 個 ) を対象としたクラスターごとの遺伝子数 png ファイルで眺めた non- DEG パターンに相当する cluster 7 の遺伝子数が最多 (7,43 個 ) となっており妥当 340

クラスタごとの遺伝子数 5% FDR 閾値を満たす遺伝子 (7,47 個 ) に限定してクラスターごとの遺伝子数を再分類 non-deg パターンに相当する cluster 7 の遺伝子数が最少 (67 個 ) となっており妥当 34

共同研究者謝辞所属敬称略他にもバグレポートやプログラム提供をいただいた諸氏 R 本体および有用なパッケージ開発者諸氏に御礼申し上げます m( )m アグリバイオ本体 TCC パッケージおよび手法比較清水謙多郎寺田透三浦文孫建強西山智明湯敏 DDBJ Pipeline Platanus および日本乳酸菌学会誌谷澤靖洋神沼英里中村保一遠野雅徳有田正規伊藤武彦鈴木チセ坂本光央 HPCI 人材養成プログラム杉原稔寺田朋子グラント基盤研究 (C)(H7-9 年度 ): ロングリード時代に対応したトランスクリプトームデータ解析ガイドラインの構築 ( 代表 ) 基盤研究 (C)(H4-6 年度 ): シークエンスに基づく比較トランスクリプトーム解析のためのガイドライン構築 ( 代表 ) 新学術領域研究 ( 研究領域提案型 )(H-6 年度 ): 非モデル生物におけるゲノム解析法の確立 ( 分担 ; 研究代表者 : 西山智明 ) NBDC との共同研究 (H6-7 年度 ):NGS 講習会関連 34