ゲノム情報解析基礎 ~ R で塩基配列解析 ~ 大学院農学生命科学研究科アグリバイオインフォマティクス教育研究プログラム門田幸二 ( かどたこうじ ) kadota@iu.a.u-tokyo.ac.jp http://www.iu.a.u-tokyo.ac.jp/~kadota/
多くのヒトが感想を述べられて 感想やコメント へのコメントいました ありがとうございます コピペではなく位置から自分が入力するのは無理そう 私もからの入力は無理です そのための多数の例題であり テンプレートを基本として必要最小限の変更で実行するのが基本です 難易度の観点から バイオスタティスティクス基礎論 回目 よりも ゲノム情報解析基礎 回目 が先のほうがありがたい 他の先生方のご都合などで決まりますので難しいところですが 私も個人的にはそれがいいと思います 前向きに検討します 講義全般 Negative( 基礎科目とはいえ 講義が止まることが多すぎ ): 名 気持ちは非常によくわかりますが できるヒト向けの講義ではありません要望 ( もっと応用編の時間を増やしてほしい ): 名 NGSハンズオン講習会を受けましょうw Positive( ちょうどよい わかりやすい ): 多数 ごっつぁんです PythonやMatlabなどに比べてGUIの使いにくさやヘルプの不足が気になった Rを利用するメリットがあれば教えて 個人利用としてはRStudioというソフトがGUIやヘルプの充実という観点からいいようです 講義では取り扱いづらいため 私は使ったことはありませんが スクリプト上で色分けする手段はあるか? あれば教えて 高機能なエディタをおススメ (Windows な私はEmEditor) Linux 上で作業をする人の多くは viやemacsというエディタを使っています 課題 のcontig_8はNを含むが Nは数えるべきなのだろうか? 私も正確なところはよくわかりませんが 多分 Nを除外して考えるのが正解 ではないだろうかと思います contigごとのgc 含量を調べると何がわかるのか気になった contigごとに違いがあるかどうかがわかる とか hoge って何? 特に意味はありません その筋のヒト が何気なしに使う用語です 農学太郎 花子みたいなものです
講義予定 4 月 日月曜日 (7:5-0:30)PC 使用 嶋田透 : ゲノムからの遺伝子予測 門田幸二 : バイオインフォマティクス基礎知識 R のイントロダクション 4 月 8 日月曜日 (7:5-0:30)PC 使用 門田幸二 :R で塩基配列解析 multi-fasta ファイルの各種解析 4 月 5 日月曜日 (7:5-0:30)PC 使用 嶋田透 : ゲノムアノテーション 遺伝子の機能推定 RNA-seq などによる発現解析 比較ゲノム解析 門田幸二 :R で塩基配列解析 R パッケージ k-mer 解析の基礎 5 月 0 日月曜日 (7:5-9:00 頃 )PC 使用 勝間進 : 非コード RNA 小分子 RNA エピジェネティクス 講義後 小テスト 全て PC 使用予定です 3
Contents パッケージ CRAN と Bioconductor 推奨パッケージインストール手順のおさらい ゲノム情報パッケージ BSgenome の概観 ヒトゲノム情報パッケージの解析 連続塩基出現頻度解析 (CpG 解析 ) k-mer 解析 仮想データ 実データ ( 課題 ) 作図 4
パッケージ R 起動直後に? 関数名 と打ち込んでも 使用法を記したウェブページが開かずにエラーが出ることがあります 5
パッケージ??alphabetFrequency と打ち込むように勧められているので打ってみる 検索結果のウェブページが表示されるので それっぽい関数名のところをクリック 6
パッケージ alphabetfrequency 関数は Biostrings というパッケージから提供されているものだと読み解く?? 関数名 は 関数名は既知だがどのパッケージから提供されているものかを知りたい場合などに利用する 7
パッケージ multi-fasta ファイルを読み込んで様々な解析ができるのは Biostrings や seqinr などの塩基配列解析用パッケージのおかげです 3citation( パッケージ名 ) で引用すべき論文がわかります 3 8
パッケージ や の部分でパッケージをロードしている これで ロードしたパッケージが提供する関数群を利用可能になる 9
パッケージ 例えば Biostrings というパッケージを library 関数を用いて読み込むことによって alphabetfrequency のような Biostrings が提供する関数群を利用できるのです ここでは 意図的に library(biostrings) を 回実行して 回目は何も表示されないということを思い出させています 実際には 回のみで大丈夫です 3?alp まで打ってから Tab キーを押すなどして タブ補完 テクを有効利用 3 0
R 本体とパッケージの関係 R 本体とパッケージ の関係は パソコンとソフト Microsoft EXCEL とアドイン Cytoscape とプラグイン のようなものという理解でよい パソコンを購入しただけの状態では できることが限られています 通常は Office やウイルス撃退ソフトなどをインストールして利用します Linux をインストールしただけの状態では できることが限られています 通常は マッピングなど各種プログラムをインストールして利用します R 本体をインストールしただけの状態では できることが限られています 各種解析を行うパッケージ ( またはライブラリ ) をインストールして利用します
CRAN と Bioconductor R パッケージの 大リポジトリ ( 貯蔵庫 ) CRAN:8,000 パッケージ以上 Bioconductor:,04 パッケージ 06 年 04 月 日現在 CRAN (The Comprehensive R Archive Network) 提供パッケージは 生命科学を含む様々な分野で利用される NGS 解析は 3 主に Bioconductor 提供パッケージを利用 3
定期的にバージョンアップ 近年のリリース頻度 R 本体 (http://www.r-project.org/) 06-04-4 に ver. 3..5 をリリース 05-06-8 に ver. 3.. をリリース 04-0-3 に ver. 3.. をリリース 0-03-30 に ver..5.0 をリリース Bioconductor (http://bioconductor.org/) は半年ごとにリリース 05-0 に ver. 3. をリリース (R ver. 3.. で動作確認 ) 提供パッケージ数 :,04 05-04 に ver. 3. をリリース (R ver. 3.. で動作確認 ) 提供パッケージ数 :,04 04-0 に ver. 3.0 をリリース (R ver. 3.. で動作確認 ) 提供パッケージ数 :934 04-04 に ver..4 をリリース (R ver. 3..0 で動作確認 ) 提供パッケージ数 :84 03-0 に ver..3 をリリース (R ver. 3.0 で動作確認 ) 提供パッケージ数 :750 03-04 に ver.. をリリース (R ver. 3.0 で動作確認 ) 提供パッケージ数 :67 0-0 に ver.. をリリース (R ver..5. で動作確認 ) 提供パッケージ数 :608 0-04 に ver..0 をリリース (R ver..5.0 で動作確認 ) 提供パッケージ数 :553 バグの修正や新たな機能がどんどん追加されている 最新版の利用をお勧め 毎年 5 月と 月ごろにバージョンアップするとよいだろう 3
Bioconductor Bioconductor に関する総説 (Review) ゲノム配列やアノテーションパッケージも Bioconductor から提供されており それらに関する言及もあり 4
パッケージのインストール 必要最小限プラスアルファ の推奨インストール手順を行えば (R で ) 塩基配列解析 で利用する多くのパッケージがインストールされます 5
パッケージのインストール これらは CRAN から提供されているものたち バイオスタティスティクス基礎論 で利用予定のパッケージは ここに書き込んでいる 6
パッケージのインストール ゲノム情報のパッケージ (BSgenome ) はBioconductor から提供されています ここでは計 6パッケージをインストールしている 例えばは マウスのmm0というバージョンのゲノム配列情報を含むパッケージの名前 (BSgenome.Mmusculus.UCSC.mm0) に相当する 3 biocliteという関数を用いて該当パッケージをインストールしています 3 7
Contents パッケージ CRAN と Bioconductor 推奨パッケージインストール手順のおさらい ゲノム情報パッケージ BSgenome の概観 ヒトゲノム情報パッケージの解析 連続塩基出現頻度解析 (CpG 解析 ) k-mer 解析 仮想データ 実データ ( 課題 ) 作図 8
BSgenome 利用の意義 ゲノム配列情報は UCSC Ensembl Illumina igenomes などのウェブサイトから取得するのが一般的ではあるが R の生物種ごとに提供されている BSgenome で取得 あるいは取り扱うことも可能 ChIP-seq 用パッケージの MEDIPS は BSgenome を利用 9
BSgenome ゲノム配列 BSgenome 0
BSgenome 黒枠部分のコードをコピペ R ver. 3..3 (Bioconductor ver. 3.) で利用可能な生物種のパッケージ名をリストアップ 83 個あることが分かる R のバージョンが古いとパッケージ数は少なくなる
BSgenome 04 年 4 月リリースのゼブラフィッシュ (Danio rerio; danrer0) のパッケージもある ヒトゲノムはこのあたり 様々なバージョン (hg7, hg8, hg9, hg38) のゲノム配列が提供されていることがわかる
BSgenome 実際にインストール済みのものを調べる この PC 環境では 7 パッケージであることがわかる 3 植物のシロイヌナズナ (Arabidopsis thaliana) のパッケージは 推奨手順通りにインストール作業をしたヒトは存在するはずです 私もインストールされてなかったりしますので なければ個別インストールで対応してください 3 3
個別インストール パッケージの個別インストール方法 パッケージ名部分を変更すれば 基本どのパッケージのインストールにも対応可能 例 : BSgenome.Athaliana.TAIR.TAIR9 4
Contents パッケージ CRAN と Bioconductor 推奨パッケージインストール手順のおさらい ゲノム情報パッケージ BSgenome の概観 ヒトゲノム情報パッケージの解析 連続塩基出現頻度解析 (CpG 解析 ) k-mer 解析 仮想データ 実データ ( 課題 ) 作図 5
BSgenome 例題 9 ヒトゲノム (GRCh38) の R パッケージを入力 3multi-FASTA ファイルを出力として得る 作業ディレクトリはどこでもよいが基本はデスクトップ上の hoge 数分かかるが 約 3.3GB のファイルが生成される ( 動作が遅くなるので ) テキストエディタで開かないで! 3 6
BSgenome 出力ファイルの内容は fasta オブジェクトに格納されている 慣れれば fasta オブジェクトの中身を R 上で直接眺めるほうが全体像をつかみやすい 7
BSgenome ~ 番染色体のみ取扱いたい場合 染色体番号の数が大きくなるほど配列長が短くなっている傾向が一目瞭然 8
BSgenome X, Y, およびミトコンドリア配列も含めたい場合 配列の並びの確認は試行錯誤 3 最初から 5 番目の要素が MT( ミトコンドリア ) だとわかっていたわけではありません 3 9
BSgenome X, Y, およびミトコンドリア配列までのサブセットを hoge0.fasta で保存したい場合 上矢印キーを何回か押してファイルに保存するためのコマンドを出し 3 水色下線部分の か所を変更すればよい 3 3 30
BSgenome 参考 こんな感じで変更して実行 やらなくてよい 実行後に hoge9.fasta よりも若干ファイルサイズの小さい hoge0.fasta が生成されていることが確認できます 決してテキストエディタで開かないで! 3
BSgenome 参考 例題 0 様々な記述形式があります やらなくてよい 3
BSgenome 6 番目以降の配列は ヒトゲノムの一部ではあるものの おそらく割り当てられる染色体が定まっていないものなどです メタゲノム解析などでヒトゲノムにマップされないリードのみ取扱いたい場合には 利用可能な全配列をマッピング時のリファレンスとして用いるのが自然だと思います 33
Contents パッケージ CRAN と Bioconductor 推奨パッケージインストール手順のおさらい ゲノム情報パッケージ BSgenome の概観 ヒトゲノム情報パッケージの解析 連続塩基出現頻度解析 (CpG 解析 ) k-mer 解析 仮想データ 実データ ( 課題 ) 作図 34
R で調べることができます ヒトゲノム中の CpG 出現確率は低い 全部で 6 通りの 連続塩基の出現頻度分布を調べると CG となる確率の実測値 (0.986%) は期待値 (4.%) よりもかなり低い 期待値 ゲノム中の GC 含量を考慮した場合 : 約 4%(A:0.95, C:0.05, G: 0.05, T:0.95) なので 0.05 0.05= 4.% ゲノム中の GC 含量を考慮しない場合 : 50%(A:0.5, C:0.5, G: 0.5, T:0.5) なので 0.5 0.5= 6.5% Lander et al., Nature, 409: 860-9, 00 35
連続塩基の出現頻度 例題 全貌を把握可能な hoge4.fa を作業ディレクトリにダウンロードして実行 36
連続塩基の出現頻度 右クリックでダウンロードし 作業ディレクトリ中に hoge4.fa があることを確認 Mac のヒトは.txt が付与されてしまう拡張子問題の解決も忘れずに! 37
連続塩基の出現頻度 Internet Explorer のヒトは CTRL と ALT キーを押しながらコードの枠内で左クリックすると全選択できます 基本はコピペ 出力ファイルの中身は tmp オブジェクトの中身と同じ 38
連続塩基の出現頻度 出力ファイルは 配列ごと ( この場合コンティグごと ) に 6 種類の 連続塩基の出現頻度をカウントしたものです 出力 :hoge.txt 39
連続塩基の出現確率 出力ファイルは 配列ごと ( この場合コンティグごと ) に 6 種類の 連続塩基の出現確率をカウントしたものです as.prob オプションを TRUE にしているだけ 出力 :hoge.txt 40
Contents パッケージ CRAN と Bioconductor 推奨パッケージインストール手順のおさらい ゲノム情報パッケージ BSgenome の概観 ヒトゲノム情報パッケージの解析 連続塩基出現頻度解析 (CpG 解析 ) k-mer 解析 仮想データ 実データ ( 課題 ) 作図 4
連続塩基の出現確率 例題 7 ヒトゲノムRパッケージを入力とすることもできます 一見ややこしいですが 3fasta オブジェクトの作成までを お約束の手順 だと思えばいいのです ( 孫建強氏提供情報 ) 3 4
連続塩基の出現確率 例題 9 は 例題 7 の記述が気になるヒト用 パッケージ名をベタで書いています 3 の tmp の中身は BSgenome.Hsapiens.NCBI.GRCh38 中で利用可能なオブジェクト名です 3 43
連続塩基の出現確率 出力 :hoge7.txt 例題 7 実行結果ファイル 約 3 分 CG の連続塩基が他に比べて確かに低いことがわかる 44
参考 連続塩基の出現頻度と確率 例題 8 染色体ごとではなく 全てをひとまとめにするやり方です 連続塩基の出現頻度順にソートして CG が少ないことを確かめています 45
k 連続塩基解析 連続塩基の解析は k= のときの k 連続塩基の解析 (k-mer 解析 ) と同じです 比較ゲノム解析 k=3 or 4 付近の値を用いてゲノムごとの頻度情報を取得し 類似性尺度として利用 アセンブル ( ゲノムやトランスクリプトーム ) k=5~00 付近の値を用いて de Bruijn グラフを作成 k-mer 頻度グラフを作成して眺め Heterozygosity の有無などを調査 モチーフ解析 転写開始点の上流配列解析 古細菌の上流 50 塩基に絞って k=4 で出現頻度解析すると おそらく TATA が上位にランクイン 発現量推定 RNA-seq 解析で リファレンスにリードをマップしてリード数をカウントするのが主流だが マッピング作業をすっ飛ばして k-mer に基づく方法で定量 Sailfish (Patro et al., Nat Biotechnol., 04) や RNA-Skim (Zhang and Wang, Bioinformatics, 04) 46
課題 任意の生物種のパッケージについて 連続塩基の出現確率を調べ 得られた結果について簡単に考察せよ ( 例題 7 のヒトゲノムや hoge4.fa を除く ) どのパッケージ ( あるいは生物種 ) を解析し どういう結果 ( 期待値と実測値 ) が得られ 例えばヒトゲノムの場合と比べてどうだったか という程度でよい 47
課題の基本的な考え方 解析する生物種の GC 含量を把握し 期待値からの差分に関する議論が重要 目的 : 連続塩基の出現頻度 (or 確率 ) を調べ 偏りの有無を調査 ヒトゲノムは CG という連続塩基の出現頻度が他 ( 特に CC, GC, GG) に比べて少ないと言われており 大まかにその傾向は確認済み 他の生物種ではどういう傾向にあるのか? ということに興味をもち調べようとしている 注意点 : 生物種ごとに GC 含量が異なる GC 含量が高いということは C と G の出現頻度が高いことを意味する それは A と T の出現頻度の相対的な低下を意味する GC 含量 50% の生物種の場合 A, C, G, T の出現確率は等しい (0.5, 0.5, 0.5, 0.5) それゆえ 計 6 種類の 連続塩基の出現確率の期待値は全て 0.5 0.5 = /6 (AA, AC, AG, AT, CA, CC, CG, CT, GA, GC, GG, GT, TA, TC, TG, TT) (/6, /6, /6, /6, /6, /6, /6, /6, /6, /6, /6, /6, /6, /6, /6, /6) 極端な例として 全て C または G のみからなる GC 含量 00% の生物種の場合 (A, C, G, T) の出現確率は (0.0, 0.5, 0.5, 0.0) となる この 連続塩基出現確率の期待値 : (AA, AC, AG, AT, CA, CC, CG, CT, GA, GC, GG, GT, TA, TC, TG, TT) (0.00, 0.00, 0.00, 0.00, 0.00, 0.5, 0.5, 0.00, 0.00, 0.5, 0.5, 0.00, 0.00, 0.00, 0.00, 0.00) 48
課題の基本的な考え方 目的 : 連続塩基の出現頻度 (or 確率 ) を調べ 偏りの有無を調査 ヒトゲノムはCGという連続塩基の出現頻度が他生物にとって意味のあることなのだろう これ ( 特にGG, CC, GC) に比べて少ないと言われており 大まかにその傾向は確認済み 他の生物種ではどういう傾向にあるのか? ということに興味をもち調べようとしている 注意点 : 生物種ごとに GC 含量が異なる GC 含量 00% の場合は CとGの出現確率はそれぞれ0.5 よって CC, CG, GC, GG の出現確率は全て0.5 0.5 = 0.5となる これが期待値 もし出現確率の実測値が例えばCCのみ高い (or 低い ) だったら 何かその GC 含量が高いということは C と G の出現頻度が高いことを意味する それは A と T の出現頻度の相対的な低下を意味する GC 含量 50% の生物種の場合 A, C, G, T の出現確率は等しい (0.5, 0.5, 0.5, 0.5) それゆえ 計 6 種類の 連続塩基の出現確率の期待値は全て 0.5 0.5 = /6 (AA, AC, AG, AT, CA, CC, CG, CT, GA, GC, GG, GT, TA, TC, TG, TT) (/6, /6, /6, /6, /6, /6, /6, /6, /6, /6, /6, /6, /6, /6, /6, /6) 極端な例として 全て C または G のみからなる GC 含量 00% の生物種の場合 (A, C, G, T) の出現確率は (0.0, 0.5, 0.5, 0.0) となる この 連続塩基出現確率の期待値 : (AA, AC, AG, AT, CA, CC, CG, CT, GA, GC, GG, GT, TA, TC, TG, TT) (0.00, 0.00, 0.00, 0.00, 0.00, 0.5, 0.5, 0.00, 0.00, 0.5, 0.5, 0.00, 0.00, 0.00, 0.00, 0.00) が 差分に関する議論が重要 という意味です 49
GC 含量情報を把握 入力が BSgenome の R パッケージでも GC 含量を計算することができる 例題 50
GC 含量情報を把握 参考 入力が BSgenome のパッケージの場合 入力ファイルというものはない コピペ実行後 ( 約 分 ) に 出力ファイルの中身に相当する tmp を実行し 3GC 含量が約 4% という情報を得る 3 5
ヒトゲノムの結果 ヒトゲノム (BSgenome.Hsapiens.NCBI.GRCh38) の GC 含量は 0.40 だった これは C と G の出現確率の合計が 0.40 ということを意味する それゆえ 各々の確率に分割すると 0.40/ = 0.05 となる 解析したパッケージ名 :BSgenome.Hsapiens.NCBI.GRCh38 ヒトゲノムの全体のGC 含量 : 約 4% 各塩基 (A, C, G, T) の出現確率 : (0.95, 0.05, 0.05, 0.95) 5
ヒトゲノムの結果 解析したパッケージ名 :BSgenome.Hsapiens.NCBI.GRCh38 ヒトゲノムの全体の GC 含量 : 約 4% AとTの出現確率の合計は GC 含量 (0.40) から 0.40 = 0.590となる それゆえ AとT 各々の確率に分割すると 0.590/ = 0.95となる 3 連続塩基の出現確率は 各塩基の出現確率の掛け算で計算可能 AとTの出現確率はともに 0.95 AA, AT, TA, TTの4 種類については その出現確率の期待値 (expected) は どれも0.95 0.95 = 0.08705 ( 約 8.7%) 各塩基 (A, C, G, T) の出現確率 : (0.95, 0.05, 0.05, 0.95) 3 AA, AT, TA, TT の出現確率の期待値 = 0.95 0.95 = 8.7% 53
ヒトゲノムの結果 再び 例題 7 の 3 連続塩基の出現確率計算結果ファイル (hoge7.txt) の考察に戻る 3 54
ヒトゲノムの結果 赤枠の数値が 出力ファイル (hoge7.txt) 中の AA, AT, TA, TT の出現確率の実測値 (observed) 概ね期待値 (8.7%) 周辺の値になっていることがわかる 考察 (discussion) としては 同一種類の連続塩基 (AA and TT) のほうが 異なる種類の連続塩基 (AT and TA) に比べて出現確率が高めである が言えるのでは 55
ヒトゲノムの結果 同じノリで 4 や 5 の残りの連続塩基の出現確率の期待値を計算することができ 実測値と比較可能 解析したパッケージ名 :BSgenome.Hsapiens.NCBI.GRCh38 ヒトゲノムの全体のGC 含量 : 約 4% 各塩基 (A, C, G, T) の出現確率 : (0.95, 0.05, 0.05, 0.95) 3 AA, AT, TA, TTの出現確率の期待値 = 0.95 0.95 = 8.7% 4 CC, CG, GC, GGの出現確率の期待値 = 0.05 0.05 = 4.% 5 AC, AG, CA, CT, GA, GT, TC, TGの出現確率の期待値 = 0.05 0.95 = 6.0% 56
ヒトゲノムの結果 解析したパッケージ名 :BSgenome.Hsapiens.NCBI.GRCh38 ヒトゲノムの全体の GC 含量 : 約 4% 各塩基 (A, C, G, T) の出現確率 : (0.95, 0.05, 0.05, 0.95) 3 AA, AT, TA, TT の出現確率の期待値 = 0.95 0.95 = 8.7% 4 CC, CG, GC, GG の出現確率の期待値 = 0.05 0.05 = 4.% 例えば CC, CG, GC, GGの出現確率の期待値は 4.% 考察 : 同一種類の連続塩基 (CC and GG) のほうが 異なる種類の連続塩基 (CG and GC) に比べて出現確率が高めである という傾向は確かにありそうだ 考察 :CGの出現確率の実測値 ( 約.0%) は 期待値 ( 約 4.%) よりもかなり低い 染色体ごとに分かれている場合は 例えばbox plotで全体像を眺める 5 AC, AG, CA, CT, GA, GT, TC, TG の出現確率の期待値 = 0.05 0.95 = 6.0% 57
Contents パッケージ CRANとBioconductor 推奨パッケージインストール手順のおさらい ゲノム情報パッケージBSgenomeの概観 ヒトゲノム情報パッケージの解析 連続塩基出現頻度解析 (CpG 解析 ) k-mer 解析 仮想データ 実データ ( 課題 ) 作図 58
作図 (box plot): 基本形 例題 0 box plot を PNG 形式ファイルで出力するやり方の基本形 59
400 pixels 作図 (box plot): 基本形 PNG ファイルのサイズを指定するところ hoge0.png 700 pixels 60
作図 (box plot): 色づけ 3 例題 color という列名のところに 連続塩基の種類ごとに色を指定した 3 タブ区切りファイル (human_mer.txt) を与えて利用 4 このファイルの情報を利用しているのは コードの下のほう 4 6
作図 (box plot): 色づけ boxplot 関数実行時のcolオプション部分で3 color 列の情報を利用していることがわかる 4 expected 列情報は 例題 では利用していない 3 4 6
400 pixels 作図 (box plot): 色づけ CG の出現確率が期待値 (4.%) より少ないのは CC, 3GC, 4 GG との相対的な関係からも明白 hoge.png 4 3 700 pixels 63
作図 (box plot): 発展形 期待値との差分を評価すべく 縦軸を log( 観測値 / 期待値 ) としてプロット 0 付近にある 連続塩基は 観測値 ( 実測された出現確率 ) が期待値とほぼ同じことを意味する この縦軸のような表現方法は一般的です hoge.png 64