3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Cross-posted on biostars: biostars.org/p/394866 $\endgroup$ – Ram RS Aug 19 at 13:45
  • 1
    $\begingroup$ 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 $\endgroup$ – Ram RS Aug 19 at 13:49
  • $\begingroup$ @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. $\endgroup$ – Nils Aug 19 at 14:25
4
$\begingroup$

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

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(
          reference=record.chrom,
          start=record.pos - 1,
          end=record.pos
      )
[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 -
$\endgroup$
  • $\begingroup$ +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. $\endgroup$ – Daniel Standage Aug 20 at 14:59
  • $\begingroup$ Done by adding type hints and doc strings :) $\endgroup$ – finswimmer Aug 20 at 16:59
  • $\begingroup$ Sorry, I should have been more clear. My edit is more of what I had in mind. :-) $\endgroup$ – Daniel Standage Aug 20 at 17:16
0
$\begingroup$

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 -
$\endgroup$
  • $\begingroup$ Thanks but my real VCF is several GB in size so I don't think this will work. $\endgroup$ – Nils Aug 19 at 11:59
  • $\begingroup$ I have edited my answer so as to include a solution making use of awk. $\endgroup$ – haci Aug 19 at 13:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.