pythonでバイオインフォマティクス - VCFファイルを扱う

GATKやSAMtoolsなどの変異解析ツールを実行した際に、含まれていた一塩基変異（SNP）、挿入（Insert）、欠失（Deletion）等の情報はVCFファイルというフォーマットで出力されます。
VCFファイルのフォーマットに関しては、GATKの説明ページを参照してください。

このVCFファイルは、フォーマットとしてあまりイケていないため、新たなフォーマットが多数提案されてきました（以下一例）。

www.kimoton.com

とはいえ、現状各種ツールが出力する変異情報に関するフォーマットはVCFファイルのことが大半なため、VCFファイルを制する者は変異情報を制するといっても過言ではない訳です。
今回はそんなVCFファイルをPython経由で扱う例を見ていきます。

VCFファイルをとってくる

毎度お馴染み1000人ゲノムプロジェクトのデータを用います。
22番染色体についてVCFファイルをダウンロードし、bgzipコマンドで.vcf.gz形式に圧縮します。 tabixでは、tabix {url} {chr}:{start}-{end}で指定したURLからVCFファイルの一部を取得することができます。因みにtabixもhtslibの一部です。

$ tabix -fh ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/20130502/supporting/vcf_with_sample_level_annotation/ALL.chr22.phase3_shapeit2_mvncall_integrated_v5_extra_anno.20130502.genotypes.vcf.gz 22:1-17000000|bgzip -c > genotypes.vcf.gz

続いて、インデックスファイルを作成します。このファイルを作成することで、ファイルの読み込み速度を高速化することができます。

$ tabix -p vcf genotypes.vcf.gz

実行すると、genotypes.vcf.gz.tbiというファイルが作成されます。これがVCFファイルのインデックスファイルになります。

上記で使用しているtabixのオプションに関しては下記を確認してください。

オプション	説明
-f, --force	既にファイルが存在していた場合、強制的に上書きする
-h, --print-header	ヘッダー行も含めて出力する
-p, --preset	インデックス作成の入力形式を指定する。有効な値は、gff、bed、sam、vcf。

pyvcfでVCFファイルを扱う

pyvcfを使うと、VCFファイルをPythonのオブジェクトとして扱うことができます。

例の如くインストールはconda若しくはpipから可能です。

pyvcf でVCFファイルを読み込むにはvcf.Readerクラスを使用します。

import vcf
v = vcf.Reader(filename='genotypes.vcf.gz')

生成されるvcf.parser.ReaderオブジェクトはVCFファイルに含まれる変異情報のジェネレータになっています。

> print(next(v))
Record(CHROM=22, POS=16050075, REF=A, ALT=[G])
> print(next(v))
Record(CHROM=22, POS=16050115, REF=G, ALT=[A])

INFO列（リード全体の情報）

INFO列では、その変異箇所について、リード全体の統計情報などが記載されます。 vcf.parser.Readerオブジェクトでは、INFO列の項目の種類は.infos属性に格納されています。

infos = v.infos
for info in infos:
    print(info)
    
# CIEND
# CIPOS
# CS
# END
# IMPRECISE
# MC
# MEINFO
# MEND
# MLEN
# MSTART
# SVLEN
# SVTYPE
# TSD
# AC
# AF
# NS
# AN
# ASN_AF
# EUR_AF
# AFR_AF
# AMR_AF
# SAN_AF
# DP

CIENDや1000人ゲノムプロジェクトに特有の項目なのであまり気にしなくてよいです。
一般的な項目としては、それぞれ下記を表しています。

INFO列の項目	表している情報
NS	サンプル数
AC	検出された変異におけるALTアレルの総数
AF	推定アレル頻度
AN	検出された変異における
DP	リード全体のリードデプス

FORMAT列（サンプルレベルの情報）

INFO列と異なり、FORMAT列にはサンプルレベルの情報が格納されています。 vcf.Readerで読みこんだオブジェクトでは、FORMAT列の項目の種類は.formats属性に格納されています。

fmts = v.formats
for fmt in fmts:
    print(fmt)

# GT
# DP

このVCFファイルではサンプルの情報として、GTとDPを格納しているようです。
それぞれ下記を表しています。

INFO列の項目	表している情報
GT	遺伝子型
DP	サンプルごとのリードデプス

レコードの情報を確認

あるレコードについて、VCFファイルに格納された変異情報を取得する例を見ていきます。

v = vcf.Reader(filename='genotypes.vcf.gz')
rec = next(v)

染色体、ポジション、ID、変異前の塩基、変異後の塩基、クオリティスコア（Phreadスコア）、フィルターの状態を取得

> print(rec.CHROM, rec.POS, rec.ID, rec.REF, rec.ALT, rec.QUAL, rec.FILTER)
22 16050075 None A [G] 100 []

INFO列、FORMAT列の値を辞書型で取得

# INFO列の値を取得
> print(rec.INFO)
{'AC': [1], 'AF': [0.000199681], 'AN': 5008, 'NS': 2504, 'DP': [8012], 'ASN_AF': [0.0], 'AMR_AF': [0.0], 'SAS_AF': ['0.0010'], 'EUR_AF': [0.0], 'EAS_AF': [''], 'AFR_AF': [0.0], 'SAN_AF': [0.0]}

# FORMAT列の値を取得
> print(rec.FORMAT)
GT:DP

変異を持つサンプル情報を取得

> samples = rec.samples

# 変異を持つサンプルがいくつあるか
> len(samples)
2504

更にあるサンプルについて、対応する変異はREF・ALT列の何番目のモノか、ヘテロ接合か否か、フェージングの状態を取得します。

> sample = samples[0]
> print(sample.gt_alleles, sample.is_het, sample.is_variant, sample.phased)
> print(int(sample['DP']))
['0', '0'] False False True 1

変異の種類ごとの個数、マルチアレルを持つ変異数を取得します。
カウントの際はdefaultdictを使用して整数で初期化された辞書型を用いていることに注意してください。

import vcf
from collections import defaultdict
f = vcf.Reader(filename='genotypes.vcf.gz')

# defaultdictで整数初期化
my_type = defaultdict(int)
num_alts = defaultdict(int)

for rec in f:
    my_type[rec.var_type, rec.var_subtype] += 1
    if rec.is_snp:
        num_alts[len(rec.ALT)] += 1
print(my_type)
print(num_alts)

## defaultdict(<class 'int'>, {('sv', 'DEL'): 6, ('indel', 'del'): 273, ('snp', 'ts'): 10054, ('indel', 'ins'): 127, ('snp', 'tv'): 5917, ('sv', 'CNV'): 2, ('snp', 'unknown'): 79, ('sv', 'SVA'): 1, ('indel', 'unknown'): 13})
## defaultdict(<class 'int'>, {1: 15971, 2: 79})

結果から、転位（transition）が最も多く10054箇所、続いて転置（transversion）が5917箇所含まれていることが分かります。また、マルチアレルは79変異含まれているようです。

リードデプスを可視化

import vcf
from collections import defaultdict
f = vcf.Reader(filename='genotypes.vcf.gz')

# defaultdictで整数初期化
sample_dp = defaultdict(int)
for rec in f:
    # SNPでない、若しくはマルチアレルの場合はスキップ
    if not rec.is_snp or len(rec.ALT) != 1:
        continue
    # リードデプスのサンプル全体における分布を集計
    for sample in rec.samples:
        dp = sample['DP']
        if dp is None:
            dp = 0
        dp = int(dp)
        sample_dp[dp] += 1

# リードデプス昇順にソート
dps = list(sample_dp.keys())
dps.sort()
# 折れ線プロットを描画
dp_dist = [sample_dp[x] for x in dps]
fig, ax = plt.subplots(figsize=(16, 9))
ax.plot(dp_dist[:50], 'r')
# 頻度最大となるリードデプスに縦線を引く
ax.axvline(dp_dist.index(max(dp_dist)))

リードデプスの分布は右側に裾を引く形となっており、最頻値は6となっているようです。

その他のVCFパーサ

pyvcfのラッパーとして、vcf_parserというOSSがあります。

github.com

下記の用途では、vcf_parserを使用すると幸せになれるかもしれません（特に一番目）。

• 複数の対立遺伝子を持つバリアントラインを分割したい
• 読み込んだ生のvcfファイル形式での出力したい
• VEPによるアノテーションを上手にパースしたい