さて、データにアノテーションを付けていく作業はvcfでやるのが効率的なのだけど、これを視覚的にわかりやすく提示するにはちょっとごちゃごちゃしすぎているので、これをtext抽出してSQLに入れてやることにする。
INFOやsampleカラムには複数の情報が詰まっていて、それも行毎にすべて揃っているわけではないので単純なawkによる抽出はできない（できなくもないが結構面倒くさい）のでpythonのモジュールを使うことにする。
pyvcfというまんまなモジュールがあったのでそれを使うことにする。
https://pypi.python.org/pypi/PyVCF
先人の情報があんまりないのだけれど、基本的には csvモジュールなんかと同じようにファイルをリードでオープンしてforで行を順番に取り出す、みたいな感じで処理できるようなのでなんとかなるでしょう。
https://bi.biopapyrus.jp/python/module/pyvcf.html
こちらに出ているスクリプトを拝借して処理をしてみると
1 28499 A 0/1 1 608335 T 0/1 1 874461 GAA 1/1 1 895287 GTCA 0/1 1 911789 TA 1/1 1 1671823 A 0/1 1 1765789 A 0/1 1 2152975 A 0/1 1 2243386 T 1/1 1 2498920 T 0/1
こんな感じに出力される。
CHROM, POS, REF, ALT, GTが出力されるはずなんだけど、どうもALTが出力されてないような。
>>> import vcf
>>> vcf_reader = vcf.Reader(open('vcf/test/example-4.0.vcf', 'r'))
>>> for record in vcf_reader:
... print record
Record(CHROM=1, POS=28499, REF=A, ALT=[C])
Record(CHROM=1, POS=608335, REF=T, ALT=[G])
Record(CHROM=1, POS=874461, REF=GAA, ALT=[GAAA])
Record(CHROM=1, POS=895287, REF=GTCA, ALT=[G])
Record(CHROM=1, POS=911789, REF=TA, ALT=[TAA])
Record(CHROM=1, POS=1671823, REF=A, ALT=[T])
Record(CHROM=1, POS=1765789, REF=A, ALT=[G])
Record(CHROM=1, POS=2152975, REF=A, ALT=[T])
Record(CHROM=1, POS=2243386, REF=T, ALT=[TG])
Record(CHROM=1, POS=2498920, REF=T, ALT=[A])こんなふうにALTはlistになっているから上のスクリプトではうまく出力できてないのだな。
なのでlistをjoinで文字列に置き換えてから引き渡せばいいのだが、listの中身がstrになってないのでまずはmapで置き換え、更にjoinすると良さそうだ。
import vcf
def format_str(x):
'''
Change numbers to a string, and change list to a string.
'''
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'])
ALT = map(str, snp_record.ALT)
ALT = ",".join(ALT)
new_record = [format_str(snp_record.CHROM),
format_str(snp_record.POS),
format_str(snp_record.REF),
format_str(ALT)]
new_record.extend(gt_list)
outfh.write('\t'.join(new_record) + '\n')
parse_vcf('sample.vcf', 'gtdata.txt')これで
1 28499 A C 0/1 1 608335 T G 0/1 1 874461 GAA GAAA 1/1 1 895287 GTCA G 0/1 1 911789 TA TAA 1/1 1 1671823 A T 0/1 1 1765789 A G 0/1 1 2152975 A T 0/1 1 2243386 T TG 1/1 1 2498920 T A 0/1
こうなった。
とここまでやっておきながら別の方法を発見。
snp_record.heterozygosity
を読み出すとheteroなら0.5 homoなら0が出力される。これでいいんじゃね？
更に
snp_record.num_het
だとheteroで１が出力されるし。
つぎはINFOからも抽出。ついでに入出力の引き渡しも合わせて最終的に出来上がったプログラムは以下の通り。
import vcf, sys
def format_str(x):
'''
Change numbers to a string, and change list to a string.
'''
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'])
ALT = map(str, snp_record.ALT)
ALT = ",".join(ALT)
GENE = ""
if 'GENE' in snp_record.INFO:
GENE = snp_record.INFO['GENE']
SNP = ""
if 'SNP_ID' in snp_record.INFO:
SNP = snp_record.INFO['SNP_ID']
new_record = [format_str(snp_record.CHROM),
format_str(snp_record.POS),
format_str(snp_record.REF),
format_str(ALT),
format_str(SNP),
format_str(GENE),
format_str(snp_record.num_het)]
# new_record.extend(gt_list)
outfh.write('\t'.join(new_record) + '\n')
args = sys.argv
parse_vcf(args[1], args[2])
