ゲノム情報解析基礎

(1)

ゲノム情報解析基礎

～ Rで塩基配列解析1 ～

1

大学院農学生命科学研究科

アグリバイオインフォマティクス教育研究プログラム

2

微生物科学イノベーション連携研究機構

門田幸二（かどたこうじ）

[email protected]

http://www.iu.a.u-tokyo.ac.jp/~kadota/

講義資料PDFが講義のページからダウンロード可能です。印刷物はありません

(2)

講義予定

2 Apr 23 2018 全てPC使用予定です 

04月16日月曜日（17:15-20:30）



嶋田透：ゲノムからの遺伝子予測



門田幸二：バイオインフォマティクス基礎知識、Rのイントロダクション



04月23日月曜日（17:15-20:30）



門田幸二：Rで塩基配列解析1、multi-FASTAファイルの各種解析



05月07日月曜日（17:15-20:30）



嶋田透：ゲノムアノテーション、遺伝子の機能推定、RNA-seqなどによ

る発現解析、比較ゲノム解析



門田幸二：Rで塩基配列解析2、Rパッケージ、k-mer解析の基礎



05月14日月曜日（17:15-19:00頃）



勝間進：非コードRNA、小分子RNA、エピジェネティクス



講義後、小テスト

(3)

各種ソフトの場所：Excel

4 Apr 23 2018 Windows7の場合。①②③Excel は行列データファイルの確認用 ① ③ ②

(5)

各種ソフトの場所：Excel

Windows10の場合。②よく使うアプリのところが③アルファベット順になっているので、EからはじまるExcelを探して行列データファイルの確認を行ってください ① ② ③ ③ ③

(6)

各種ソフトの場所：エディタ

6 Apr 23 2018 Rコマンドを利用する際のエディタは、 R付属のものを推奨。主目的は二重クォーテーション問題の回避。私は EmEditorというテキストエディタを使用 ②

(7)

各講義科目へのアクセス

①教育プログラム、②各講義のページ、③「ゲノム情報解析基礎」の場合

① _②

(8)

(Rで)塩基配列解析

8 Apr 23 2018 ①「（Rで）塩基配列解析」を用いて講義を進めます。レポート課題で用いるhoge7.fa は、②からダウンロードしてください ② ①

(9)

(遺伝子)アノテーション

遺伝子ごとに、どの染色体のどの座標上に存在するのかなどの情報を含むタブ区切りテキストファイル。①GFF/GTF形式ファイルが有名 ①

(10)

アノテーションツール

10

Apr 23 2018

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

(11)

ツールの一例

①

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

(12)

ツールの一例

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

(13)

GFF/GTF形式ファイルの例

GFF3形式（シロイヌナズナ; TAIR10_GFF3_genes.gff）

GTF形式（ゼブラフィッシュ; Danio_rerio.Zv9.75.gtf）

他にrefFlat形式など様々なファイル形式が存在します。このようなファイルを入力として、任意の（染色体上にある遺伝子の）サブセットを抽出することができます

(14)

解析基礎2

14 Apr 23 2018 目的：アノテーションファイル（annotation.txt）中の第1列目に対して、リストファイル（genelist1.txt）中の文字列と一致する行を抜き出して、hoge1.txtというファイル名で出力したい入力1：アノテーションファイル（annotation.txt）入力2：リストファイル（genelist1.txt）出力：hoge1.txt

(15)

解析基礎2

目的：アノテーションファイル（annotation.txt）中の第1列目に対して、リストファイル（genelist1.txt）中の文字列と一致する行を抜き出して、hoge1.txtというファイル名で出力したい ① ②

(16)

解析基礎2

16 Apr 23 2018 ①作業ディレクトリは「デスクトップ – hoge」。hoge フォルダ中にannotation.txtとgenelist1.txtが存在するという前提。②貸与PCはkadotaではなくstudent ② ① ①

(17)

基本はコピペ

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

(18)

実行結果

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

(19)

実行結果

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

(20)

色の説明

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

(21)

応用

このサンプルコードは①1列目でキーワード検索する場合。別のリストファイルを読み込んで

4列目で検索したい場合のやり方を示します

(22)

解答例

1.

目的のキーワードリストを含むファイルを作成し（例：

list.txt

）

2.

該当箇所を変更し、Rコンソール画面上でコピペ

22 Apr 23 2018 ①メモ帳など任意のエディタでリストファイル(list.txt) を作成し、②「デスクトップ - hoge」フォルダ上で保存 ① ②

(23)

解答例

1.

目的のキーワードリストを含むファイルを作成し（例：

list.txt

）

2.

該当箇所を変更し、Rコンソール画面上でコピペ

一連の作業手順を記述したスクリプトを1 つのファイルとして保存することをお勧め

(24)

警告メッセージ

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

(26)

ありがちなミス1

26 Apr 23 2018 ①作業ディレクトリの変更を忘れているため、②in_f1で指定した最初のファイルの、 ③読み込み段階でエラーが出る。つまり、実際に行ったフォルダ中にはannotation.txt というファイルは存在しないということ ① ② ③

(27)

ありがちなミス2

必要な入力ファイルが作業ディレクトリ中に存在しない。この場合、①in_f2で指定したgenelist1.txtが存在しないため、②それの読み込み段階でエラーが出ている ① ②

(28)

ありがちなミス3

28 Apr 23 2018 ①出力予定のファイル名(hoge1.txt)と同じものをExcel など別のプログラムで開いているため、②最後の write.table関数のところでエラーが出る。対処法は、出力ファイル名を変更するか、開いている別のプログラムを閉じる。MacはExcelで開いていても保存できる ① ②

(29)

ありがちなミス4

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

(30)

コード内部の説明

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

(32)

読み込み

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

(33)

行列data

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

(34)

dimで行数と列数を表示

34 Apr 23 2018 ①オブジェクトdataの行数と列数は11と4。 ②ウェブページ中の表記が灰色なのは、特にやらなくてもいいコマンドだから ① ②

(35)

行列の要素へのアクセス

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

(36)

Tips：上下左右の矢印キー

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

(37)

行列の要素へのアクセス

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

(38)

行列の要素へのアクセス

38 Apr 23 2018 行列dataの要素へのアクセスは [行, 列]。①2列目の情報のみ抽出 ①

(39)

行列の要素へのアクセス

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

(40)

Tips：関数とオプション

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

(41)

Tips：タブ補完

参考列番号を指定する以外にも特定の列を表示するやり方がある。head=TRUEで入力ファイルを読み込むと、列の名前を利用することができる。①subcellular_location列の情報を抽出したい場合は、②「data$su」くらいまで打ち込んでから③Tabキーを押すと、 ① ② ③

(42)

Tips：タブ補完

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

(43)

Tips：table関数

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

(44)

Tips：ソート

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

(45)

Tips：is.element関数

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

(46)

Tips：“二重クォーテーション”

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

(47)

目的をおさらい

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

(48)

コード内部の説明

48 Apr 23 2018 ①論理値ベクトルobjを用いて、② TRUEの要素に対応する行を抽出入力2：リストファイル（genelist1.txt） ① ②

(49)

Tips：昔話

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

(50)

50 Apr 23 2018 genelist1.txt 例題1と12は手順が異なるだけで実質的に同じ。「①～④」と「⑤～⑥」を見比べるとよい ① ② ③ ④ ⑤ ⑥

(51)

入力: annotation.txt

出力: hoge12.txt

ヘッダー行が①ある場合（例題12）

(52)

52 Apr 23 2018 入力: annotation2.txt 出力: hoge13.txt ヘッダー行が②ない場合（例題13）。③ 出力時にヘッダー部分を表示させない ② ③

(53)

ヘッダー行が①ある場合（例題12）と②ない場合（例題 13）の主な違い。③ヘッダー行がない場合は出力ファイルにもヘッダー行は必要ないのでそう明記しているだけ。逆に言えば、デフォルトはcol.names=Tということ ② ①

(54)

FASTA形式

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

(56)

ゲノム配列

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

(57)

基本情報取得

出力: hoge1.txt 入力: hoge4.fa multi-FASTAファイルを読み込んで、トータルの配列長、染色体数（コンティグ数）、配列長の平均、中央値、最大値、最小値、N50、GC含量を計算した結果を返すコードを実行してみよう

(58)

基本情報取得

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

(59)

ダウンロードと確認

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

(60)

コピペ

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

(61)

実行結果

コピペ後に、①list.files()で、②出力ファイルとして指定したファイル名の③hoge1.txtが作成されていることを確認 ① ② ③

(62)

実行結果

62 Apr 23 2018 ①の出力ファイルをテキストエディタやExcelで眺めてもよいが、②オブジェクトtmpの中身を③write.table関数を用いて出力しているだけなので、④R上で眺めている ① ② ② ④ ③

(63)

入出力の関係確認

出力: hoge1.txt 入力: hoge4.fa ①contig_1が最短、contig_2が最長。②N50 の値は65 bpであり、③contig_3の長さと同じ ① ① ② ③

(64)

N

50



アセンブル結果の評価基準の1つ



長いコンティグから足していってTotal_length

の

50 %に達したときの

コンティグの長さ



一般に数値が大きいほどよい

64 Apr 23 2018 contig_2 (103 bp) contig_3 (65 bp) contig_4 (49 bp) contig_1 (24 bp) Total_length (241 bp) Total_length ×0.5 (120.5 bp) 出力: hoge1.txt averageだと外れ値の影響を受けやすく、medianだと短いコンティグが多くを占める場合に不都合らしい

(65)

課題1

左記のコードをhoge4.faの代わりにhoge7.fa を入力として実行し、1. 全配列長（配列長の総和）、2. N50の値、および3. GC含量を示せ hoge7.fa ①

(66)

課題1

66 Apr 23 2018 ①作業ディレクトリの確認、②解析したいファイル(hoge7.fa)の存在確認、③推奨エディタの起動 ③ ① ②

(67)

課題1

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

(68)

課題1

68 Apr 23 2018 変更後のコードを全選択して①コピーし、R Console画面上で②ペースト。全選択後に③を押すやり方でもよい ① ② ③

(69)

コード内部の説明

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

(70)

コード内部の説明

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

(71)

コード内部の説明

width(fasta)に①sum関数を適用すれば、トータルの配列長（配列長の総和）になる。そのものずばりを②の部分で利用 ① ②

(72)

コード内部の説明

72 Apr 23 2018 ①length関数は要素数を返す。この場合、fastaオブジェクトの要素数（つまりコンティグ数）を返す ①

(73)

Tips：条件判定

①50 bp以上のコンティグからなるサブセットの抽出ができそうだ!

(74)

Tips：条件判定

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

(75)

任意の領域の切り出し

76 Apr 23 2018 入力：sample1.fasta 出力：hoge1.fasta >kadota AGTGACGGTCTT >kadota TGACGGT ①subseq関数を用いて、任意の領域の配列を切り出すことができます。②例題1は、 3-9塩基目の配列を切り出すやり方です ① ②

(77)

コピペ

入力：sample1.fasta 出力：hoge1.fasta >kadota AGTGACGGTCTT >kadota TGACGGT ①入力ファイル読み込み直後は、②12 bp。③ そのfastaオブジェクトが、④任意の範囲を抽出するsubseq関数の入力として使われている ③ ① ① ② ④ ③ ④

(78)

コピペ

78 Apr 23 2018 入力：sample1.fasta 出力：hoge1.fasta >kadota AGTGACGGTCTT >kadota TGACGGT ①部分配列の抽出後のfastaオブジェクトは、② 7bp。③そのfastaオブジェクトが、④出力用関数であるwriteXStringSetの入力として使われている ③ ① ① ② ③ ④ ④

(79)

Tips：関数のオプション

入力：sample1.fasta 出力：hoge1.fasta >kadota AGTGACGGTCTT >kadota TGACGGT subseq関数実行時に、①数値を直接指定してもいいし、②オプション名を明記してもよい ① ②

(80)

Tips：関数のオプション

80 Apr 23 2018 入力：sample1.fasta >kadota AGTGACGGTCTT ①原因既知状態でエラーを出す。 ②「3番目の位置から5塩基分抽出」というwidthオプションを利用出力：hoge1.fasta >kadota TGACGGT ① ②

(81)

Tips：関数の使用法

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

(82)

任意の領域の切り出し

82 Apr 23 2018 ① ② ②例題4。入力がmulti-FASTAファイル（hoge4.fa）で、リストファイル（list_sub2.txt）で指定した複数領域を切り出したい場合

(83)

GC含量計算部分の説明

84 Apr 23 2018 ① ②右のサイドバーを下に移動させるとGC含量計算部分を見られる ②

(85)

GC含量計算部分の説明

①fastaオブジェクトを入力として、②配列全体のGC含量（57.68%）を得る部分を解説します

①

(86)

GC含量計算部分の説明

86 Apr 23 2018 ①行頭に#がついている箇所はコメントアウトされている箇所に相当する。実際に消すのではなく、一時的に機能しないようにしている。②が機能していないコメント部分だったことを思い出せば納得できるはず ① ②

(87)

GC含量計算部分の説明

黒枠部分を再度コピペしたのち、①fasta オブジェクトの中身を表示させたところ

(88)

GC含量計算部分の説明

88 Apr 23 2018 ①alphabetFrequency関数は、塩基ごとの出現回数を返す ① ①

(89)

GC含量計算部分の説明

DNA配列上の①Mは「A or C」、②Rは「A or G」などというルールがあるようです ① ② ① ②

(90)

GC含量計算部分の説明

90 Apr 23 2018 ①dim関数は行列の行数と列数を返す。 alphabetFrequency関数出力結果は、4 行×18列からなることが分かる。②最初の4列分のみ抽出するやり方。キーボードの上下キーを上手に利用して最小限の労力でキータイプ（あるいはコピペ）すべし! ① ②

(91)

GC含量計算部分の説明

任意のサブセットを取得可能。 ①2:3やc(1,4)などをうまく利用

(92)

GC含量計算部分の説明

92 Apr 23 2018 黒丸中の数値はcontig_1中のAの数が4個、赤丸中の数値は、contig_4中のTの数が10個であるということ。 ①rowSums関数は行ごとの和を返す。②の記述法でも同じ結果となる ① ① ②

(93)

GC含量計算部分の説明

①rowSums関数の入力として、ACGTのみのカウント数を与えているが、その結果（返り値）は、配列中にNなどを含まない場合は実質的に配列ごとの配列長と同じ。Nを含むhoge7.fa でやってみると違いがわかる。尚、②applyという関数でもrowSumsと同じ結果が得られる ② ② ①

(94)

GC含量計算部分の説明

94 Apr 23 2018 ① ①オブジェクトCG中には、配列ごとのCとGのカウント数が格納されている。②オブジェクトACGT中には、配列ごとのA, C, G, Tのカウント数が格納されている ① ① ②

(95)

GC含量計算部分の説明

例えば49塩基からなるcontig_4中に、①ACGTのいずれかの塩基が49個、②CGの数は25個あることを意味する。③sum関数は、ベクトルの要素の和を返す ① ② ③

(96)

GC含量計算部分の説明

96 Apr 23 2018 ①ここではsum関数を用いて配列全体の総和で GC含量計算をしているが、②CG/ACGTとやると、配列ごとのGC含量を得られる。例えば、 contig_1はCGの数が16個で、ACGTの数が24 個。それゆえGC含量は16/24 = 0.6666667 ① ① ②

(97)

配列ごとのGC含量計算

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

(98)

配列ごとのGC含量計算

98 Apr 23 2018 ①出力ファイル(hoge1.txt)中の②1 番右側の列が配列ごとのGC含量 ① ②

(99)

配列ごとのGC含量計算

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

(100)

課題2

100

Apr 23 2018

hoge7.faの配列ごとのGC含量を示せ

ゲノム情報解析基礎