text: user146: I am interested in converting a VCF file into a FASTA file given a reference sequence with Python or R. 
  Samtools/BCFtools (Heng Li) provides a Perl script vcfutils.pl which does this, the function vcf2fq (lines 469-528)
  This script has been modified by others to convert InDels as well, e.g. this by David Eccles 
  ./vcf2fq.pl -f <input.fasta> <all-site.vcf> > <output.fastq>
  
  Given that converting a VCF into a FASTA given some reference is a somewhat standard procedure in bioinformatics, are there standard R or Python implementations used by the community?
  
  user776: Assuming you have samtools-indexed FASTA files for your reference genome, and have BEDOPS and htslib installed, you could use bed2faidxsta.pl to convert BED to FASTA:
  $ vcf2bed < variants.vcf | bed2faidxsta.pl --fastaDir=/path/to/fasta > variants.fa
  
  To index your FASTA files, you could do the following:
  $ cd /path/to/fasta
  $ bgzip -c my_genome.fa > my_genome.fa.bgz
  $ samtools faidx my_genome.fa.bgz
  
  
  user492: I tried the pyfaidx approach reported by Jared, but I got unhandled exceptions that I could not handle easily within the framework. I had bcftools lying around, so I tried bcftools consensus and it worked like a charm. However, you do need an indexed VCF.
  bcftools consensus all-site.vcf.gz < input.fasta > output.fasta
  
  NB: bcftools consensus has a few options specified with the --haplotype argument for choosing which alleles should be incorporated into the FASTA file.
  
  user1269: You can whip up something quite easily in Python using pyfaidx, which allows you to access a FASTA file by genomic coordinates. So you'd just have to pull the coordinates from the VCF to get the sequence, then output in FASTA format with whatever unique identifier for each sequence you want to use.
  Edit: The FastaVariant class is new in pyfaidx, I hadn't seen it before. It makes this even easier, as you won't have to manually mess with the VCF to get the positions. Something like the following should work:
  consensus = FastaVariant('yourgenome.fasta', 'variants.vcf.gz', het=True, hom=True)
  
  out = open("variants.fasta", "w")
  
  for chrom in consensus.keys:
      for var in consensus[chrom].variants_sites:
          record = consensus[chrom][var-1:var]
          print(record, file=out)
  
  out.close()
  
  
original_text: 146: <p>I am interested in converting a VCF file into a FASTA file given a reference sequence with Python or R. </p>
  
  <p>Samtools/BCFtools (Heng Li) provides a Perl script <a href="https://github.com/lh3/samtools/blob/master/bcftools/vcfutils.pl" rel="noreferrer"><code>vcfutils.pl</code></a> which does this, the function <code>vcf2fq</code> (lines 469-528)</p>
  
  <p>This script has been modified by others to convert InDels as well, e.g. <a href="https://github.com/gringer/bioinfscripts/blob/master/vcf2fq.pl" rel="noreferrer">this</a> by David Eccles </p>
  
  <pre><code>./vcf2fq.pl -f &lt;input.fasta&gt; &lt;all-site.vcf&gt; &gt; &lt;output.fastq&gt;
  </code></pre>
  
  <p>Given that converting a VCF into a FASTA given some reference is a somewhat standard procedure in bioinformatics, are there standard R or Python implementations used by the community?</p>
   776: <p>Assuming you have <code>samtools</code>-indexed FASTA files for your reference genome, and have <a href="http://bedops.readthedocs.io" rel="nofollow noreferrer">BEDOPS</a> and <a href="https://github.com/samtools/htslib" rel="nofollow noreferrer">htslib</a> installed, you could use <a href="https://gist.github.com/alexpreynolds/fa9b0f90e181e3b4f640" rel="nofollow noreferrer">bed2faidxsta.pl</a> to convert BED to FASTA:</p>
  
  <pre><code>$ vcf2bed &lt; variants.vcf | bed2faidxsta.pl --fastaDir=/path/to/fasta &gt; variants.fa
  </code></pre>
  
  <p>To index your FASTA files, you could do the following:</p>
  
  <pre><code>$ cd /path/to/fasta
  $ bgzip -c my_genome.fa &gt; my_genome.fa.bgz
  $ samtools faidx my_genome.fa.bgz
  </code></pre>
   492: <p>I tried the <code>pyfaidx</code> approach reported by Jared, but I got unhandled exceptions that I could not handle easily within the framework. I had <code>bcftools</code> lying around, so I tried <a href="https://samtools.github.io/bcftools/howtos/consensus-sequence.html" rel="nofollow noreferrer"><code>bcftools consensus</code></a> and it worked like a charm. However, you do need an indexed VCF.</p>
  <pre class="lang-bash prettyprint-override"><code>bcftools consensus all-site.vcf.gz &lt; input.fasta &gt; output.fasta
  </code></pre>
  <p>NB: <code>bcftools consensus</code> has <a href="https://samtools.github.io/bcftools/bcftools.html#consensus" rel="nofollow noreferrer">a few options</a> specified with the <code>--haplotype</code> argument for choosing which alleles should be incorporated into the FASTA file.</p>
   1269: <p>You can whip up something quite easily in Python using <a href="https://github.com/mdshw5/pyfaidx" rel="nofollow noreferrer">pyfaidx</a>, which allows you to access a FASTA file by genomic coordinates. So you'd just have to pull the coordinates from the VCF to get the sequence, then output in FASTA format with whatever unique identifier for each sequence you want to use.</p>
  
  <p>Edit: The <code>FastaVariant</code> class is new in <code>pyfaidx</code>, I hadn't seen it before. It makes this even easier, as you won't have to manually mess with the VCF to get the positions. Something like the following should work:</p>
  
  <pre><code>consensus = FastaVariant('yourgenome.fasta', 'variants.vcf.gz', het=True, hom=True)
  
  out = open("variants.fasta", "w")
  
  for chrom in consensus.keys:
      for var in consensus[chrom].variants_sites:
          record = consensus[chrom][var-1:var]
          print(record, file=out)
  
  out.close()
  </code></pre>
  