0から始めるエクソームデータ解析 v1.0
目次
○はじめに
○解析環境を整えよう
○解析用ツールのインストールと参照ファイルの入手
○エクソーム解析
○サンプルデータを用いた解析
○おわりに
下の図は、DNAシーケンサー1ランあたりで生成される遺伝子配列データ量の年次変化を示しています。
2005年から爆発的にデータ量が増えていることがひと目で分かると思いますが、これは、次世代シーケンサー技術の開発によるものです。
一方で、大量に生成されたデータを解析するための技術は、専門性が高まっており、必要な情報を集めようと思っても前提知識を要するものが多かったり、日本語の情報が少なかったり、という状況のため、やや敷居が高くなっている感があります。
本コンテンツでは、次世代シーケンサーのアプリケーションの一つである「エクソーム解析」の方法について、私自身、コマンドラインを使った解析についての知識が殆どない中から学んだ経験を生かして、できるだけ分かりやすく、解説したいと思います。
次世代シーケンサーの利用方法には様々なものがありますが、その王道は、やはりゲノムの配列解析でしょう。ただ、全ゲノム（ヒトの場合およそ30億塩基対）を解析するためには、まだまだ時間もお金もかかる、というのが現状です。
その中で、非常に効率がいい方法として流行しているのが、「エクソーム解析」です。この方法では、先にWetな実験でエクソン（Exon: 転写→翻訳を経て、最終的にタンパク質の配列を決定するゲノム中の領域）のみを網羅的にキャプチャーし、選択されてきた分子のみを対象として次世代シーケンサー実験を行います。エクソン全体＝エクソーム（Exon+ギリシャ語の「すべて・完全｣などを意味する接尾辞のomeでExome）のサイズは、全ゲノムの1-1.5%程度なので、これによって、コストと時間は大幅に削減され、一方で、タンパク質配列に影響を及ぼす重要な部分の配列情報は網羅的に得ることができます。エクソーム解析の主な目的は、変異と疾患などの表現型の関連を調べることで、実際に、この方法で、稀なメンデル型遺伝疾患の原因遺伝子を同定した、という論文が次々と報告されています。また、遺伝的影響を受けるものの、その形式はより複雑な疾患（生活習慣病、自己免疫疾患、精神神経疾患など）のリスク変異同定にも、現在進行形で活用されつつあります。
さて、解析には当然、それ用のソフトが必要になってきますが、WINDOWS OSでは上手く動かない、動かそうとすると手間がかかるものが多いので、別のOSを使う必要があります。ここでは、LINUX OSを使うことにします。LINUXをインストールする方法としては、すっかりOSを入れ替える、WINDOWSなど他のOSと両方動くようにする、などの方法もありますが、一番お手軽なのは、仮想PCを動かすという方法ではないかと思います（WINDOWS7にあるXpモードもこの方法です）。
仮想PCを動かすソフトウェアもいろいろあるのですが、ここでは、わりと人気、無償、64ビットの仮想マシンが動かせる（基本的に、次世代シーケンサー解析は64ビットマシンでないと厳しいです）、などの理由から、VMware Playerというソフトを使ってみます。VMware Playerのダウンロードはこちら（http://www.vmware.com/jp/products/desktop_virtualization/player/overview）から。登録が必要になります。インストーラのダウンロードが終わったら、ウィザードに従って、インストールを完了させましょう（デフォルトの設定で特に問題ないと思います）。
インストールが終わったら、PCの再起動の後、VMware Playerを起動してみます。下のようなウィンドウが出てくると思うので、「新規仮想マシンの作成」をクリックします。
すると、ゲストOS（ここでは、LINUX OSを入れたい）をどのようにインストールするかを聞かれるのですが、このままではインストールできません。インストーラディスクイメージファイル（インストール用の光ディスクの中身をまとめたファイル、というイメージで良いかと思います）というのを入手する必要があります。
ひとことでLINUX OSといっても、いろいろな種類があるのですが、有名などころとしては、Ubuntu、Fedoraとかいう名前のものがあります。ここでは、Ubuntuベースで、バイオインフォマティクス解析用のツールが最初からある程度入っている、BioLinuxというパッケージのイメージファイル（isoファイル）を使ってみます。ファイルのダウンロードは、http://nebc.nerc.ac.uk/tools/bio-linux/bio-linux-6.0から、Download Nowに進んで、distro.ibiblio.org/bio-linux/iso/bio-linux-6-latest.iso.というところをクリックすれば、特に登録しなくても入手できます。
さて、BioLinuxのisoファイルが手に入ったら、VMwareの新規仮想マシン作成画面に戻って、インストールディスクイメージファイルのところで、先ほど入手したファイルを選択して、次に進みます。ゲストOSはLinux、バージョンはUbuntu64ビットを選択し、続いて格納場所などを指定します。ディスクの容量は、少し余裕を持たせておくといいでしょう。メモリなどは後からでも変更できますが、「ハードウェアをカスタマイズ」のボタンから、この時点で増やしておくこともできます（スムーズに解析を進めるには、8Gぐらいは欲しいところ）。完了をクリックすると、言語、時刻、キーボードレイアウト、ユーザー名、パスワードなどを選択、入力した後、自動でインストールが進みます。
これでBioLinuxのインストールは完了ですが、いろいろとVMware Playerを使いやすくするために、VMware Toolsというのがあったほうが便利です。上部、「仮想マシン」のタブに、「VMware-Toolsのインストール」という項目があり、ここを押すと、フォルダが開くので、まずは中のファイルをデスクトップに移動してみましょう。
続いてインストールにとりかかるわけですが、残念ながら、WINDOWSでよくやるように、インストーラをダブルクリックし、ウィザードに従って「次へ」を順番にクリックすればOK、という訳にはいきません。
インストールは、ターミナル（WINDOWSのコマンドプロンプトのようなもの）を、上部の黒い四角のアイコンをクリックして起動し、コマンドを入力して行う必要があります。この辺が、経験がない人からすると敷居が高いのですが、殆どの次世代シーケンサー解析ツールは、ターミナルからのコマンド入力で動かす必要があるので、少しずつ慣れていきましょう。
ターミナルを起動すると、下のような画面が出てきますが、これは、自分が誰で今どこの場所にいるかを示しています。この場所を「ディレクトリ」と表現します（WINDOWSのフォルダのようなものです）が、ターミナルを使ったコマンド入力では、何も指定しないと今いるディレクトリ（カレントディレクトリ）にあるファイルを扱うことになるので、頻繁にディレクトリを移動する必要があります。そのためのコマンドがcd（Change Directoryの意）です。ただ、目的地にたどり着くためには、自分がどこなら行けるのかが分からないと不便です。そこで、まずはls（LiStの意）と打って、Enterを押してみてください。カレントディレクトリにあるファイルとディレクトリのリストが出てきたかと思います。青い文字はディレクトリを意味しているので、これらの場所なら、直に移動することができます。先ほど、VMware-Toolsのインストール用ファイルはデスクトップに置いたので、
cd Desktop
と入力して、デスクトップに移動します。
ここで、下の方に何か出ている文字を見ると、Tarを使用してインストーラを展開しろ、と書いてあるので、それに従います。
tarというのは、コマンド名でもありますが、ファイルをひとまとめにしたものの拡張子でもあります（とりあえずはzipみたいなものと考えてください）。
いろいろな拡張子のファイルを展開（解凍）するためのコマンドは、こちら（http://uguisu.skr.jp/Windows/tar.html）によくまとまっています。記載されているTAR形式の解凍法に従って、
tar xvf [VMware Toolsのファイル名.tar]
と入力して解凍します（[]の中には、実際のファイル名を入れてください）。tarとxvf、xvfとファイル名の間には、半角スペースを空ける必要があります。小文字大文字が違っても認識されなくなります。ちなみに、xvfと半角スペースの後で、タブを2回押すと、カレントディレクトリにあるファイル一覧が表示され、さらにタブを押すと、順番にファイルが選択されます。また、xvf Vと打ってタブを押すと、Vから始まるファイルの、続きの名前が自動入力されるので、楽することができます。
さて、ファイルが解凍されたところで下の方に出ている文字を再び見ると、vmware-install.plを実行しろ、と書いてあります。
このファイルは、デスクトップではなく、解凍したフォルダの中にあるので、そちらに移動します。
cd vmware-tools-distrib
次にvmware-install.plの実行ですが、そのままvmware-install.plと打っても上手くいきません。前に何も書かずに文字列を入力すると、パスが通っている（http://www.obenri.com/_operation/path_through.html）ディレクトリにある実行ファイルは動かすことができますが（cdとかtarを入力した時が、これに当てはまります）、そうでないコマンドやファイルなどを実行することはできないのです。
vmware-install.plを実行するには、今ここにあるvmware-install.plを実行しろ、とターミナルに命令する必要があります。そのためには、
./vmware-install.pl
と入力します。"."はカレントディレクトリを示す記号で、カレントディレクトリにあるvmware-install.plというファイルを実行しろ、という意味になります。
さて、上記のように入力すると…上手くいきませんね。
何やら、re-run this program as the super userと表示されています。これは、あなたにはこのファイルを実行してインストールする権限がないので、スーパーユーザーとして実行しなさい、ということを意味しています。スーパーユーザーというのは、まぁ名前を見ただけで何か凄そうなのは分かると思いますが、アドミニストレーター権限のようなものです。スーパーユーザーとして実行するためには、
sudo ./vmware-install.pl
と入力して下さい。その後、パスワードを要求されるので、アカウントパスワードを入力すると、無事実行することができます。
また、常にスーパーユーザーとしてコマンドを実行するには、
sudo su
と入力して、パスワードを入れてください。
ユーザー名も変更されたかと思います。
さて、VMware-Toolsのインストールに戻りますが、ファイルを実行後、勝手に進むのに従ってEnterを押していけば大丈夫です。
これで、ゲストOSとホストOS（元から入っているOS）の行き来がスムーズになったり、コピペが共通してできたりするようになります（上手くいかないときは、ターミナルでvmware-userと入力した後、何度かゲストとホストの間でコピペを試してみてください）。
また、重要な機能として、共有フォルダを設定することができるようになります。
これによって、ゲスト、ホストの両方のOSで、スムーズに解析したり、ファイルを移動したりすることが可能になります。
ここでは、必要最低限のことしか記しませんでしたが、UNIX/LINUXの基本的なコマンド、ディレクトリ構造、絶対パスと相対パスなどの用語について、ある程度理解があった方が、後々の解析がスムーズになると思いますので、適宜、調べていただけますと幸いです。
また、だいたいのコマンドは、
[コマンド名] -h
とか、
[コマンド名] --help
で、ヘルプを表示させることができます。
共有フォルダの設定と、仮想マシンの設定は、いずれも「仮想マシン」タブの「仮想マシンの設定」で可能です。繰り返しになりますが、次世代シーケンサー解析には、かなりのマシン性能が必要となります。ハードウェアタブのところで、メモリとプロセッサの設定をしておきましょう（メモリの方は、変更するにはPlayerをいったん終了する必要があります）。メモリは、本体のメモリよりちょっと少なめぐらいに設定するのがいいと思いますが、8Gぐらいは欲しいところです。プロセッサは、4スレッド以上あるのなら、最大の4にしておくのが良いと思います。共有フォルダは、オプションタブの方で、「常に有効」にして、「追加」ボタンから、ウィザードを使って設定します。設定が済んだ共有フォルダは、ゲスト側では、”/mnt/hgfs/[指定した共有フォルダ名]”のディレクトリに存在しています。
さて、LINUX OSのインストールが済んだところで、次は解析ツールのインストールになります。
その前に、エクソーム解析の基本的な流れは、
①Wetな実験と一次配列解析による、ショートリードの配列と各塩基のクオリティを示したファイル（最も一般的な形式はFastqファイル）の生成
②ショートリードの参照配列へのマッピング
③参照配列との比較による変異（一塩基変異=Single-Nucleotide Variant: SNV、欠失/挿入変異=Insertion/Deletion: Indel）のコール
④変異のアノテーション（関連情報の付加）
という感じになるのですが、ここでは、②-④についてとりあげます。
今回インストールするツールのリストは、
・BWA
・GATK
・SAMTools
・Picard
・snpEff
の5種類で、
各ステップにおいて、
②→BWA
③→GATK
④→GATKとsnpEff
を基本的に用います。
また、SAMToolsとGATKを様々なファイルの操作に、Picardを重複リードの除去用います。
Burrows-Wheeler Alignerの略。マッピングの際ギャップを許すので、Indelの同定に向いているようです。解析時間はわりとかかります。多分現時点では最も使われているマッピング用のツール。
インストール方法：
http://sourceforge.net/projects/bio-bwa/files/
から、最新のファイル（2011/10月時点でbwa-0.5.9.tar.bz2）をダウンロード。
ダウンロードしたファイルがあるディレクトリに移動。
http://uguisu.skr.jp/Windows/tar.htmlを参考に、
tar -jxf bwa-0.5.9.tar.bz2
で解凍。
cd bwa-0.5.9
で移動。
ここで、
make
と入力するとインストールが勝手に始まります。
cp bwa /usr/bin/bwa
で、実行ファイルであるbwaというバイナリファイルを、パスが通っているディレクトリにコピー
bwa
と入力して、バージョンや使い方が出てきたらインストール成功です。
The Genome Analysis Toolkitの略。いろんなツールの集合体で、一口には説明できないので、http://www.broadinstitute.org/gsa/wiki/index.php/The_Genome_Analysis_Toolkitをご参照下さい。今回は、変異のコール、アノテーションやマッピングのキャリブレーションなどに使います。
インストール方法：
Broad Instituteの中の人でなければ、
ftp://ftp.broadinstitute.org/pub/gsa/GenomeAnalysisTK/GenomeAnalysisTK-latest.tar.bz2からダウンロード。
ダウンロードしたファイルがあるディレクトリに移動。
tar -jxf GenomeAnalysisTK-latest.tar.bz2
で解凍。
以上でインストールは完了。
このツールは、JAVAという言語で書かれていて、
毎回
java -jar Path/To/GenomeAnalysisTK.jar
で起動することになります。（Path/Toというのは、GenomeAnalysisTK.jarがあるディレクトリへの絶対パスを示しています。参照コマンドを書いてあるサイトなどでは、こういう表記法がされていることが多いのですが、そのままPath/Toと入力しても動きません。）
SAMToolsのSAMはSequence Alignment/Mapの意味。SAMフォーマットは、ショートリードを参照配列にマッピングした結果のファイル形式として、最も一般的なものです。SAMをデータ圧縮のためにバイナリ化した形式がBAMフォーマットで、SAM→BAMの変換、BAMのソートによる容量の削減、ソートしたBAMのインデックス作成などもSAMToolsで行います。また、変異のコールもできますが、最近はGATKのGenotyperでやった方が精度が高い、とのことで、そちらでやることが増えています。
SAMToolsはもとからBioLinuxには入っていますが、いちおう最新版に置き換えておきましょう。
http://sourceforge.net/projects/samtools/files/
から、最新のファイル（2011/10月時点でsamtools-0.1.18.tar.bz2）をダウンロード。
ダウンロードしたファイルがあるディレクトリに移動。
tar -jxf samtools-0.1.18.tar.bz2
で解凍。
cd samtools-0.1.18
で移動。
make
でインストール。
cp samtools /usr/bin/samtools
でパスが通っているディレクトリにコピー。
samtools
で、バージョンが0.1.18になっていればOKです。
ぴかーると読みます。なぜPicardかはよく知りません。SAM、BAMファイルを操作するためのJAVA言語で書かれたツール群で、ここでは重複リード（Pair endの両リードとも配列が完全に一致しているもの。特定のPCR産物が偏って増幅されることなどによって生じ、カバレッジのばらつきや、変異コールのエラーにつながる。）を除去する目的で使います。
http://sourceforge.net/projects/picard/files
から、最新のファイル（2011/10月時点でpicard-tools-1.53.zip）をダウンロード。
ダウンロードしたファイルがあるディレクトリに移動。
unzip picard-tools-1.53.zip
で解凍（Lhaplusとかでもいいが）。
これも、
java -jar Path/To/XXX.jar
で起動するので、以上でインストール完了。
SNPのEffectをアノテーションとしてつけるためのツール。イントロンかエクソンか、アミノ酸置換を起こすか起こさないか、などの情報が付加できます。
http://snpeff.sourceforge.net/download.htmlから、
コアプログラム（http://sourceforge.net/projects/snpeff/files/snpEff_v2_0_3_core.zip/download）をダウンロードして解凍。
ヒトのデータベース（2011/10月時点で最新はGRCh37.64 http://sourceforge.net/projects/snpeff/files/databases/v2_0_3/snpEff_v2_0_3_GRCh37.64.zip）をダウンロードして解凍。
dataというディレクトリができるので、snpEff_v2_0_3ディレクトリの中に格納します。
GATKなど同様、
java -jar Path/To/snpEff.jar
で動かします。
ツールが揃ったところで、さぁ解析、と行かないところがもどかしいですが、その前に、様々な参照ファイル（ヒトゲノムの参照配列、変異のリストなど）が必要になってくるので、それらをひと通り入手しましょう。
このファイルを入手できるサイトは複数ありますが、Broad InstituteのFTPサイトを例としてあげておきます。
パスワードを要求されますが、空欄のままでOKです。
こちらから、bundle、1.2、b37と進んだところにある、
human_g1k_v37.fasta.gz
が1000ゲノムプロジェクトで使われたヒトの参照配列
human_g1k_v37.fasta.fai.gz
が、そのインデックスファイルです。
fastaというのが、一般的なゲノム参照配列の書式になっています。
また、このファイルでは、各染色体の番号を
1, 2, 3…という風に表記していますが、
chr1, chr2, chr3…という書式になっている参照ファイルもあります。
この表記法、染色体の順番ともに揃っていないと、解析が上手くいかないところがでてくるので、注意が必要です。
（http://www.broadinstitute.org/gsa/wiki/index.php/Contig_orderingをご参考下さい）
上記サイトには、古いバージョンのヒトゲノム参照配列などもあります。
変異のリストもいろいろなバージョンがあります。
最新のものは、dbSNP（http://www.ncbi.nlm.nih.gov/projects/SNP/）のFTPサイトからダウンロードできます（ここのところ、1000ゲノムのデータが出てきたことで頻回に更新されています）。ファイルの詳細については、各ディレクトリのReadmeファイルをご参照下さい。
また、ファイルのフォーマットとしては、今はvcf形式（Eメールクライアントなどで使われている電子名刺の形式と同じ名前ですが、別のものです）が一般的になってきています。
少し古いもので良ければ、上述のBroad InstituteのFTPサイトにも、
dbsnp_132.b37.vcf.gz
dbsnp_132.b37.vcf.idx.gz
という変異リストのファイルと、そのインデックスファイルがあります。
今回はこのファイルを使ってみることにします。
また、短いIndelがある部分は、ショートリードを上手くマッピングすることが難しくなりやすいのですが、その問題点を解決する方法として、既知のIndelの周辺で、再アラインメントを行う、という方法があります。そのために、Indelだけのリストが必要になります。
これも、Broad InstituteのFTPサイトから
1000G_biallelic.indels.b37.vcf.gz
1000G_biallelic.indels.b37.vcf.idx.gz
をダウンロードすれば、1000ゲノムプロジェクトで見つかったIndelのリストが入手できます。
また、dbSNPや1000ゲノムのvcfファイルには、たいていINFOのフィールドにVC (Variant classification) の項目があり、こちらにSNPかINDELかが書かれています。
grepというUNIXコマンド（特定の文字列がある行、ない行だけを抜き出すことなどができる）を使った下のようなコマンド
zcat [vcfファイルの名前.vcf.gz] | grep -e ^# -e INDEL > [Indelだけのファイルの名前.vcf]
で、ダウンロードしたgzファイルから、INDELのある行だけを抽出したファイルを作成することもできます。（"|"はパイプ、">"は出力のリダイレクトと呼ばれるもので、入出力の受け渡しに便利な表記法です。）
エクソーム解析の場合、ターゲット領域が各キットごとに決まっています。
データのクオリティコントロール（平均カバレッジ、5X以上で読めている領域の割合など）のために、このファイルも入手しておきましょう。
ファイルは、エクソームキャプチャー用のキットの、各メーカーサイトから入手できると思いますが、
Illumina TruSeqは https://icom.illumina.com/download/summary/Ue6ccqdHiEmCuHC3cBLpxQ
Agilent SureSelectは https://earray.chem.agilent.com/earray/ のeArrayのサイトから、"exon"のキーワードで探してもらえれば、ダウンロードできます（いずれも登録が必要です）。
eArrayのサイトからダウンロードできるSureSelectのファイルの場合は、chr1, chr2, chr3…という書式になっているので、上でダウンロードしたFastaファイルなどと合わせる必要があります。
ここでは、sedというコマンドを使って修正してみます。
sed s/chr// [エクソームターゲット領域のファイル.bed] > [修正後のエクソームターゲット領域のファイル.bed]
sedは、sで置換機能を使うことを指定し、"sed s/対象文字列/置換文字列/ ファイル名"という形で入力すると、各行の最初に出てくる対象文字列（ここではchr）が、置換文字列（ここでは、何も文字がない）で置き換えられます。
また、最初のヘッダ行も削除しておく必要がありますが、これは、普通にテキストエディタを使ってもらえばいいと思います。
この項は、http://www.vlsci.org.au/lscc/exome-pipelineの解析パイプラインがよくまとまっていたので、そちらを参考に作成させていただきました。
いよいよ解析にとりかかりますが、その前に、解析をスムーズに進めるためのディレクトリ構造を作成し、また、何度も使うディレクトリ、ファイルの場所を事前に指定しておきます。
ディレクトリ構造はあくまでも一例ですが、下の図のような形にしておくと、操作がしやすいかと思います。
共有フォルダ下にdata、exon_outputのディレクトリをつくり、dataの中にnextGen_exome、その中に、各ツールの解凍後のディレクトリ、各参照ファイルのディレクトリが入ったReference_files、fastqファイルのディレクトリが入ったexomeのディレクトリを、それぞれ移動、作成しています。exon_outputの中には、各プロジェクトのディレクトリを作成しています。
続いて、何度も使うディレクトリ、ファイルの場所を、下のように入力して事前に指定しておきます。（各ツールのバージョンが違う場合は、それぞれ修正してください。）
fastq_name=sequence.txt
i=6
#Fastqファイルの指定。Illuminaの次世代シーケンサーのアウトプットファイルは、デフォルト名がs_[レーン名]_[Pair endの場合、"1"か"2"]_sequence.txtになります。i=[レーン名]（ここでは、"6"にしてある）と指定すれば、ファイル名を変更しなくていいようにしています。
fastq_dir=/mnt/hgfs/share/data/nextGen_exome/exome/fastq/
#Fastqファイルがあるディレクトリを指定。
outdir=/mnt/hgfs/share/exon_outputs/s_6
#ファイルを出力するディレクトリを指定
reference=/mnt/hgfs/share/data/nextGen_exome/Reference_files/g1k_v37/human_g1k_v37.fasta
#参照配列のファイルを指定。
dbSNP_file=/mnt/hgfs/share/data/nextGen_exome/Reference_files/snp/dbsnp_132.b37.vcf
#変異リストのファイルを指定。
indel_ref_file=/mnt/hgfs/share/data/nextGen_exome/Reference_files/indel/1000G_biallelic.indels.b37.vcf
#Indelのファイルを指定。
exon_capture_file=/mnt/hgfs/share/data/nextGen_exome/Reference_files/exome_target/SureSelect_All_Exon_G3362_with_annotation.hg19.bed
#エクソームターゲットのファイルを指定。
bwa_dir=/mnt/hgfs/share/data/nextGen_exome/bwa-0.5.9/
#BWAのバイナリファイルがあるディレクトリを指定。
samtools_dir=/mnt/hgfs/share/data/nextGen_exome/samtools-0.1.18
#samtoolsのバイナリファイルがあるディレクトリを指定。
gatk_dir=/mnt/hgfs/share/data/nextGen_exome/GenomeAnalysisTK-1.2-38-g0016c70/
#GATKのバイナリファイルがあるディレクトリを指定。
picard_dir=/mnt/hgfs/share/data/nextGen_exome/picard-tools-1.53/
#Picardのバイナリファイルがあるディレクトリを指定。
snpEff_dir=/mnt/hgfs/share/data/nextGen_exome/snpEff_v2_0_3/
#snpEffのバイナリファイルがあるディレクトリを指定。
#$の後に、先ほど入力した左辺の名前を入れると、右辺のディレクトリを指定することができます。
#下のコマンドで、outdirの下に、途中で出力されるファイルを格納するためのディレクトリを作成し、outdir/intermediateに移動します。
mkdir $outdir/intermediate
mkdir $outdir/intermediate/tmp
mkdir $outdir/intermediate/coverage/
cd $outdir/intermediate
#BWAを使って、マッピングを行います。
#参照FastaファイルのIndexを、indexコマンドを使って作成します。
$bwa_dir/bwa index -p $reference -a bwtsw $reference
#参照ファイルごとに、一度だけIndexを行う必要があります。-pで接頭辞（Prefix）を指定。参照Fastaファイルと同じでOK。-aでアルゴリズムを指定。アルゴリズムはbwtswとisがありますが、ヒト全ゲノムのような2GB以上のファイルの場合はbwtsw。オプションの後に、半角スペースを空けて参照Fastaファイルを指定。
#それぞれのFastqファイルのショートリードを、alnコマンドで参照配列にアラインメントします。
$bwa_dir/bwa aln -t 4 $reference $fastq_dir"s_"$i"_1_"$fastq_name > "s_"$i"_1_.sai"
$bwa_dir/bwa aln -t 4 $reference $fastq_dir"s_"$i"_2_"$fastq_name > "s_"$i"_2_.sai"
#-tは、処理に使うスレッド数。オプションの後、PrexixとFastqファイルを指定。Fastqファイルを出力したIlluminaパイプラインソフトウェアのバージョンが1.3より新しい時は、-Iのオプションをつけます。次のステップで入力するファイルは拡張子が.saiでないと通らないので、そのように出力します。
#sampeで、Pair EndのリードからSAMファイルを作成します（Single endの場合はsamse）。
$bwa_dir/bwa sampe -P $reference "s_"$i"_1_.sai" "s_"$i"_2_.sai" $fastq_dir"s_"$i"_1_"$fastq_name $fastq_dir"s_"$i"_2_"$fastq_name -r "@RG\tID:01\tSM:s6\tPL:Illumina" > "s_"$i".sam"
#-Pで先ほどalnしたsaiファイルを指定し、sai1 sai2 fastq1 fastq2の順番で並べます。-rは、Read groupのヘッダーをつけるためのオプションで、SM（Sample）とPL（Platform）がないとGATKが動きません。ID、SMは適当でOK。PLは、今のところ454, LS454, Illumina, Solid, ABI_Solid, CGのいずれかにする必要があるようです。出力は.samで。
#SamtoolsのviewコマンドでSAMを、容量が少ないBAMに変換。
$samtools_dir/samtools view -bS "s_"$i".sam" > "s_"$i".bam"
#-bで出力がBAMであること、-Sで入力がSAMであることを指定。
rm "s_"$i".sam"
#いらなくなったSAMファイルをrm (remove) します。
#Samtoolsのsortコマンドを使って、BAMファイルをソート。これをすると、さらに容量が減ります。
$samtools_dir/samtools sort "s_"$i".bam" "s_"$i"_sorted"
#ソートが終わったら、BAMにインデックスを付ける必要があります。
$samtools_dir/samtools index "s_"$i"_sorted.bam"
#PicardのMarkDuplicatesで、重複リードを除去します。
java -Xmx2G -jar $picard_dir/MarkDuplicates.jar ASSUME_SORTED=true REMOVE_DUPLICATES=true INPUT="s_"$i"_sorted.bam" OUTPUT="s_"$i"_sorted_rem.bam" METRICS_FILE=duplicate VALIDATION_STRINGENCY=LENIENT
#重複リード除去後のBAMにインデックス作成。
$samtools_dir/samtools index "s_"$i"_sorted_rem.bam"
#PicardのCreateSequenceDictionaryで、Intarvalファイルの作成に必要な、参照配列のdictファイルを作成。
java -Xmx2G -jar $picard_dir/CreateSequenceDictionary.jar REFERENCE=$reference OUTPUT=/mnt/hgfs/share/data/nextGen_exome/Reference_files/g1k_v37/human_g1k_v37.dict
#GATKのRealignerTargetCreatorを使って、s_"$i"_sorted.bamから再アラインメントした方がよさそうなIntervalのファイルを作成。
java -jar $gatk_dir/GenomeAnalysisTK.jar -T RealignerTargetCreator -R $reference -I "s_"$i"_sorted_rem.bam" --known $indel_ref_file -log intervals.log -L $exon_capture_file -o "s_"$i".intervals"
#GATKのIndelRealignerで、先ほど調べたIntervalの再アラインメントを施行。
java -Djava.io.tmpdir=$outdir/intermediate/tmp -jar $gatk_dir/GenomeAnalysisTK.jar -T IndelRealigner -R $reference -I "s_"$i"_sorted_rem.bam" -log realigned.log -targetIntervals "s_"$i".intervals" -o "s_"$i"_realigned.bam"
#PicardのFixMateInformationで、Pair endの情報を修正。
java -Xmx2G -jar $picard_dir/FixMateInformation.jar INPUT="s_"$i"_realigned.bam" SO=coordinate VALIDATION_STRINGENCY=SILENT TMP_DIR=$outdir/intermediate/tmp
#再アラインメント後のBAMファイルにインデックス作成。
$samtools_dir/samtools index "s_"$i"_realigned.bam"
#IndelRealignerの時にできていたインデックスファイルを削除（これをしないと、次のステップがうまくいかない）。
rm "s_"$i"_realigned.bai"
#gatkのCountCovariatesとTableRecalibrationを使って、再アラインメントしたファイルの塩基クオリティスコアを再計算。
java -jar $gatk_dir/GenomeAnalysisTK.jar -I "s_"$i"_realigned.bam" -R $reference --knownSites $dbSNP_file -nt 4 -T CountCovariates -l INFO -log recalc.log -recalFile recal_data.csv -cov ReadGroupCovariate -cov QualityScoreCovariate -cov CycleCovariate -cov DinucCovariate
java -jar $gatk_dir/GenomeAnalysisTK.jar -I "s_"$i"_realigned.bam" -R $reference -T TableRecalibration -l INFO -recalFile recal_data.csv -log new_qual.log --out "s_"$i"_realigned_alnRecal.bam"
#再アラインメント、再計算後のBAMファイルにインデックス作成。
$samtools_dir/samtools index "s_"$i"_realigned_alnRecal.bam"
#GATKのDepthOfCoverageでカバレッジの計算。
java -jar $gatk_dir/GenomeAnalysisTK.jar -T DepthOfCoverage -R $reference -I "s_"$i"_realigned_alnRecal.bam" -o $outdir/intermediate/coverage/"s_"$i"_depth" -L $exon_capture_file -omitBaseOutput -ct 5 -ct 10 -ct 20 -ct 30
#-ctで、5X以上カバレッジの割合、10Xカバレッジ以上の割合などを表示させることができます。-geneListで、エクソンごと、遺伝子ごとなどのカバレッジを調べることもできます。
#BedtoolsのcoverageBedを使う方法は割愛。
#GATKのUnifiedGenotyperを使って、SNVのコールをします。
java -jar $gatk_dir/GenomeAnalysisTK.jar -T UnifiedGenotyper -R $reference -I "s_"$i"_realigned_alnRecal.bam" -nt 4 -L $exon_capture_file -D $dbSNP_file -A AlleleBalance -A DepthOfCoverage -o $outdir/"s_"$i"_SNV.vcf" -stand_call_conf 50.0 -stand_emit_conf 10.0 -dcov 200 -out_mode EMIT_VARIANTS_ONLY
#-Lで領域を指定してやると、そこにあるSNVだけ出力できます。-stand_call_confと-stand_emit_confで、コールする変異、ファイルに吐き出す変異の最小Phred Scoreを指定。-dcovは、これ以上のカバレッジがあるときに、入力した数のリードをランダムにサンプリングして、処理にかかる時間を短縮します。-out_modeで、変異の部分だけ出力するか、全ての塩基について出力するか、などが設定できます。
#続いて、GATKのVariantFiltrationを使って、SNVのフィルタリングを行います。
#まずoutdirに移動。
cd $outdir
java -jar $gatk_dir/GenomeAnalysisTK.jar -T VariantFiltration -R $reference -V "s_"$i"_SNV.vcf" -o "s_"$i"_SNV_filtered.vcf" --clusterWindowSize 10 --filterExpression "QUAL < 30.0 || QD < 5.0 || HRun > 5 || SB > -0.10" --filterName "HARD_TO_VALIDATE"
#--clusterWindowSizeは、SNVがクラスターとなって見つかる（=本当にたくさんSNVがあるのかもしれないけど、False positiveの可能性も高い）場所を調べるためのウィンドウサイズを決めます（ここでは10塩基づつウィンドウをずらして調べる）。--filterExpressionで、vcfファイルのINFOフィールドにある項目のうち、フィルタリングに使うものを指定。複数指定するときは、||で区切ります。--filterNameで、filterExpressionごとに、フィルターをパスしなかったSNVにつけるタグの名前を指定できます。
#まずoutdir/intermediateに戻ります。
cd $outdir/intermediate
#GATKのUnifiedGenotyperを使って、Indelのコールをします。
java -Xmx2G -jar $gatk_dir/GenomeAnalysisTK.jar -T UnifiedGenotyper -glm INDEL -R $reference -I "s_"$i"_realigned_alnRecal.bam" -nt 4 -L $exon_capture_file -D $dbSNP_file -A AlleleBalance -A DepthOfCoverage -o $outdir/"s_"$i"_indel.vcf" -stand_call_conf 50.0 -stand_emit_conf 10.0 -dcov 200
#-glm（genotype_likelihoods_model）で、INDELコールであることを指定します。
#続いて、Indelのフィルタリングを行います。
#まずoutdirに移動。
cd $outdir
#GATKのVariantFiltrationを使って、フィルタリングを行います。
java -jar $gatk_dir/GenomeAnalysisTK.jar -T VariantFiltration -R $reference -V $"s_"$i"_indel.vcf" -o $outdir/"s_"$i"_indel_filtered.vcf" --filterExpression "MQ0 >= 4 && ((MQ0 / (1.0 * DP)) > 0.1)" --filterName "HARD_TO_VALIDATE" --filterExpression "SB >= -1.0" --filterName "StrandBiasFilter" --filterExpression "QUAL < 10" --filterName "QualFilter"
#SNVのときと同じようにパラメーター設定。詳細は、ヘルプなどもご参照ください。
#CombineVariantで、SNVのファイルとIndelのファイルをまとめます。
java -Xmx2g -jar $gatk_dir/GenomeAnalysisTK.jar -R $reference -T CombineVariants -V:SNV,vcf $outdir/"s_"$i"_SNV_filtered.vcf" -V:INDEL,vcf $outdir/"s_"$i"_indel_filtered.vcf" -o "s_"$i"_SNV_indel_filtered.vcf"
#snpEffを使って、変異のアノテーションを行います。
#カレントディレクトリにsnpEffのconfigファイルが存在している必要があるので移動
cd $snpEff_dir
#snpEffによるアノテーションを実行
java -Xmx4G -jar ./snpEff.jar eff -v -i vcf -o vcf GRCh37.63 $outdir/"s_"$i"_SNV_indel_filtered.vcf" > $outdir/"s_"$i"_SNV_indel_filtered.snpEff.vcf"
#以上で、アノテーションつきの、SNVとIndelをコールしたvcfファイルが作成できます。
ということで、一通り方法を解説させていただきましたが、実際に生データを扱ってみないと実感がわかない、というのが正直なところだと思います。
ですので、サンプルデータを用いた解析を行なってみましょう。
生データを公開しているサイトとしては、NCBIのSRAが有名ですが、こちらでは、SRA形式のファイルしか配布されておらず、各自でFastqファイルに変換する必要があります。
変換は、SRA Toolkitを使えば可能ですが、面倒な場合は、DDBJ（DNA Data Bank of Japan）のDRAからなら、Fastqファイルを直にダウンロードできます。
今回は、劣性遺伝のメンデル型遺伝病の原因遺伝子を、両親が近親婚の症例から得たDNAのエクソーム解析によって同定することを目的とした研究の、公開されている生データをとりあげてみます。
研究の論文はこちら
http://www.ncbi.nlm.nih.gov/pubmed/21936905
生データはこちら
http://trace.ddbj.nig.ac.jp/DRASearch/study?acc=SRP007801
から、右側のNavigation部分のExperiment欄、Fastqのところから、ファイルがダウンロードできるFTPサイトに移動できます。
材料が揃ったところで、上述のコマンドをコピペで入力していくと、こんな感じのファイルができると思います。
ここまでダウンサイズできると、あとはExcelなどでの解析でも十分対応できるのですが、ここでは、annovarというツールを使って、簡単に疾患の原因遺伝子候補を絞り込む方法を紹介したいと思います。
#例のように、最新バージョンのものをダウンロードして、解凍します。
#annovarでの解析を始める前に、まず、上の解析で出来たvcfファイルから、PASSのフラグがついている変異だけを抽出します。
grep -e ^# -e PASS $outdir/"s_"$i"_SNV_indel_filtered.snpEff.vcf" > $outdir/"s_"$i"_SNV_indel_filtered.snpEff.PASS.vcf"
#annovarのディレクトリに移動します。
cd annovar
#./convert2annovar.plを使って、vcf形式から、annovarで解析可能な形式に変換します。
./convert2annovar.pl $outdir/"s_"$i"_SNV_indel_filtered.snpEff.PASS.vcf" -format vcf4 > $outdir/"s_"$i"_SNV_indel_filtered.snpEff.PASS.txt"
#ここでは、自動で劣性遺伝のメンデル型遺伝病と仮定した場合の絞り込みを行なってくれる、./auto_annovar.plを使ってみます。
#フィルタリングのステップの詳細は、こちらをご参照下さい。また、annovarは、これ以外にも、いろいろとカスタマイズ可能なフィルタリング法を持っていますので、適宜、試していただければと思います。
#まず、解析に必要なデータベースファイルをダウンロードします。
./annotate_variation.pl -downdb -buildver hg19 gene humandb
./annotate_variation.pl -downdb -buildver hg19 knownGene humandb
./annotate_variation.pl -downdb -buildver hg19 segdup humandb
./annotate_variation.pl -downdb -buildver hg19 phastConsElements46way humandb
./annotate_variation.pl -downdb -buildver hg19 -webfrom 1000g2010nov humandb
./annotate_variation.pl -downdb -buildver hg19 -webfrom annovar snp130 humandb
#./auto_annovar.plを使って、変異の絞り込みを行います。
./auto_annovar.pl --buildver hg19 -model recessive $outdir/"s_"$i"_SNV_indel_filtered.snpEff.PASS.txt humandb
#最終的に生成してきたファイルを見ると、20個の遺伝子名が残っており、その中には、原因遺伝子と考えられたSHROOM3も残っているかと思います。
という感じで、ひと通り、エクソームデータ解析の流れを紹介してみました。
ただ、今回とりあげたデータ解析法で、必ず上手くいくとは限らないのが現実です。
例えば、サンプルデータでとりあげた研究では、SNPチップを用いたabsence-of-heterozygosity region（ホモの変異が連続して続く領域、近親婚の両親から共通して遺伝したと考えられる）の解析を組み合わせて、更なる絞り込みを行なっているのですが、こういった方法も有効です。
また、原因変異が同定できない理由としては、
・調べた疾患が、浸透率100%のメンデル型遺伝でない。
・表現型の定義が間違っている。
・原因変異がコールされていない（原因変異がExome外にある場合を含む）。
・LowQualやHARD_TO_VARIDATEに真の原因変異が潜んでいる。
・dbSNPに登録されいるが、疾患の原因になる変異を誤って消している。
・種間の保存の程度や予測プログラムによる判定では重要でなさそうだが、疾患の原因になる変異を誤って消している。
などが考えられると思います。
さらに、dbSNPについてのTipsとして、
・dbSNPには、既知のメンデル型遺伝病の原因変異も登録されている。
・dbSNPは、130から1000ゲノムのデータが入ってきている。
・1000ゲノム以降のデータには、劣性で病気の原因となる、頻度の低い変異が入っている可能性がある。
・1000ゲノム以降のデータには、FalsePositiveも少なくないかもしれない。
といったことも考慮する必要があると考えられます。
dbSNP135では、いずれかの人種で頻度が1%以上、などの基準を満たした変異のみを含むファイルが作成されるようなので、稀な疾患の原因遺伝子を探す場合、今後はこのファイルを使うのが良いかもしれません。
というわけで、まだまだ一筋縄ではいきませんし、複雑な遺伝様式をとるありふれた疾患の研究では、さらに解析法を工夫する必要がありますが、エクソーム解析が、遺伝的影響を受ける疾患の研究において、非常に有望な手段であることは間違いないところでしょう。
最後に、エクソーム解析を用いた疾患研究についての、良いReviewを紹介しておきますので、ご参照下さい。
Exome sequencing as a tool for Mendelian disease gene discovery
（Mendelian disease gene discoveryとありますが、今後のcommon diseaseへの応用についても述べてあります。）
Written by 高田篤（mail: atakata@brain.riken.jp）