This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Vcf-Consensus To Generate Consensus File

I have a VCF file and a reference fasta file, and am trying to create a consensus sequence using vcf-consensus.

First I create a tabix file, and then I run

cat NC_002516.2.fa | vcf-consensus 1_e5_d30.vcf.gz > out.fa

and I get an error

The fasta sequence does not match the REF at NC_002516.2:837550. G() in .fa, GT in .vcf at /group/RDI/gatk/vcftools/vcftools_0.1.9/perl/vcf-consensus line 18 main::error('The fasta sequence does not match the REF at NC_002516.2:8375...') called at /group/RDI/gatk/vcftools/vcftools_0.1.9/perl/vcf-consensus line 129 main::apply_variant('HASH(0x13a8b2a0)', 'ARRAY(0x13e1ac70)') called at /group/RDI/gatk/vcftools/vcftools_0.1.9/perl/vcf-consensus line 87 main::do_consensus('HASH(0x13a8b2a0)') called at /group/RDI/gatk/vcftools/vcftools_0.1.9/perl/vcf-consensus line 9

Has anyone else encountered this? Looking at the VCF file for location 837550,

NC_002516.2 521975 . AC A 3149.48 PASS AC=2;AF=1.00;AN=2;DP=79;FS=0.000;HRun=1;HaplotypeScore=315.8867;MQ=59.23;MQ0=0;QD=39.87;SB=-1477.05 GT:AD:DP:GQ:PL 1/1:1,78:79:99:3192,238,0 NC_002516.2 574043 . A AGGG 5241.60 PASS AC=2;AF=1.00;AN=2;DP=78;FS=0.000;HRun=1;HaplotypeScore=335.5748;MQ=42.53;MQ0=0;QD=67.20;SB=-2258.53 GT:AD:DP:GQ:PL 1/1:14,63:78:99:5242,235,0 NC_002516.2 675693 . G A 1670.23 PASS AC=2;AF=1.00;AN=2;DP=69;Dels=0.00;FS=0.000;HRun=0;HaplotypeScore=16.7733;MQ=57.48;MQ0=0;QD=24.21;SB=-839.28 GT:AD:DP:GQ:PL 1/1:0,69:69:99:1703,207,0 AC=2;AF=1.00;AN=2;DP=82;Dels=0.00;FS=0.000;HRun=1;HaplotypeScore=10.8533;MQ=59.44;MQ0=0;QD=24.28;SB=-814.59 GT:AD:DP:GQ:PL 1/1:0,82:82:99:2024,246,0 NC_002516.2 837550 . GT G 2892.26 PASS AC=2;AF=1.00;AN=2;DP=72;FS=0.000;HRun=1;HaplotypeScore=278.7822;MQ=59.96;MQ0=0;QD=40.17;SB=-1333.27 GT:AD:DP:GQ:PL 1/1:0

837550 is an indel. The thing is, there have been indels before, so that does not look to be the problem. Thanks

vcf consensus next-gen
• 7,228 views • link
updated by Ram Ram • written by Nupur G • Nupur G

are your sure the assembly version your vcf file is based on and the assembly version of your fasta file are the same

• 123 views • link
written by Laura Laura

2 answers

Yes I am...I looked in the code, and it seems that it was a QC check that it was failing it. So I commented it. But not sure if that is a good thing.

• 0 views • link
written by Nupur G • Nupur G

I got the same error of The fasta sequence does not match the REF, and fixed this problem with dos2unix command

dos2unix NC_002516.2.fa
cat NC_002516.2.fa | vcf-consensus 1_e5_d30.vcf.gz > out.fa
• 0 views • link
updated by Ram Ram • written by liangjiao.xue liangjiao.xue

Log in to answer this question.

More posts like this
  • Clinical significance explanation
    written by abedkurdi10 • abedkurdi10
    Hi all, I have a question about CLNSIG field in clinvar database: 5 73981270 rs820878 T C 1777.96 PASS AC=2;AF=1;AN=2;DP=68;Dels=0;FS=0;HRun=0;HaplotypeScore=9.785;MQ=48.79;MQ0=0;QD=26.15;SB=-1169.48;TI=NM_000521;GI=HEXB;FC=Missense_L62S;EXON;CLNSIG=5,2;CLNDBN=Sandhoff_disease\x2c_infantile_type,not_specified;CLNORIGIN=1,1 GT:AD:DP:GQ:PL:VF:GQX 1/1:0,68:68:99:1811,132,0:1:99 5 73992881 rs10805890 …
  • GATK DP filter for multiple samples
    written by sun.nation sun.nation
    I used GATK to call SNPs and filter based on DP < 10. I did get SNPs but having problems to get filtered SNPs. Command …
  • Haploid data in VCFTools
    written by whiteeagle115 • whiteeagle115
    Hi everyone! I'm having problems on using haploid data (from VCF format) in VCFtools. First, I have a merged VCF with 9 samples and I …
  • Minor Allele Frequency calculation from VCF file
    written by a.james a.james
    Hello All I am using [GATK RNA-seq variant pipeline][1] for finding mutation/variant calling on the list of gene given in the following command line java-1.7 …
  • How to obtain chromosome number from given scaffold number
    written by Parimala Devi Parimala Devi
    Hi, I am working on *Saccharomyces cerevisiae*, Y55 strain. I obtained my reference sequence from [here][1]. And this is how the `Y55_Stanford_2014_JRIF00000000.fsa` sequence looks. It …
  • To remove genotype "./." from vcf file with awk
    written by ivivek_ngs ivivek_ngs
    Hi, I need help with awk commands, I have 4 samples in my vcf files, so field $10,$11,$12,$13 are the fields which have the genotype …
  • How Can I Make A Multi Alignment From The Vcf File Of Snp'S.
    written by satish.nair • satish.nair
    Hi Everyone, I have a VCF file of SNP,s that I called using GATK and now I need to make a multi alignment with this …
  • Missing Allele Balance Tag In Vcf File
    written by Johan Johan
    <p>I have a problem with my vcf files.</p> <p>I have called SNPs using the GATK UnifiedGenotyper with the AlleleBalance annotator. For some of the SNPs …
  • Unifiedgenotyper Bug?
    written by Leandro Batista • Leandro Batista
    <p>Hi everyone,</p> <p>I'm calling SNPs with UnifiedGenotyper:</p> <p><em>java -Xmx2g -jar /Volumes/DroboLeandro/Analysis/Software/GenomeAnalysisTK/GenomeAnalysisTK.jar -T UnifiedGenotyper -R /Volumes/DroboLeandro/Analysis/Genomes/Reference/NCBIM37.chr.fa -I BALBcByJ.recal.bam --dbsnp /Volumes/DroboLeandro/Analysis/dbSNP/dbSNP.vcf --genotype_likelihoods_model SNP --metrics_file BALBcByJ.metrics.snp.txt -stand_call_conf 30 …
  • Calling Of Mutil-Allelic Snp Using Gatk
    written by Bioscientist Bioscientist
    <p>I have four samples in a trio (and actually they are all patients) I tried using GATK-UnifiedGenotyper to call SNP/indel independently for each of them; …