生命情報科学

VCF

ファイルフォーマット

VCF (variant call format) はゲノムの塩基の変異データの保存する際に利用するファイル形式である。SNPs 解析などによく用いられる。シーケンシングデータをリファレンス配列にマッピングしたとき、リファレンス配列上の塩基とそこにマッピングされたシーケンシングデータ上の塩基などの情報が記載される。

次には A. thaliana の 1001 genomes でダウンロードした VCF の例である。この例にあるにように、VCF ファイルの先頭には ## から始まる部分があり、アノテーション情報などが記載されている。それから、# から始まるヘッダー行に続いて、実際のデータが記載される。

##fileformat=VCFv4.1
##FILTER=<ID=q25,Description="Quality below 25"g>
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
##FILTER=<ID=q30,Description="Quality below 30">
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	10001
1	135	.	C	T	28	PASS	DP=5	GT:GQ:DP	1|1:28:5
1	138	.	T	C	36	PASS	DP=5	GT:GQ:DP	1|1:36:5
1	346	.	C	T	32	PASS	DP=3	GT:GQ:DP	1|1:32:3
1	508	.	T	C	27	PASS	DP=4	GT:GQ:DP	1|1:27:4
1	657	.	C	T	32	PASS	DP=3	GT:GQ:DP	1|1:32:3
1	698	.	G	A	32	PASS	DP=6	GT:GQ:DP	1|1:32:6
1	730	.	G	A	28	PASS	DP=3	GT:GQ:DP	1|1:28:3
1	833	.	C	T	36	PASS	DP=4	GT:GQ:DP	1|1:36:4
1	892	.	T	G	36	PASS	DP=4	GT:GQ:DP	1|1:36:4
1	953	.	T	C	32	PASS	DP=5	GT:GQ:DP	1|1:32:5
1	1260	.	T	A	30	PASS	DP=3	GT:GQ:DP	1|1:30:3
1	2309	.	T	A	36	PASS	DP=4	GT:GQ:DP	1|1:36:4

VCF ヘッダー行

ヘッダー行はタブ区切で、最初の 8 列が固定されている。9 列目には 10 列目以降に記載されているデータのフォーマットを書かれ、10 列目以降には実際データが書かれている。サンプルが複数ある場合は、11 列目、12 列目と列を増やして対応する。

  1. CHROM: 染色体番号。
  2. POS: 塩基の位置。
  3. ID: 塩基（変異）の ID。該当変異が dbSNP に登録されていれば、rs から始まるアクセッション番号を記載する。
  4. REF: リファレンス配列上の塩基。
  5. ALT: シーケンシングデータ上の塩基。
  6. QUAL: シーケンシングデータのスコア。
  7. FILTER: VCF を作成する際に利用したフィルタリング条件をすべて満たす場合は PASS と記載される。そうでない場合は、満たされない条件がセミコロン区切りで記載される。例えば q10;s50 ならば、クオリティが 10 未満かつ 50% 以下のサンプルがこの変異を持つ、という意味になる。
  8. INFO: key:data の形式で追加情報を記載する。AA, AC, AF, AN, BQ, CIGAR, DB, DP, END, H2, MQ, MQ0, NS, SB, SOMATIC, VALIDATED がすでに定義されているが、これ以外のキーを独自に定義できる。

塩基の変異情報

VCF ファイルの REF 列と ALT 列を見比べることによって、変異なのか、挿入あるいは欠損なのかを知ることができる。例えば、次の例にある 793657 の塩基はリファレンス配列では C であるのに対して、シーケンシングデータ上では A である。そのため、793657 の塩基は C から A の変異であるとわかる。850756 ではリファレンス上に 2 塩基あるのに対してシーケンシングデータ上では 1 塩基となっているため、これは A が欠損したとわかる。854093 も同様に考えて、T が挿入されているとわかる。

#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	10001
1	793657	.	C	A	40	PASS	DP=13	GT:GQ:DP	1|1:40:13
1	850756	.	TA	T	34	PASS	DP=4	GT:GQ:DP	1|1:34:4
1	854093	.	C	CT	40	PASS	DP=12	GT:GQ:DP	1|1:40:12

遺伝子型の情報

遺伝子型の情報は VCF ファイルの第 10 列目以降に記載する。次の例では、FORMAT 列に GT:GO:DP と書かれ、10001 列では 1|1:40:13 のように記載されている。これは GT = 1|1, GO = 40, DP = 13 を表している。GT で記載されているデータが遺伝子型である。

#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	10001
2	448314	.	C	T	38	PASS	DP=6	GT:GQ:DP	1|1:38:6
2	448314	.	G	C	27	PASS	DP=2	GT:GQ:DP	0|0:27:2
2	793657	.	C	A	40	PASS	DP=13	GT:GQ:DP	1|0:40:13
2	3237192	.	C	T	40	PASS	DP=8	GT:GQ:DP	0|1:40:8
2	3899408	.	G	C	32	PASS	DP=8	GT:GQ:DP	.|.:32:8

GT に記載されるデータは二倍体であれば「|」または「/」を真ん中に起き、左右に 0, 1, 2, … などの数値を記載する。上の例では 448314 位置にある塩基は T をホモに持つ。448314 位置にある塩基は G をホモに持つ。これに対して、793657 位置にある塩基は C と A をヘテロに持つ。データが観測されない場合は . と記載する。

VCFtools

VCFtools は VCF ファイルを扱うためのツール一群である。

複数の VCF ファイルを一つの VCF ファイルにマージする例。

tabix -p vcf intersection_9998.vcf.gz
tabix -p vcf intersection_9999.vcf.gz

vcf-merge intersection_9998.vcf.gz intersection_9999.vcf.gz | bgzip -c > out.vcf.gz

gzip -dc out.vcf.gz | head -n 20
## ##fileformat=VCFv4.2
## ##FILTER=<ID=q25,Description="Quality below 25">
## ##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
## ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
## ##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
## ##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
## ##FILTER=<ID=q30,Description="Quality below 30">
## ##source_20160702.1=vcf-merge(v0.1.14-12-gcdb80b8) intersection_9998.vcf.gz intersection_9999.vcf.gz
## ##sourceFiles_20160702.1=0:intersection_9998.vcf.gz,1:intersection_9999.vcf.gz
## ##INFO=<ID=SF,Number=.,Type=String,Description="Source File (index to sourceFiles, f when filtered)">
## ##INFO=<ID=AC,Number=.,Type=Integer,Description="Allele count in genotypes">
## ##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
## #CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	9998	9999
## 1	110	.	G	T	34.00	PASS	AC=2;AN=2;DP=4;SF=1	GT:GQ:DP	.	1|1:34:4
## 1	125	.	G	T	32.00	PASS	AC=2;AN=2;DP=9;SF=0	GT:GQ:DP	1|1:32:9	.
## 1	214	.	G	A	40.00	PASS	AC=2;AN=2;DP=7;SF=1	GT:GQ:DP	.	1|1:40:7
## 1	502	.	T	C	40.00	PASS	AC=2;AN=2;DP=12;SF=1	GT:DP:GQ	.	1|1:12:40
## 1	508	.	T	C	36.00	PASS	AC=4;AN=4;DP=15;SF=0,1	GT:GQ:DP	1|1:32:6	1|1:40:9
## 1	645	.	C	A	40.00	PASS	AC=2;AN=2;DP=12;SF=1	GT:GQ:DP	.	1|1:40:12
## 1	698	.	G	A	38.00	PASS	AC=4;AN=4;DP=14;SF=0,1	GT:DP:GQ	1|1:7:36	1|1:7:40

PyVCF

PyVCF は VCF ファイルをパースするための Python ライブラリーである。PyVCF を利用することで VCF ファイルから必要なデータを簡単に取得できる。

インストール

PyVCF は pip コマンドでインストールできる。

pip install pyvcf --user

使い方

VCF ファイルから各サンプルの genotype (GT) 値を取得してタブ区切りのファイルに保存する例。

import os
import sys
import argparse
import vcf

'''
Usage:
    python getGT.py -i sample.vcf -o output.txt
'''

def format_str(x):
    if x is None:
        y = ''
    else:
        if isinstance(x, list):
            if not x:
                y = ';'.join(x)
            else:
                y = ''
        else:
            y = str(x)
    return y

def parse_vcf(i, o):
    with open(o, 'w') as outfh:
        vcf_fh = vcf.Reader(open(i, 'r'))
        for snp_record in vcf_fh:
            gt_list = []
            for sample in snp_record.samples:
                gt_list.append(sample['GT'])
            
            new_record = [format_str(snp_record.CHROM),
                          format_str(snp_record.POS),
                          format_str(snp_record.ID),
                          format_str(snp_record.REF),
                          format_str(snp_record.ALT)]
            new_record.extend(gt_list)
            outfh.write('\t'.join(new_record) + '\n')


if __name__ == '__main__':
    parser = argparse.ArgumentParser(description = 'Get Sample information.')
    parser.add_argument('-i', '--input', required = True)
    parser.add_argument('-o', '--output', required = True)
    args = parser.parse_args()
    parse_vcf(args.input, args.output)