I am trying to generate a consensus sequence from a sam file generated with bowtie2.
I tried samtools:
samtools mpileup -uf ref.fa aln.bam | bcftools view -cg - | vcfutils.pl vcf2fq > cns.fq
but when I visualise the sam file the consensus is calling the reference seq, not the sequenced base even though the coverage is ~ 20x.
Does anyone have any suggestions?
Thanks in advance,
Fiona
You want to use the "call" command of bcftools rather than the "view" command.
The command is:
bcftools call --output-type v
the "--output-type v" parameter is to output uncompressed VCF
Hope this helps!
Kind regards,
Lee
Thanks. I see that vcftools has been updated and I have since changed to call -c but still get the same issue with incorrect base calls. I don't think I was getting this issue when using an older version of samtools.
Oh - sorry - I misunderstood your question.
Is it happening for all SNPs or just individual SNPs? If it's happening for individual SNPS, then perhaps the issue is with some kind of default filter like minimum base quality or alignment quality.
The older versions of samtools and bcftools are still available for download at:
http://sourceforge.net/projects/samtools/files/samtools/
If you'd like to upload the bam file for a small region around one of the incorrectly called bases, I'd be happy to download it and take a look. I'd also need to know which reference file to use (if publicly available) or need the corresponding section of the reference file (if it's not publicly available)
You have to call SNP and Indels. Use freebayes or Varscan to call SNPs and Indels with parameters what you need. Then use bcftools consensus or vcf-consensus.
If you have gaps in your assembly, mask the regions with N.
https://samtools.github.io/bcftools/howtos/consensus-sequence.html
- Badges
- Science method