【VCFPyの使い方】VCFファイルの変異情報を分析する

VCF形式はゲノム解析において重要なファイルフォーマットであり、VCFPyはそれを扱うための効率的なPythonライブラリです。

本記事では、VCFファイルフォーマットの構造と特徴、およびVCFファイルをpythonで取り扱えうことのできるVCFPyの使い方を解説します。

VCFとは

VCF (Variant Call Format) はゲノムの塩基変異データを保存するために広く使用されているファイル形式です。シーケンシングデータをリファレンス配列にマッピングした際の情報を記録します。具体的には、リファレンス配列上の塩基とそこにマッピングされたシーケンシングデータ上の塩基などの情報が記載されています。

VCFファイルは、次世代シーケンサーから得られる膨大な変異情報を構造化された形式で保存し、効率的に扱うことを可能にします。

VCFファイルの例

以下は、最小限の情報を持つVCFファイルの例です。このようなファイルをテキストエディタで作成し、.vcf拡張子で保存できます。

4行目の1と12345の間など、間隔が空く部分はタブ区切りとなります。

##fileformat=VCFv4.2
##source=Example
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Sample1
1 12345 rs123 A G 50 PASS . GT 0/1
1 67890 . T C 99 PASS . GT 1/1

VCFファイルのフォーマット

VCFファイルは、大きく分けて以下の部分から成り立っています。

1. メタ情報セクション：ファイルの先頭にある、「##」から始まる行で、アノテーション情報などが記載されています。

2. ヘッダー行：「#」から始まる行で、データ列の名前が記載されています。

3. データ行：実際の変異データが記載されています。それぞれの行は一つの変異に対応し、タブ区切りのフィールドで構成されています。

VCFの主要フィールド

VCFファイルのデータ行には、以下の情報が含まれています。

• CHROM：染色体番号

• POS：塩基の位置（1開始）

• ID：変異の識別子（例: dbSNP ID rs123）。.の場合、識別子なし

• REF：リファレンス配列上の塩基

• ALT：変異した塩基配列（例: G,Tなど複数も可）。

• QUAL：変異コールの品質スコア（Phredスケール、例: 99）。

• FILTER：フィルタリング結果（例: PASSは合格、または特定条件名）。

• INFO：追加情報。セミコロン区切りで記載（例: DP=100;AF=0.5など）。

• FORMAT：サンプルデータの形式を示す（例: GT:DP:AD）。
• Sample1, Sample2…：各サンプルごとのデータ（FORMAT列に従って記載される）。

QUALの補足

• QUALは、変異コールの信頼性を示すスコアで、Phredスケールで表現されます。

• このスコアは、ALT列に記載された変異が誤りである確率を示します。

• 計算式：

  $$\text{QUAL}=-10\cdot {\log }_{10}(\text{誤り確率})$$ 

  • 例: QUAL=30の場合、誤り確率は 

    ${10}^{-3}$10−3 (0.1%)。

高いQUALスコアの意味

• QUAL値が高いほど、その変異が正しい（リファレンス配列と異なる）とする信頼性が高いことを意味します。

• 一般的な基準：

  • QUAL ≥ 30: 高信頼度の変異コール。

  • QUAL < 20: 信頼性が低く、フィルタリング対象となる可能性が高い。

QUAL値が「.」の場合

• QUALが「.」の場合、その位置について品質スコアが計算されていないことを意味します。

注意点

• QUALスコアは、リード深度や他の要因に依存するため、単独で使用するよりも他のフィルタリング基準（例: INFO列のDPやMQ）と組み合わせて評価することが推奨されます。

FILTERの補足

• FILTER列には、各変異が特定のフィルタリング基準を満たしているかどうかが記載されています。

• 主な値：

  • PASS: 全てのフィルタリング基準を通過したことを示します。

  • フィルタ名: 特定のフィルタ基準を満たさなかった場合、そのフィルタ名（例: low_qual, strand_bias）が記載されます。

  • “.”: フィルタリングが適用されていないことを示します

PASSの意味

• 「PASS」は、変異が全ての事前定義されたフィルタリング基準（例えば、品質スコアやリード深度など）を満たしていることを意味します

• 一般的なフィルタリング基準：

  • QUAL ≥ 30

  • リード深度 (DP) が十分に高い

  • 対立遺伝子頻度 (AF) が適切な範囲内

  • マッピング品質 (MQ) が高い

特定フィルタ名の例

以下はFILTER列に記載される可能性があるフィルタ名とその意味です

• low_qual: QUAL値が低すぎる。

• strand_bias: 対立遺伝子が片方向からのみ検出された。

• low_allele_frac: 対立遺伝子頻度が低すぎる。

• contamination: 異常な混入による変異と判断された。

INFOの補足

INFO列は、変異に関する追加情報を記載するフィールドです。複数の情報がセミコロン (;) 区切りで記載されており、それぞれのキーと値がペアになっています。

以下は、よく使われるINFOフィールドのキーとその意味です。

FORMATの補足

FORMAT列は、各サンプルに関するデータ形式を定義するフィールドです。複数のフィールドがコロン (:) 区切りで記載されており、これに従ってサンプルデータが記録されます。

以下は、よく使われるFORMATフィールドのキーとその意味です。

VCFPyの使い方

VCFPyとは

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

下記のリンクからVCFPyの公式ドキュメントを確認できます。

VCFPyのインストール

VCFPyはpipコマンドを使用して簡単にインストールすることができます。

pip install vcfpy

ファイルの読み込み

ここで使用するサンプルとなるVCFファイルはVCFPyのGithubリポジトリから、input.vcfをダウンロードします。

import vcfpy

# VCFファイルを読み込む
reader = vcfpy.Reader.from_path('input.vcf') 

# ヘッダー行を作成して表示する
header = ['#CHROM', 'POS', 'REF', 'ALT'] + reader.header.samples.names 

# データ行を1行ずつ取り出し、加工して表示する
for record in reader:
    if not record.is_snv():
        continue
    line = [record.CHROM, record.POS, record.REF]
    line += [alt.value for alt in record.ALT]
    line += [call.data.get('GT') or './.' for call in record.calls]
    print('\t'.join(map(str, line)))

#CHROM POS REF ALT BLANK NA12878 NA12891 NA12892 NA19238 NA19239 NA19240chr22 42522392 G A 0/0 0/1 0/1 0/0 0/0 0/0 0/0 
chr22 42522597 C T 0/1 0/0 0/0 0/0 0/0 0/0 0/0 
chr22 42522613 G C 0/1 0/1 0/0 0/1 0/1 0/1 0/1 
(以下省略)

まとめ

本記事では、VCFファイルフォーマットの構造と特徴を解説し、さらにVCFファイルをpythonで取り扱うことのできる、VCFPyの使い方を解説しました。

VCFファイルはテキスト形式ですが、非常に大きなファイルとなるため、なかなかそのまま開くことが難しいです。今回のVCFPyのようなツールを利用して、ぜひ中身を色々確認してみてください！