お家でできるMacBookでやる次世代シーケンスデータ解析 pdf

全文

(1)緒方法親. お家でできるMac Bookでやる次世代シーケンスデータ解析.

(2) まえがき. ご連絡は @n0rr もしくは [email protected] までお願いします。 2012/03/25 2012/03/30 追記.

(3) 次世代シーケンサーが市場に出てから４年経った。日本を. 在することは事実であるが、全てではない。実際に私もカイ. 含む世界の多くの研究機関に導入されて活躍している。その. コのトランスクリプトームを実験条件ごとに比較して発現変. 一方で、次世代シーケンサーの出力する膨大なテキストデー. 動遺伝子を抽出し、それらの遺伝子のアノテーションを付け. タの解析を行う人材の不足が指摘されてきた。私自身、修士. る作業をMac Book Pro１台で行ってきた。Mac にはもともと. 課程までは無料統計ソフトパッケージのRを使うことはある. ターミナルという端末が付いているし、それを通してUNIXを. ものの、”ウェットの”研究者として活動しており、博士後期. 使ったり簡単なプログラムを作成することができる。プログ. 課程に進んで次世代シーケンサーの解析を行う必要が出てき. ラムと言っても恐れることはない。わずかに１０個のコマン. た時には困った。幸いにも当時つくばの研究所にいらっ. ドを覚えさえすれば次世代シーケンスデータ解析には足りて. しゃった古崎さんにUNIXの使い方の手ほどきをしてもら. しまうし、分からないことがあってもグーグル先生が教えて. い、さらに当時農工大の学部生だった西村くんにUNIX/. くれる。どうしようもなかったらNGS現場の会や. Linuxについて習うことが出来たが、もし彼らのような人が周. SEQanswersなどのコミュニティーで質問すればいいのだ。. りにいなかったら絶望的であった。. これは、次世代シーケンサーを使いたい研究者、配属され. パソコンに向かう時間の長いデータ解析を主とする研究者. た研究室にたまたま塩漬けになってる次世代シーケンスデー. とを”ドライな”研究者と呼ぶことがあるが、単に次世代シー. タがあった学部生だけでなく、生物が好きで研究したいけど. ケンサーに対応することだけでなく今後さらに複雑で多量の. 事情があって研究機関に所属できない人も対象に、Mac. データを出力する実験機器が登場するであろうことを考える. Book１台を使って追加予算なしに自由な研究が出来るレベル. と、データ解析を”ドライな”研究者に任せてしまおうと思う. まで達して欲しいと思って書いている。今や最も貴重な資源. ことは解釈の視点を放棄してしまうことになるのではない. は試薬でも試料でもデータでもなく、データを解析できるわ. か。そこで、もっと多くの”ウェットな”研究者が次世代シー. ずかな技術とそれを解釈できる柔軟なアイデアなのだ。みん. ケンサーのデータ解析を行うようになる未来を期待して、こ. なが参加してくれれば日本の生物学はもの凄い早さで進むと. れを書くことにした。大型のサーバーが必要となる解析が存. 思う。 ii.

(4) チャプタ 1. ターミナルで UNIXを使ってやるところ. bowtieの使い方発現量推定多型解析.

(5) セクション 1. 多くの場合で何かを学ぶときまずは基本からとなるだろう. Bowtieでマッピングする. し、私もターミナルやUNIXを使う前にLinux標準教科書の読解を勧められたが、プログラミング自体に興味のない生物学者にとって基本の習得は苦痛であるし何の意味があるのか分からない。それよりも、とにかく何かを動かして生物学的な意味のある結果を得られた方が楽しいしわかると思う。このセクションではもっとも広く用いられているマッピングツー. できるようになること 1. 無料のマッピングツールであるBowtieを使う 2. Bowtieの出力した結果の意味が分かる. ルであるBowtieを使って、大腸菌のシーケンスデータを大腸菌のゲノムデータにマッピングする演習を行ってもらう。さっそく手元のMacのターミナルを起動して欲しい。私の Mac Book Proは２０１１年の３月に札幌のApple Storeで１３. 3. ”ディレクトリ”と”ファイル”の概念が分かる. 万円くらいで購入したものを途中でLionにアップデートしたもので、Mac OS X バージョン10.7.2 プロセッサ 2.3GHz Intel. ※ 文中のnorは私のユーザー名なので、各々自分の macのユーザー名に合わせて読んで欲しい。. Core i5 メモリ8G 1333MHz DDR3 機動ディスク Macintosh HD となっている。画面は13-inch。HDDは320GBでこれは２つ程度のプロジェクトを管理するのに足りる容量である。たくさんのプロジェクトがある場合には外付けHDDで対応している。ターミナルの起動だが、アプリケーションの中のユーティリティーの中に黒い四角いアイコンのターミナルがあるのでクリックすればよい。私の環境ではユーティリティの左から4番目の一番下のところに表示されている。よく使うので 4.

(6) Dockに入れてしまうのがよいだろう。. レクトリに今いることを示している。ここで、ディレクトリの中をリストアップするlsコマンドを使おう。すると、 Last login: Sat Jan 21 18:07:52 on ttys000 dhcp186-162:~ nor$ ls Desktop Documents Downloads. クリックすると、ターミナルが起動して以下のウィンドウが. Library Movies. 出てくる。. Music Pictures Public R analtools dhcp186-162:~ nor$. と出てくる。これをFinderで表現すると、これは格好が悪いので、[ターミナル]→[環境設定]→[起動] の[起動時に開く]で[Pro]を選択しておくと黒くて透けてるターミナルを使うことが出来る。起動すると、 Last login: Sat Jan 21 18:07:52 on ttys000 dhcp186-162:~ nor$. こういうことになる。このlsコマンドにオプションを付ける. と書いてあり、今起動したことと、ユーザーとしてnorがログ. とさらに多くの情報を得ることが出来る。たとえば-a オプ. インしていることが分かる。norは私のMac Book Proのユー. ションを使って隠しファイルも表示させると、. ザー名である。これは、/Users ディレクトリの中のnorディ. 5.

(7) dhcp186-162:~ nor$ ls -a. _64.zip/download いきなりこのリンクを入力するとすぐにダ. .. ウンロードが始まってしまうので、sourceforgeのリンクから. ... http://sourceforge.net/projects/bowtie-bio/files/bowtie/0.12. 7/ 自分の欲しいファイルをクリックするのが良い。今回は. .CFUserTextEncoding .DS_Store .RData. bowtie-0.12.7-macos-10.5-x86_64.zip を手に入れた。ダウン. .Rapp.history. という具合に隠しファイルも出てくし、-lオプションでは. ロードされたファイルは通常ダウンロードディレクトリに入っているので、これをドラッグドロップして先ほど作った. dhcp186-162:~ nor$ ls -l. analtoolsディレクトリに移動させて欲しい。ちゃんと出来て. total 2944 drwx------+ 16 nor staff. 544 1 21 18:16 Desktop. drwx------+ 204 nor staff 6936 1 20 20:53 Documents. いれば、ターミナルに戻って、cd コマンドを使ってanaltools ディレクトリに移動し、lsコマンドで確認できる。. drwx------+ 96 nor staff 3264 1 20 01:14 Downloads. とそれぞれの大きさや権限が表示される。. dhcp186-162:~ nor$ cd analtools dhcp186-162:analtools nor$ ls. では、Bowtieをしまうディレクトリを作ろう。私の環境では. とした時にbowtie-0.12.7-macos-10.5-x86_64.zipが出てくる. 解析ツールをanaltoolsというディレクトリにしまっているの. はずだ。次にbowtie-0.12.7-macos-10.5-x86_64.zipの解凍を. で、. 行う。解凍には、unzipコマンドを使う。. dhcp186-162:~ nor$ mkdir analtools. dhcp186-162:analtools nor$ unzip bowtie-0.12.7-macos-10.5-x86_64.zip. とした。これで、/Users/nor/ディレクトリの中に/Users/. これで bowtie-0.12.7-macos-10.5-x86_64.zipが解凍されて/. nor/analtoolsディレクトリが出来たことになる。. Users/nor/analtoolsディレクトリ内にbowtie-0.12.7ディレクトリが出来ているはずなのでlsコマンドで確認してほしい。. 次はBowtieを手に入れよう。Bowtieは以下のリンクから手に入れることが出来る。http://sourceforge.net/projects/. また、cdコマンドを使ってbowtie-0.12.7ディレクトリへ移動する。. bowtie-bio/files/bowtie/0.12.7/bowtie-0.12.7-macos-10.5-x86 6.

(8) dhcp186-162:Downloads nor$ cd bowtie-0.12.7 dhcp186-162:bowtie-0.12.7 nor$. これでBowtieを使う準備が整った。ふつうはmake作業があるのだが、Bowtieは解凍すればすぐに使えるのだ。試しに. e_coli_1000.fqがそれだ。more コマンドを使ってのぞいてみることにする。 dhcp186-162:reads nor$ more e_coli_1000.fq @r0. dhcp186-162:bowtie-0.12.7 nor$ /Users/nor/analtools/bowtie-0.12.7/bowtie. GAACGATACCCACCCAACTATCGCCATTCCAGCAT. とすれば、. +. No index, query, or output file specified!. EDCCCBAAAA@@@@?>===<;;9:99987776554. Usage:. @r1. bowtie [options]* <ebwt> {-1 <m1> -2 <m2> | --12 <r> | <s>} [<hit>]. の後に数十行に渡ってBowtieのマニュアルを読むことが出来. CCGAACTGGATGTCTCATGGGATAAAAATCATCCG + EDCCCBAAAA@@@@?>===<;;9:99987776554. る。ターミナルでの入力中にtabキーを押すとファイル名や. @r2. ディレクトリ名を補完できて楽である。次にシーケンスファ. TCAAAATTGTTATAGTATAACACTGTTGCTTTATG. イルの準備をする。といってもデモデータが準備されてお. + EDCCCBAAAA@@@@?>===<;;9:99987776554. り、/Users/nor/analtools/bowtie-0.12.7/reads ディレクトリ. more コマンドで出力されたはじめの12行はこうなっている。. に入っているのでそれを確認するだけである。. この12行は3つのシーケンス（ショートリードともいう）を表す。4行で１単位になっているのがfastqファイルの特徴であ. dhcp186-162:bowtie-0.12.7 nor$ cd reads/. る。. dhcp186-162:reads nor$ ls e_coli_1000.fa . e_coli_10000snp.fa. e_coli_1000_1.fq. e_coli_1000.fq . e_coli_10000snp.fq. e_coli_1000_2.fa. @r0. e_coli_1000.raw . e_coli_1000_1.fa. e_coli_1000_2.fq. GAACGATACCCACCCAACTATCGCCATTCCAGCAT. いろいろ入ってる。今回使うのはイルミナ社の次世代シーケンサーGenome Analyzer Ⅱxの出力形式であるfastqのファイルなので、拡張子が.fqになっているファイルをみる。. + EDCCCBAAAA@@@@?>===<;;9:99987776554. これが１単位で、１行目がシーケンスのID、２行目がシーケンス、３行目は+、４行目はシーケンスのクオリティになっている。３行目の+以降にシーケンスのIDが入っていること 7.

(9) もあるが、ファイルサイズを小さくするために+だけになって. 渡している。wc コマンドは文字数を数えるコマンドで、-l オ. いることも多い。@ID、シーケンス、+ID、シーケンスクオ. プションを使って行を数えるようにしている。この”|”は大事. リティーの４行になっていることを知っていることは非常に. なので覚えておくこと。４行で１シーケンスなので、この. 重要なので覚えておくこと。more コマンドを終了させたい時. e_coli_1000.fqには1000本のシーケンスが含まれていること. は、[control]を押しながら[c]を押す。テキストエディタなど. が分かった。次にリファレンスゲノムデータを見よう。ゲノムデータは /Users/nor/analtools/bowtie-0.12.7/genomes ディレクトリに入っているので確認する。. ターミナル内でのコントロール. dhcp186-162:reads nor$ cd .. dhcp186-162:bowtie-0.12.7 nor$ cd genomes/ dhcp186-162:genomes nor$ ls ふだんのコントロール. NC_008253.fna. cd コマンドのcd .. では、ひとつ上のディレクトリへの移動を他のアプリケーションでコピーなどのショートカットキーでつかうコントロールとは異なるので注意すること。. 示す。ここでは /Users/nor/analtools/bowtie-0.12.7/genome ディレクトリからcd ..を使って一度 /Users/nor/analtools/bowtie-0.12.7/ ディレクトリに移動. このfastqファイルに何本のシーケンスが入っているのか数え. し、cd genomes/ で. てみよう。. /Users/nor/analtools/bowtie-0.12.7/genomes ディレクトリに. dhcp186-162:reads nor$ more e_coli_1000.fq | wc -l 4000. 移動しているが、これをcd ../genomes/と省略しても構わない。その後lsコマンドで内容を確認している。ゲノムデータと. ここではパイプ”|”とwc コマンドを使った。”|”はUNIX上で処. してNC_008253.fnaが見つかった。拡張子 .fna はfas-. 理を続けて行うことを示すものであり、ここではmore コマン. taフォーマットを示している。またmore を使って覗いたり、. ドの出力である”e_coli_1000.fq”の内容をwcコマンドに引き. wc を使って数えたりしてみよう。 8.

(10) dhcp186-162:genomes nor$ more NC_008253.fna >gi|110640213|ref|NC_008253.1| Escherichia coli 536, complete genome AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGA GTGTCTGATAGCAGCTTCTGAACTGGTTACCTGCCGTGAGTAAATTAAAATTTTA TTGACTTAGGTCACTAAATACTTTAACCAATATAGGCATAGCGCACAGACAGATA AAAATTACAGAGTACACAACATCCATGAAACGCATTAGCACCACC dhcp186-162:genomes nor$ more NC_008253.fna | wc 70557 70562 5009545. 70kの長さがあるようだ。このリファレンスファイルに先ほど. installしないで使うことをお勧めする。もし何度も使って面倒なようなら、その時だけPathを書かなくても良いように、 dhcp186-162:genomes nor$ PATH=$PATH:/Users/nor/analtools/bowtie-0.12.7/. と入力してしまっても良い。さて、bowtiebuildが完了し、新しいファイルが出来ていれば成功だ。ls コマンドで確認すると、. のシーケンスをマッピングするために最後の準備が必要である。ここから本格的にbowtieを使って解析を進めることになる。以下のように行う。. dhcp186-162:genomes nor$ ls NC_008253.1.ebwt NC_008253.4.ebwt NC_008253.rev.2.ebwt NC_008253.2.ebwt NC_008253.fna NC_008253.3.ebwt NC_008253.rev.1.ebwt. dhcp186-162:genomes nor$ /Users/nor/analtools/bowtie-0.12.7/bowtie-build NC_008253.fna NC_008253. いつくかのファイルが生成されている。これで準備が整った. bowtie-buildという解析プログラムを使うためにbowtie-build. ので、さっそく先のシーケンスをマッピングしたい。コマン. までのPathを記述し入力ファイルとしてNC_008253.fnaを指. ドは以下の通りである。. 定、さらに出力ファイル名としてNC_008253を指定している。Pathという概念ははじめてでてきたが、これまで使って. dhcp186-162:genomes nor$ /Users/nor/analtools/bowtie-0.12.7/bowtie -l 32 -n 2 --best -p 1 /Users/nor/analtools/bowtie-0.12.7/genomes/NC_008253 /Users/nor/analtools/bowtie-0.12.7/reads/e_coli_1000.fq | more. きたls やcd、wc などは/usr/bin/というディレクトリに入っているためにPathなしでも動くが、bowtie-buildはそこに入っていないためにPathが必要となる。/Users/nor/analtools/ bowtie-0.12.7/までがPathである。いちいち解析ツールを/ use/binに置いても構わないが（そこに置くことをinstallという）、次世代シーケンスデータ解析でつかう解析ツールたちは更新が盛んなのでバージョン管理を楽に行うためにあえて. 実際には改行を含まない。文字と文字の間には半角スペースが入っている。ひとつずつ説明すると、最初の/Users/nor/ analtools/bowtie-0.12.7/bowtie でbowtieまでのPathとbowtieを動かす命令をしている。次の -l 32 はbowtieのオプションでシーケンスの最初から32塩基までをマッピングに使うことを意味する。続く-n 2 ではミスマッチの許容数を示し、これはシーケンスエラーやSNPを考慮するためであ 9.

(11) る。--bestはマルチヒットオプションで、リファレンス配列上. 回の場合には/Users/nor/analtools/bowtie-0.12.7/genomes. の複数の部位にマッピングできた場合、最もきれいに納まっ. ディレクトリに保存されるが、他の場所、例えば/Users/nor/. たものを選ぶように指示している。最後のオプション-p 1 は. Desktop に保存したい場合には、. 使うプロセッサ数を指定するもので、私の環境では最大4と. ファレンスファイルを指定し、シーケンスファイルまでの. dhcp186-162:genomes nor$ /Users/nor/analtools/bowtie-0.12.7/bowtie -l 32 -n 2 --best -p 1 /Users/nor/analtools/bowtie-0.12.7/genomes/NC_008253 /Users/nor/analtools/bowtie-0.12.7/reads/e_coli_1000.fq > /Users/nor/Desktop/bowtie_mapped_e_coli_1000_2_NC_008253.txt. Pathとシーケンスファイルを指定することでbowtieは動く。. とすれば良い。書き出して保存すれば、ターミナルの方では. なっている。その後にリファレンスファイルまでのPathとリ. 結果をmore でターミナル上に表示するようにしている。. # reads processed: 1000 # reads with at least one reported alignment: 697 (69.70%) # reads that failed to align: 303 (30.30%) Reported 697 alignments to 1 output stream(s). 上のような出力は得られただろうか。左のカラムから説明すると、シーケンスID、ワトソン鎖かクリック鎖か、張り付いた先のリファレンス配列の名前、張り付いた場所、シーケンス、シーケンスクオリティ、マルチヒット、ミスマッチの場所、となっている。これを保存する場合には、 dhcp186-162:genomes nor$ /Users/nor/analtools/bowtie-0.12.7/bowtie -l 32 -n 2 --best -p 1 /Users/nor/analtools/bowtie-0.12.7/genomes/NC_008253 /Users/nor/analtools/bowtie-0.12.7/reads/e_coli_1000.fq > bowtie_mapped_e_coli_1000_2_NC_008253.txt. という具合にマッピングの統計が表示される。今回は、1000本のシーケンスを解析し、697本がリファレンス配列にヒット、303本がどこにもヒットしなかった。ここまでがbowiteの基本的な使い方である。いつくかのコマンドも出てきた。今回の作業の要約は以下の通りである。. 1. ターミナルの起動. とすればよい。”>” は書き出しを命令するもので、マッピングの結果をbowtie_mapped_e_coli_1000_2_NC_008253.txt. 2. analtools ディレクトリの作成(ls, mkdir). という名前のテキストファイルにして保存するようになる。. 3. Bowtie の入手とanaltools ディレクトリへの保存. 特に保存場所を指定しない場合は現在いるディレクトリ、今 10.

(12) 4. analtools ディレクリへの移動(cd) 5. Bowtie圧縮ファイルの解凍(unzip) 6. シーケンスデータの確認(ls, cd, more, wc) 7. リファレンスデータの確認(ls, cd, more, wc) 8. リファレンスデータのインデックス化(bowtie-build) 9. マッピング(bowtie). 11.

(13) セクション 2. 前回のセクションではBowtieのパッケージに元から入って. マッピングして発現量を推定する. いるシーケンスファイルとリファレンスファイルをBowtieで解析した。それらのファイルはよく準備されていたために、そのまま利用することができたが、通常自分の解析したいデータを解析する場合には、ファイルの準備が必要でこれに多くの時間が割かれることも多い。地味な作業であるし、パソコンが得意な人に取っては常識以前の問題と捉えられてい. できるようになること 1. Bowtieの結果を解析して発現量を調べる. るためにこの部分に紙面を割いたテキストを見たことがない。しかし先に述べたように非常に重要なので少々冗長にはなるが我慢していただきたい。. 2. DDBJシーケンスアーカイブから欲しいシーケンスデータを手に入れる. さっそくシーケンスデータを手に入れることにする。シーケンスデータはDDBJ Sequence Read Archive. 3. リファレンスファイルを手に入れる. http://trace.ddbj.nig.ac.jp/dra/index.shtml からダウンロードすることが出来る。[データの検索・ダウンロード]を選択してDRASearch http://trace.ddbj.nig.ac.jp/DRASearch/ を見る。今回は、キイロショウジョウバエKc167セルラインの mRNAシーケンスデータ SRR029023 を例として利用することにする。DRASearchのAccession の項目にSRR029023と入力し、Searchをクリックする。. 12.

(14) SRR029023がキイロショウジョウバエのKc167セルライン由来であること、Single リードであること（Pair end でない)、 Trizolを使ってサンプル調製したこと、Illumina のGAを使って36bpシーケンスしたことが分かる。ダウンロードした SRR029023.fastq.bz2は通常ダウンロードディレクトリ返ってくる結果がこちら。. /Users/nor/Download/ 内に入っているのでFinder で確認しよう。このシーケンスデータ解析を管理するために新しくディレクトリを作ることにする。ターミナルを起動して、 dhcp186-162:~ nor$ mkdir drkc167 dhcp186-162:~ nor$ cd drkc167. /Users/nor/ ディレクトリ内に /Users/nor/drkc167/ ディレクトリを作成し、drkc167 ディレクトリに移動する。Finderで FASTQをクリックするとSRR029023.fastq.bz2のFTPサイト. SRR020923.fastq.bz2を/Users/nor/drkc167 ディレクトリに. に繋がるので、SRR029023.fastq.bz2をクリックしてダウン. 入れても構わないが、ターミナルを使うと. ロードする。先のDRASearchの結果のページのStudyをクリックすると、どんなデータであるかの説明が見れる。. dhcp186-162:~ drkc167$ cp /Users/nor/Download/SRR029023.fastq.bz2 /Users/nor/drkc167/SRR029023.fastq.bz2. とする。改行は入っていない。半角スペースが入っている。ここでls コマンドを使えばSRR029023.fastq.bz2 が/Users/ nor/drkc167/ ディレクトリに入っていることを確かめられる。拡張子が.fastq.bz2になっていることに注目して欲しい。これはbzip2形式で圧縮されたfastqファイルであることを示し 13.

(15) ている。解凍しなければ使うことは出来ないので、解凍する。コマンドは、以下の通り。. dhcp186-162:drkc167 nor$ wc -l SRR029023.fastq 25079000 SRR029023.fastq. 今回のSRR029023.fastqではIDが SRR029023.1. dhcp186-162:drkc167 nor$ bzip2 -dc SRR029023.fastq.bz2 > SRR029023.fastq. HWI-EAS440_42_2_1_N_Nとなっている。また、シーケン. “>” 以下を書かないとターミナル上にfastqファイルが展開さ. スの本数は25079000/4で、6269750本あることもわかった。. れるので試してみても良いでしょう。ls で確認すると、. 次に、リファレンスファイルを準備する。生物種によって、. dhcp186-162:drkc167 nor$ ls -l. 公開されているリファレンスファイルの状態はさまざまであ. total 1738520. るため、この作業は特に注意を必要とする。今回はキイロ. -rw-r--r-- 1 nor staff 768425977 1 22 19:36 SRR029023.fastq. ショウジョウバエを使うので、比較的容易である。Ensemble. -rw-r--r--@ 1 nor staff 121695777 1 22 19:34 SRR029023.fastq.bz2. ファイルサイズが大きくなっていることが分かる。次にmore コマンドでfastqファイルの中身を、wcコマンドでシーケンス. で、cDNAリファレンスファイルを探す。 http://asia.ensembl.org/info/data/ftp/index.html. の本数を数える。 dhcp186-162:drkc167 nor$ more SRR029023.fastq @SRR029023.1 HWI-EAS440_42_2_1_1600_2041 ACCGAAGCCATACGAGAGGCCCGCGGACTTCATCGA + IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII2I @SRR029023.2 HWI-EAS440_42_2_1_785_608. 図のように多くの生物種のリファレンス配列が公開されてい. ACGAGTTCCAGAACTTGCACGCCGTTGCCAAGGAGA. る。下の方を探して行くと、. + IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIAIG @SRR029023.3 HWI-EAS440_42_2_1_1408_814 ATAAAGTTCATGGCATTTGAAGAAAAACTGTGTTTA + IIIIII+IIIIIIIIIIIIII1IIIIIIIIIIIIII. キイロショウジョウバエ Drosophila melanogaster を見つけた。cDNA配列が欲しいので、左から２番目の[FASTA]をクリック。３つのファイルが保存されたFTPサイトに繋がる。 14.

(16) とする。前のセクションでインデックス化を行った時と比べて、長い時間がかかるはずだ。genomic DNAのリファレンスだともっと長くかかる。インデックス化が終わったらマッピこういうところの[README]は良く読んでおいた方が良い。. ングを行う。コマンドは以下のように。. 今回は、allかab initioのどちらかを選ぶわけだが、README によると、allは知られている全部のcDNAが載っていて、ab initioの方にはゲノムから遺伝子構造を予測するツールを使って作成した予測cDNAが載っていると書いてある。どちらを選んでも今回は問題がないが、allの方を使うことにする。リンクをクリックするとダウンロードが始まり、通常/Users/ nor/Download/ ディレクトリ内に保存される。これを/ Users/nor/drkc167/ ディレクトリに移動させ、以下のコマンドで解凍する。. n0rr:drkc167 nor$ /Users/nor/analtools/bowtie-0.12.7/bowtie -l 36 -n 2 -p 4 Drosophila_melanogaster.BDGP5.65.cdna.all SRR029023.fastq > bowmapped_l36n2_SRR029023.txt. 念のため説明すると、まず/Users/nor/analtools/ bowtie-0.12.7/bowtie でbowtieまでのpathを記述し、bowtie を使う。次に-l 36 でシーケンスの何塩基をマッピングで使うのかを宣言する。今回はシーケンスの全長を示す。さらに-n 2 で１本のシーケンスあたりいくつのミスマッチを許容するか宣言する。最後の-p 4 のオプションで、CPUを４つ使うことを宣言している。続いて、リファンレスを指定し、シーケ. n0rr:drkc167 nor$ gunzip Drosophila_melanogaster.BDGP5.65.cdna.all.fa.gz. ンスを指定する。こうしてマッピングされた結果を”>” を使ってテキストファイル(bowmapped_l36n2_SRR029023.txt)に. ここまでで、シーケンスデータとリファレンスデータが. っ. たのでここからは前回のセッションと同様にbowtieを使った. 保存するようになっている。さて、more などで結果を見てみると、. リファレンスのインデックス化と、マッピングを行って行く。まずインデックス化を行う。bowtieは前のセクションで解凍したままになっているはずなので、コマンドは n0rr:drkc167 nor$ /Users/nor/analtools/bowtie-0.12.7/bowtie-build Drosophila_melanogaster.BDGP5.65.cdna.all.fa Drosophila_melanogaster.BDGP5.65.cdna.all. ７つのカラムができている。ひとつめのカラムにはスペース 15.

(17) が入っているようだ。カラムを順番に見て行くと、. に何本のショートリードがマッピングされたかを数えること. SRR029023.12 HWI-EAS440_42_2_1_1326_1394. ができ、すなわち発現量の推定となる。その前にマッピング. シーケンスのIDが書いてある。. の統計を見ておこう。bowtieが終わると、 # reads processed: 6269750. +. # reads with at least one reported alignment: 5303371 (84.59%). ワトソン鎖にマップされたかクリック鎖にマップされたか。. # reads that failed to align: 966379 (15.41%). 今回はサンプル調整の際にそれを見分けていないので特に意. Reported 5303371 alignments to 1 output stream(s). 味がない。. という結果が表示される。マップ率が余りにも低い場合には. FBtr0089424 当たったリファレンスの名前。今回はここが最重要。. ナのプロトコルでコントロールとして行われるphiXというファージのゲノムリシーケンシングでは98.5%が通常マップさ. 444 当たったリファレンス上の位置。 AGCGCATCCGCGTCAAGCTCGACGGCTCCCAGCTGG シーケンス IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII8 シーケンスのクオリティ 3. シーケンスに問題がなかったか調査する必要がある。イルミ. マッピングのクオリティ（0が一番いい）. れ、ヒトゲノムのリシーケンシングでは80%程度というから、ゲノムの複雑さにも依存するようである。また、カイコのmRNAでは50-60%という報告があり、リファレンスの整備状況にも依存していることがうかがえる。さて、マッピングデータの解釈に戻るが、発現量の推定には３つ目のカラムに含まれるリファレンス情報を集めて数えれば良いと書いた。さっそくunixを使って行う。コマンドを以下に示す。. 19:G>C リファレンスはGだけどシーケンスはCだった。さて、これらの情報の中で、今回最も大事なのは３つめのカラムである。どのトランスクリプトにヒットしたかが記述されているので、これを数えるとそれぞれのトランスクリプト. n0rr:drkc167 nor$ more bowmapped_l36n2_SRR029023.txt | awk -F "\t" '{print $3}' | sort | uniq -c | awk '{print $2"\t"$1}' > count_bowmapped_l36n2_SRR029023.txt. 紙面の都合上改行されているが、実際には改行は入っていない。半角スペースで区切られているだけだ。コマンドを順に 16.

(18) 説明すると、最初のmore で見る。これを終了させるには[con-. A. trol]と[c]を押す。今回は終わらせずに”|” で次の処理に繋げ. A. る。次の処理はawk だ。awk はunixに入っているスクリプト. B B. 言語。-F はフィールドの区切りを示していて、今回は””に囲まれた\t ごとに区切ることを宣言している。\t はタブを示す. C. 重複をまとめる. 表現である。続く’{print }’ で出力を示す。今回は、’{print $3}’ となっているので、タブ区切りで３つ目のカラムを出力するようになる。これで、どのトランスクリプトに当たった. n0rr:drkc167 nor$ more ABC.txt | sort | uniq A B. かの情報だけを集めてくることができる。続くsort では、行を順番に並び替える処理をしている。並び変わったところでu-. C. いくつのをまとめたのか数える. niq コマンドに-c オプションを付けて重複を数える。ここまでで、カウント数、リファレンス名となっているので、次の awkで$1と$2の順番を入れ替えて保存している。A, B, Cの並んでいるファイルを例にすると、もとのファイルを n0rr:drkc167 nor$ more ABC.txt. n0rr:drkc167 nor$ more ABC.txt | sort | uniq -c 2A 2B 1C. awkで整える. A n0rr:drkc167 nor$ more ABC.txt | sort | uniq -c | awk '{print $2"\t"$1}' B A. 2. B. 2. C. 1. C A B. 並び替えて n0rr:drkc167 nor$ more ABC.txt | sort | more. 17.

(19) 以上の操作で発現量の推定が出来た。結果のファイルを見て. データをcDNAにマッピングして発現量を推定できるように. みよう。１列目に遺伝子名が、２列目に発現量（その遺伝子. なった。今回の要約は以下の通りである。. にマッピングされたショートリードの数）が示されている。 n0rr:drkc167 nor$ more count_bowmapped_l36n2_SRR029023.txt 1 FBtr0005009. 2. FBtr0005088. 642. FBtr0005673. 100. FBtr0005674. 218. FBtr0006151. 16. FBtr0070003. 452. FBtr0070007. 474. 1. 解析用のディレクトリを作る (mkdir) 2. シーケンスファイルのダウンロード 3. シーケンスファイルの解析ディレクトリへの移動 (cp) 4. シーケンスファイルの解凍 (bzip2 -dc) 5. シーケンスデータの確認 (more, wc -l). 結果のファイルの一行目の 6. リファレンスファイルのダウンロード 1. は要らないのでテキストエディタで開いて削除してしまお. 7. リファレンスファイルの解析ディレクトリへの移動. う。実際のデータでは下に発現していた遺伝子の種類の数だ. 8. リファレンスファイルの解凍 (gunzip). け続いている。 wc コマンドで数えると、18774行あり、つまり、18774個の遺伝子の発現がこの実験で推定されたことに. 9. リファレンスファイルのインデックス化 (bowtie-build). なる。リファレンスのcDNA数は. 10.マッピング (bowtie). n0rr:drkc167 nor$ more Drosophila_melanogaster.BDGP5.65.cdna.all.fa | grep ">" | wc -l. 11.マッピング結果の確認 (more). 25415. より、25415個であるので、18774/25415=74%の遺伝子が発. 12.発現量の推定(more, awk, sort, uniq -c). 現していたことになる。ここまでで、mRNAのシーケンス 18.

(20) セクション 3. 前回のセクションでは、bowtieのマッピングデータを解析. 多型解析. してそれぞれのcDNAの発現量を推定した。リファレンスデータとシーケンスデータの準備はできるようになっていただけただろうか。このセッションではマッピングデータを解析して多型解析を行ってもらう。bowtieでマッピングする際に-n オプションで１つのリードあたりのミスマッチ塩基数許容量を指定できる。前回のセッションはこれを２としたので、１. できるようになること 1. 解析ツールをダウンロードしてコンパイルすることが出来る。 2. シーケンスデータをリファレンス配列にマッ. 本のリードの中でリファレンス配列との非一致塩基が２つ以内のリードだけがマッピングされている。通常、シーケンスしたサンプルのゲノムがリファレンス配列を作成する時に用いられたサンプルのゲノムと完全に一致することはなく、 SNPやIndel などの多型を含む部分のマッピングはミスマッチ. ピングして、SNPやIndel などの多型候補の抽. 許容によってうまく行われている。これらのミスマッチ部位. 出ができる。. の情報を集めることで、シーケンスしたサンプルのもつSNP やIndel などの多型について知ることが出来る。bowtieはミスマッチを許容することが出来るためSNPを持つ配列のマッピングが可能である一方で、Indelを持つ配列をマッピングすることは出来ない。Indelを含む多型の解析にはbwaなど他のマッピングツールの利用が必要であるので、今回のセッションでは、ヒトゲノムのエキソームデータを用い、bwaを用いたマッピングと多型抽出を行う。. 19.

(21) 前回のセッションと同様に、シーケンスデータとリファレン. する。ペアの場合には、シングルリードと同様に片側をシー. スデータの準備を行い、続いてマッピング、その後にマッピ. ケンスした後に、反対側のアダプターからシーケンスを行う. ングデータを解析して多型の抽出を行ってもらう。. ので、結果として、フローセル上の１つのクラスターから前. 今回扱うデータはHapMapにも登録されている日本人女性のエキソームデータだ。IDはH_IJ-NA18973となっている。 DRAで探すと、 Experiment: http://trace.ddbj.nig.ac.jp/DRASearch/experiment?acc=SRX 033211. 後２種類のシーケンスが得られる。今回は、 SRR077861_1.fastq.bz2 と SRR077861_2.fastq.bz2 の２つのファイルとして保存されているので、それぞれをダウンロードしよう。SRR077861.fastq.bz2という小さなファイルも置いてあるが、これはペアエンドの対応がうまくつかないはぐれシーケンスなので今回は無視する。ちなみにエキソームシーケンスとは、DNAの中からエキソンのところをプローブで集. Study: http://trace.ddbj.nig.ac.jp/DRASearch/study?acc=SRP00407 6. めて来てシーケンスした、偏らせたDNAシーケンシングの手法のことである。ゲノムのほとんどはcDNAでない部分である. Run: http://trace.ddbj.nig.ac.jp/DRASearch/run?acc=SRR077861. ので、多型の影響がより大きいと考えられるエキソーム部分. で見つかった。Runのリンクから、前回のセッションと同様. 行うことを目指している。ゲノムリシーケンスはこちらにあ. に.fastq.bz2形式のシーケンスデータをダウンロードする。前. る：http://trace.ddbj.nig.ac.jp/DRASearch/experiment?. 回と大きく異なっている点は、今回のシーケンスデータがペ. acc=ERX000157. アエンドのものであることだ。ペアエンドの詳細については. まずはターミナルを起動し、/Users/nor/. イルミナ社HPなどに詳しいのでそちらを参考にしていただき. ディレクトリ以下に今回の解析ディレクトリを作成する。. たい。シングルリードでは、前後をアダプター配列に挟まれ. を濃縮することで、重要な部位についての深度の高い解析を. dhcp186-162:~ nor$ mkdir homo_SNP_SRR077861. たサンプル由来のDNAを片方のアダプター側からシーケンス 20.

(22) 先ほどダウンロードしたシーケンスファイルを解析ディレクトリに移動させ、bzip2 コマンド（オプションはdc）を用いて解凍する。 dhcp186-162:~ homo_SNP_SRR077861$ bzip2 -dc SRR077861_1.fastq.bz2 > SRR07786_1.fastq dhcp186-162:~ homo_SNP_SRR077861$ bzip2 -dc SRR077861_2.fastq.bz2 > SRR07786_2.fastq. ここまでで、シーケンスファイルとリファレンスファイルがった。次に、今回はじめて用いるマッピングツールbwa のダウンロードとコンパイルを行う。bwaのsource forge のページ http://sourceforge.net/projects/bio-bwa/files/ からダウンロード出来る。現在最新版は0.6.1であるが、ここで例示する. ここまででうまくいかない場合は、セクション２の13ページ. のは0.5.9 である。最新版を使ってみて、同じコマンドが通ら. を参照して欲しい。続いてリファレンスファイルを取得す. ない場合にバージョンの落ちたものを使うと良いのではない. る。以下のリンクから. だろうか。さっそくダウンロードし、/Users/nor/analtools/. ftp://genome-ftp.cse.ucsc.edu/goldenPath/hg19/bigZips/. ディレクトリに移動させよう。次世代シーケンス解析ツール. chromFa.tar.gz をダウンロードする。シーケンスファイルと同様に解析ディレクトリに移動させ、解凍する。. は更新が早いので、自分が使っている間にも更新されて行く。バージョンを分けて管理するために、/Users/nor/ analtools/ ディレクトリ内に、bwa ディレクトリを作成し、/. dhcp186-162:~ homo_SNP_SRR077861$ tar zxvf chromFa.tar.gz. 解凍出来たら、cat コマンドでまとめる。. Users/nor/analtools/bwa/ ディレクトリ以下にダウンロードしたbwa-0.6.1.tar.bz2を置いても良い。私の環境ではそう. dhcp186-162:~ homo_SNP_SRR077861$ cat *.fa > chr_all.fa. なっている。次のコマンドで解凍する。. として、ファイル名が.faで終わる全てのファイルをひとつの. n0rr:nakagawa nor$ cd /Users/nor/analtools/bwa/. マルチファスタファイルにまとめる。今回の作業では多くのファイルが生成されて煩雑になるので、リファレンス用のディレクトリを作ってしまってしまおう。. n0rr:bwa nor$ tar zxvf bwa-0.6.1.tar.bz2. 解凍されると、/Users/nor/analtools/bwa/bwa-0.6.1/ ディレクトリができる。そこに移動して、コンパイルを行う。. dhcp186-162:~ homo_SNP_SRR077861 nor $ mkdir reference n0rr:bwa nor$ cd bwa-0.6.1 dhcp186-162:~ homo_SNP_SRR077861 nor $ mv chr_all.fa ./reference/chr_all.fa n0rr:bwa-0.6.1 nor$ make. 21.

(23) このmakeがコンパイルである。これをするとb-. n0rr:bwa-0.6.1 nor$ /Users/nor/analtools/bwa/bwa-0.6.1/bwa. waが使えるようになる。引き続いてmake install すること. Program: bwa (alignment via Burrows-Wheeler transformation). で、bwaを/usr/bin/ ディレクトリに置いてpathなしで動かす. Version: 0.6.1-r104 Contact: Heng Li <[email protected]>. ことも出来るが、前述の通りバージョン管理が重要なので、あえてこれを行わない。. Usage: bwa <command> [options]. Command: index. さて、コンパイルはうまく行っただろうか。. aln. index sequences in the FASTA format. gapped/ungapped alignment. samse. generate alignment (single ended). sampe. generate alignment (paired ended). bwasw. BWA-SW for long queries. -bash: make: command not found. fastmap. identify super-maximal exact matches. と表示されて止まっていないだろうか。xcode が入っていな. fa2pac. norr:bwa-0.6.1 nor$ make. い場合にはうまくいかない。また、私の環境ではSnow Leopard からLion にアップグレードした時にどこかに行ってしまったので再インストールが必要であった。xcodeは以下のリンク https://developer.apple.com/xcode/ から入手することができる。指示に従ってインストールして欲しい。ファイルサイズが3G程度あり、私の環境ではダウンロードに30分かかった。また、そのままインストールすると10.5sdkが消されてMacOSx10.5用のアプリ開発が不可能になるので、他所にコ. convert FASTA to PAC format. pac2bwt. generate BWT from PAC. pac2bwtgen alternative algorithm for generating BWT bwtupdate bwt2sa pac2cspac stdsw. update .bwt to the new format generate SA from BWT and Occ convert PAC to color-space PAC. standard SW/NW alignment. 以上のようにオプションを表示させるはずだ。ここからbwa を使ってマッピングを進める。bwaはbowtieと比べてコマンドを入力する機会が多いので、bwaまでのpathを先に入れてしまっても良いかもしれない。. ピーした。アプリ開発の予定はないので消えてしまっても問. n0rr:bwa-0.6.1 nor$ PATH=$PATH:/Users/nor/analtools/bwa/bwa-0.6.1/. 題ない。. 入れておくと、bwaと打っただけで動くようになる。. うまく行っていれば、pathを書けばbwaが動いて. n0rr:bwa-0.6.1 nor$ bwa. 22.

(24) 今回のマッピングでも、まずはリファレンスファイルのインデックス化から始まる。リファレンスファイルは /Users/. dhcp186-162:homo_SNP_SRR077861 nor$ /Users/nor/analtools/bwa/bwa-0.5.9/bwa aln ./reference/chr_all.fa SRR077861_2.fastq > aln_sa2.sai. ２つのマッピング結果をマージして結果ファイルを得る。. nor/homo_SNP_SRR077861/reference ディレクトリ内にあり、今は、/Users/nor/homo_SNP_SRR077861 ディレクトリにいるので、コマンドは以下のようになる。ここの例では bwa 0.5.9を使っていることに注意して欲しい。 dhcp186-162:homo_SNP_SRR077861 nor$ /Users/nor/analtools/bwa/bwa-0.5.9/bwa index -a bwtsw ./reference/chr_all.fa. dhcp186-162:homo_SNP_SRR077861 nor$ /Users/nor/analtools/bwa/bwa-0.5.9/bwa sampe ./reference/chr_all.fa aln_sa1.sai aln_sa2.sai SRR077861_1.fastq SRR077861_2.fastq > bwa_aln.sam. 一時間半ほどで完了する。ここまででマッピングは終わりだ。結果ファイルを見てみよう。結果ファイルの拡張子に注目して欲しい。.sam はSequence Alignment/Map (SAM)形式. リファレンスファイルの指定に新しい表現を使っている。./. と呼ばれ、次世代シーケンス解析で頻繁に使われるフォー. reference/chr_all.fa によって、今いるディレクトリ内にある. マットである。詳細な説明が以下のリンクに示すsamtools の. reference ディレクトリ内の chr_all.fa を指定しているのだ。. ページに記載されているので、参照していただきたい。. もちろんここを. http://samtools.sourceforge.net/samtools.shtml 実はbowtie. /Users/nor/homo_SNP_SRR077861/reference/chr_all.fa と. もマッピングの結果をSAM形式で出力することができる。や. パスを全て書いても構わない。私の環境では１時間程度でイ. り方はいたって簡単で、-S オプションをつけるだけである。. ンデックス化が終わった。続いてシーケンスの片側のマッピ. 引き続いて多型抽出を行う。多型抽出にはsamtools. ングである。bwaではペアのリードを別々にマッピングして. http://samtools.sourceforge.net/ を用いる。これもsource-. 後でマージするようになっている。. forge のリンクからダウンロードすることができる。 http://sourceforge.net/projects/samtools/files/samtools/. dhcp186-162:homo_SNP_SRR077861 nor$ /Users/nor/analtools/bwa/bwa-0.5.9/bwa aln ./reference/chr_all.fa SRR077861_1.fastq > aln_sa1.sai. samtools もいくつかのバージョンを使い分けることがあるの. これも私の環境では２時間半かかった。続いて反対側のシー. で、/Users/nor/analtools ディレクトリにsamtools のディレ. ケンスのマッピング。同じだけ時間がかかる。. クトリを作って、中に各バージョンのsamtools を格納しておくのが良い。しつこいようだがコマンドを書くと、 23.

(25) dhcp186-162:homo_SNP_SRR077861 nor$ cd /Users/nor/analtools/ dhcp186-162:analtools nor$ mkdir samtools dhcp186-162:analtools nor$ cd samtools. ここにダウンロードしたsamtools-0.1.18.tar.bz2を移し、 dhcp186-162:analtools nor$ tar zxvf samtools samtools-0.1.18.tar.bz2. と解凍する。そして、samtools-0.1.18のディレクトリで dhcp186-162:analtools nor$ cd samtools samtools-0.1.18 dhcp186-162:analtools-0.1.18 nor$ make. としてコンパイルを行う。本セッション２度目のコンパイルだが、慣れて来ただろうか。今後必要なツールが現れる度に. bamファイルを並び替えている。この状態のファイルを入力ファイルとして必要とするツールもあるので、この、ソートされたバイナリーファイルの存在を良く覚えておくこと。次はいよいよ多型抽出である。 dhcp186-162:homo_SNP_SRR077861 nor$ /Users/nor/analtools/samtools/samtools-0.1.18/samtools mpileup -uf ./reference/chr_all.fa bwa_aln.sam.bam.sorted.bam | /Users/nor/analtools/samtools/samtools-0.1.18/bcftools/bcftools view -bvcg - > bwa_aln.var.raw.bcf. samtools のmpileup を使って多型抽出を行う。（-f オプションだけにすることで１塩基単位のカバレッジをみることもできる）多型ファイルをbcf形式で保存する。そのうち、信頼性の高いものを抽出するためにvarFilterを使う。. 行う必要があり、場合によってはコンパイルがうまく行かないものを解決しなければいけないこともある。では、さっそ. dhcp186-162:homo_SNP_SRR077861 nor$ /Users/nor/analtools/samtools/samtools-0.1.18/bcftools/bcftools view bwa_aln.var.raw.bcf | . くsamtoolsを使って、多型解析を進めていこう。. /Users/nor/analtools/samtools/samtools-0.1.18/bcftools/vcfutils.pl varFilter -D 500> bwa_aln.var.vcf. dhcp186-162:homo_SNP_SRR077861 nor$ /Users/nor/analtools/samtools/samtools-0.1.18/samtools view -b -o bwa_aln.sam.bam -S -t ./reference/chr_all.fa.fai bwa_aln.sam. 以上で、vcf 形式の多型ファイルが得られた。vcf形式の詳細についても、samtools のHP. ここでは、マッピング結果のsamファイルをsamtools view を. http://samtools.sourceforge.net/samtools.shtml に記載され. 使ってbamファイルに変換している。bam はバイナリーファ. ているので参考にしていただきたい。. イルなので、more コマンド等で中身を覗いても読めない。 dhcp186-162:homo_SNP_SRR077861 nor$ /Users/nor/analtools/samtools/samtools-0.1.18/samtools sort bwa_aln.sam.bam bwa_aln.sam.bam.sorted. ここまでをまとめると、 1. ファイルの解凍(bzip2 -dc/ tar zxvf) 2. ファイルの移動(mv) 24.

(26) 3. コンパイル(make) 4. マッピング(bwa) 5. ファイルの変換(samtools) 6. 多型抽出(samtools) Mac Book Air を使っている場合に、ファイルサイズが大きくなりすぎて入りきらない問題が確認されている。ファイルの出力先を外付けでデバイスなどにして解決することもできるが、それではおっくうなので、もう少し小さなファイルで完了できるE. coli. の多型抽出を付記する。シーケンス. 今回は/Users/nor/otherPJ/macde/ ディレクトリ内にecoli ディレクトリを作ってそこで作業する予定なので、もしotherPJ ディレクトリがなければ $nor cd /Users/nor/ $nor mkdir otherPJ. ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/fastq/SRA008/S RA008271/SRX003267 以下のSRR014475.fastq.bz2. $nor cd otherPJ $nor mkdir macde. リファレンス $nor cd macde http://www.ncbi.nlm.nih.gov/nuccore/NC_000913 画面右のSend より”File” “FASTA” を選んで”Create File” をクリックするとsequence.fasta というファイルが/Users/nor/ Downloas/ ディレクトリ以下にダウンロードされる。. $nor mkdir ecoli $nor cd ecoli とする。/Users/nor/Download/ ディレクトリ内のシーケンスファイル(SRR014475.fastq.bz2)とリファレンスファイル 25.

(27) (sequence.fasta)をFinder 上で/Users/nor/otherPJ/macde/ ecoli/ ディレクトリまで移動させても構わないが、今回はmv コマンドを用いる。 nor$ mv /Users/nor/Downloads/SRR014475.fastq.bz2 ./. 思ったよりマップ率低い。ここからsamtoolsである。. 最後の./ は今いるディレクトリの中に、という意味である。. fai作り. 同じようにリファレンスも取ってくる。 nor$ mv /Users/nor/Downloads/sequence.fasta ./. nor$ /Users/nor/analtools/samtools/samtools-0.1.18/samtools faidx sequence.fasta bam化. 今回のマッピングはbowtie でやることにする。まずはリファレンスのインデックスを作る。 nor$ /Users/nor/analtools/bowtie-0.12.7/bowtie-build sequence.fasta sequence 続いてシーケンスの解凍とマッピング nor$ bzip2 -cd SRR014475.fastq.bz2 > SRR014475.fastq nor$ /Users/nor/analtools/bowtie-0.12.7/bowtie -n 2 -l 36 -p 4 -S sequence SRR014475.fastq > bowmap_SRR014475.sam. nor$ /Users/nor/analtools/samtools/samtools-0.1.18/samtools view -b -o aln.sam.bam -S -t sequence.fasta.fai bowmap_SRR014475.sam sort nor$ /Users/nor/analtools/samtools/samtools-0.1.18/samtools sort aln.sam.bam aln.sam.bam.sorted 多型抽出 nor$ /Users/nor/analtools/samtools/samtools-0.1.18/samtools mpileup -uf sequence.fasta aln.sam.bam.sorted.bam | 26.

(28) /Users/nor/analtools/samtools/samtools-0.1.18/bcftools/bcft ools view -bvcg - > aln.var.raw.bcf varfilter nor$ /Users/nor/analtools/samtools/samtools-0.1.18/bcftools/bcft ools view aln.var.raw.bcf | /Users/nor/analtools/samtools/samtools-0.1.18/bcftools/vcf utils.pl varFilter > aln.var.vcf. PCR duplicate を除いたbamから多型抽出 nor$ /Users/nor/analtools/samtools/samtools-0.1.18/samtools mpileup -uf sequence.fasta aln.uniq.bam | /Users/nor/analtools/samtools/samtools-0.1.18/bcftools/bcft ools view -bvcg - > aln.uniq.var.raw.bcf つづいてvarfilter. このaln.var.vcf が多型ファイルになる。. $nor /Users/nor/analtools/samtools/samtools-0.1.18/bcftools/bcft ools view aln.uniq.var.raw.bcf | /Users/nor/analtools/samtools/samtools-0.1.18/bcftools/vcf utils.pl varFilter > aln.uniq.var.vcf. さらに１塩基ごとのカバレッジを求める. あまりかわらなかった。. nor$ /Users/nor/analtools/samtools/samtools-0.1.18/samtools mpileup -f sequence.fasta aln.sam.bam.sorted.bam > aln.cov.txt PCR duplicate を除いたbamからカバレッジ求める samtools rmdup でPCR duplicate を除いてみる nor$ /Users/nor/analtools/samtools/samtools-0.1.18/samtools rmdup -s aln.sam.bam.sorted.bam aln.uniq.bam. nor$ /Users/nor/analtools/samtools/samtools-0.1.18/samtools mpileup -f sequence.fasta aln.uniq.bam > aln.uniq.cov.txt 参考までに。 27.

(29) チャプタ 2. Rをつかってやるところ. 発現変動遺伝子抽出カバレッジの統計をとる.

(30) セクション 1. ターミナル上でUNIXを使っておこなってきたのはテキスト. 発現変動遺伝子抽出. データの解析である。次世代シーケンサー使うことの大きな利点として、定量的であることが挙げられる。これを生かすためには、数の解析が必要である。ここではRを使って数の解析を行う。Rには非常に多くの統計パッケージがそろっていて便利である。. できるようになること. Rのホームページは http://www.r-project.org/ である。インストール方法については、このサイト. 1. カウントデータを整形してRに読み込ませることができる 2. Rを使って簡単なグラフを書くことができる 3. Rのパッケージを使って発現変動遺伝子抽出ができる. http://cse.naro.affrc.go.jp/takezawa/r-tips/r/01.html で非常に詳しく日本語で書かれている。筑波大学のミラーサイト http://cran.md.tsukuba.ac.jp/bin/macosx/ よりダウンロードする2012/03/24現在での最新バージョンは2.14.2であるようだ。私の環境では2.13.0を使っている。 Rを用いた次世代シーケンスデータ解析については、門田先生のホームページが丁寧でわかりやすい。http:// www.iu.a.u-tokyo.ac.jp/~kadota/r_seq.html かなり勉強させて頂いたし、パッケージのインストールの不具合については直接教えてもいただいた。ここではデータの整形を中心に一通りの発現変動遺伝子抽出をデモンストレーションする。今回でもデータとするのはSRAに登録されているマウスの肝臓 RNA-seq データ 29.

(31) http://trace.ddbj.nig.ac.jp/DRASearch/run?acc=SRR001358 と、同じくマウスの骨格筋データ http://trace.ddbj.nig.ac.jp/DRASearch/run?acc=SRR001361 である。SRR001358, SRR001361 のそれぞれのアクセッション番号を検索していただくとすぐに出てくる。これをbowtie. これら２つのファイルを合体、整形してRの読み込めるファイルを作ろう。今のファイルには、どちらか一方でしか発現していない遺伝子の記述に問題がある。ここで使うのはjoin というコマンドである。早速作業を始める。. を用いてcDNA(. まずはどちらか一方でも発現した遺伝子についてのリストを. ftp://ftp.ensembl.org/pub/release-66/fasta/mus_musculus/cdna/ 内の. つくる。. Mus_musculus.NCBIM37.66.cdna.all.fa.gz). にマッピングし、-S オプ. ションで.sam 形式の出力を得てそれを awk -F "\t" '{print $3}' | sort | uniq -c の一連のコマンドラインでカウントデータに変換した。作業は/Users/nor/ ディレクトリ以下に/ Users/nor/otherPJ/macde/ ディレクトリを作って行った。. nor$ cat gene_counts_liver.txt gene_counts_skemuscle.txt | awk '{print $2}' | sort | uniq > genename.txt たくさんのデータがある場合には、リファレンスファイルからこのリストを作ってしまった方が楽だ。その場合には $ more Mus_musculus.NCBIM37.66.cdna.all.fa | grep ">" | awk -F " " '{print $1}' | tr -d ">" | sort > genename.txt とすればよい。次に、カウントデータファイルを整形してリストファイルとのjoin に備える。 $ more gene_counts_liver.txt | awk ‘{print $2”\t”$1}’ > temp1. このようになっている。骨格筋側もあるので、. 30.

(32) カラムの順序を入れ替えただけである。join は一列目を見るので、このように整形しておくことが重要である。また、 sort されていることも重要で、今回の場合には発現量を推定する操作で遺伝子名がsort された状態になっているのでこの. コマンドは nor$ join -1 1 -2 1 -a 1 -a 1 genename.txt temp1 | awk '{print $1"\t"$2+0}' | more これをtemp11に書き出して. まま使える。tophat やCASAVA など他の方法で得た発現量データを扱う際にはこの点をよく注意して欲しい。もしsort が必要な場合には more filename | sort > filename2 とすれば. nor$ join -1 1 -2 1 -a 1 -a 1 genename.txt temp1 | awk '{print $1"\t"$2+0}' > temp11. よい。引き続き骨格筋のファイルを整形する。. とする。同様に骨格筋の方も. nor$ more gene_counts_skemuscle.txt | awk '{print $2"\t"$1}' > temp2. nor$ join -1 1 -2 1 -a 1 -a 1 genename.txt temp2 | awk '{print $1"\t"$2+0}' > temp21. 整形したそれぞれの発現量ファイルを遺伝子リストとjoin さ. ここで、wc コマンドを使って行数が一致していることを確か. せる。. めよう。. temp1 とtemp2 は行数が違うのにtemp11 とtemp21の行数が一致していることがわかる。UNIX上のテキスト操作では、狙い通りにものが進んでいないことがあるので、こうやって処 31.

(33) 理の度に処理の完了を確認していくことが重要である。最後. ３列目はいらないのでawk で削除。これをtemp3として書き. にこれらのファイルを一つにまとめてヘッダーをつけよう。. 出す。. まとめるのにはpaste コマンドを用いる. nor$ paste temp11 temp21 | awk '{print $1"\t"$2"\t"$4}' > temp3 次にヘッダーを作ります。 $ vi header.txt とするとvimというテキストエディタが起動するのでi を押してINSERTモードにして. 左側（肝臓）と右側（骨格筋）の遺伝子名が一致していることを確かめてほしい。片方でだけ発現している遺伝子の項目には、もう片方で０が入っているのがポイントだ。. genenameタブliverタブskmus、と書き込んで[esc]ボタンを押し、[:][w][q]と入力しエンター。vim の使い方はここに詳しい。http://archiva.jp/web/tool/vim_basic.html もちろんテキストエディタを開いて同様の入力をしたファイルをheader.txtという名前で保存しても同じことであるが、改行コードがCRになってしまったりするので、mi http://www.mimikaki.net/ を使って改行コードをLFに直す必 32.

(34) 要が出てくるかもしれない。. では、これを先のtemp3と合体させる。 . これが正解。これに好きな前をつけて書き出して完成だ。 nor$ cat header.txt temp3 > for_R_liv_vs_skm.txt 次にこれをRに読み込ませよう。Rを起動すると、. こういった感じである。データを読み込ませるにはデータまでのパスを書けばよいので、 data<-read.table("/Users/nor/otherPJ/macde/for_R_liv_vs_ skm.txt",header=T,row.names=1) とする。簡単に説明すると、”data”として（名前は好きなもので何でもよい）/Users/nor/otherPJ/macde/ ディレクトリ 33.

(35) に入ってるfor_R_liv_vs_skm.txt を読み込め、ヘッダーがつ. この図では一つの遺伝子をひとつの点として、liver での発現. いていて、それぞれのデータの名前は１列目に書いてある。. 量を水平軸に、skelmでの発現量を垂直軸にプロットしてい. という意味がある。続いて、data をattachしてliver とskelm. る。左上にあるのと、右下にあるのがなんとなくの発現変動. を読み込ませる。summary を使うとデータの確認ができるの. 遺伝子(Differentially Expressed Genes/ DEG) である。これに. で習慣としてここまでを一気に行うとよい。. DESeqという解析パッケージを用いてDEG抽出を行う。 DESeqのインストールには門田先生のHP http://www.iu.a.u-tokyo.ac.jp/~kadota/r_seq.html#install の通りにやればうまくいく。 source("http://www.bioconductor.org/biocLite.R") biocLite("DESeq", dependencies=TRUE). 試しにプロットしてみよう。 plot(log(liver),log(skelm),pch=20). としてDESeqを手に入れる。DESeqを選んだ理由は、これが反復のないデータ(n=1)でも解析できるからである。以下にコマンドラインを示す。 out_f <- "DESeq_out_liv_vs_skelm.txt" param1 <- 1 param2 <- 1 library(DESeq) data.cl <- c(rep(1, param1), rep(2, param2)) cds <- newCountDataSet(data, data.cl) 34.

(36) cds <- estimateSizeFactors(cds) sizeFactors(cds) cds <- estimateVarianceFunctions(cds, method="blind") out <- nbinomTest(cds, 1, 2) logratio <- out$log2FoldChange FDR <- out$padj rank_FDR <- rank(FDR, ties.method="min") tmp <- cbind(rownames(data), data, logratio, FDR, rank_FDR). 今回の比較ではFDRが十分に小さいものが得られなかった。私の. write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F). 手元のデータ（同じ方法でDEG55個(FDR<0.05)）のプロッ. 以上である。. える。ここまででDEG抽出は終了である。. トと比較すると発現量の小さい遺伝子がかなり多いように見. 結果は、/Users/nor/ディレクトリ以下に DESeq_out_liv_vs_skelm.txtという名前で格納されている。これをエクセルなどに貼り、フィルタを使ってFDR順に並び替えるとDEGが抽出できる。FDRはFalse Discovery Rate の略で、ここではp値と同じような判断の指標として扱っている。. 35.

(37) 続いて、もっとおおまかにトランスクリプトームをつかむ方. 左が肝臓で、右が骨格筋のものである。それぞれの遺伝子の. 法を紹介する。. 発現量（リード数）を多い順に並び替えてプロットした。も. par(mfrow=c(1,2)). うひとつ、ハエの脂肪体のデータを載せる。こちらの図は先のふたつよりもかなりかたよっているようだ。. barplot(matrix(sort(liver/sum(liver)))) barplot(matrix(sort(skelm/sum(skelm))),ylim=c(0,1)) とすると以下の図が描ける。. 36.

(38) セクション 2. ここではgenomic DNA のコピー数変異解析を行う。ショウ. カバレッジの統計をとる. ジョウバエのセルラインkc167のDNA リシーケンシングデータをショウジョウバエゲノムにマッピングし、リファレンス配列上の１塩基ごとのカバレッジ（マッピングされたショートリードの厚み）を数えてgenomic DNA のコピー数変異を解析する。シーケンスデータはSRR034211を、DRAよりダウンロードして用いた。リファレンスデータには、ハエゲノムを. できるようになること 1. カバレッジの統計をグラフを描いて眺めることができる. 用い、maqというマッピングツールを用いてマッピングした。maqは以下のサイトから入手する。 http://maq.sourceforge.net/ maq user’s manual http://maq.sourceforge.net/maq-man.shtml のMapass2 Work Flow は次世代シーケンス解析の概念が分かりやすくまとめられているので参考にするとよい。maqを解凍し、コンパイルする。PATHを通して、以下のように maq easyrunと引き続くpileupを実行する。 /Users/nor/analtools/maq-0.7.1/scripts/maq.pl easyrun ./Dmel_genome/no_space_dmg.fa srr034211.fastq && cd easyrun && maq pileup ref.bfa all.map > all_pileup.txt この結果をRを使って可視化したいので、以下のコマンドでR 向けにデータを整える。 awk '($1=="chr2L"){print $2"\t"$4}' all_pileup.txt > chr2L_pileup_r.txt && 37.

(39) awk '($1=="chr3L"){print $2"\t"$4}' all_pileup.txt > chr3L_pileup_r.txt && awk '($1=="chr2R"){print $2"\t"$4}' all_pileup.txt > chr2R_pileup_r.txt && awk '($1=="chr3R"){print $2"\t"$4}' all_pileup.txt > chr3R_pileup_r.txt && awk '($1=="chr4"){print $2"\t"$4}' all_pileup.txt > chr4_pileup_r.txt && awk '($1=="chrX"){print $2"\t"$4}' all_pileup.txt > chrX_pileup_r.txt && awk '($1=="chrMT"){print $2"\t"$4}' all_pileup.txt > chrMT_pileup_r.txt これをRに読み込ませてグラフを描くだけである。 chr3R_pileup_r.txt は１列目がポジション、２列目がカバ. 同様にしてすべての染色体について行った結果がこれだ。. レッジになっているので、前のセクションと同様にヘッダーをつけてRが読み込めるようにする。今回は、pos と count とつけた。あとは plot(pos,count,pch=20,col=ifelse(count < 1, "red", "black")) とすれば描ける。. 38.

(40) 簡単な統計も取ってみた。0rate というのはカバレッジが０だったものの数を染色体の長さで割ったものである。ミトコンドリアで0率が低いことがわかる。ただし、ミトコンドリアは平均カバレッジが13倍なのでこれによる可能性もある。平均カバレッジのゲノムvs ミトコンドリアの比較によって、１ゲノムあたりのミトコンドリアゲノム数、すなわち一細胞あたりのミトコンドリア数が10個弱であることも推定される。chr3Lに非常にカバレッジの高い部位があった。このデータはハエのkc167セルラインのものである。昆虫のセルラインが樹立される過程で染色体の異数化が起こること MTはミトコンドリアで、カバレッジが高いように見えるのは. が知られている。また、ほ乳類のセルラインではテロメア伸長酵素が発現していることが知られているが、昆虫細胞の場. 極端に高いのがないためにグラフの縦軸が小さいのだ。. 合はというと、発現していないことが明らかになっている。 chr. 0rate. Cov.. 逆に、正常に発生した個体の全身のさまざまな場所でテロメ. length. Zero. Max. chr3L. 24543557. 9480442. 0.3862701. 1.393538. 33323. chrX. 22422827. 9185534. 0.4096510. 1.298916. 5129. chrMT. 19517. 4297. 0.2201670. 13.291797. 91. chr2L. 23011544. 9162103. 0.3981525. 1.257954. 1650. chr2R. 21146708. 7369697. 0.3485033. 1.475968. 613. ラインでテロメア伸長は起きているのか調べたくてこれを. chr3R. 27905053. 10216401. 0.3661129. 1.243935. 509. やった。結果は、特定の部位のコピー数変異を見つけてい. chr4. 1351857. 776293. 0.5742420. 1.609679. 629. ア伸長酵素と考えられる遺伝子が発現している。そもそも、ハエのテロメアはトランスポゾンをsmall RNAで制御することで伸長しており、テロメア伸長酵素は必要ない。昆虫セル. る。これらの部位は反復配列ではなかった。. 39.

(41) チャプタ 3. その他の問題と対処. bcl変換とデマルチプレックスクオリティチェック＆コントロール記憶媒体のフォーマット.

(42) セクション 1. illumina社のGenome Analyzer Ⅱxは、bclファイルとstats. bcl変換とデマルチプレックス. ファイルとしてシーケンス結果を出力する。しかしながら、通常外注したシーケンスはインデックス毎に分かれたfastq ファイルで納品されるし、自前で読んだ場合でもCASAVAのマニュアルに従えばデマルチプレックス済のfastqが手に入るのでほとんど気にすることはない。だが、まれにここが上手く行かない場合があるらしいので、やり方を書いておく。. できるようになることまず、bcl変換だが、これにはIllumina社のCASAVAを使. 1. CASAVAを使ってillumina Genome Analyzer. う。CASAVAは有償ソフトウェアなので、illuminaユーザーで. Ⅱxの出力ファイルを変換し、fastqファイルを. ないと手に入れることが出来ない。手に入れたらインストー. 得る。. ル手順に従って、以下のように行う。CASAVAのバージョン. 2. 得られたfastqがマルチプレックスのインデッ. が異なる場合には、コマンドのv1.8.0のところを適宜変更すること。. クスを含む場合にデマルチプレックスをUNIX で行い、コマンドラインをまとめてbashプログラムにする。. $ tar -jxf CASAVA_v1.8.0.tar.bz2 $ mkdir CASAVA_v1.8.0-build $ cd CASAVA_v1.8.0-build $ ../CASAVA_v1.8.0/configure. 3. bashプログラムを実行する. $ make $ make install. 私の環境では、ダウンロードしたCASAVA_v1.8.0.tar.bz2 を/ Users/nor/analtools/ディレクトリに格納後、/analtools ディレクトリ内にてこれを行った。make installするのでどこで 41.

(43) 行っても同じである。これを使ってbcl変換を行うために、必. --output-dir /Volumes/NOGEBLANC/120114_5 --use-bases-mask y72. 要な入力ファイルの確認を行う。基本的には、Genome Ana-. 改行は入っていない。最後の--use-bases-mask y72 では、. lyzer の出力ファイルで 120114_HWUSI-. 65bp(シーケンス)+7bp(インデックスシーケンス)を合わせて. EAS1748_00049_FC こういう名前の物があればよい。ファ. １つの配列として出力するように書いた。これがうまく行く. イルの内容物についてはCASAVAのマニュアルに記載がある. と、. ので別途参照すること。外付けHDD内にこのファイルを置いてあってもbcl変換が可能であるが、外付けHDDの名前にスペースが入っているとうまくいかないことがあるのでアンダーバー等に変更しておくこと。次に、出力フォルダを決める。fastqは圧縮されたfastq.gzで出力されるので、70bp前後. [2012-01-22 23:09:49] [configureBclToFastq.pl]. INFO: Creating. directory '/Volumes/NOGEBLANC/120114_5' [2012-01-22 23:09:49] [configureBclToFastq.pl]. INFO: Basecalling software: RTA. [2012-01-22 23:09:49] [configureBclToFastq.pl]. INFO:. version: 1.09 (build 35) [2012-01-22 23:09:49] [configureBclToFastq.pl]. WARNING:. 'LocationFileType' element not found in /Volumes/LOGITEC_HD/test/120114_HWUSI-EAS1748_00049_FC/Data/Intensities/BaseCalls/../RTAConfiguration.xml. のシーケンスの場合には１フローセル辺り10GB程度になる。 fastqを解凍することを考えると、100GB程度の空き容量のあるメディア内で行うことが望ましい。今回は、外付けHDDで. [2012-01-22 23:09:49] [configureBclToFastq.pl]. INFO: Original. use-bases mask: y72 [2012-01-22 23:09:49] [configureBclToFastq.pl]. INFO: Guessed use-bases. mask: yyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyy,IIIIIIn Creating directory '/Volumes/NOGEBLANC/120114_5/Project_FC/Sample_lane1'. あるLOGITEC_HD内のtestというディレクトリに格納した出力ファイル /Volumes/LOGITEC_HD/test/120114_HWUSIEAS1748_00049_FC/ を入力ファイルにし、これまた外付け. Creating directory '/Volumes/NOGEBLANC/120114_5/Project_FC/Sample_lane2' Creating directory '/Volumes/NOGEBLANC/120114_5/Project_FC/Sample_lane3' Creating directory '/Volumes/NOGEBLANC/120114_5/Project_FC/Sample_lane4' Creating directory '/Volumes/NOGEBLANC/120114_5/Project_FC/Sample_lane5' Creating directory '/Volumes/NOGEBLANC/120114_5/Project_FC/Sample_lane6'. SSDであるNOGEBLANC内に /Volumes/NOGEBLANC/ 120114_5 という出力ファイルを作成して行うことにした。コマンドは、以下のようである。. Creating directory '/Volumes/NOGEBLANC/120114_5/Project_FC/Sample_lane7' Creating directory '/Volumes/NOGEBLANC/120114_5/Project_FC/Sample_lane8' [2012-01-22 23:09:50] [configureBclToFastq.pl]. WARNING: Couldn't. determine control software version due to missing /Volumes/LOGITEC_HD/test/120114_HWUSI-EAS1748_00049_FC/Data/Intensities/BaseCalls/../../../runParameters.x ml. nor$ configureBclToFastq.pl --input-dir. [2012-01-22 23:09:50] [configureBclToFastq.pl]. WARNING: use-base-mask. length is 72 for read 1: expected length from config file is 65. /Volumes/LOGITEC_HD/test/120114_HWUSI-EAS1748_00049_FC/Data/Inten sities/BaseCalls/. [2012-01-22 23:09:50] [configureBclToFastq.pl]. INFO: Read 1: length =. 72: mask = yyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyy. 42.