How to manipulate a reference FASTA or bam to include variants from a VCF? 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.


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)? 
What about if I wanted to do the same with a bedpe of germline structural variants?
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 \