ゲノム情報解析基礎

Transcription

1 講義資料 PDF が講義のページからダウンロード可能です 印刷物はありません ゲノム情報解析基礎 ~ R で塩基配列解析 ~ 大学院農学生命科学研究科アグリバイオインフォマティクス教育研究プログラム 微生物科学イノベーション連携研究機構門田幸二 ( かどたこうじ ) kadota@iu.a.u-tokyo.ac.jp

2 講義予定 04 月 6 日月曜日 (7:5-0:30) 嶋田透 : ゲノムからの遺伝子予測 門田幸二 : バイオインフォマティクス基礎知識 R のイントロダクション 04 月 3 日月曜日 (7:5-0:30) 門田幸二 :R で塩基配列解析 multi-fasta ファイルの各種解析 05 月 07 日月曜日 (7:5-0:30) 嶋田透 : ゲノムアノテーション 遺伝子の機能推定 RNA-seq などによる発現解析 比較ゲノム解析 門田幸二 :R で塩基配列解析 R パッケージ k-mer 解析の基礎 05 月 4 日月曜日 (7:5-9:00 頃 ) 勝間進 : 非コード RNA 小分子 RNA エピジェネティクス 講義後 小テスト 全て PC 使用予定です

3 Contents 行列形式ファイルの解析基礎 ( アノテーションファイルを例に ) 例題をテンプレートとして任意の解析を行う基本手順 入力ファイルの最後の改行の有無 ありがちなミスとエラーメッセージ コード内部の説明 ( 行列演算の基礎 ) multi-fasta ファイルからの各種情報抽出 基本情報取得 ( コンティグ数 配列長 N50 GC 含量 ) 任意の領域の切り出し GC 含量計算部分の説明 3

4 各種ソフトの場所 :Excel Windows7 の場合 3Excel は行列データファイルの確認用 3 4

5 各種ソフトの場所 :Excel Windows0 の場合 よく使うアプリのところが 3 アルファベット順になっているので E からはじまる Excel を探して行列データファイルの確認を行ってください

6 各種ソフトの場所 : エディタ R コマンドを利用する際のエディタは R 付属のものを推奨 主目的は二重クォーテーション問題の回避 私は EmEditor というテキストエディタを使用 6

7 各講義科目へのアクセス 教育プログラム 各講義のページ 3 ゲノム情報解析基礎 の場合 3 7

8 (R で ) 塩基配列解析 (R で ) 塩基配列解析 を用いて講義を進めます レポート課題で用いる hoge7.fa は からダウンロードしてください 8

9 ( 遺伝子 ) アノテーション 遺伝子ごとに どの染色体のどの座標上に存在するのかなどの情報を含むタブ区切りテキストファイル GFF/GTF 形式ファイルが有名 9

10 アノテーションツール PubMed で annotation tool で検索した結果 308 年 3 月 30 日のときは,898 件ヒット これらの多くは自動アノテーションプログラムであり 入力がゲノム配列 ( 通常は後述する FASTA 形式ファイル ) 出力がアノテーション結果 (GFF 形式ファイルを含む ) である 3 0

11 ツールの一例 DFAST は 乳酸菌など原核生物 (prokaryotes) 用のアノテーションツール ゲノム配列決定とアノテーションは密接に関連している 原核生物 = バクテリア + 古細菌

12 ツールの一例 ここでも書かれているが オリジナルの DFAST はウェブツールであり その日本語解説もある DFAST を利用すると 公共塩基配列データベースの つである DDBJ への登録も便利です

13 GFF/GTF 形式ファイルの例 GFF3 形式 ( シロイヌナズナ ; TAIR0_GFF3_genes.gff) 他に refflat 形式など様々なファイル形式が存在します このようなファイルを入力として 任意の ( 染色体上にある遺伝子の ) サブセットを抽出することができます GTF 形式 ( ゼブラフィッシュ ; Danio_rerio.Zv9.75.gtf) 3

14 解析基礎 目的 : アノテーションファイル (annotation.txt) 中の第 列目に対して リストファイル (genelist.txt) 中の文字列と一致する行を抜き出して hoge.txt というファイル名で出力したい 入力 : アノテーションファイル (annotation.txt) 出力 :hoge.txt 入力 : リストファイル (genelist.txt) 4

15 解析基礎 目的 : アノテーションファイル (annotation.txt) 中の第 列目に対して リストファイル (genelist.txt) 中の文字列と一致する行を抜き出して hoge.txt というファイル名で出力したい 5

16 解析基礎 作業ディレクトリは デスクトップ hoge hoge フォルダ中に annotation.txt と genelist.txt が存在するという前提 貸与 PC は kadota ではなく student 6

17 基本はコピペ 一連のコマンド群をコピーして作業ディレクトリは デスクトッ R Console 画面上でペースト ブラウザがInternet プ Explorer hoge hoge の場合は フォルダ中 CTRL とALTキーを押しながらコードの枠内で左クリックすると にannotation.txtとgenelist.txt 全選択できます トリプルクリックでも全選択可能が存在するという前提 7

18 実行結果 コピペ実行後に list.files() 出力ファイル名として指定した hoge.txt が生成されているのがわかる list.files() で表示される結果 と 実行後の hoge フォルダの中身 は当然同じ 実行前の hoge フォルダ 実行後の hoge フォルダ 8

19 実行結果 out というオブジェクトの中身を write.table という関数でファイルに出力しています この場合 出力ファイル hoge.txt の中身は R コンソール画面中で out と打ち込むことで見られる 実行後の hoge フォルダ 9

20 色の説明 R コード中の色の使い分けについて説明します はじめに のところに の情報があります 0

21 応用 このサンプルコードは 列目でキーワード検索する場合 別のリストファイルを読み込んで 4 列目で検索したい場合のやり方を示します

22 解答例. 目的のキーワードリストを含むファイルを作成し ( 例 :list.txt). 該当箇所を変更し Rコンソール画面上でコピペ メモ帳など任意のエディタでリストファイル (list.txt) を作成し デスクトップ - hoge フォルダ上で保存

23 解答例. 目的のキーワードリストを含むファイルを作成し ( 例 :list.txt). 該当箇所を変更し Rコンソール画面上でコピペ 一連の作業手順を記述したスクリプトを つのファイルとして保存することをお勧め 3

24 Contents 行列形式ファイルの解析基礎 ( アノテーションファイルを例に ) 例題をテンプレートとして任意の解析を行う基本手順 入力ファイルの最後の改行の有無 ありがちなミスとエラーメッセージ コード内部の説明 ( 行列演算の基礎 ) multi-fasta ファイルからの各種情報抽出 基本情報取得 ( コンティグ数 配列長 N50 GC 含量 ) 任意の領域の切り出し GC 含量計算部分の説明 4

25 警告メッセージ list.txt ファイル作成時に membrane と打った後に改行を 入れた場合と 入れない場合の挙動の違いを把握し 後学のために警告メッセージの意味を理解しておくとよい この場合は結果には影響していないことがわかる R は警告メッセージの記述内容が比較的分かりやすいのでよく読むべし 5

26 ありがちなミス 作業ディレクトリの変更を忘れているため in_f で指定した最初のファイルの 3 読み込み段階でエラーが出る つまり 実際に行ったフォルダ中には annotation.txt というファイルは存在しないということ 3 6

27 ありがちなミス 必要な入力ファイルが作業ディレクトリ中に存在しない この場合 in_f で指定した genelist.txt が存在しないため それの読み込み段階でエラーが出ている 7

28 ありがちなミス 3 出力予定のファイル名 (hoge.txt) と同じものを Excel など別のプログラムで開いているため 最後の write.table 関数のところでエラーが出る 対処法は 出力ファイル名を変更するか 開いている別のプログラムを閉じる Mac は Excel で開いていても保存できる 8

29 ありがちなミス 4 実行スクリプトをコピーする際 最後の行のところで改行を含ませずに 3R Console 画面上でペーストしたため 最後のコマンドが実行されない ( 出力ファイルが生成されない ) これも比較的ありがちなパターンです コピペ後に無意識にリターンキーを押すことを心がけるだけでもよいでしょう 3 9

30 Contents 行列形式ファイルの解析基礎 ( アノテーションファイルを例に ) 例題をテンプレートとして任意の解析を行う基本手順 入力ファイルの最後の改行の有無 ありがちなミスとエラーメッセージ コード内部の説明 ( 行列演算の基礎 ) multi-fasta ファイルからの各種情報抽出 基本情報取得 ( コンティグ数 配列長 N50 GC 含量 ) 任意の領域の切り出し GC 含量計算部分の説明 30

31 コード内部の説明 コードの中身を説明します 黒枠部分を再度コピペ 3

32 読み込み in_f で指定したファイルを読み込め 読み込むファイルの最初の行はヘッダー部分 3 ファイルの区切り文字はタブです 4 読み込んだ結果を data という名前で取り扱う 4 3 3

33 行列 data data と打ってリターン 入力ファイルの中身を正しく読み込めていることがわかる header=true としているので 3 このように見えて列名として認識される 3 33

34 dim で行数と列数を表示 オブジェクト data の行数と列数は と 4 ウェブページ中の表記が灰色なのは 特にやらなくてもいいコマンドだから 34

35 行列の要素へのアクセス 行列 data の要素へのアクセスは [ 行, 列 ] humei は 読み込み元ファイルの annotation.txt 中では 7 行 4 列目だが 行目をヘッダー行としているので 3 6 行 4 列目とする必要がある 利用例は ファイル読み込み時に x 行 y 列目に不具合がある のようなエラーが出た時のトラブルシューティングなど 3 35

36 Tips: 上下左右の矢印キー 上矢印キーを押すと 直前に打ったコマンドが表示される 最初から全部打ち直すのではなく 上下左右の矢印キーを有効に利用し最小限の労力で打つべし! 36

37 行列の要素へのアクセス 行列 data の要素へのアクセスは [ 行, 列 ] 行目の情報のみ抽出 読み込み時に head=true としていたので ヘッダー行がついている 37

38 行列の要素へのアクセス 行列 data の要素へのアクセスは [ 行, 列 ] 列目の情報のみ抽出 38

39 行列の要素へのアクセス 行列 dataの要素へのアクセスは [ 行, 列 ] param 列目の情報のみ抽出 paramには という数値が代入されていたのでこうなる 39

40 Tips: 関数とオプション 参考 行列 data の最初の数行を表示したい場合は head 関数を利用 n=3 というオプションを利用すると 最初の 3 行分のみ表示 関数ごとに様々なオプションを利用可能 40

41 Tips: タブ補完 参考 列番号を指定する以外にも特定の列を表示するやり方がある head=true で入力ファイルを読み込むと 列の名前を利用することができる subcellular_location 列の情報を抽出したい場合は data$su くらいまで打ち込んでから 3Tab キーを押すと 3 4

42 Tips: タブ補完 参考 列名中の su からはじまる文字列を補完して表示してくれる Tab キーを用いた補完機能 という意味で タブ補完 という このテクニックは Linux でも利用可能 4

43 Tips:table 関数 参考 hoge の中身は subcellular_location 列の情報からなる 3table 関数は ベクトル中の要素ごとの出現回数を返す NGS データ中の特定のリードの出現回数 や アノテーションファイル中の染色体ごとの遺伝子数 など 様々な局面で利用可能 3 43

44 Tips: ソート 参考 sort 関数と併用することで全体像を俯瞰可能 例えば nuclear に局在する遺伝子数が最も多く 4 個であった などが簡単にわかる 44

45 Tips:is.element 関数 is.element 関数は hoge ベクトルに対して nuclear の文字が存在する場所を TRUE 存在しない場所を FALSE として返す as.character 関数は 文字列ベクトルとして取り扱いたい場合に利用 45

46 Tips: 二重クォーテーション 二重クォーテーションが自動で変更されるエディタは非推奨です 日本語の二重クォーテーションもだめです Microsoft Word や PDF ファイル中のコードのコピペ時によくハマります 46

47 目的をおさらい 目的は 数万 ~ 数百万行からなるファイルを読み込んで特定のキーワードを含む行のみ取り出すテクニックの習得 is.element や as.character はここで利用 47

48 コード内部の説明 論理値ベクトル obj を用いて TRUE の要素に対応する行を抽出 入力 : リストファイル (genelist.txt) 48

49 Tips: 昔話 コード作成当時は as.character 関数を用いてデータの型を文字列ベクトルに揃えていた 少なくとも現在 (R ver 以降 ) は この関数がなくても大丈夫なようだ 同じ関数でもバージョンによって挙動が異なる ( バージョンの違いの一例 ) 49

50 例題 と は手順が異なるだけで実質的に同じ ~4 と 5~6 を見比べるとよい genelist.txt

51 ヘッダー行が ある場合 ( 例題 ) 入力 : annotation.txt 出力 : hoge.txt 5

52 ヘッダー行が ない場合 ( 例題 3) 3 出力時にヘッダー部分を表示させない 入力 : annotation.txt 3 出力 : hoge3.txt 5

53 ヘッダー行が ある場合 ( 例題 ) と ない場合 ( 例題 3) の主な違い 3 ヘッダー行がない場合は出力ファイルにもヘッダー行は必要ないのでそう明記しているだけ 逆に言えば デフォルトは col.names=t ということ 3 53

54 Contents 行列形式ファイルの解析基礎 ( アノテーションファイルを例に ) 例題をテンプレートとして任意の解析を行う基本手順 入力ファイルの最後の改行の有無 ありがちなミスとエラーメッセージ コード内部の説明 ( 行列演算の基礎 ) multi-fasta ファイルからの各種情報抽出 基本情報取得 ( コンティグ数 配列長 N50 GC 含量 ) 任意の領域の切り出し GC 含量計算部分の説明 54

55 FASTA 形式 R で multi-fasta ファイルを読み込んで自在に解析できます ゲノム配列解析 FASTA 形式ファイルの解析 ここでは全体像を完全に把握すべく hoge4.fa ファイル ( ダウンロードは後程 ) を仮想ゲノム配列ファイルとして取り扱う 55

56 ゲノム配列 実際のゲノム配列はここからも取得可能 R で染色体ごとの配列長や GC 含量の計算ができる 56

57 基本情報取得 multi-fasta ファイルを読み込んで トータルの配列長 染色体数 ( コンティグ数 ) 配列長の平均 中央値 最大値 最小値 N50 GC 含量を計算した結果を返すコードを実行してみよう 入力 : hoge4.fa 出力 : hoge.txt 57

58 基本情報取得 FASTA 形式 基本情報を取得 コードの最初のほうに入力ファイルと出力ファイルを記述するので コピペ実行結果としてどういう名前のファイルが出力されるべきかわかる 58

59 ダウンロードと確認 例題の入力ファイル (hoge4.fa) をダウンロード R 上で作業ディレクトリの確認 3 貸与 PC は ここが student 4 作業ディレクトリに解析したい入力ファイルがあることを確認 他のファイルがあってもよい

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

61 実行結果 コピペ後に list.files() で 出力ファイルとして指定したファイル名の 3hoge.txt が作成されていることを確認 3 6

62 実行結果 の出力ファイルをテキストエディタや Excel で眺めてもよいが オブジェクト tmp の中身を 3write.table 関数を用いて出力しているだけなので 4R 上で眺めている 3 4 6

63 入出力の関係確認 contig_ が最短 contig_ が最長 N50 の値は 65 bp であり 3contig_3 の長さと同じ 入力 : hoge4.fa 3 出力 : hoge.txt 63

64 N50 average だと外れ値の影響を受けやすく median だと短いコンティグが多くを占める場合に不都合らしい アセンブル結果の評価基準の つ 長いコンティグから足していって Total_length の 50% に達したときのコンティグの長さ 一般に数値が大きいほどよい 出力 : hoge.txt contig_ (03 bp) Total_length 0.5 (0.5 bp) contig_3 (65 bp) contig_4 (49 bp) contig_ (4 bp) Total_length (4 bp) 64

65 課題 左記のコードを hoge4.fa の代わりに hoge7.fa を入力として実行し. 全配列長 ( 配列長の総和 ). N50 の値 および 3. GC 含量を示せ hoge7.fa 65

66 課題 作業ディレクトリの確認 解析したいファイル (hoge7.fa) の存在確認 3 推奨エディタの起動 3 66

67 課題 テンプレートコードをコピペし 入力ファイル名部分を変更 67

68 課題 変更後のコードを全選択して コピーし R Console 画面上で ペースト 全選択後に 3 を押すやり方でもよい 3 68

69 コード内部の説明 コードの中身を説明します 黒枠部分を再度コピペ 69

70 コード内部の説明 入力ファイル情報を格納したものが (DNAStringSet という形式で保持されている ) fasta オブジェクト width の位置にあるのがコンティグごとの配列長情報 配列長情報は width(fasta) で数値ベクトルとして抽出可能 70

71 コード内部の説明 width(fasta) に sum 関数を適用すれば トータルの配列長 ( 配列長の総和 ) になる そのものずばりを の部分で利用 7

72 コード内部の説明 length 関数は要素数を返す この場合 fasta オブジェクトの要素数 ( つまりコンティグ数 ) を返す 7

73 Tips: 条件判定 50 bp 以上のコンティグからなるサブセットの抽出ができそうだ! 73

74 Tips: 条件判定 参考 サブセット抽出テクニックはここで使っている コードの中身が分かると応用範囲が飛躍的に増大 一定以上のスキルをもつバイオインフォマティシャンは 例題を探すよりも自分で作るヒトのほうが多いかも 74

75 Contents 行列形式ファイルの解析基礎 ( アノテーションファイルを例に ) 例題をテンプレートとして任意の解析を行う基本手順 入力ファイルの最後の改行の有無 ありがちなミスとエラーメッセージ コード内部の説明 ( 行列演算の基礎 ) multi-fasta ファイルからの各種情報抽出 基本情報取得 ( コンティグ数 配列長 N50 GC 含量 ) 任意の領域の切り出し GC 含量計算部分の説明 75

76 任意の領域の切り出し subseq 関数を用いて 任意の領域の配列を切り出すことができます 例題 は 3-9 塩基目の配列を切り出すやり方です 入力 :sample.fasta >kadota AGTGACGGTCTT 出力 :hoge.fasta >kadota TGACGGT 76

77 コピペ 入力ファイル読み込み直後は bp 3 その fasta オブジェクトが 4 任意の範囲を抽出する subseq 関数の入力として使われている 入力 :sample.fasta >kadota AGTGACGGTCTT 出力 :hoge.fasta >kadota TGACGGT 77

78 コピペ 部分配列の抽出後の fasta オブジェクトは 7bp 3 その fasta オブジェクトが 4 出力用関数である writexstringset の入力として使われている 3 4 入力 :sample.fasta >kadota AGTGACGGTCTT 出力 :hoge.fasta >kadota TGACGGT

79 Tips: 関数のオプション subseq 関数実行時に 数値を直接指定してもいいし オプション名を明記してもよい 入力 :sample.fasta >kadota AGTGACGGTCTT 出力 :hoge.fasta >kadota TGACGGT 79

80 Tips: 関数のオプション 原因既知状態でエラーを出す 3 番目の位置から 5 塩基分抽出 という width オプションを利用 入力 :sample.fasta >kadota AGTGACGGTCTT 出力 :hoge.fasta >kadota TGACGGT 80

81 Tips: 関数の使用法? 関数名 で使用法を記したウェブページが開く ページの下のほうに 大抵の場合使用例が掲載されている 使用法既知の関数のマニュアルをいくつか読んで慣れておく 8

82 任意の領域の切り出し 例題 4 入力が multi-fasta ファイル (hoge4.fa) で リストファイル (list_sub.txt) で指定した複数領域を切り出したい場合 8

83 Contents 行列形式ファイルの解析基礎 ( アノテーションファイルを例に ) 例題をテンプレートとして任意の解析を行う基本手順 入力ファイルの最後の改行の有無 ありがちなミスとエラーメッセージ コード内部の説明 ( 行列演算の基礎 ) multi-fasta ファイルからの各種情報抽出 基本情報取得 ( コンティグ数 配列長 N50 GC 含量 ) 任意の領域の切り出し GC 含量計算部分の説明 83

84 GC 含量計算部分の説明 右のサイドバーを下に移動させると GC 含量計算部分を見られる 84

85 GC 含量計算部分の説明 fastaオブジェクトを入力として 配列全体のGC 含量 (57.68%) を得る部分を解説します 85

86 GC 含量計算部分の説明 行頭に # がついている箇所はコメントアウトされている箇所に相当する 実際に消すのではなく 一時的に機能しないようにしている が機能していないコメント部分だったことを思い出せば納得できるはず 86

87 GC 含量計算部分の説明 黒枠部分を再度コピペしたのち fasta オブジェクトの中身を表示させたところ 87

88 GC 含量計算部分の説明 alphabetfrequency 関数は 塩基ごとの出現回数を返す 88

89 GC 含量計算部分の説明 DNA 配列上の M は A or C R は A or G などというルールがあるようです 89

90 GC 含量計算部分の説明 dim 関数は行列の行数と列数を返す alphabetfrequency 関数出力結果は 4 行 8 列からなることが分かる 最初の4 列分のみ抽出するやり方 キーボードの上下キーを上手に利用して最小限の労力でキータイプ ( あるいはコピペ ) すべし! 90

91 GC 含量計算部分の説明 任意のサブセットを取得可能 :3 や c(,4) などをうまく利用 9

92 GC 含量計算部分の説明 黒丸中の数値はcontig_ 中のAの数が4 個 赤丸中の数値は contig_4 中のTの数が0 個であるということ rowsums 関数は行ごとの和を返す の記述法でも同じ結果となる 9

93 GC 含量計算部分の説明 rowsums 関数の入力として ACGTのみのカウント数を与えているが その結果 ( 返り値 ) は 配列中にNなどを含まない場合は実質的に配列ごとの配列長と同じ Nを含むhoge7.fa でやってみると違いがわかる 尚 applyという関数でもrowsumsと同じ結果が得られる 93

94 GC 含量計算部分の説明 オブジェクトCG 中には 配列ごとのCとGのカウント数が格納されている オブジェクトACGT 中には 配列ごとのA, C, G, Tのカウント数が格納されている 94

95 GC 含量計算部分の説明 例えば49 塩基からなるcontig_4 中に ACGTのいずれかの塩基が49 個 CGの数は5 個あることを意味する 3sum 関数は ベクトルの要素の和を返す 3 95

96 GC 含量計算部分の説明 ここではsum 関数を用いて配列全体の総和で GC 含量計算をしているが CG/ACGTとやると 配列ごとのGC 含量を得られる 例えば contig_はcgの数が6 個で ACGTの数が4 個 それゆえGC 含量は6/4 =

97 配列ごとの GC 含量計算 sum 関数を用いずに CG/ACGT とやって 配列ごとの GC 含量を得るための項目 記述内容がほぼ同じことが分かる 97

98 配列ごとの GC 含量計算 出力ファイル (hoge.txt) 中の 番右側の列が配列ごとの GC 含量 98

99 配列ごとの GC 含量計算 ACGT 列は4 種類の塩基のみの出現数 Length 列は配列長情報を表す 配列長は ACGT 以外の全てを含むので その差分 (Length - ACGT) がNなどのACGT 以外の塩基のトータルの出現回数ということになる Length ACGTという関係 99

100 課題 hoge7.fa の配列ごとの GC 含量を示せ hoge7.fa 00