How to manipulate a reference FASTA or bam to include variants from a VCF?

Question

I have some software which takes fastas as the input. I need to include SNVs and InDels from a VCF into the reference hg38 and then use this.

The problem is, I don't know of an algorithmically sound way to do this.

1. Are there any existing software packages which could do this efficiently? Is it easier to output a FASTA, or a bam (and then convert to a FASTA)?

2. What about if I wanted to do the same with a bedpe of germline structural variants?

Answer 1

You could convert VCF to BED via vcf2bed --snvs, vcf2bed --insertions, and vcf2bed --deletions, and then use samtools faidx by way of a wrapper script to convert BED to FASTA, e.g.:

$ vcf2bed --snvs < variants.vcf | bed2faidxsta.pl > snvs.fa
$ vcf2bed --insertions < variants.vcf | bed2faidxsta.pl > insertions.fa
$ vcf2bed --deletions < variants.vcf | bed2faidxsta.pl > deletions.fa

You need FASTA files for your reference genome, which have been indexed with samtools faidx, e.g., for hg38:

$ cd /foo/bar/baz
$ wget ftp://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/*.fa.gz
$ for fn in `ls *.fa.gz`; do gunzip $fn; done
$ for fn in `ls *.fa`; do samtools faidx $fn; done

Once you have indexed FASTA files somewhere on your file system, you can pipe BED to the bed2faidxsta.pl script, to get out FASTA sequences.

Answer 2

GATK has a solution that might work for you: FastaAlternateReferenceMaker, which : "Given a variant callset, this tool replaces the reference bases at variation sites with the bases supplied in the corresponding callset records."

Input

The reference, requested intervals, and any number of variant ROD files.

Output

A FASTA file representing the requested intervals.

Usage example

java -jar GenomeAnalysisTK.jar \
   -T FastaAlternateReferenceMaker \
   -R reference.fasta \
   -o output.fasta \
   -L input.intervals \
   -V input.vcf \

Answer 3

There's a vcf2fq sub-program that was written as part of vcfutils to convert a VCF file into a fastq file given a reference sequence. Unfortunately this doesn't work properly with INDELs (it will just mask them, rather than actually converting them), so I wrote a modification to implement INDEL correction as well:

./vcf2fq.pl -f <input.fasta> <all-site.vcf> > <output.fastq>

It adds and removes sequence from the reference, as specified by the INDEL information in the VCF file. The code keeps a record of where in the reference sequence the INDELs should go, then inserts (or deletes) the INDELs, updating the reference position as necessary. It should work with things like deletions right next to insertions, but hasn't been extensively tested.

Answer 4

Here's an update how this work in 2019. (The answer is quite similar to my update in this thread):

bgzip and index your vcf file.

$ bgzip -c input.vcf > input.vcf.gz
$ tabix input.vcf.gz

Now use bcftools consensus:

$ bcftools consensus -f genome.fa input.vcf.gz > consensus.fa