This is a test version of Biostars. For the public version, visit https://www.biostars.org.
using bcftools/samtools to generate consensus using reference + variants; output is a truncated reference

I am using samtools/bcftools to generate a 'de novo' consensus assembly for use with rc454 to clean my sequencing data. The outputted consensus sequence is identical to the reference (see the BLAST alignment graphic at the bottom of the post) till it is truncated at ~5000 of ~7000 bp.

I don't understand why it's happening or what to do about it.

I've used

  1. samtools mpileup to generate .vcf
  2. bcftools consensus to (try to) generate a consensus assembly.

According to this guy, the consensus function is buggy. Does anyone else have similar issues? Any ideas of what could be going wrong or a viable workaround? I have to process ~150 files w/ 10,000-60,000 reads of lengths 450-720 bp.​

Details of commands, options, and outputs

### command to generate Variant calls
samtools mpileup -vf mnv1.fa -r mnv1:5473-5922 APL000000034.34.sorted.bam > APL000000034.34.amp2.vcf

[fai_load] build FASTA index.
[mpileup] 1 samples in 1 input files
<mpileup> Set max per-file depth to 8000

### consensus command and error message
bcftools consensus -f mnv1.fa APL000000034.34.amp2.vcf > scratch.fa
Symbolic alleles other than <DEL> are currently not supported: <X> at mnv1:5474

### first three lines of purported "consensus" sequence
head -3 scratch.fa
>mnv1 gi|113478394|ref|NC_008311.1| Murine norovirus 1, complete genome
GTGAAATGAGGATGGCAACGCCATCTTCTGCGCCCTCTGTGCGCAACACAGAGAAACGCA
AAAACAAGAAGGCTTCGTCTAAAGCTAGTGTCTCCTTTGGAGCACCTAGCCCCCTCTCTT

### first three lines of reference sequence
head -3 mnv1.fa
>mnv1 gi|113478394|ref|NC_008311.1| Murine norovirus 1, complete genome
GTGAAATGAGGATGGCAACGCCATCTTCTGCGCCCTCTGTGCGCAACACAGAGAAACGCAAAAACAAGAA
GGCTTCGTCTAAAGCTAGTGTCTCCTTTGGAGCACCTAGCCCCCTCTCTTCGGAGAGCGAAGACGAAATT

< image not found >

samtools bcftools mosaik consensus rc454
• 8,645 views • link
updated by Ram Ram • written by nirvan • nirvan

Can you show the variant at mnv1:5474? I assume it's some sort of duplication or other large structural variant. The consensus command can't handle them, so you might need to filter them out. However, I think GATK has a similar tool and perhaps it can handle this better.

• 0 views • link
written by Devon Ryan Devon Ryan

Variant at mnv1:5474

grep 5474 APL000000034.34.amp2.vcf

mnv1    5474    .    A    <X>    0    .    DP=1321;I16=671,645,0,0,39501,1.24369e+06,0,0,44237,1.49124e+06,0,0,1358,1456,0,0;QS=1,0;MQSB=2.39844e-05;MQ0F=0    PL    0,255,255
• 0 views • link
updated by Ram Ram • written by nirvan • nirvan

Interesting, what version of bcftools is this? I would have expected the most recent ones to hand gVCFs.

• 0 views • link
written by Devon Ryan Devon Ryan

This is bcftools1.2. I don't believe the earlier versions have "consensus" function.

• 0 views • link
updated by Ram Ram • written by nirvan • nirvan

Ah, I'd forgotten that I'd written that, good catch.

• 0 views • link
written by Devon Ryan Devon Ryan

0 answers

No answers yet.

Log in to answer this question.

More posts like this