NGSデータ解析まとめ

サカナ研究者の手探りNGS解析(おもに進化生物学)

GenomeやRNA-Seq assemblyの結果を確認する: assembly-statsについて

de novo assembleで得られたゲノムやRNA-Seqのデータについての統計値(contig長の平均・中央値・最大値、またN50-N100の値など)を確認したいときに便利なツールとして、assembly-statsを紹介します。

以下のサイトからソースコードをダウンロードしてコンパイルします。

github.com

コンパイルの方法はReadme.mdおよびサイトの指示に従って下さい。cmakeのインストールは必要かもしれません。

使用方法は

assembly-stats Contig.fasta > Contig.stats

のような感じで、outputはこんな感じになります。

stats for Namazu_OE.fasta
sum = 366946480, n = 324791, ave = 1129.79, largest = 25875
N50 = 2488, n = 45468
N60 = 1933, n = 62163
N70 = 1358, n = 84637
N80 = 773, n = 120174
N90 = 382, n = 189190
N100 = 201, n = 324791
N_count = 0
Gaps = 0

非モデル生物ゲノムのde novo assembly(その2):DBG2OLCを使用する(2): 解析の流れ

(2016-03-04 暫定版です)

前回の続きです。Illumina NGS (HiSeq, MiSeqなど)のデータとPacBioのデータの両方を使ってゲノムのde novo assembleを行う (hybrid assembly)方法について書きます。

ここではDGB2OLCのManualに示されているサンプルデータを用いて解析を行います。以下は解析の手順(ほぼマニュアルの翻訳)になります。

メモ:Step 1(Illumina readsのassemble)とStep 2(PacBioとIllumina contigのhybrid assemble)では、blasrとSparcは不要です。blasrがインストールできなくても、hybrid assembleのcontig(backbone_raw.fasta)を得ることはできます。Step 3(consensus配列の作成)にはblasrが必要になります。

サンプルデータのダウンロード(FTPサイト)
Illumina reads (fastq): ftp://qb.cshl.edu/schatz/ectools/w303/Illumina_500bp_2x300_R1.fastq
PacBio reads (fasta): ftp://qb.cshl.edu/schatz/ectools/w303/Pacbio.fasta.gz

Step 0: 下準備

まず、解析を行う適当なディレクトリを用意して、そのディレクトリに移動します。サンプルデータ(Illumina, PacBio)はそのディレクトリに置きます。
次に、実際に解析を行うディレクトリを1つ作成します。解析の結果、たくさんのファイルが出力されるので、適当なディレクトリを作っておかないと混乱します。

# 解析を行うディレクトリ // 日付 + sample
mkdir 160304_sample
# カレントディレクトリを移動
cd 160304_sample

Step 1: Illumina readsをassembleする

最初に、Illumina readsをassembleしてcontigを作成します。Assembleはどのようなソフトで行っても良いようですが(ようするに、contigのfastaファイルがあれば良い)、DBG2OLCが用意しているサンプルデータでは、SparseAssemblerを使用しています。

以下のコマンドを実行します。詳しくはManualを参照。

# SparseAssemblerによるIllumina readsのde novo assemble
SparseAssembler LD 0 k 51 g 15 NodeCovTh 1 EdgeCovTh 0 GS 12000000 f ../Illumina_500bp_2x300_R1.fastq

解析がうまくいけば、いくつかのファイルができます。Contigのファイルは"Contigs.txt"になります。

Step 2: DBG2OLCを実行する

PacBio readsとIllumina readsから作られたcontigsを組み合わせたhybrid assemblyを行います。DBG2OLCを使用します。いくつかの変数(パラメーター)を設定する必要があります。

以下のようなコマンドで実行します。

# step2 // 
DBG2OLC k 17 AdaptiveTh 0.0001 KmerCovTh 2 MinOverlap 20 RemoveChimera 1 Contigs Contigs.txt f ../Pacbio.fasta

サンプルの例ではPacBioのデータが1つのファイルしかないですが、もし複数のsubreadsのファイルがある、また複数cell分のFASTA fileが存在する場合、catなどを使って1つのファイルに合体して使えば問題ないようです。

# たとえば、3つのFASTA fileを1つに合体する
cat subread1.fasta subread2.fasta subread3.fasta > PacBio.fasta

以下、変数の説明です。

# 3つのメジャーな変数(あんまりわかってません。。すいません)

AdaptiveTh adaptive k-mer matchingのしきい値。もし"M < AdaptiveTh × Contig_Length"ならば、 そのcontigはlong readのanchorとしては使えない。(Mはcontig (Illumina)とlong read (PacBio)の間でmatchするk-mers)
KmerCovTh fixed k-mer matchingのしきい値。あるcontigについて、もし"M < KmerCovTh"ならば、そのcontigはlong readのanchorとしては使われない。
MinOverlap Long readsのpair間で(許容する)最小のoverlap scoreの値。

Manualの[Miscellaneous]のところには、これらのパラメーター設定について、以下のような記述があります。

[Miscellaneous]
At this point, the parameters may be fine-tuned to get better performance. As with SparseAssembler, LD 1 can be used to load the compressed reads/anchored reads.
Suggested tuning range is provided here:
For 10x/20x PacBio data: KmerCovTh 2-5, MinOverlap 10-30, AdaptiveTh 0.001~0.01.
For 50x-100x PacBio data: KmerCovTh 2-10, MinOverlap 50-150, AdaptiveTh 0.01-0.02.

(和訳)この点について、これらの(3つの)パラメーターを最適化することで、より良いパフォーマンスが得られる可能性がある。Compressed read/anchored readsを読み込むため、SparseAssemblerではLD 1にする設定が利用できる(<-- よくわかりません。すいません)
10x/10x PacBio data: KmerCovTh 2-5, MinOverlap 10-30, AdaptiveTh 0.001~0.01にする
50x-100x PacBio data: KmerCovTh 2-10, MinOverlap 50-150, AdaptiveTh 0.01-0.02にする

また、実際に変数をどのように設定するかは「いくつかのパターンで実行してみて、より良い結果を出すものを探すべし(意訳)」とのことです。

その他の、あまりflexibleでない(重要性は相対的に低い)パラメーターについて、

k-mer size 17が良い
Contigs Illumina readsのassmbleにより得られたcontigのファイル (FASTA)
MinLen 最小のread length(しきい値?)
RemoveChimera (たぶんPCR artefactで生じた)キメラのリードを除く。データのcoverageが10x以上の場合、1に設定する

Step 3: ゲノムの「コンセンサス配列」を作成する

ここでblasrとSparcが必要になります。blasrおよびSparcのインストールについては前のエントリを参照して下さい。

(1) Sparcのプログラムおよびいくつかのスクリプトを、解析を行うディレクトリにコピーする

コピーでもシンボリックリンクでもいいですが、Sparcフォルダの中にある以下のものを解析を行うディレクトリに置きます。

compiled/Sparc
utility/RunSparcConsensus.txt
utility/SeqIO.py
utility/split_and_run_sparc.sh
utility/split_reads_by_backbone.py

(2) ctg_pb.fastaを作成する

"ctg_pb.fasta"はIllumina contigのファイルとPacBioのファイルを合体したファイルで、次の作業のinput fileになります。

cat Contig.txt ../PacBio.fasta > ctg_pb.fasta

(3) シェルスクリプトsplit_and_run_sparc.shを実行する

以下のコマンドを実行して、解析開始です。Consensus配列の結果は、"results"というフォルダ内にできます。

sh split_and_run_sparc.sh backbone_raw.fasta DBG2OLC_Consensus_info.txt ctg_pb.fasta results 2>cns_log.txt

ここでのinput fileは以下の4つです

(1) backbone_raw.fasta (by DBG2OLC)
(2) DBG2OLC_Consensus_info.txt (by DBG2OLC)
(3) DBG contigs (in fasta format)
(4) PacBio reads (in fasta format)

3と4を合体したのが"ctg_pb.fasta"になります。

非モデル生物ゲノムのde novo assembly(その2):DBG2OLCを使用する(1): 下準備とインストール

ここでは、Illumina NGS (HiSeq, MiSeqなど)のデータとPacBioのデータの両方を使ってゲノムのde novo assembleを行う (hybrid assembly)方法について書きます。

具体的には、DBG2OLCを使った方法について書いていきます。私の使用している環境は、以下になります。

OS: Ubuntu 12.04.4 LTS
CPU: 2 x Xeon E5-2620v2 2.10GHz
RAM: 256 GB DDR3-1600 (16GB x 16)

1. DBG2OLCのダウンロード
以下のサイトからダウンロードします。
https://sites.google.com/site/dbg2olc/

コンパイルは必要ありません(パッケージに実行可能ファイルが付属しています。PATHを通しておくと良いです。

2. 解析のための下準備:blasrのインストール準備
DBG2OLCで解析を行うためには、blasrをインストールする必要がありますが、これが結構大変です。以下は、blasrのインストール手順と、引っかかったところのメモです。

以下の作業、HDF5およびblasrの展開・コンパイルなどを行ったPATHは、/home/user/Tools/になります。

(1) HDF5のインストール

blasrをインストールするためには、先にHDF5というのをインストールしておく必要があります。

# 参考URL
www.hdfgroup.org

ソースコードをURLからダウンロードしてコンパイルします。

# installation manual
https://www.hdfgroup.org/ftp/HDF5/current/src/unpacked/release_docs/INSTALL

# ./configure --prefix=のPATHについて
prefixのディレクトリは適当な場所にする (/usr/local/hdf5など、あるいはlocalなどこか)
今回、私のLinuxでは/home/user/Toolsに設定しました <--注意:あまりよくないのでいずれ変更すること
HDF5をインストールすると、include, libなどのディレクトリがここで指定したPATH (e.g. /home/user/Tools)に作られます。

(2) updating gcc and g++ (2.6.3 -> 2.8)

私のLinuxに入っていたg++(2.6.3)ではblasrがコンパイルできなかったので(2016.3.2現在)、gcc g++をupdateする必要がありました。
# 参考URL
金星☆ちゃんねる: C++11のためにGCCの最新版をインストールする

# commands
sudo add-apt-repository ppa:ubuntu-toolchain-r/test
sudo apt-get update
sudo apt-get install gcc-4.8 g++-4.8
sudo update-alternatives --install /usr/bin/gcc gcc /usr/bin/gcc-4.8 20
sudo update-alternatives --install /usr/bin/g++ g++ /usr/bin/g++-4.8 20
sudo update-alternatives --install /usr/bin/gcc gcc /usr/bin/gcc-4.6 10
sudo update-alternatives --install /usr/bin/g++ g++ /usr/bin/g++-4.6 10
# 元に戻す方法
sudo update-alternatives --config gcc
sudo update-alternatives --config g++

3. blasrのインストール

注意:以下については、blasrのパッケージに付属している"README.INSTALL.md"ファイルや、blasr/README.md at master · mchaisso/blasr · GitHubも参照して下さい。

(1) blasrのダウンロード

以下のサイトからダウンロード、展開します。
github.com

(2) get libcpp

以下を実行して、libcppの内容をupdateします。

git pull --rebase origin master && git submodule update --init
make update-submodule

うまくいかない場合:以下URLからlibcpp folderの内容を入手して、既存のlibcppと置き換えます(ちなみに、私はうまくいきませんでした)。
github.com

(3) configure.py実行

configure.pyにHDF5_INCLUDE, HDF5_LIBのPATHを与えて実行します。以下実行例です。ただし、HDF5_INCLUDE, HDF5_LIBのPATHの指定がうまくいかないことがあります。その場合の対処法については後述します。<--やり方がそもそも間違えている可能性もあり。あとで調べること 

# configure.pyを実行 // 160302
./configure.py --shared --sub --no-pbbam HDF5_INCLUDE=/home/user/Tools/include HDF5_LIB=/home/user/Tools/lib

(4) make configure-submodule実行

以下のコマンドを実行します。

make configure-submodule

(5) Build // README.INSTALL.mdも参照

以下、コマンドを実行します。

# libcpp libraryをmake
make build-submodule
# make
make

# エラーが出て止まらなければ動作確認
./blasr

動かない場合は、エラーメッセージを見て適宜対応します。
以下のようなエラーが出る場合は、LD_LIBRARY_PATHが通っていないことが原因です。

./blasr
blasr: error while loading shared libraries: libpbihdf.so: cannot open shared object file: No such file or directory

このような場合は、以下を実行します。

## export LD_LIBRARY_PATHの追加
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/home/user/Tools/blasr/libcpp/hdf:/home/user/Tools/blasr/libcpp/alignment:/home/user/Tools/blasr/libcpp/pbdata:/home/user/Tools/lib

このexportの命令は、.profileとかの設定ファイルに書いておくと良いのかもしれません。<--ちょっとよくわからない(要勉強!)

以上、コンパイルがうまくいけば

./blasr    # blasr実行
blasr - a program to map reads to a genome
 usage: blasr reads genome 
 Run with -h for a list of commands 
          -help for verbose discussion of how to run blasr.

In release v3.0.1 of BLASR, command-line options will use the 
single dash/double dash convention: 
Character options are preceded by a single dash. (Example: -v) 
Word options are preceded by a double dash. (Example: --verbose) 
Please modify your scripts accordingly when BLASR v3.0.1 is released. 

のように表示されます。

(6) おまけ
最後に、blasrのPATHを通しておくと良いです(方法はいろいろあると思います)。

# /usr/local/bin/にコピー
sudo cp blasr /usr/local/bin

4. DBC2OLCの実行に必要なプログラムのダウンロード

それぞれ、以下のサイトからダウンロードして展開します。いずれもコンパイルは必要ありません(パッケージにすでに付属しています)。

(1) Sparcのダウンロード
github.com
PATHを通す必要はないです。

(2) SperseAssemblerのダウンロード
github.com
こちらは、PATHを通しておくと良いです。

次回は、DBG2OLCの具体的な使い方について説明します。

非モデル生物ゲノムのde novo assembly(その1): はじめに、いくつかの方法

2016年現在、次世代シーケンサーを使って、非モデル生物のdraft genomeを個人レベルで決定することもだいぶ現実的になってきました。ここでは、最近試みているNGSで得られたゲノム配列データのde novo assembleについて考えていきます。
(現在、いろいろ試しているので、ちょっとずつ更新する予定です // 2016-03-02)

非モデル生物の全ゲノムを決定する方法として、ここでは以下の2つを取り上げます。

1. Illumina HiSeqなどによるペアエンドリード & HiSeqによるメイトペアリード

Platanus assemblerを使用する方法
以下のエントリを参照(未作成 160302)

2. Illumina HiSeqリード & PacBioリードによるhybrid assembly

(1) PlatanusによるIllumina short readsのde novo assemble, scaffolding & FGAPによるgap closing
以下のエントリを参照(未作成 160302)

(2) DBG2OLCを使用する方法
インストール方法
使用方法
PacBioによるde novo assembly およびhybrid assemblyについての参考サイト
github.com

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で新規書類を開き、
[データ] --> [外部データの取り込み] --> [テキストファイルのインポート]
を選択します(下図)。

f:id:hashiyuki:20160125142746j:plain

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

f:id:hashiyuki:20160125142933j:plain

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

f:id:hashiyuki:20160125143236j:plain

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