Question: Getting A Vcf File From A Fasta Alignment
4
gravatar for Bioch'Ti
4.4 years ago by
Bioch'Ti1000
France (Avignon)
Bioch'Ti1000 wrote:

Hi All,

I'm wondering if anybody is aware of a tool that converts a fasta alignment (that includes fasta sequences of different individuals) into a vcf file? I know that the opposite is feasible (vcf to fasta) using the GATK dedicated tool, but I want to perform the opposite !

Thanks for your help, C.

vcf fasta alignment conversion • 11k views
ADD COMMENTlink modified 4.4 years ago • written 4.4 years ago by Bioch'Ti1000

what's your input ? a CLUSTALW aln file ?

ADD REPLYlink written 4.4 years ago by Pierre Lindenbaum110k

The input file is a txt file (.fas for example) that contains only the aligned sequences (two lines per individual : one of comment, the second one for the sequence itself). Here is an example:

>Ind1
ACGTGGCTAGATCA
>Ind2
ACGTGGCTAGATCA
>Ind3
ACGTGCCTAGATCA
ADD REPLYlink modified 4.4 years ago by Pierre Lindenbaum110k • written 4.4 years ago by Bioch'Ti1000

how are managed the indels ?

ADD REPLYlink written 4.4 years ago by Pierre Lindenbaum110k

It will be nice to manage the indels too to make the script more general. But in my case, I do not need to consider these markers.

ADD REPLYlink written 4.4 years ago by Bioch'Ti1000

I added a support to MSA/Fasta. See below.

ADD REPLYlink written 4.4 years ago by Pierre Lindenbaum110k
8
gravatar for Pierre Lindenbaum
4.4 years ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum110k wrote:

I quickly wrote something for a CLUSTALW alignment (not fasta). See: https://github.com/lindenb/jvarkit/wiki/Biostar94573

It seems to work with the following clustalw input:

$ curl https://raw.github.com/biopython/biopython/master/Tests/Clustalw/opuntia.aln

CLUSTAL W (1.81) multiple sequence alignment


gi|6273285|gb|AF191659.1|AF191      TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA
gi|6273284|gb|AF191658.1|AF191      TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA
gi|6273287|gb|AF191661.1|AF191      TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA
gi|6273286|gb|AF191660.1|AF191      TATACATAAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA
gi|6273290|gb|AF191664.1|AF191      TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA
gi|6273289|gb|AF191663.1|AF191      TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA
gi|6273291|gb|AF191665.1|AF191      TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA
                                    ******* **** *************************************

gi|6273285|gb|AF191659.1|AF191      TATATA----------ATATATTTCAAATTTCCTTATATACCCAAATATA
gi|6273284|gb|AF191658.1|AF191      TATATATA--------ATATATTTCAAATTTCCTTATATACCCAAATATA
gi|6273287|gb|AF191661.1|AF191      TATATA----------ATATATTTCAAATTTCCTTATATATCCAAATATA
gi|6273286|gb|AF191660.1|AF191      TATATA----------ATATATTTATAATTTCCTTATATATCCAAATATA
gi|6273290|gb|AF191664.1|AF191      TATATATATA------ATATATTTCAAATTCCCTTATATATCCAAATATA
gi|6273289|gb|AF191663.1|AF191      TATATATATA------ATATATTTCAAATTCCCTTATATATCCAAATATA
gi|6273291|gb|AF191665.1|AF191      TATATATATATATATAATATATTTCAAATTCCCTTATATATCCAAATATA
                                    ******          ********  **** ********* *********

gi|6273285|gb|AF191659.1|AF191      AAAATATCTAATAAATTAGATGAATATCAAAGAATCCATTGATTTAGTGT
gi|6273284|gb|AF191658.1|AF191      AAAATATCTAATAAATTAGATGAATATCAAAGAATCTATTGATTTAGTGT
gi|6273287|gb|AF191661.1|AF191      AAAATATCTAATAAATTAGATGAATATCAAAGAATCTATTGATTTAGTGT
gi|6273286|gb|AF191660.1|AF191      AAAATATCTAATAAATTAGATGAATATCAAAGAATCTATTGATTTAGTGT
gi|6273290|gb|AF191664.1|AF191      AAAATATCTAATAAATTAGATGAATATCAAAGAATCTATTGATTTAGTGT
gi|6273289|gb|AF191663.1|AF191      AAAATATCTAATAAATTAGATGAATATCAAAGAATCTATTGATTTAGTAT
gi|6273291|gb|AF191665.1|AF191      AAAATATCTAATAAATTAGATGAATATCAAAGAATCTATTGATTTAGTGT
                                    ************************************ *********** *

gi|6273285|gb|AF191659.1|AF191      ACCAGA
gi|6273284|gb|AF191658.1|AF191      ACCAGA
gi|6273287|gb|AF191661.1|AF191      ACCAGA
gi|6273286|gb|AF191660.1|AF191      ACCAGA
gi|6273290|gb|AF191664.1|AF191      ACCAGA
gi|6273289|gb|AF191663.1|AF191      ACCAGA
gi|6273291|gb|AF191665.1|AF191      ACCAGA
                                    ******



$ curl https://raw.github.com/biopython/biopython/master/Tests/Clustalw/opuntia.aln" |\
  java -jar dist/biostar94573.jar

##fileformat=VCFv4.1
##Biostar94573CmdLine=
##Biostar94573Version=ca765415946f3ed0827af0773128178bc6aa2f62
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read="" depth"="">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Approximate read="" depth."="">
##contig=<ID=chrUn,length=156>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  gi|6273284|gb|AF191658.1|AF191  gi|6273285|gb|AF191659.1|AF191  gi|6273286|gb|AF191660.1|AF191  gi|6273287|gb|AF191661.1|AF191  gi|6273289|gb|AF191663.1|AF191  gi|6273290|gb|AF191664.1|AF191  gi|6273291|gb|AF191665.1|AF191
chrUn   8   .   T   A   .   .   DP=7    GT:DP   0:1 0:1 1:1 0:1 0:1 0:1 0:1
chrUn   13  .   A   G   .   .   DP=7    GT:DP   0:1 0:1 0:1 0:1 1:1 1:1 1:1
chrUn   56  .   ATATATATATA ATA,A,ATATA .   .   DP=7    GT:DP   1:1 2:1 2:1 2:1 3:1 3:1 0:1
chrUn   74  .   TCA TAT .   .   DP=7    GT:DP   0:1 0:1 1:1 0:1 0:1 0:1 0:1
chrUn   81  .   T   C   .   .   DP=7    GT:DP   0:1 0:1 0:1 0:1 1:1 1:1 1:1
chrUn   91  .   T   C   .   .   DP=7    GT:DP   1:1 1:1 0:1 0:1 0:1 0:1 0:1
chrUn   137 .   T   C   .   .   DP=7    GT:DP   0:1 1:1 0:1 0:1 0:1 0:1 0:1
chrUn   149 .   G   A   .   .   DP=7    GT:DP   0:1 0:1 0:1 0:1 1:1 0:1 0:1

Edit

I added a support for fasta/MSA

$ cat input.msa 
>Ind1
ACGTGGCTAGATCA
>Ind2
ACGTGGCTAGATCA
>Ind3
ACGTGCCTAGATCA

get the output:

$ cat input.msa | java -jar dist/biostar94573.jar 
##fileformat=VCFv4.1
##Biostar94573CmdLine=
##Biostar94573Version=4094eaed3dd2c9364309c20f88c6f79ad54a7450
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth.">
##contig=<ID=chrUn,length=14>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Ind1    Ind2    Ind3
chrUn   6   .   G   C   .   .   DP=3    GT:DP   0:1 0:1 1
ADD COMMENTlink modified 4.4 years ago • written 4.4 years ago by Pierre Lindenbaum110k

Hi Pierre,

I am looking for a program for converting fasta -> vcf and I found this discussion and your code very helpful. However, I have a couple of questions regarding this output - 

#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Ind1    Ind2    Ind3
chrUn   6   .   G   C   .   .   DP=3    GT:DP   0:1 0:1 1

Shouldnt the genotype for Ind1 = 0/0 and Ind2 = 0/0 and Ind3 = 1/1? Also does this program allow input as a consensus sequence - for instance, there could be sites such Y and R in the sequence and the vcf could translate these as a A/T or C/G?

ADD REPLYlink written 3.7 years ago by diviya.smith50

you're right:I've fixed the bug !  thank you

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  Ind1    Ind2    Ind3 chrUn   6       .       G       C       .       .       DP=3    GT:DP   0/0:1   0/0:1   1/1:1

 

consensus sequence: no (I don't have time for this)

 

ADD REPLYlink written 3.7 years ago by Pierre Lindenbaum110k

Thanks! This is very helpful.

ADD REPLYlink written 3.7 years ago by diviya.smith50

Hi Pierre, I have many vcf files now, generated courtesy of your script. I only need SNP counts for the purposes of my analyses. Is there a way to get SNP counts (only) for all my vcf files? I thought to use grep, but I noticed I was getting reports with indels, too.

 

ADD REPLYlink written 3.4 years ago by modernsynthesis0
1

use bcftools filter with TYPE="snp"  http://samtools.github.io/bcftools/bcftools.html#expressions

ADD REPLYlink written 3.4 years ago by Pierre Lindenbaum110k

[filter.c:1298 filters_init1] Error: the tag "INFO/snp" is not defined in the VCF header


error in vcf file generated by your script, when using bcftools filter with TYPE="snp".

ADD REPLYlink written 2.9 years ago by he.zhengshan0

https://vcftools.github.io/man_latest.html

vcftools --vcf input_file.vcf --remove-indels --recode --recode-INFO-all --out SNPs_only

I finally fix it!

ADD REPLYlink written 2.9 years ago by he.zhengshan0

hrUn    14005    .    TCG    TCA,TAG ##actually 2 snp

chrUn    28819    .    TTCA    TGAA,TTCC ##actually 3snp

Could you seperate them, not together?

 

ADD REPLYlink written 2.9 years ago by he.zhengshan0

Pierre, many thanks for providing this script, works very well. I used this script for one fasta alignment, but I would like to use it for thousands of fasta alignments. Is there a way to do this separately for each locus all at once?

 

ADD REPLYlink modified 3.4 years ago • written 3.4 years ago by modernsynthesis0

just concatenate

* the files  with cat ?

* or merge the vcf later ?

ADD REPLYlink written 3.4 years ago by Pierre Lindenbaum110k

Right, but I think that would be violating homology, no? Desired workflow: fasta > vcf > snp count (for each locus/alignment). Many thanks in advance.
 

ADD REPLYlink modified 3.4 years ago • written 3.4 years ago by modernsynthesis0

Apologies if I am not making any sense. Either way, thank you for providing script.

ADD REPLYlink written 3.4 years ago by modernsynthesis0

Should I use a for loop?

ADD REPLYlink written 3.4 years ago by modernsynthesis0

Hi Pierre,

 

Thanks for this! Would it be possible to output EVERY position on the reference as VCF format? (no coverage = ./., and if same as reference = 0/0)? Otherwise I can't merge VCF files of different samples together....

ADD REPLYlink modified 2.9 years ago • written 2.9 years ago by Coryza360

ask this as an issue at https://github.com/lindenb/jvarkit/issues . Anyway, I away from my sources for a few days.

ADD REPLYlink written 2.9 years ago by Pierre Lindenbaum110k
1
gravatar for Bioch'Ti
4.4 years ago by
Bioch'Ti1000
France (Avignon)
Bioch'Ti1000 wrote:

Hi,

I also found that PGDSpider (http://www.cmpg.unibe.ch/software/PGDSpider/) is able to convert fasta alignment into a vcf file. The problem is that it does not handle concatened alignements in a single fasta file (eg sequences of different sizes).

hope this helps Best, C.

ADD COMMENTlink written 4.4 years ago by Bioch'Ti1000
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1101 users visited in the last hour