This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Require suggestion for python tools that can implement VCF file into FASTA files

I am working on 1001 genomes from the 1001 genomes project. I wish to implement the VCF file into FASTA files incorporating the indels, etc to better work on the genomes. I know about GATK but is there any other tool that can help me achieve this, preferably in python environment?

fasta vcf
• 3,114 views • link
updated by Alban Nabla • Alban Nabla • written by Rachayita • Rachayita

First do you really mean ‘1001’ genomes or ‘1000’ genomes? Then, what do you mean by ‘implement’, do you really mean convert? Finally, it’s not clear how introducing indels into a fasta makes any sense. Why would having the information in fasta format help?

• 0 views • link
written by 4galaxy77 4galaxy77

Thank you for your response. 1001 genome project is a catalog of Arabidopsis thaliana genetic variation. I want to incorporate the SNPs of the VCF into the reference genome so that I have separate fasta files for all the individual genomes corresponding to all the individual VCFs.

• 0 views • link
written by Rachayita • Rachayita

2 answers

Have a look at bcftools consensus:

Create consensus sequence by applying VCF variants to a reference fasta file.

(Not a Python tool, though.)

• 0 views • link
written by Jukka Matilainen • Jukka Matilainen

Could you please give me an idea on how to do that? Thank you!

• 0 views • link
written by Rachayita • Rachayita

There's an examples section in the link.

• 0 views • link
written by 4galaxy77 4galaxy77

I think vcf-consensus-builder may be what you are after.

• 0 views • link
written by Alban Nabla • Alban Nabla

Thank you for your response. I looked into vcf consensus builder, but looks like it needs a depth file to work. I have just the vcf files and the reference genome fasta. Any other alternatives are highly appreciated.

• 0 views • link
written by Rachayita • Rachayita

I am afraid if you don't have access to the .bam files or a depth file, it will be hard to build a consensus.

I drafted a (clumsy) script to integrate the vcfs for one sequence/chromosome:

from Bio.Seq import Seq, MutableSeq
from cyvcf2 import VCF
from Bio import SeqIO

# parse the files and extract one chromosome
genome = SeqIO.parse("genome.fasta", "fasta")
chrom = next(genome)
chrom1A = chrom
chrom1B = chrom
vcf = VCF('variations.vcf')
chromvars = vcf('1')

def integrate(vcf,chrom):
    chromA = MutableSeq(chrom.seq)
    chromB = MutableSeq(chrom.seq)
    shiftA = 0
    shiftB = 0
    for var in vcf:
        if var.var_type == 'snp':
            if len(var.ALT) == 1:
                chromA[var.POS+shiftA-1] = var.ALT[0]
                chromB[var.POS+shiftB-1] = var.ALT[0]
            if len(var.ALT) == 2:
                chromA[var.POS+shiftA-1] = var.ALT[0]
                chromB[var.POS+shiftB-1] = var.ALT[1]
        if var.var_type == 'indel':
            reflen = len(var.REF)
            if len(var.ALT) == 1:
                chromA = chromA[:var.POS+shiftA-1] + var.ALT[0] + chromA[var.POS+shiftA+reflen:]
                chromB = chromB[:var.POS+shiftB-1] + var.ALT[0] + chromB[var.POS+shiftB+reflen:]
                shiftA = shiftA + len(var.ALT[0]) - reflen
                shiftB = shiftB + len(var.ALT[0]) - reflen
            if len(var.ALT) == 2:
                chromA = chromA[:var.POS+shiftA-1] + var.ALT[0] + chromA[var.POS+shiftA+reflen:]
                chromB = chromB[:var.POS+shiftB-1] + var.ALT[1] + chromB[var.POS+shiftB+reflen:]
                shiftA = shiftA + len(var.ALT[0]) - reflen
                shiftB = shiftB + len(var.ALT[1]) - reflen
    seq1 = Seq(chromA)
    seq2 = Seq(chromB)
    return seq1, seq2

# integrate using the above function
chromA, chromB = integrate(chromvars,chrom)
chrom1A.seq = chromA
chrom1B.seq = chromB

# output to fasta
SeqIO.write(chrom1A, 'modchrom1A.fsa', 'fasta')
SeqIO.write(chrom1B, 'modchro1B.fsa', 'fasta')

This is for one chromosome, but you can iterate through a full genome and automate the whole process. I am sure there are better ways but I hope this helps.

• 0 views • link
written by Alban Nabla • Alban Nabla

Log in to answer this question.

More posts like this