Fix dash/empty ALT alleles for deletions in VCF

Question

I have a VCF file that has deletions specified like this, with the REF allele containing only the deleted portion and the ALT allele just saying -:

chr1 101 . AAA -

I would like to combine this with another VCF in which deletions are specified like this, with REF and ALT both extended by 1bp to the left to ensure neither is empty:

chr1 100 . GAAA G

Is there a tool to convert from one of these representations to the other? I have already tried vt normalize and bcftools norm but neither of them seem to do this.

I don't really mind which representation I end up with although I think the second is the one more consistent with the spec.

Answer 1

According to the specs the ALT is not allowed to contain -. While technically supported, 'N' as the alternative allele is sometimes used unnecessarily and this script can also fix those cases if one exchanges the two "-" within the fix_record function for "N".

Here is a little python (>=3.6) script that is able to fix it. It makes use of the pysam module. You need a reference file and this fasta file must be indexed by samtools faidx reference.fa. Please note that it doesn't support multiple values in the ALT column for records that contain a -.

Usage:

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

The script:

#!/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(region=f"{record.chrom}:{record.pos - 1}-{record.pos - 1}")    
            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()

A lot of the code above is boilerplate for running a script on the command line. The critical code for the transformation is in the fix_record function. The relevant code is below, expanded, with line numbers and a description of each command.

[1]   ref = self.fasta.fetch(
          region=f"{record.chrom}:{record.pos - 1}-{record.pos - 1}"
      )
[2]   record.pos = record.pos - 1
[3]   record.ref = f"{ref}{record.ref}"
[4]   record.alts = tuple(
          map(
              lambda alt: ref if alt == "-" else f"{ref}{alt}",
              record.alts
          )
      )

1. This grabs the base preceding the inserted/deleted sequence.
2. This corrects the position of the variant to account for prepending another nucleotide.
3. This defines the new corrected REF allele: concatenate the base grabbed in [1] with the original REF allele
4. Apply an anonymous (lambda) function to each ALT allele; the function will replace - with the base grabbed in [1], or preprend the base to the allele if it is not -

Answer 2

Probably it is better to use dedicated tools, however, here is a quick R workaround. It is relatively straightforward to remove the extra base (first letter of column 4) in the second type of output:

# just a workaround to encolose the second file in a data frame
# you can just read your file as a table using fread(), read.csv(), ...
test <- data.frame(
  col1 = "chr1",
  col2 = 100,
  col3 = ".",
  col4 = "GAAA",
  col5 = "G",
  stringsAsFactors = FALSE
)

> test
  col1 col2 col3 col4 col5
1 chr1  100    . GAAA    G

# remove the first letter in column 4
test$col4 <- substring(test$col4, 2)
# converting column 4 to "-"
test$col5 <- "-"

> test
  col1 col2 col3 col4 col5
1 chr1  100    .  AAA    -

EDIT:

Here is how it can be done with awk for speed and memory concerns: The logic is the same as above, remove the first letter in column 4 and then change column 5 to "-".

awk '{$4 = substr($4,2); print}' infile > outfile
awk '{$5 = "-"; print}' outfile > outfile2

infile

chr1 100 . GAAA G

outfile

chr1 100 . AAA G

outfile2

chr1 100 . AAA -

Answer 3

One more thing. If empty ALT allele is not "-", you can change it such as "*", as follows.

    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))

Answer 4

This is my first post and I will add some explanations. When I tried to use this script, print(*map(self.fix_record, self.vcf), sep="", end=""), it will increase memory usage because it will load everything into memory and then start to apply fix_record function. So, I just used two lines for line in self.vcf: print(self.fix_record(line), sep="", end=""). Because pysam uses generator by default, if we process one by one by pysam's VariantFile function, it will consume very little memory. In my case, the VCF file was more than 200GB, so I had to fix this script for my VCF file.

Second my VCF file used "*" for empty ALT, so I fixed two lines (1) if "*" in record.alts: (2) record.alts = tuple(map(lambda alt: ref if alt == "*" else f"{ref}{alt}", record.alts))

#!/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="")
        for line in self.vcf:
            print(self.fix_record(line), sep="", end="")

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

Answer 5

Edited version of Namshin's code this can handle if REF or ALT is '-' ( these records in vcf file are observed when there is insertion/deletion). Added documentation in the script to explain code block I have modified. Hope this useful. Tested example with GRCh37

python fix_vcf.py -f genome.fa -o out_fixed.vcf input.vcf

example input.vcf:

##fileformat=VCFv4.2
#CHROM  POS ID  REF ALT QUAL  FILTER  INFO
7 6027157 . - T 1 PASS  NS=1  GQ  10
9 98218616  . - C 1 PASS  NS=1  GQ  10
17  41245200  . AT  - 1 PASS  NS=1  GQ  10
13  32929213  . C - 1 PASS  NS=1  GQ  10

out_fixed.vcf:

##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
#CHROM  POS ID  REF ALT QUAL  FILTER  INFO
7 6027157 . T TT  1 PASS  NS=1
9 98218616  . A AC  1 PASS  NS=1
17  41245199  . GAT G 1 PASS  NS=1
13  32929212  . CC  C 1 PASS  NS=1

# -*- 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.

        If the REF column of the records contains a `-`, the variant position will change.
        For reaching this the reference base for the position is determined and prepended to the current ALT bases.
        The REF value `-` is replaced with the determined ref base on position. 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 - 2, end=record.pos -1)

            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))

        if "-" in record.ref:
            ref = self.fasta.fetch(reference=record.chrom, start=record.pos - 1, end=record.pos)

            record.pos = record.pos
            record.ref = f"{ref}"
            record.alts = tuple(map(lambda alt: 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="")
        for line in self.vcf:
            print(self.fix_record(line), sep="", end="")

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