bcftools consensus command

Question

[Begonia_pavonina]

Hello everyone, I try to use the bcftools consensus command, but it seems something is missing. I have designed my commands on the two sources below, and the manual.

Error while trying to get consensus Fastq from BAM https://github.com/samtools/bcftools/wiki/HOWTOs#consensus-calling

The first command tried:

bcftools consensus -s sample.vcf.gz -f reference.fa > output.fa

And additionally:

cat reference.fa | bcftools consensus -s sample.vcf.gz  > output.fa

But it does not seem to work in both cases. Would anyone have an idea of what could be wrong here?

[GenoMax]

Have you checked this tutorial by @finswimmer: Generating consensus sequence from bam file

[Begonia_pavonina]

It looks a very interesting pipeline, but I keep having problems using tabix on .vcf.gz.

It requires files zipped with bgzip, but even though does not process the files.

[GenoMax]

but I keep having problems using tabix on .vcf.gz.

What kind of problems? Do you see an error? Please post along with your exact commands.

[Begonia_pavonina]

I have not changed the command used in the tutorial first, then tried with/without different flags. Here is the original version.

bcftools mpileup -Ou -f ref.fa input.bam | bcftools call -Ou -mv | bcftools norm -f ref.fa -Oz -o output.vcf.gz

tabix file.vcf.gz 
>tbx_index_build failed: file.vcf.gz

bcftools index file.vcf.gz
>index: failed to open file.vcf.gz

The two indexing steps do not work.

[GenoMax]

You need to use your own file names. Names in the tutorial are for example. Surely your input BAM file is not called input.bam unless you renamed your original file to have that name.

What exactly did you call your output VCF files? Please use that name in your index command.

[Begonia_pavonina]

Hello genomax, thanks for your answer. I have renamed the files in my post up here for commodity purpose, as one would use user@host in place of the actual address of a server. They are indeed not labelled as such in my working space nor in my command.

The command is fine, I just find out the wrong module was loaded in the script including bcftools. Thanks for your help!

[GenoMax]

I see. So things are working again?

[Begonia_pavonina]

Yes it is now. Thanks for your help!

[mhe]

great tool - really