Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-04-16T13:32:08.047Z
Answers.LastActivityDate: 2021-07-03T15:12:45.313Z
Answers.LastEditDate: 2021-07-03T15:12:45.313Z
Answers.UserDisplayName: winni2k
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-11-13T22:36:19.737Z
Answers.Text: Thank you. This wouldn't be too difficult to convert into Python/R, but I'm asking for scripts which have this already completed.
Answers.UserDisplayName: ShanZhengYang
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-11-13T22:46:42.487Z
Answers.Text: If you have Python available, you very probably already have Perl. In any case, you can use `system` in R to call external scripts (written in any language). See `?system` for more details.
Answers.UserDisplayName: Alex Reynolds
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-11-13T21:45:19.470Z
Answers.LastActivityDate: 2017-11-13T21:45:19.470Z
Answers.UserDisplayName: Alex Reynolds
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-11-14T08:22:15.503Z
Answers.Text: Hi Jared, could you expand your answer, with the code to use to perform this? It is my the first time with pyfaidx and I have problems finding how to use it. It seems something like `consensus = FastaVariant('tests/data/chr22.fasta', 'tests/data/chr22.vcf.gz', het=True, hom=True)` should work?
Answers.UserDisplayName: llrs
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-11-14T18:52:30.900Z
Answers.Text: @Llopis yes, it should be something like that. The `FastaVariant` class is actually new, I wasn't aware it existed, so thanks! I've updated my answer with code I believe should do as intended.
Answers.UserDisplayName: Jared Andrews
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-04-19T10:18:27.690Z
Answers.Text: It seems like `FastaVariant` [ignores indels](https://github.com/mdshw5/pyfaidx/blob/a088141d8c18d0420eb74a1ce1578c3c1ea7ef8f/pyfaidx/__init__.py#L1111).
Answers.UserDisplayName: Eli Korvigo
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-11-13T18:51:09.507Z
Answers.LastActivityDate: 2017-11-14T18:51:42.240Z
Answers.LastEditDate: 2017-11-14T18:51:42.240Z
Answers.UserDisplayName: Jared Andrews
Body: <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>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-11-12T14:06:10.273Z
Comments.Text: I use the pybedtools package which is a wrapper around bedtools. Bedtools will work on VCF files in addition to BED and GTF. Why do you want a Python or R implementation instead of using Heng's solution?
Comments.UserDisplayName: Bioathlete
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-11-12T18:41:20.553Z
Comments.Text: "Why do you want a Python or R implementation instead of using Heng's solution?" Sometimes it's easier to integrate R code with R code, and Python code with Python
Comments.UserDisplayName: ShanZhengYang
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-11-13T10:51:11.670Z
Comments.Text: Based on the name and a quick look at the code, the function you mention produces fastq and not fasta. Do you want fasta or fastq?
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-11-13T13:56:55.697Z
Comments.Text: @terdon Either works
Comments.UserDisplayName: ShanZhengYang
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-11-13T15:05:22.150Z
Comments.Text: Typically I find using the best tool is more important than what language it was written in.
Comments.UserDisplayName: Bioathlete
Tags: <r><python><fasta><vcf><samtools>
Title: Converting a VCF into a FASTA given a reference with Python, R
ThreadText: # Converting a VCF into a FASTA given a reference with Python, R
  
  (Asked by: username_0 on 2017-11-12)
  
  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`](https://github.com/lh3/samtools/blob/master/bcftools/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](https://github.com/gringer/bioinfscripts/blob/master/vcf2fq.pl) 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?
  
  **Tags:** r, python, fasta, vcf, samtools
  
  > **Comments:**
  > **username_1**: I use the pybedtools package which is a wrapper around bedtools. Bedtools will work on VCF files in addition to BED and GTF. Why do you want a Python or R implementation instead of using Heng's solution?  
  > **username_0**: "Why do you want a Python or R implementation instead of using Heng's solution?" Sometimes it's easier to integrate R code with R code, and Python code with Python  
  > **username_2**: Based on the name and a quick look at the code, the function you mention produces fastq and not fasta. Do you want fasta or fastq?  
  > **username_0**: @username_2 Either works  
  > **username_1**: Typically I find using the best tool is more important than what language it was written in.  
  
  ---
  
  ## 3 answers
  
  ### Answer by username_5 on 2017-11-13. Score: 5 | Accepted answer
  
  You can whip up something quite easily in Python using [pyfaidx](https://github.com/mdshw5/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()
  
  ```
  
  > **Comments:**
  > **username_6**: Hi username_5, could you expand your answer, with the code to use to perform this? It is my the first time with pyfaidx and I have problems finding how to use it. It seems something like `consensus = FastaVariant('tests/data/chr22.fasta', 'tests/data/chr22.vcf.gz', het=True, hom=True)` should work?  
  > **username_5**: @Llopis yes, it should be something like that. The `FastaVariant` class is actually new, I wasn't aware it existed, so thanks! I've updated my answer with code I believe should do as intended.  
  > **username_7**: It seems like `FastaVariant` [ignores indels](https://github.com/mdshw5/pyfaidx/blob/a088141d8c18d0420eb74a1ce1578c3c1ea7ef8f/pyfaidx/__init__.py#L1111).  
  
  ---
  
  ### Answer by username_3 on 2019-04-16. Score: 3
  
  I tried the `pyfaidx` approach reported by username_5, but I got unhandled exceptions that I could not handle easily within the framework. I had `bcftools` lying around, so I tried [`bcftools consensus`](https://samtools.github.io/bcftools/howtos/consensus-sequence.html) and it worked like a charm. However, you do need an indexed VCF.
  
  ```bash
  bcftools consensus all-site.vcf.gz < input.fasta > output.fasta
  
  ```
  
  NB: `bcftools consensus` has [a few options](https://samtools.github.io/bcftools/bcftools.html#consensus) specified with the `--haplotype` argument for choosing which alleles should be incorporated into the FASTA file.
  
  
  ---
  
  ### Answer by username_4 on 2017-11-13. Score: 2
  
  Assuming you have `samtools`-indexed FASTA files for your reference genome, and have [BEDOPS](http://bedops.readthedocs.io) and [htslib](https://github.com/samtools/htslib) installed, you could use [bed2faidxsta.pl](https://gist.github.com/alexpreynolds/fa9b0f90e181e3b4f640) 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
  
  ```
  
  > **Comments:**
  > **username_0**: Thank you. This wouldn't be too difficult to convert into Python/R, but I'm asking for scripts which have this already completed.  
  > **username_4**: If you have Python available, you very probably already have Perl. In any case, you can use `system` in R to call external scripts (written in any language). See `?system` for more details.  
  
  ---
  
  