This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Consensus sequence with characters apart from ATGC

Dear All,

I have aligned the reads to the reference genome using Bowtie2. The generated sam was then converted to bam and sorted. The sorted file using the following command i tried to generate consensus sequence:

samtools mpileup -uf ref.fa aligned_sorted.bam | bcftools call -c | vcfutils.pl vcf2fq > aligned.fasta

However, a fasta file is generated. Viewing the tail of the file i get the follwoing:

HEHHHHHEHEEHEBBEBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
BBBBBBBEEEEEEHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHEHHHH
HHHHHHHHHHHHHHHEEEHHHEHHCCC@CCC@CEEEEEEEEEEECEEEEEEEEEEHHHHH
EHHHHHHHHHHHHHHHHHHHHHHHEHHKKKKKKKKKKKHEEEHECEEEEEEEECCEEEEE
EEEEEEEEEEEEEEEEHHHHHHHEKHHKKHKNNNNNQQQQNQQQNQNQQQQQQQTQTTQQ
TQQTQQQQQQNQQHNTTTTTTTTTTTQQQQNQQQQQQQQQQQQQQQQNQQQQQQQQQTQQ
QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQNNNNQQQQKKKKK
KKKKKKKHHHHHHHHHHHHHHHHEEEEEEEEEEEEEEEEEEEEEEEBBBBBBBBBBBBEE
EEEEEEEEEEEEEEEEEEEEEEEEEEEECCCCCCCCCFFFFFFFFFFFFFFFFFFFFFFF
FFFFFFFFFFFFFFFFFFFFFFCCCCCCCC

What does this mean. I expected the file to have charcaters 'ATGC'

Thanks in advance..

Regards, sam

next-gen
• 2,979 views • link
updated by finswimmer finswimmer • written by sambioinfo2018 • sambioinfo2018

Is this base quality ?

Your last command is vcfutils.pl vcf2fq, so you get a fastq file

• 1 views • link
written by Bastien Hervé Bastien Hervé

Thanks Bastien for your reply. But why is base quality printed in a consensus fasta file. I would like to use this fasta file for further analysis.

• 1 views • link
written by sambioinfo2018 • sambioinfo2018

1 answer

Hello,

if you like to have a consensu fasta you have to use bcftools consensus.

Furthermore which version of bcftools/samtools are you using. Since v1.9 samtools mpileup is deprecated and bcftools mpileup should be used instead, see:

Note that using "samtools mpileup" to generate BCF or VCF files is now
deprecated.  To output these formats, please use "bcftools mpileup" instead

fin swimmer

• 1 views • link
written by finswimmer finswimmer

Thanks finswimmer. i am using samtools v1.8 and bcftools v1.8. So should i use the following command to generate the consensus sequence:

bcftools mpileup -Ou -f reference.fa alignments.bam | bcftools call -mv -Oz -o calls.vcf.gz
tabix calls.vcf.gz

cat reference.fa | bcftools consensus calls.vcf.gz > consensus.fa

Thanks!

• 1 views • link
written by sambioinfo2018 • sambioinfo2018

Looks fine to me.

You could think about to output bcftools call to compressed bcf and indexing with bcftools index.

I prefer to use the -f option instead of piping the reference to bcftools consensus.

fin swimmer

• 1 views • link
written by finswimmer finswimmer

Log in to answer this question.

More posts like this