10
$\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?

$\endgroup$
  • $\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 Nov 12 '17 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 Nov 12 '17 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 Nov 13 '17 at 10:51
  • $\begingroup$ @terdon Either works $\endgroup$ – ShanZhengYang Nov 13 '17 at 13:56
  • 3
    $\begingroup$ Typically I find using the best tool is more important than what language it was written in. $\endgroup$ – Bioathlete Nov 13 '17 at 15:05
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$
  • 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 Nov 14 '17 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 Nov 14 '17 at 18:52
  • 1
    $\begingroup$ It seems like FastaVariant ignores indels. $\endgroup$ – Eli Korvigo Apr 19 '18 at 10:18
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$
  • $\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 Nov 13 '17 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 Nov 13 '17 at 22:46
2
$\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 --haplotypes argument for choosing which alleles should be incorporated into the FASTA file.

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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