VCFファイルの表記揺れについて

VCF（Variant Caller Format）は、一塩基変異（SNP）、挿入・欠損（INDEL）、および構造変異（SV）を表すために使用される標準化されたテキストファイルフォーマットです。詳細に知りたい方はGATKの説明ページを参照してください。

このフォーマットを用いると、変異に関する情報を格納することができるのですが、VCFの生成方法に応じてREF（変異前のアレル）及びALT（変異後のアレル）の記載方法が異なる場合が多々あります。

VCF同士をマージしたい場合であったり、特定のVCFフォーマットしか受け付けないツールを使いたい場合とかに困るわけです。今回はこれについて書きます。

なんのこっちゃ

いやよくわからんという方は下記の質問で取り上げられていた例を転記します。
下記の①、②はいずれも同じ変異（chr1, pos101-103のAAA欠損）を表しているわけですが、表記の仕方が違います。①では欠損した部分のみが記載されており、②では欠損した部分 + 手前の塩基が記載されています。
こういった違いが使用するツールによって異なってくるわけです。

①

#CHROM    POS    ID    REF    ALT
chr1      101    .    AAA     -

②

#CHROM    POS    ID    REF    ALT
chr1      100    .    GAAA    G

参考：Fix dash/empty ALT alleles for deletions in VCF

対処例

遺伝的変異データのオープンアクセスデータベースであるEVA (European Variation Archive) への登録に使用されるパイプライン（EVA pipeline）では、DB登録前にREF、ALT表現の正規化を行っているみたいです。
一部抜粋するとこんな感じ。Deletionに関しては共通項を削除してしまうようです。

参考：Allele normalisation - EVA

VCFの定義的にどうなん？

VCFの定義に立ち返ってみましょう。厳密なVCFの定義に関しては、 Samtools/Hts-specs のレポジトリにドキュメント化されているのでこれを参照します。

REF

各塩基配列はA,C,G,T,Nのいずれかを用いて表されます（大文字小文字区別なし）。複数の塩基配列が記載されることもあります。

以下、大事な部分なのでドキュメントから引用します。

The value in the POS field refers to the position of the first base in the String. For simple insertions and deletions in which either the REF or one of the ALT alleles would otherwise be null/empty, the REF and ALT Strings must include the base before the event (which must be reflected in the POS field), unless the event occurs at position 1 on the contig in which case it must include the base after the event
引用：The Variant Call Format Specification - 5 Representing variation in VCF records

POSはREFの最初の塩基を指しています。REFもしくはALTアレルがnull/emptyとなる場合、VCFではイベントが発生した塩基の1つ前の塩基を含めなくてはなりません（イベントがコンティグのポジション1で発生した場合を除く）。

なお、リファレンス配列がIUPAC ambiguity codesを含んでいる場合、定義上の最初の塩基に変換されます（R = A/GならA）。

ALT

各塩基配列はA,C,G,T,N,*のいずれか（大文字小文字区別なし）、変異が存在しないことを示す .、カッコつきの文字列（<INS>のような形式）、breakend replacement（t[p[のような形式）で表されます。マルチアレルの場合はカンマで区切られて表されます。
*はoverlapping deletionを示す際に用います。
参考：Spanning or overlapping deletions (* allele) - GATK Forum

使用例

おまけでアレルパターンごとのVCFの表記例を見ていきましょう。

例1

例	アレル配列	変異
Ref	`a t C g a`	Cがリファレンスと同じ
1	`a t G g a`	CがGに置換
2	`a t - g a`	Cが欠損
3	`a t CAg a`	AがCの後に挿入

VCFで表現するとこうなります。
注意しなければいけないのが、欠損変異の場合、欠損発生位置の1つ前のポジションから見た変異がVCFに記録されることになります（ポジション3で起きたＣの欠損変異はPOS:2, TC→T）。また、VCFはPOSでソートされます。

#CHROM    POS    ID    REF    ALT         QUAL    FILTER    INFO
20        2      .     TC     T           .       PASS      DP=100
20        3      .     C      G           .       PASS      DP=100
20        3      .     C      CA          .       PASS      DP=100

例2

例	アレル配列	変異
Ref	`a t C g a`	Cがリファレンスと同じ
1	`a t G g a`	CがGに置換
2	`a t - g a`	Cが欠損

この場合、アレルは tC, tG, t となります。
VCFで表現するとこうなります。

#CHROM    POS    ID    REF    ALT         QUAL    FILTER    INFO
20        2      .     TC     TC,TG,T     .       PASS      DP=100

例3

例	アレル配列	変異
Ref	`a t C g a`	Cがリファレンスと同じ
1	`a t - g a`	Cが欠損
2	`a t - - a`	CとGが欠損
3	`a t CAg a`	AがCの後に挿入

この場合、アレルは　tCg, tg, t, tCAg となります。
VCFで表現するとこうなります。

#CHROM    POS    ID    REF    ALT         QUAL    FILTER    INFO
20        2      .     TCG    TG,T,TCAG   .       PASS      DP=100

参考：The Variant Call Format Specification - 5 Representing variation in VCF records

つまり？

少し関係ないところも見ていきましたが、大事なのは太字にしたところです。
REF/ALTいずれかのアレル配列が空となる場合、VCFではイベントが発生した塩基の1つ前の塩基を含めなくてはならないのです。 -なんて記号使って空を表現してはならないのです。

正しい対処法

先ほどの質問の回答でfinswimmer氏がご丁寧にスクリプトを記載してくれております。
内容を理解した上で、ありがたく使わせて頂きましょう。

$ samtools faidx reference.fa

でリファレンス配列のインデックスを作成した後、下記のように使えるようです。

$ python fixvcf.py -f reference.fasta -o out.vcf in.vcf

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import argparse
import sys

import pysam


class Application:
    def __init__(self):
        self.args = Application.get_args()
        self.vcf = pysam.VariantFile(sys.stdin)
        self.fasta = pysam.FastaFile(self.args.fasta)

    @staticmethod
    def get_args() -> argparse.Namespace:
        """
        Read command line arguments and set input and output sources.

        :return: parsed command line arguments
        """
        parser = argparse.ArgumentParser(prog="fixvcf.py", description="Fixes deletion represented by -")
        parser.add_argument('--version', action='version', version='%(prog)s 0.1')

        parser.add_argument("input", help="input file, use - to read from stdin")
        parser.add_argument("-f", "--fasta", help="fasta reference file", required=True)
        parser.add_argument("-o", "--output", help="output file")

        if len(sys.argv) == 1:
            parser.print_help(sys.stderr)
            sys.exit(1)

        args = parser.parse_args()

        if args.output:
            sys.stdout = open(args.output, "w")

        if args.input != "-":
            sys.stdin = open(args.input, "r")

        return args

    def fix_record(self, record: pysam.libcbcf.VariantRecord) -> pysam.libcbcf.VariantRecord:
        """
        If the ALT column of the records contains a `-`, the variant position is shift one position to the left.
        For reaching this the reference base for the position is determined and prepended to the current REF bases.
        The ALT value `-` is replaced with the determined ref base on position before. In case there are multiple ALT
        values the determined ref base is prepended to them.

        :param record: vcf file record
        :return: fixed vcf record
        """
        if "-" in record.alts:
            ref = self.fasta.fetch(reference=record.chrom, start=record.pos - 1, end=record.pos)

            record.pos = record.pos - 1
            record.ref = f"{ref}{record.ref}"
            record.alts = tuple(map(lambda alt: ref if alt == "-" else f"{ref}{alt}", record.alts))

        return record

    def start(self):
        """
        Print the header of the vcf file. Afterwards the fix_record method is applied to each record in the vcf file.

        For more information about python's map() see: https://docs.python.org/3.6/library/functions.html#map

        For more information about the unpacking mechanism using `*`see:
        https://docs.python.org/3/tutorial/controlflow.html#tut-unpacking-arguments
        """
        print(self.vcf.header, end="")
        print(*map(self.fix_record, self.vcf), sep="", end="")


if __name__ == "__main__":
    app = Application()
    app.start()