This is a test version of Biostars. For the public version, visit https://www.biostars.org.
bcftools consensus command

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?

vcf bcftools
• 11,738 views • link
updated by Ram Ram • written by Begonia_pavonina • Begonia_pavonina

great tool - really

• 0 views • link
written by mhe • mhe

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.

• 0 views • link
written by Begonia_pavonina • Begonia_pavonina

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.

• 0 views • link
written by GenoMax GenoMax

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.

• 0 views • link
updated by GenoMax GenoMax • written by Begonia_pavonina • Begonia_pavonina

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.

• 0 views • link
written by GenoMax GenoMax

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!

• 0 views • link
written by Begonia_pavonina • Begonia_pavonina

I see. So things are working again?

• 0 views • link
written by GenoMax GenoMax

0 answers

No answers yet.

Log in to answer this question.

More posts like this