ゲノム情報解析基礎 ～ Rで塩基配列解析1 ～ 大学院農学生命科学研究科 アグリバイオインフォマティクス教育研究プログラム 門田幸二（かどた こうじ） [email protected] http://www.iu.a.u-tokyo.ac.jp/~kadota/ Apr 14 2015 1 全てPC使用予定です 講義日程 4月 7日火曜日（17:15-20:30）PC使用 嶋田透：ゲノムからの遺伝子予測 門田幸二：バイオインフォマティクス基礎知識 4月14日火曜日（17:15-20:30）PC使用 門田幸二：Rで塩基配列解析1 4月21日火曜日（17:15-20:30）PC使用 嶋田透：ゲノムアノテーション、遺伝子の機能推定、EST、 RNA-seqなどによる発現解析 門田幸二：Rで塩基配列解析2 4月28日火曜日（17:15-19:00頃）PC使用 Apr 14 2015 勝間進：非コードRNA、小分子RNA、エピジェネティクス 講義後、小テスト 2 Contents 行列形式ファイルの解析基礎（アノテーションファイルを例に） 例題をテンプレートとして任意の解析を行う基本手順 ありがちなミスとエラーメッセージ 入力ファイルの最後の改行の有無 コード内部の説明（行列演算の基礎） multi-FASTAファイルからの各種情報抽出 Apr 14 2015 基本情報取得（コンティグ数、配列長、N50、GC含量） 任意の領域の切り出し GC含量計算部分の説明 3 各種ソフトの場所 Apr 14 2015 Excelは行列データファイルの確認用。 エディタはR付属のものを推奨。二重 クォーテーション問題を回避するのが主 な目的。門田はEmEditorを使っています。 4 「ゲノム情報解析基礎」の場合 各講義科目へのアクセス ① ② ③ Apr 14 2015 5 (Rで)塩基配列解析 Apr 14 2015 「（Rで）塩基配列解析」を 用いて講義を進めます。 レポート課題で用いる sample4.fastaはここから ダウンロードしてください。 6 (遺伝子)アノテーション Apr 14 2015 遺伝子ごとに、どの染色体のどの座 標上に存在するのかなどの情報を含 むタブ区切りテキストファイル。 GFF/GTF形式ファイルが有名です。 7 GFF/GTF形式ファイルの例 GFF3形式（シロイヌナズナ; TAIR10_GFF3_genes.gff） 他にrefFlat形式など様々な ファイル形式が存在します。 このようなファイルを入力と して、任意の（染色体上にあ る遺伝子の）サブセットを抽 出することができます。 GTF形式（ゼブラフィッシュ; Danio_rerio.Zv9.75.gtf） Apr 14 2015 8 解析基礎2 目的：アノテーションファイル （annotation.txt）中の第1列目に対して、 リストファイル（genelist1.txt）中の文字 列と一致する行を抜き出して、 hoge1.txtというファイル名で出力したい 入力1：アノテーションファイル（annotation.txt） 出力：hoge1.txt 入力2：リストファイル（genelist1.txt） Apr 14 2015 9 解析基礎2 Apr 14 2015 目的：アノテーションファイル （annotation.txt）中の第1列目に対して、 リストファイル（genelist1.txt）中の文字 列と一致する行を抜き出して、 hoge1.txtというファイル名で出力したい 10 解析基礎2 Apr 14 2015 作業ディレクトリは「デスクトッ プ – hoge」。hogeフォルダ中 にannotation.txtとgenelist1.txt が存在するという前提。貸与 PCは黒矢印部分が「kadota」 ではなく「iu」。 11 基本はコピペ ①一連のコマンド群をコピーして 作業ディレクトリは「デスクトッ ②RプConsole画面上でペースト。 – hoge」。hogeフォルダ中 ブラウザがInternet にannotation.txtとgenelist1.txt Explorerの場 合は、CTRLとALTキーを押しな が存在するという前提。 がらコードの枠内で左クリックす ると、全選択できます。 ① ② Apr 14 2015 12 実行結果 「list.files()で表示される結 果」と「実行後のhogeフォル ダの中身」は当然同じです 実行前のhogeフォルダ 実行後のhogeフォルダ Apr 14 2015 13 実行結果 outというオブジェクトの中 身をwrite.tableという関数で ファイルに出力しています。 この場合、出力ファイル (hoge1.txt)の中身は、Rコン ソール画面中で「out」と打 ち込むことでも見られます。 実行後のhogeフォルダ Apr 14 2015 14 色の説明 Apr 14 2015 Rコード中の色の使い分け について説明します。 15 応用 Apr 14 2015 このサンプルコードは1列 目でキーワード検索する場 合。別のリストファイルを読 み込んで4列目で検索した い場合のやり方を示します。 16 解答例 1. 2. 「メモ帳」など任意のエディタで リストファイル(list.txt)を作成 目的のキーワードリストを含むファイルを作成し（例：list.txt） 該当箇所を変更し、Rコンソール画面上でコピペ Apr 14 2015 17 解答例 1. 2. 一連の作業手順を記述した スクリプトを1つのファイルと して保存することをお勧め 目的のキーワードリストを含むファイルを作成し（例：list.txt） 該当箇所を変更し、Rコンソール画面上でコピペ Apr 14 2015 18 ありがちなミス1 Apr 14 2015 作業ディレクトリの変更を忘 れているため、in_f1で指定し た最初のファイルの読み込み 段階でエラーが出る。つまり、 実際に行ったフォルダ中には annotation.txtというファイル は存在しないということ。 19 ありがちなミス2 Apr 14 2015 必要な入力ファイルが作業 ディレクトリ中に存在しない。 この場合、in_f2で指定した genelist1.txtが存在しないた め、それの読み込み段階でエ ラーが出ている。それゆえ、 その情報を用いているコマン ド部分でエラーが出ている。 20 ありがちなミス3 Apr 14 2015 出力予定のファイル名と同じものを エクセルなど別のプログラムで開い ているため、最後のwrite.table関数 のところでエラーが出る。対処法は、 出力ファイル名を変更するか、開い ている別のプログラムを閉じる。 21 ありがちなミス4 Apr 14 2015 実行スクリプトをコピーする際、最後 の行のところで改行を含ませずにR Console画面上でペーストしたため、 最後のコマンドが実行されない（出 力ファイルが生成されない）。これも 比較的ありがちなパターンです。コピ ペ後に無意識にリターンキーを押す ことを心がけるだけでもよいでしょう。 22 警告メッセージ Apr 14 2015 list.txtファイル作成時に、membraneと打った後に改 行を入れた場合（左）と入れない場合（右）の挙動の 違いを把握し、後学のために警告メッセージの意味 を理解しておくとよい。この場合は結果には影響して いないことがわかる。Rは警告メッセージの記述内 容が比較的分かりやすいのでよく読むべし。 23 Contents 行列形式ファイルの解析基礎（アノテーションファイルを例に） 例題をテンプレートとして任意の解析を行う基本手順 ありがちなミスとエラーメッセージ 入力ファイルの最後の改行の有無 コード内部の説明（行列演算の基礎） multi-FASTAファイルからの各種情報抽出 Apr 14 2015 基本情報取得（コンティグ数、配列長、N50、GC含量） 任意の領域の切り出し GC含量計算部分の説明 24 コード内部の説明 Apr 14 2015 コードの中身を説明します。 黒枠部分を再度コピペ。 25 ①in_f1で指定したファイルを読み込め ②読み込むファイルの最初の行はヘッダー部分 ③ファイルの区切り文字はタブです ④読み込んだ結果をdataという名前で取り扱う 読み込み ④ Apr 14 2015 ① ② ③ 26 行列data ② Apr 14 2015 dataと打ってリターン。入力 ファイルの中身を正しく読 み込めていることがわかる。 header=TRUEとしているの でこのように見える。 27 行数と列数を表示 Apr 14 2015 オブジェクトdataの行数と列数は11と4。 ウェブページ中の表記が灰色なのは、 特にやらなくてもいいコマンドだから。 28 行列の要素へのアクセス Apr 14 2015 行列dataの要素へのアクセスは[行, 列]。 読み込んだ元ファイルである annotation.txt中では7行×4列目だが、 1行目をヘッダー行としているので6行 ×4列目とする必要がある。利用例は、 ファイル読み込み時に「x行×y列目に 不具合がある」のようなエラーが出た時 のトラブルシューティングなど。 29 行列の要素へのアクセス Apr 14 2015 上矢印キーを押すと、直前に打っ たコマンドが表示される。最初か ら全部打ち直すのではなく、上下 左右の矢印キーを有効に利用し 最小限の労力で打つべし! 30 行列の要素へのアクセス Apr 14 2015 行列dataの要素へのアクセスは [行, 列]。2行目の情報のみ抽出。 読み込み時にhead=TRUEとして いたので、ヘッダー行がついて いることが分かる。 31 行列の要素へのアクセス Apr 14 2015 行列dataの要素へのアクセスは [行, 列]。2列目の情報のみ抽出。 32 行列の要素へのアクセス Apr 14 2015 行列dataの要素へのアクセスは [行, 列]。param列目の情報の み抽出。paramには1という数値 が代入されている。 33 参考 Tips：関数とオプション Apr 14 2015 行列dataの最初の数行を表示 したい場合は、head関数を利用。 n=3というオプションを利用する と最初の3行分のみ表示。関数 ごとに様々なオプションを利用 可能です。 34 参考 Tips：タブ補完 Apr 14 2015 列番号を指定する以外にも 特定の列を表示するやり方が 存在します。head=TRUEで入力 ファイルを読み込むと、列の名 前を利用することができます。 subcellular_location列の情報を 抽出したい場合は、「data$su」く らいまで打ち込んでからTab キーを押す。 35 参考 Tips：タブ補完 Apr 14 2015 列名中のsuからはじまる文字列 を補完して表示してくれる。「Tab キーを用いた補完機能」という 意味で「タブ補完」という。このテ クニックはLinuxでも利用可能。 36 参考 Tips：table関数 Apr 14 2015 table関数は、ベクトル中の要素 ごとの出現回数を返す。「NGS データ中の特定のリードの出現 回数」や、「アノテーションファイ ル中の染色体ごとの遺伝子数」 など、様々な局面で利用可能。 37 参考 Tips：ソート Apr 14 2015 sort関数と併用すること で全体像を俯瞰可能 38 Tips：is.element関数 Apr 14 2015 hogeベクトルに対して、”nuclear”の 文字が存在する場所をTRUE、存在 しない場所をFALSEとして返す。 as.character関数は、文字列ベクト ルとして取り扱いたい場合に利用。 39 Tips：“二重クォーテーション” Apr 14 2015 二重クォーテーションが自動で 変更されるエディタは非推奨で す。日本語の二重クォーテー ションもだめです。Microsoft WordやPDFファイル中のコード のコピペ時によくハマります。 40 目的をおさらい Apr 14 2015 目的は、数万～数百万行からなるファ イルを読み込んで特定のキーワードを 含む行のみ取り出すテクニックを習得。 41 目的をおさらい Apr 14 2015 論理値ベクトルobjを用いてTRUEの 要素に対応する行を抽出している 42 目的をおさらい Apr 14 2015 コード作成当時はas.character関数を用 いてデータの型を文字列ベクトルに揃え ていたが、少なくとも現在(R ver. 3.1.3)は この関数がなくても大丈夫なようだ。 43 Tips：論理値ベクトル ファイル名：rcode_20140908.txt Apr 14 2015 参考 is.element関数は、 1番目の引数で 指定したベクトル(この場合x)の要素 が2番目の引数で指定したベクトル (この場合y)の要素として含まれるか 否かをTRUEまたはFALSEで返す。 44 参考 Tips：論理値ベクトル 全遺伝子から特定の遺伝 子群の位置情報を取得した い場合などによく用います ファイル名：rcode_20140908.txt Apr 14 2015 45 1と12は手順が異なるだ けで実質的に同じです genelist1.txt Apr 14 2015 46 このコードはヘッダー行 がある場合のものです 入力: annotation.txt 出力: hoge12.txt Apr 14 2015 47 このコードはヘッダー行 がない場合のものです 入力: annotation2.txt 出力: hoge13.txt Apr 14 2015 48 ヘッダー行がある場合 ヘッダー行がない場合 Apr 14 2015 49 Contents 行列形式ファイルの解析基礎（アノテーションファイルを例に） 例題をテンプレートとして任意の解析を行う基本手順 ありがちなミスとエラーメッセージ 入力ファイルの最後の改行の有無 コード内部の説明（行列演算の基礎） multi-FASTAファイルからの各種情報抽出 Apr 14 2015 基本情報取得（コンティグ数、配列長、N50、GC含量） 任意の領域の切り出し GC含量計算部分の説明 50 FASTA形式 Apr 14 2015 Rでmulti-FASTAファイルを読 み込んで自在に解析できます。 ゲノム配列解析≒FASTA形式 ファイルの解析。ここでは全体 像を完全に把握すべくhoge4.fa ファイルを仮想ゲノム配列ファ イルとして取り扱う。 51 ゲノム配列 Apr 14 2015 実際のゲノム配列はここからも 取得可能。Rで染色体ごとの配 列長やGC含量の計算ができる。 52 基本情報取得 multi-FASTAファイルを読み込んで、 トータルの配列長、染色体数（コンティ グ数）、配列長の平均、中央値、最大 値、最小値、N50、GC含量を計算した 結果を返すコードを実行してみよう。 入力: hoge4.fa 出力: hoge1.txt Apr 14 2015 53 基本情報取得 Apr 14 2015 コードの最初のほうに入力ファ イルと出力ファイルを記述する ので、コピペで実行した結果とし てどういう名前のファイルが出 力されるべきかわかる。 54 ①例題の入力ファイル （hoge4.fa)をダウンロード。②R 上で作業ディレクトリの確認。③ ① 作業ディレクトリに解析したい入 力ファイルがあることを確認。 基本情報取得 ② ③ Apr 14 2015 55 基本情報取得 ① ①一連のコマンド群をコピーして ②R Console画面上でペースト。 ブラウザがInternet Explorerの場 合は、CTRLとALTキーを押しな がらコードの枠内で左クリックす ると、全選択できます。 ② Apr 14 2015 56 実行結果 Apr 14 2015 コピペ後にlist.files()で、出力ファイル として指定したファイル名のhoge1.txt が作成されていることを確認。 57 実行結果 Apr 14 2015 出力ファイルをテキストエディタやExcel で眺めてもよいが、オブジェクトtmpの 中身を出力しているだけなのでR上で 眺めている。 58 実行結果 contig_1の配列が最短、contig_2 の配列が最長であることがわか る。入力と出力の関係を確認。 入力: hoge4.fa 出力: hoge1.txt Apr 14 2015 59 averageだと外れ値の影響を受けや すく、medianだと短いコンティグが多 くを占める場合に不都合らしい。 N50 アセンブル結果の評価基準の一つ 出力: hoge1.txt 長いコンティグから足していってTotal_length の50%に達したときのコンティグの長さ 一般に数値が大きいほどよい contig_2 (103 bp) contig_3 (65 bp) contig_4 (49 bp) contig_1 (24 bp) Total_length ×0.5 (120.5 bp) Total_length (241 bp) Apr 14 2015 60 課題1 アセンブル結果の評価基準の一つ Apr 14 2015 左記のコードをhoge4.faの代わりに sample4.fastaを入力として実行し、 ①全配列長（配列長の総和）、② N50の値、および③GC含量を示せ。 長いコンティグから足していってTotal_length の50%に達したときのコンティグの長さ 一般に数値が大きいほどよい 61 ①作業ディレクトリの確認、 ②解析したいファイルの存在 確認、③推奨エディタの起動。 課題1 ③ ① ② Apr 14 2015 62 課題1 Apr 14 2015 テンプレートコードをコピペし、 入力ファイル名部分を変更。 63 実行結果 課題1の実行結果 Apr 14 2015 64 コード内部の説明 Apr 14 2015 コードの中身を説明します。 黒枠部分を再度コピペ。 65 コード内部の説明 Apr 14 2015 入力ファイル情報を格納し たものがfastaオブジェクト。 widthの位置にあるのがコ ンティグごとの配列長情報。 66 コード内部の説明 Apr 14 2015 width(fasta)にsum関数を適用すれば、 トータルの配列長（配列長の総和）になる。 67 コード内部の説明 Apr 14 2015 length関数は要素数を返す。こ の場合、fastaオブジェクトの要素 数（つまりコンティグ数）を返す。 68 Tips：条件判定 Apr 14 2015 50 bp以上のコンティグからなる サブセットの抽出ができそうだ! 69 Tips：条件判定 Apr 14 2015 コードの中身が分かると応用範囲 が飛躍的に増大。一定以上のスキ ルをもつバイオインフォマティシャ ンは「例題を探す」よりも「自分で作 る」ヒトのほうが多いかも…。 70 Contents 行列形式ファイルの解析基礎（アノテーションファイルを例に） 例題をテンプレートとして任意の解析を行う基本手順 ありがちなミスとエラーメッセージ 入力ファイルの最後の改行の有無 コード内部の説明（行列演算の基礎） multi-FASTAファイルからの各種情報抽出 Apr 14 2015 基本情報取得（コンティグ数、配列長、N50、GC含量） 任意の領域の切り出し GC含量計算部分の説明 71 任意の領域の切り出し subseq関数を用いて、任意の領域 の配列を切り出すことができます。 入力：sample1.fasta >kadota AGTGACGGTCTT 出力：hoge1.fasta >kadota TGACGGT Apr 14 2015 72 subseq関数実行時に、①数 値を直接指定してもいいし、② オプション名を明記してもよい。 Tips：関数のオプション ① 入力：sample1.fasta >kadota AGTGACGGTCTT Apr 14 2015 出力：hoge1.fasta >kadota TGACGGT ② 73 Tips：関数のオプション ①原因既知状態でエラーを出す。 ②「3番目の位置から5塩基分抽 出」という他のオプションを利用。 ① ② 入力：sample1.fasta >kadota AGTGACGGTCTT Apr 14 2015 出力：hoge1.fasta >kadota TGACGGT 74 Tips：関数の使用法 「?関数名」で使用法を記した ウェブページが開く。ページ の下のほうに、大抵の場合使 用例が掲載されている。使用 法既知の関数のマニュアルを いくつか読んで慣れておく。 … Apr 14 2015 75 任意の領域の切り出し Apr 14 2015 入力がmulti-FASTAファイル （hoge4.fa）で、リストファイル （list_sub2.txt）で指定した複 数領域を切り出したい場合 76 Contents 行列形式ファイルの解析基礎（アノテーションファイルを例に） 例題をテンプレートとして任意の解析を行う基本手順 ありがちなミスとエラーメッセージ 入力ファイルの最後の改行の有無 コード内部の説明（行列演算の基礎） multi-FASTAファイルからの各種情報抽出 Apr 14 2015 基本情報取得（コンティグ数、配列長、N50、GC含量） 任意の領域の切り出し GC含量計算部分の説明 77 GC含量計算部分の説明 右のサイドバーを下に移動させ るとGC含量計算部分を見られる。 ① ② Apr 14 2015 78 GC含量計算部分の説明 Apr 14 2015 fastaオブジェクトを出 発点として、配列全体 のGC含量（57.68%）を 得るところの説明です。 79 GC含量計算部分の説明 Apr 14 2015 黒枠部分を再度コピペした のち、fastaオブジェクトの 中身を表示させたところ。 80 GC含量計算部分の説明 Apr 14 2015 alphabetFrequency関数は、 塩基ごとの出現回数を返す。 81 GC含量計算部分の説明 DNA配列上のMは「A or C」、 Rは「A or G」などという ルールがあるようです。 http://en.wikipedia.org/wiki/FASTA_format Apr 14 2015 82 GC含量計算部分の説明 Apr 14 2015 dim関数は行列の行数と列数 を返す。alphabetFrequency 関数出力結果は、4行×18列 からなることが分かる。キー ボードの上下キーを上手に利 用して最小限の労力でキータ イプ（あるいはコピペ）すべし! 83 GC含量計算部分の説明 Apr 14 2015 任意のサブセットを取得可能。 2:3やc(1,4)などをうまく利用。 84 GC含量計算部分の説明 Apr 14 2015 黒丸中の数値はcontig_1中の Aの数が4個、赤丸中の数値は、 contig_4中のTの数が10個で あるということ。 rowSums関数 は行ごとの和を返す。 85 GC含量計算部分の説明 Apr 14 2015 rowSums関数の入力として、 ACGTのみのカウント数を 与えているが、その結果（返 り値）は、配列中にNなどを 含まない場合は実質的にコ ンティグごとの配列長と同じ。 86 GC含量計算部分の説明 Apr 14 2015 オブジェクトCG中には、配列ごとの CとGのカウント数が格納されてい る。オブジェクトACGT中には、配列 ごとのA, C, G, Tのカウント数が格 納されている。例えば49塩基から なるcontig_4中に、ACGTの4種類 の塩基が49個、CGの数は25個あ ることを意味する。sum関数は、ベ クトルの要素の和を返す。 87 GC含量計算部分の説明 Apr 14 2015 ここではsum関数を用いて配 列全体の総和でGC含量計算 をしているが、sum関数を用い ずに「CG/ACGT」とやると、コ ンティグごとのGC含量を得ら れる。例えば、contig_1はCG の数が16個で、ACGTの数が 24個。それゆえGC含量は 16/24 = 0.6666667となる。 88 配列ごとのGC含量計算 Apr 14 2015 sum関数を用いずに「CG/ACGT」 とやって、コンティグごとのGC含量 を得るための項目。記述内容がほ ぼ同じであることが分かる。 89 配列ごとのGC含量計算 Apr 14 2015 出力ファイル(hoge1.txt) 中の一番右側の列が配 列ごとのGC含量です。 90 配列ごとのGC含量計算 Apr 14 2015 ACGT列は4種類の塩基のみの 出現数、Length列は配列長情 報を表す。配列長は、ACGT以 外の全てを含むので、2つの数 値の差分(Length - ACGT)がN などのACGT以外の塩基のトー タルの出現回数ということにな る。Length ≧ ACGTという関係。 91 課題2 sample4.fastaの配列 ごとのGC含量を示せ。 アセンブル Apr 14 2015 92 ①作業ディレクトリの確認、 ②解析したいファイルの存在 確認、③推奨エディタの起動。 課題2 ③ ① ② Apr 14 2015 93 課題2 Apr 14 2015 テンプレートコードをコピペし、 入力ファイル名部分を変更。 94 課題2の実行結果 ①出力ファイル(hoge1.txt)中 の一番右側の列が配列ごと のGC含量です。②gene_4の 配列中にACGT以外の文字 が1つあることがわかる。 ② Apr 14 2015 ① 95
© Copyright 2024 ExpyDoc
