bcftools consensus does not produce an output
1
1
Entering edit mode
4.4 years ago
ieie ▴ 10

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

genome bcftools • 2.3k views
ADD COMMENT
0
Entering edit mode

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

ADD REPLY
3
Entering edit mode
4.4 years ago

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

ADD COMMENT
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

I am using bcftools 1.9.

ADD REPLY

Login before adding your answer.

Traffic: 2508 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6