VCFファイルをPythonで扱う（cyvcf2の紹介）

アメリエフ2年目になりました、iijm_lです🐰

最近は、変異データと仲良しこよしなので、
勉強がてら新しいVCFパースツールを試してみたいと思います。

以前弊社ブログにて、PyVCFをご紹介しました。

staffblog.amelieff.jp

今回は、cyvcf2というツールを使って、サンプルコードを作成してみます。

cyvcf2のcyは、Cythonから取っているようです。
Pythonモジュールは計算スピードがネックになることも多いですが、cyvcf2 では htslib を活用することでパフォーマンス向上を図っています。

↓↓　以下、論文です。↓↓
cyvcf2: fast, flexible variant analysis with Python | Bioinformatics | Oxford Academic

インストール

anacondaで仮想環境の作成と同時にインストールしました。

# 仮想環境の作成およびcyvcf2のインストール
$ conda create -n cyvcf2test python=3.6 cyvcf2 -c bioconda
# 仮想環境にログイン
$ conda activate cyvcf2test

サンプルコード

cyvcf2を用いて、SNV(Single Nucleotide Variant)とIndel(insertion/deletion)の集計を行うスクリプトを書いてみました。

$ vim snv_or_indel.py

### snv_or_indel.py

import sys
from cyvcf2 import VCF

# 第1引数: VCFファイル
vcf_file = sys.argv[1]

# VCFの染色体のリストを取得する
chrom_list = VCF(vcf_file).seqnames

# 集計するためのdictを作成する
count_dict = {}
for chrom_name in chrom_list:
    count_dict[chrom_name] = {'SNV': 0, 'indel': 0}

# VCFの変異データを1行ずつ読み込む
for variant in VCF(vcf_file):
    chrom = variant.CHROM
    ref = variant.REF
    # multi-allele を除外
    if len(variant.ALT) != 1:
        continue
    alt = variant.ALT[0]
    # SNV or Indel ?
    if len(ref) == 1 and len(alt) == 1:
        # SNV
        count_dict[chrom]['SNV'] += 1
    else:
        # Indel
        count_dict[chrom]['indel'] += 1

for chrom, count in count_dict.items():
    snv = str(count['SNV'])
    indel = str(count['indel'])
    print(f'{chrom},{snv},{indel}')

真ん中のループ内部の処理をざっくり書くと、こんな感じ。

• 変異データからREF、ALTを取得
  
• マルチアレル（複数のALTを持つ）を除外
  
• REFとALTが1文字ならばSNV、それ以外はIndelとして染色体別にカウント
  

注意点として、variant.REFはstr型ですが、variant.ALTはlist型で返ってきます。
これは、VCFのALT列には複数のアレル情報が記載される場合があるからです。


cyvcf2では「.vcf」はもちろん、「.vcf.gz」も同じコードで読み込めます！
実行してみますと〜

$ python snv_or_indel.py sample_all_fil.vcf 
chr1,68769,34878
chr2,58456,30159
chr3,49001,25605
chr4,42792,19755
chr5,40020,21724
chr6,45415,21602
chr7,37937,19121
chr8,36292,18192
chr9,30670,15405
chr10,37861,18912
chr11,39330,19755
chr12,36197,18558
chr13,22115,9853
chr14,23019,11597
chr15,23583,12024
chr16,27580,13249
chr17,23792,14089
chr18,18136,8958
chr19,20656,10952
chr20,18714,10370
chr21,9935,4632
chr22,11950,6547
chrX,13953,10174
chrY,993,932
chrM,41,5

大体SNV：Indel = 2:1 くらいという結果に...🧐


バイオデータを簡単に扱えるツールは他にもたくさんありますので、お試しください😃
計算パフォーマンスを比較してみても面白いですね！