# Fix dash/empty ALT alleles for deletions in VCF

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.

• Cross-posted on biostars: biostars.org/p/394866 – Ram RS Aug 19 '19 at 13:45
• OP, the first of the two representations is so outdated that most tools might error out. You're better off insisting on converting to format #2 – Ram RS Aug 19 '19 at 13:49
• @RamRS I agree. Would definitely prefer a way to convert format 1 to format 2 than vice versa. The other way round could be done with a simple awk solution though, since it doesn't require knowledge of the reference sequence. – Nils Aug 19 '19 at 14:25

According to the specs the ALT is not allowed to contain -.

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("-f", "--fasta", help="fasta reference file", required=True)

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 the unpacking mechanism using *see:
https://docs.python.org/3/tutorial/controlflow.html#tut-unpacking-arguments
"""
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 -
• +1, thanks for sharing a working script. It might be even more helpful for innocent bystanders to emphasize and explain the code that does the actual transformation. – Daniel Standage Aug 20 '19 at 14:59
• Done by adding type hints and doc strings :) – finswimmer Aug 20 '19 at 16:59
• Sorry, I should have been more clear. My edit is more of what I had in mind. :-) – Daniel Standage Aug 20 '19 at 17:16
• Thanks finswimmer for the script. It's been very helpful to me. I would like to extend it for the case of REF alleles with "-" as in: chr1 6029086 - A to: chr1 6029086 CA C I have made the attempt to do it by extend this solution for other cases but couldn't find the right way to do it for all of them. For example in: chr1 1684347 - CCT I would change it to (REF would change to "C"(the position grabbed)+ALT and ALT would change to the position grabbed("C") which worked in the previous example): chr1 1684347 CCCT C But it's wrong cause actually should be like: chr1 1684347 C CCCT How – Dagsbio Oct 24 '19 at 15:05
• Generator version of the python script. It won't consume your memmory at all. #print(*map(self.fix_record, self.vcf), sep="", end="") for line in self.vcf: print(self.fix_record(line), sep="", end="") – Namshin Kim Feb 16 '20 at 1:15

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

• Thanks but my real VCF is several GB in size so I don't think this will work. – Nils Aug 19 '19 at 11:59
• I have edited my answer so as to include a solution making use of awk. – haci Aug 19 '19 at 13:22

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

• Hi Namshin, thanks for contributing, but this is not the way we usually answer questions on this website. You should not assume that a reader reads the whole thread. If it's an extension for a different answer, meaning which one and also explain what is the benefit of your extension. Note you can  your answer. – Kamil S Jaron Feb 20 '20 at 11:54

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("-f", "--fasta", help="fasta reference file", required=True)

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.

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


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("-f", "--fasta", help="fasta reference file", required=True)

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.

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