Consensus sequence with characters apart from ATGC

Question

[sambioinfo2018]

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

[Bastien Hervé]

Is this base quality ?

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

[sambioinfo2018]

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.

Answer 1

[finswimmer]

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

[sambioinfo2018]

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!

[sambioinfo2018]

Or do you have the better way to do the same?

[finswimmer]

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

[sambioinfo2018]

Thanks Fin for your help!