bcftools consensus does not produce an output

Question

[ieie]

Hi everybody

I am trying to get a consensus sequences with bcftools with this commands:

bcftools mpileup -Ou -f reference.fa alignments_sorted.bam | bcftools call -mv -Oz -o calls.vcf.gz
tabix calls.vcf.gz
cat reference.fa | bcftools consensus calls.vcf.gz > consensus.fa

I then get this message:

The fasta sequence does not match the REF allele at MG564755.1:8368:
.vcf: [A]
.vcf: [AAATAT] <- (ALT)
.fa:  [T]CAAATCATCGACGGAAACGGAAAGAGAGGGATT

and no output. I saw a previous post where the same question was asked, but following the solution given I could not solve the problem. thanks a lot

[adrien.cv]

Hi everybody, I almost have the same problem as before … With the command line:

bcftools consensus -I -f refseq.fa input.vcf.gz > consensus.fa

I have only the bi-allelic SNP in IUPAC code in the consensus.fa and not the tri-allelic SNP or indel. Please, have you got an idea ? or an option in the command line ? thanks a lot

Answer 1

[Kevin Blighe]

A few options to try:

First try to set / reset the reference alleles in your VCF as per the reference FASTA that you're using:

bcftools norm --check-ref s --fasta-ref reference.fa -Ov calls.vcf.gz

Then run the bcftools consensus command

---------------------------------

Another option would be to skip the sites where there is discordance between the reference FASTA and the REF allele in your VCF, with something like:

bcftools consensus --fasta-ref reference.fa --exclude GT!="ref" calls.vcf.gz

------------------------------------------------

Could also try:

bcftools consensus --fasta-ref reference.fa --haplotype R calls.vcf.gz > out.fa
bcftools consensus --fasta-ref out.fa --haplotype A calls.vcf.gz

Kevin

[finswimmer]

And in addition the standard question: What version of bcftools are you using?

[ieie]

I am using bcftools 1.9.