<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>

CC BY-SA 4.0
2019-04-16T13:32:08.047Z
2021-07-03T15:12:45.313Z
2021-07-03T15:12:45.313Z
winni2k
<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>

CC BY-SA 3.0
2017-11-13T22:36:19.737Z
Thank you. This wouldn't be too difficult to convert into Python/R, but I'm asking for scripts which have this already completed.
ShanZhengYang
CC BY-SA 3.0
2017-11-13T22:46:42.487Z
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.
Alex Reynolds
CC BY-SA 3.0
2017-11-13T21:45:19.470Z
2017-11-13T21:45:19.470Z
Alex Reynolds
<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>

CC BY-SA 3.0
2017-11-14T08:22:15.503Z
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?
llrs
CC BY-SA 3.0
2017-11-14T18:52:30.900Z
@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.
Jared Andrews
CC BY-SA 3.0
2018-04-19T10:18:27.690Z
It seems like `FastaVariant` [ignores indels](https://github.com/mdshw5/pyfaidx/blob/a088141d8c18d0420eb74a1ce1578c3c1ea7ef8f/pyfaidx/__init__.py#L1111).
Eli Korvigo
CC BY-SA 3.0
2017-11-13T18:51:09.507Z
2017-11-14T18:51:42.240Z
2017-11-14T18:51:42.240Z
Jared Andrews
<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>

CC BY-SA 3.0
2017-11-12T14:06:10.273Z
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?
Bioathlete
CC BY-SA 3.0
2017-11-12T18:41:20.553Z
"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
ShanZhengYang
CC BY-SA 3.0
2017-11-13T10:51:11.670Z
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?
terdon
CC BY-SA 3.0
2017-11-13T13:56:55.697Z
@terdon Either works
ShanZhengYang
CC BY-SA 3.0
2017-11-13T15:05:22.150Z
Typically I find using the best tool is more important than what language it was written in.
Bioathlete
<r><python><fasta><vcf><samtools>
Converting a VCF into a FASTA given a reference with Python, R
# 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 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   > 

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**: @terdon 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, 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 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?  
> **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 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.


---

### Answer by username_4 on 2017-11-13. Score: 2

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

> **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.  

---

