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
•
written
by
ieie •
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
And in addition the standard question: What version of bcftools are you using?
Log in to answer this question.
Hi everybody, I almost have the same problem as before … With the command line:
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