This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Reality check for consensus fasta generated with bcftools

Hello,

I have generated a consensus seuqwnce with bcftools with:

bcftools mpileup -Ou -f ref.fa input_AlnSrtDedSrt.bam | bcftools call -Ou -mv | bcftools norm -f ref.fa -Oz -o input.vcf.gz
tabix input.vcf.gz
bcftools consensus -f ref.fa input.vcf.gz > out.fa

The file out.fa is there, but it is 100% identical to the ref.fa, even having the same header (I assumed bcftools got the header from ref.fa).

Is this procedure OK or I simply dumped ref.fa into out.fa?

reference-genome bcftools assembly fasta
• 3,780 views • link
updated by Ram Ram • written by marongiu.luigi marongiu.luigi

Can you show a few samples lines from the VCF? How are you sure that the files are identical?

• 0 views • link
written by Ram Ram

it is 100% identical to the ref.fa

And you have confirmed that by doing an alignment?

• 0 views • link
written by GenoMax GenoMax 
$ head ref.fa
>MN855675.1 Siphoviridae sp. isolate 61, complete genome
GAAACAATGGATATTACTGCTACAGGGACCCAAGGACGGGTAAAGAGTTTGGATTAGGCAGAGACAGGCG
AATCGCAATCACTGAAGCTATACAGGCCAACATTGAGTTATTTTCAGGACACAAACACAAGCCTCTGACA
GCGAGAATCAACAGTGATAATTCCGTTACGTTACATTCATGGCTTGATCGCTACGAAAAAATCCTGGCCA
GCAGAGGAATCAAGCAGAAGACACTCATAAATTACATGAGCAAAATTAAAGCAATAAGGAGGGGTCTGCC
TGATGCTCCACTTGAAGACATCACCACAAAAGAAATTGCGGCAATGCTCAATGGATACATAGACGAGGGC
AAGGCGGCGTCAGCCAAGTTAATCAGATCAACACTGAGCGATGCATTCCGAGAGGCAATAGCTGAAGGCC
ATATAACAACAAACCATGTCGCTGCCACTCGCGCAGCAAAATCAGAGGTAAGGAGATCAAGACTTACGGC
TGACGAATACCTGAAAATTTATCAAGCAGCAGAATCATCACCATGTTGGCTCAGACTTGCAATGGAACTG
GCTGTTGTTACCGGGCAACGAGTTGGTGATTTATGCGAAATGAAGTGGTCTGATATCGTAGATGGATATC
$ head out.fa 
>MN855675.1 Siphoviridae sp. isolate 61, complete genome
GAAACAATGGATATTACTGCTACAGGGACCCAAGGACGGGTAAAGAGTTTGGATTAGGCA
GAGACAGGCGAATCGCAATCACTGAAGCTATACAGGCCAACATTGAGTTATTTTCAGGAC
ACAAACACAAGCCTCTGACAGCGAGAATCAACAGTGATAATTCCGTTACGTTACATTCAT
GGCTTGATCGCTACGAAAAAATCCTGGCCAGCAGAGGAATCAAGCAGAAGACACTCATAA
ATTACATGAGCAAAATTAAAGCAATAAGGAGGGGTCTGCCTGATGCTCCACTTGAAGACA
TCACCACAAAAGAAATTGCGGCAATGCTCAATGGATACATAGACGAGGGCAAGGCGGCGT
CAGCCAAGTTAATCAGATCAACACTGAGCGATGCATTCCGAGAGGCAATAGCTGAAGGCC
ATATAACAACAAACCATGTCGCTGCCACTCGCGCAGCAAAATCAGAGGTAAGGAGATCAA
GACTTACGGCTGACGAATACCTGAAAATTTATCAAGCAGCAGAATCATCACCATGTTGGC

The alignment gives 100% identity. The only thing that is different is the line span...

• 0 views • link
written by marongiu.luigi marongiu.luigi

Do the contig names in the VCF match fasta headers in the ref file?

• 0 views • link
written by Ram Ram

Please show us a few sample lines from the VCF. zgrep -A5 "^#CHR" input.vcf.gz should suffice.

• 0 views • link
written by Ram Ram

Following, as I'm having the same issue of bcftools consensus generating an output file identical to my input reference sequence, without having applied any variants from the .vcf.gz generated by bcftools mpileup. I concatenated the output consenseus sequence with the reference into a multifasta and aligned to confirm they are exactly the same, including the headers. Confirmed that the .vcf.gz file contains expected contents.

For reference, here's my code:

bcftools mpileup --max-depth 500 -f -A $referenceSeq $sample.plastome.rmdup.bam | bcftools call -mv -Oz -o $sample.calls.vcf.gz
bcftools index $sample.calls.vcf.gz
bcftools consensus -f $referenceSeq $sample.calls.vcf.gz -o $sample.plastome.fasta
• 0 views • link
updated by Ram Ram • written by pwchan • pwchan

0 answers

No answers yet.

Log in to answer this question.

More posts like this