11
$\begingroup$

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?

Improve this question
$\endgroup$
5
  • $\begingroup$ 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? $\endgroup$
    – Bioathlete
    Commented Nov 12, 2017 at 14:06
  • 1
    $\begingroup$ "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 $\endgroup$
    – ShanZhengYang
    Commented Nov 12, 2017 at 18:41
  • $\begingroup$ 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? $\endgroup$
    – terdon ♦
    Commented Nov 13, 2017 at 10:51
  • $\begingroup$ @terdon Either works $\endgroup$
    – ShanZhengYang
    Commented Nov 13, 2017 at 13:56
  • 3
    $\begingroup$ Typically I find using the best tool is more important than what language it was written in. $\endgroup$
    – Bioathlete
    Commented Nov 13, 2017 at 15:05

3 Answers 3

Reset to default
5
$\begingroup$

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()
Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ 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? $\endgroup$
    – llrs
    Commented Nov 14, 2017 at 8:22
  • $\begingroup$ @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. $\endgroup$
    – Jared Andrews
    Commented Nov 14, 2017 at 18:52
  • 2
    $\begingroup$ It seems like FastaVariant ignores indels. $\endgroup$
    – Eli Korvigo
    Commented Apr 19, 2018 at 10:18
3
$\begingroup$

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.

Improve this answer
$\endgroup$
2
$\begingroup$

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
Improve this answer
$\endgroup$
2
  • $\begingroup$ Thank you. This wouldn't be too difficult to convert into Python/R, but I'm asking for scripts which have this already completed. $\endgroup$
    – ShanZhengYang
    Commented Nov 13, 2017 at 22:36
  • $\begingroup$ 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. $\endgroup$
    – Alex Reynolds
    Commented Nov 13, 2017 at 22:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.