VCF fileの操作：grepを使用する方法

以前、VCFファイルから必要なデータを抽出する方法として、Rを使う方法を少し書きました。複雑な操作をする時はRが便利と思いますが、比較的簡単な操作（たとえば特定の情報を持つSNPsの行を抽出する場合など）を行う場合は、UNIXのコマンドである"grep"を利用することで比較的簡単にデータ抽出を行うことができます。grepは、指定したファイル内を検索して、マッチする文字列を含む行を返すコマンドです（詳しくはこちらなど参照）。

ここでは、grepを用いたデータ抽出の方法を紹介します。

(1) GNU grep (ggrep)のインストール

解析を行う前に、Mac OS Xの場合は、より高速に動作するGNU grepをインストールします（OS XにはBSD grepがすでに入っていますが、こちらよりGNU grepのほうが速いらしい）。

Homebrewを使用します。

brew update  #brewをupdateする
brew tap homebrew/dupes
brew install homebrew/dupes/grep

無事インストールされているかどうかを確認するため、以下のコマンドを入力します。

which ggrep

うまくインストールされていると"/usr/local/bin/ggrep"など、ggrepのPATHが表示されます。

また、versionの確認は以下のコマンドで。

ggrep --version

GNU grepのダウンロードに際して、参考にしたページは以下です。

qiita.com

(2) サンプルVCF fileのダウンロード

ここでは、Variant Effect PredictorなどでアノテーションされたVCFファイルのサンプルとして以下のサイトから"MHC_test4.var.annotated.vcf"というファイルをダウンロードして下さい。次に、適当な場所にフォルダを作成し、その中にサンプルVCFファイル（MHC_test4.var.annotated.vcf）を保存して、そのフォルダに移動して下さい。

(3) grepによるデータの抽出その1: 特定の遺伝子上にあるSNPs

例えば、HLA-AのSNPsを検索する場合、ggrepで以下のように実行します。

ggrep -w HLA-A MHC_test4.var.annotated.vcf > HLA-A.vcf

結果のファイルHLA-A.vcfの中身を見ると、以下のようになっています。

less HLA-A.vcf    # HLA-A.vcfの中身を表示

6|29941244|.|C|G|999|.|DP=138;VDB=0.0000;AF1=0.3334;AC1=2;DP4=58,11,48,20;MQ=49;FQ=403;PV4=0.068,0.36,0.37,0.38;CSQ=G|upstream_gene_variant|MODIFIER|HLA-A|ENSG00000206503|Transcript|ENST00000396634|protein_coding||||||||||rs2571420|16|1|HGNC|HGNC:4931||||G:0.0790|G:0.08|G:0.05|G:0.02||G:0.14|||||||||||GT:PL:GQ|1/1:255,205,0:99|0/0:0,33,191:33|0/0:0,175,255:99
6|29941256|.|A|G|31.9|.|DP=150;VDB=0.0001;AF1=0.3229;AC1=2;DP4=101,36,9,1;MQ=50;FQ=33.5;PV4=0.45,4.9e-34,1,1;CSQ=G|upstream_gene_variant|MODIFIER|HLA-A|ENSG00000206503|Transcript|ENST00000396634|protein_coding||||||||||rs9260078|4|1|HGNC|HGNC:4931||||A:0.2897|G:0.79|G:0.76|G:0.74||G:0.62|||||||||||GT:PL:GQ|0/0:0,220,255:99|1/1:76,18,0:12|0/0:0,187,255:99
6|29941267|.|G|C|999|.|DP=159;VDB=0.0003;AF1=0.33;AC1=2;DP4=105,41,11,1;MQ=50;FQ=999;PV4=0.18,1,1,0.49;CSQ=C|5_prime_UTR_variant|MODIFIER|HLA-A|ENSG00000206503|Transcript|ENST00000396634|protein_coding|1/10||||8|||||rs28749141||1|HGNC|HGNC:4931||||C:0.3168|C:0.30|C:0.33|C:0.28||C:0.35|||||||||||GT:PL:GQ|0/0:0,247,255:99|1/1:208,23,0:17|0/0:0,187,255:99
6|29941293|.|T|C|999|.|DP=182;VDB=0.0019;AF1=0.3333;AC1=2;DP4=118,44,14,5;MQ=51;FQ=999;PV4=1,1,1,0.5;CSQ=C|5_prime_UTR_variant|MODIFIER|HLA-A|ENSG00000206503|Transcript|ENST00000396634|protein_coding|1/10||||34|||||rs9260079||1|HGNC|HGNC:4931||||||||||||||||||||GT:PL:GQ|0/0:0,255,255:99|1/1:255,57,0:51|0/0:0,211,255:99
6|29941294|.|G|A|999|.|DP=180;VDB=0.0027;AF1=0.3333;AC1=2;DP4=71,16,61,31;MQ=51;FQ=999;PV4=0.027,7.8e-66,0.34,0.01;CSQ=A|5_prime_UTR_variant|MODIFIER|HLA-A|ENSG00000206503|Transcript|ENST00000396634|protein_coding|1/10||||35|||||rs9260080||1|HGNC|HGNC:4931||||||||||||||||||||GT:PL:GQ|1/1:255,255,0:99|0/0:0,57,255:57|0/0:0,161,255:99
6|29941376|.|G|T|999|.|DP=195;VDB=0.0063;AF1=0.3333;AC1=2;DP4=106,67,14,6;MQ=54;FQ=999;PV4=0.63,0.068,0.46,1;CSQ=T|intron_variant|MODIFIER|HLA-A|ENSG00000206503|Transcript|ENST00000396634|protein_coding||1/9||||||||rs9260081||1|HGNC|HGNC:4931||||G:0.3232|T:0.69|T:0.74|T:0.72||T:0.61|||||||||||GT:PL:GQ|0/0:0,255,255:99|1/1:255,60,0:54|0/0:0,208,255:99
6|29941404|.|A|C|999|.|DP=208;VDB=0.0134;AF1=0.6667;AC1=4;DP4=46,23,75,59;MQ=54;FQ=212;PV4=0.17,1.9e-100,1e-06,0.48;CSQ=C|intron_variant|MODIFIER|HLA-A|ENSG00000206503|Transcript|ENST00000396634|protein_coding||1/9||||||||rs1632889||1|HGNC|HGNC:4931||||A:0.1579|C:0.78|C:0.82|C:0.95||C:0.80|||||||||||GT:PL:GQ|1/1:200,255,0:99|1/1:182,72,0:72|0/0:0,147,255:99
6|29941420|.|G|C|999|.|DP=220;VDB=0.0173;AF1=0.6667;AC1=4;DP4=48,27,73,63;MQ=57;FQ=212;PV4=0.19,3.5e-81,3e-06,1;CSQ=C|intron_variant|MODIFIER|HLA-A|ENSG00000206503|Transcript|ENST00000396634|protein_coding||1/9||||||||rs1655919||1|HGNC|HGNC:4931||||G:0.0955|C:0.94|C:0.89|C:0.98||C:0.83|||||||||||GT:PL:GQ|1/1:255,255,0:99|1/1:243,72,0:72|0/0:0,199,255:99
6|29941428|.|C|T|999|.|DP=208;VDB=0.0181;AF1=0.3333;AC1=2;DP4=63,33,58,50;MQ=57;FQ=999;PV4=0.089,5.1e-68,0.0018,1;CSQ=T|intron_variant|MODIFIER|HLA-A|ENSG00000206503|Transcript|ENST00000396634|protein_coding||1/9||||||||rs1632888||1|HGNC|HGNC:4931||||T:0.1414|T:0.07|T:0.08|T:0.20||T:0.17|||||||||||GT:PL:GQ|1/1:204,255,0:99|0/0:0,87,255:87|0/0:0,179,255:99
6|29941440|.|T|G|84.2|.|DP=204;VDB=0.0177;AF1=0.3333;AC1=2;DP4=57,36,53,53;MQ=57;FQ=86;PV4=0.12,3.8e-94,0.0029,0.45;CSQ=G|intron_variant|MODIFIER|HLA-A|ENSG00000206503|Transcript|ENST00000396634|protein_coding||1/9||||||||rs1632887||1|HGNC|HGNC:4931||||||||||||||||||||GT:PL:GQ|1/1:129,255,0:99|0/0:0,99,255:99|0/0:0,161,255:99
6|29941442|.|T|C|95.2|.|DP=211;VDB=0.0171;AF1=0.3333;AC1=2;DP4=56,37,32,36;MQ=58;FQ=97;PV4=0.11,1.5e-99,0.14,0.11;CSQ=C|intron_variant|MODIFIER|HLA-A|ENSG00000206503|Transcript|ENST00000396634|protein_coding||1/9||||||||rs1632886||1|HGNC|HGNC:4931||||||||||||||||||||GT:PL:GQ|1/1:140,202,0:99|0/0:0,99,255:99|0/0:0,164,255:99

このように、"HLA-A"を含むすべてのSNPsの行が抽出されています。

(4) grepによるデータの抽出その2: 特定の性質を持つSNPs

特定の性質を持つSNPsを、consequencesなどの属性をキーにして抽出することができます（SNPごとの属性についてはこのエントリを参照）。例えば、ナンセンス変異 (stop_gained)になるSNPをすべて抽出する場合、ggrepで以下のように実行します。

ggrep stop_gained MHC_test4.var.annotated.vcf > stop_gained.vcf

あるいは、ミスセンス変異のSNPsの場合、

ggrep missense MHC_test4.var.annotated.vcf > missense.vcf

(5) grepによるデータの抽出その3: 有害な (deleterious) 変異

SIFTやPolyPhenなどにより予測された、塩基置換による個体への影響をキーにして、特定のSNPsを抽出することもできます（SIFTやPolyPhenについては、以下のエントリを参照）。

たとえば、有害 (deleterious)な変異を起こすSNPsを抽出するには、以下のようにします。

ggrep deleterious MHC_test4.var.annotated.vcf > deleterious.vcf

(6) grepによるデータの抽出その4: 複数の遺伝子上にあるSNPsを同時に抽出する

たとえば、特定の疾患に対して複数の関連遺伝子が存在していて、それらの遺伝子上にあるSNPsを同時に抽出したいような場合、以下のような方法が使えます。

まず、最初に関連遺伝子をリストアップしたファイルを作ります。たとえば、

HLA-A
HLA-B
HLA-DRB5
HLA-DQB1

これを、たとえばファイル名"candidates.txt"として保存します。

そして、以下のように実行します。

ggrep -wFf candidates.txt MHC_test4.var.annotated.vcf > candidates.vcf

4つの遺伝子上に存在するSNPsのリストになっています（確認して下さい）。

さらに、ここから絞り込むこともできます。たとえば、genotypeがref/alt allelesのheterozygoteになっているものに絞る場合は、以下のようにします。

ggrep 0\/1 candidates.vcf > candidates.hetero.vcf

(7) おまけ：抽出したSNPsのリストをExcelで開く

抽出したSNPsのリストは"|"で区切られたファイルになっています。そのままでは見辛いので、Excelに読み込んで開くと良いです。

Excelで新規書類を開き、
[データ] --> [外部データの取り込み] --> [テキストファイルのインポート]
を選択します（下図）。

選択対象を「すべてのファイル」にし、読み込むVCFファイルを選択して「インポート」をクリックします（下図）。

テキストファイルウィザードの最初の画面は、そのまま「次へ」をクリックします。2つ目の画面で、「区切り文字」のところで「その他」をチェックし、区切り文字入力場所に"|"を入力（下図）し、「次へ」をクリックします。

あとはそのまま進めば完了です。

IGVの使い方(1): ヒトゲノムannotation情報の取り込み

IGVにRefSeqなど、ヒトゲノム上の遺伝子アノテーションを取り込む方法についてです。

IGVでは、"File" メニューの"Load from Server"から、ゲノム上の位置に関連付けられたアノテーション情報を読み込むことができます。

(1) "File" -> "Load from Server"の中身について

このメニューを選択すると、以下のような画面が出ます。

このうち、Available Datasetの一番上にある"Annotations"の中に、いくつかの呼び出し可能なデータがあります。"Annotations"の左側の三角をクリックすると、下図のようにデータの各項目が出てきます。

さらに各項目の左側の三角をクリックすると、各データ内の各項目が出てきます（下図）。たとえば"Genes"の項目の中にも、"Ensembl Genes"や"UCSC Genes"など、異なる由来のアノテーションが含まれます。

(2) "Annotations"の各項目について

• Genes: 遺伝子のアノテーション。EnsemblやUCSCなど、複数の異なる機関によるアノテーションが存在する。通常はどれか一つでOK。
• Phenotype and Disease Associations: 表現型および病気との関連（直訳）、"OMIM Genes"は遺伝子、"GWAS catalog"はGenome wide association studiesの情報で、SNPsです。
• Sequence and Regulation: CpG islandの位置及び領域ごとのGC%を表示します。"GC percentage"がわりと重いです。要注意。
• Variation and Repeats: "dbSNP"はSNPsのデータベースで、"Repeat Masker"はそのまま、RepeatMaskerで予測した繰り返し配列の位置を表示します。
• Comparative Genomics: これについては別のエントリで詳細を見てみます（予定）

(3) 一般的な注意

"File" -> "Load from Server"でデータを取り込む際に、データ項目のtreeの上流の方でチェックを入れてしまうと（たとえば、"Annotations"そのものにチェックを入れる）、膨大なデータを取り込もうとするので、IGV自体が止まってしまいます。ただ、これはちょっとわかりにくい仕様のような気がします。IGVのdocumentsを見ると、ちゃんと以下のような注意書きがあるのですが。

Warning: Selecting a folder selects all of its subfolders and all of the datasets in those folders. This can potentially be a huge amount of data. To be sure you are loading only the data you want, open collapsed folders and select only the datasets of interest.
和訳- 注意：一つのフォルダを選択すると、そのフォルダ内のすべてのサブフォルダ及びそれらのフォルダ内のすべてのデータセットを選択してしまいます。この場合、潜在的に膨大な量のデータになることがあります。すべてのフォルダを開いて、必要なデータセットだけを選択して、必要なデータのみを読み込むよう、気を付けて下さい。

IGV (Integrative Genomic Viewer)を使ってみる

次に、IGVで実際にリファレンスゲノムにreadsがマップされた結果を見てみます。

まず、IGVを開きます。デスクトップ上のIGVのエイリアスをダブルクリックします。このとき、IGVでは"Loading genome"という表示が一瞬出ますが、このときリファレンスゲノムとしてhg19を読み込んでいるようです。

IGVが起動したら、自分のサンプルを読み込みます。

(1) BAM fileを読み込む

"File"メニューから"Load from File..."を選択し、開きたいBAM file（indexあり）を選択します。
Mappingされた状態の見方ですが、まず画面左上の2番目の選択する場所（図1）で、"All"から"chr1"などに選択を変更します。

図1. IGVの染色体番号を選択する部分

次に、IGV画面右上のZoomを選択する部分（図2）で、倍率を適当に上げていくと（+側に動かす）、BAMファイルに記録されている個々のreadの位置が表示されます。

図2. IGVの右上（Zoom）

ただし、最初は見ている染色体上の位置が、何もないところ（セントロメア上）になっているので、赤で示されているカーソルを適当な位置に動かします（図3）。

図3. 染色体の位置とカーソル（赤い四角）

(2) VCF file（SNP等のアノテーション）を読み込む

"File"メニューから"Load from File..."を選択し、開きたいVCF fileを選択します。VCFファイルは圧縮された状態（.gz）をそのまま開けます（解凍する必要なし）。

うまく開けると図4のような状態になります。

図4. IGV window (BAM fileとVCF fileを読み込んだ状態）

"Sample"の行に示されているのがreferenceのgenotypeの状態で、水色はreferenceと異なるalleleのhomozygote, 青はheterozygoteになっています。色の部分にカーソルを合わせると、各SNPの情報を見ることができます（おそらくVCF fileに記載されているもの）。

IGVの具体的な使い方については次回以降のエントリで説明します。

IGV (Integrative Genomic Viewer)のインストール（Macの場合）

Broad Instituteが提供しているIGV (Integrative Genomic Viewer)のインストールについての覚え書きです。

まず最初に、以下のIGVのPageを開きます。

Home | Integrative Genomics Viewer

"Download"のところをクリックします。最初は、ユーザ登録を求められるので、登録のページに移動してメールアドレスを入力します。メールアドレスを送ると、ダウンロードのページに移動できます。

IGVのインストール方法は3つあります。
(1) Mac applicationとしてIGVをインストールする
(2) "Use the Java Web Start buttons"を使用する
(3) バイナリをダウンロードして、コマンドラインからIGVを立ち上げる

このうち(1)ですが、このAppを使用するためには、MacにJava 7がインストールされている必要があります。Java 6 (JRE1.6)を使用している場合は、(2)の方法でインストールする必要があります。自分のMacでJavaのversionを確認するためには、ターミナルを起動して

java -verison

と入力します。

たとえば私のMac (OS10.8.5) <--ちょっと古い　では、こうでした。Version 6です。

java version "1.6.0_65"
Java(TM) SE Runtime Environment (build 1.6.0_65-b14-462-11M4609)
Java HotSpot(TM) 64-Bit Server VM (build 20.65-b04-462, mixed mode)

そこで、(2)の方法で行きます（もし最初からJava 7がインストールされている場合は、(1)でダウンロードしたアプリをダブルクリックするだけでIGVが開くと思います）。

まず、Dowanloadページから"Java Web Start"をインストールします。うまくいけば、ダウンロードされたファイル"igv_hm.jnlp"（注）をダブルクリックするとIGVのインストールが始まり、起動します。他のソフト（テキストエディタ等）で開いてしまう場合、「右クリック」-> 「このアプリケーションで開く」-> 「Java Web Start.app」とします。

うまく行かない場合、Java Runtime Environmentをupdateする必要があり、この場合はOSの指示に従って最新のJREをインストールします。

インストールが完了すると、IGVのエイリアスがデスクトップに生成されます。

(3) の方法については試していません（Linuxなどではこちらを使う？）

注：igv_hm.jnlpがダウンロードされる位置

/private/var/folders/nz/cgn91zkn2vg1lp8fbnvmbgvr0000gn/T/igv_hm.jnlp

BAM fileからunmapped readsを抽出する（samtoolsを使用する）

BAM file (.bam)から、reference sequencesにmapされなかったreadsを抽出する方法について。以下、ちょっとした覚え書き。

BAM fileはbwaなどでreference sequence(s)にNGS readsをマッピングした結果のoutput file であるSAM file (.sam)を圧縮したファイル。基本的には、samtools view コマンドで作成します。

samtools view -bS example.sam > example.bam

マッピング結果から多型サイトを抽出したり、contigをexportしたり、などの作業は、BAM fileをもとに行うことが多いです。また、BAM indexを作成して、samtools tviewやTabletなどでmapping結果を見る際にも使用します。

samtools sort example.bam example.sort
samtools index example.sort.bam
samtools view example.sort.bam reference.fasta

以下は、BAM fileから"mappingされなかった" readsをFASTQ formatで抽出する方法について。

具体的には、novoalignの解説サイト

Extracting Unmapped Reads from a BAM File Produced by NovoAlign | Novocraft

を参考にしています。ていうか、以下はこのサイトの（ほぼ）日本語訳です。

• • • • • -

以下ここから

Novoalignで作られたBAM fileからunmapped readsを抽出する (extracting)方法について

Introduction

Referenceにmapされなかった配列をBAM fileから抽出することがが必要となるケースはしばしば生じる。Samtoolsを使うことで、この作業は容易に行うことができる。
（中略）

この作業は以下の2つのステップで行われる。

1. Unmapped readsを、readnameでsortされたBAM fileから抽出する
2. BAM fileをFASTQ read fileに変換する

FASTQ fileの生成をskipして、(unmapped) readsを直接、BAM fileから再度アラインメントすることも可能である。BAMをFASTQに変換するために、Hydra packageの"bamToFastq"を使用する。
（中略）

Novoalignによるアラインメント

# Illumina paired-end: mean=500, sd=100
novoalign -d database.nix -f reads1.fastq -reads2.fastq -oSAM  -i 500 100 | samtools view -bS - alignments

Unmapped readsに対するフィルタリング

次に、以下3つのカテゴリーでreadsのフィルタリングを行う。（注：ここではIllumina NGSによるpaired end readsを想定している）

1. Paired endの一方はマッピングされているが、それ自体はマッピングされていない
2. それ自体はマッピングされているが、paired endのもう一方はマッピングされていない
3. Paired endの両方ともマッピングされていない

1-3の各カテゴリは、以下のsamtools filtering commandで抽出することができる。

samtools view -u  -f 4 -F264 alignments.bam  > tmps1.bam
samtools view -u -f 8 -F 260 alignments.bam  > tmps2.bam
samtools view -u -f 12 -F 256 alignments.bam > tmps3.bam

上のコマンドは、sortされたBAM fileでも、sortされていないBAM fileでも可能である。

次に、これら各カテゴリのBAM fileを結合 (merge)し、readsの名前順にsortする。ここでのsortはpaired endの相手を正しく得るために必要である。

samtools merge -u - tmps[123].bam | samtools sort -n - unmapped

ReadsをBAM fileから直接再アラインメントする

novoalign  -F BAMPE -f unmapped.bam -d …

または

ReadsをFASTQ formatで抽出する

次に、unmapped read pairsをそれぞれのFASTQ filesに抽出する（注：paired end readsを想定）。ここではHydra-SV packageのbamToFastqプログラムを利用する。

bamToFastq -bam unmapped.bam -fq1 unmapped_reads1.fastq -fq2 unmapped_reads2.fastq

Single endの場合はもっと簡単に、以下の1行コマンドで抽出できる。

samtools view -uf 4 alignments.bam | bamToFastq -bam stdin -fq1 unmapped.fastq -fq2  empty.fastq

ここまで

• • • • • -

わからないのは、3つのカテゴリをsamtoolsで抽出するところの

samtools view -u  -f 4 -F 264 alignments.bam  > tmps1.bam
samtools view -u -f 8 -F 260 alignments.bam  > tmps2.bam
samtools view -u -f 12 -F 256 alignments.bam > tmps3.bam

の部分ですが、samtoolsにおけるoption -fや-Fを理解するためには、SAM fileの形式を理解する必要がありそうです。SAM fileがどのように記述されているか、については、日本語で以下のサイトが参考になりました。

http://crusade1096.web.fc2.com/sam.html

SAM fileの2列目は、FLAGという形式で、readの状況を表わしているようです（2進数の10進法表示！？）。いずれにしても、samtoolsを使いこなすには、まずSAM fileの意味を理解する必要がありそうです・・・