bcftools consensus command
0
1
Entering edit mode
4.1 years ago

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 • 8.6k views
ADD COMMENT
1
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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.

ADD REPLY
1
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

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!

ADD REPLY
1
Entering edit mode

I see. So things are working again?

ADD REPLY
0
Entering edit mode

Yes it is now. Thanks for your help!

ADD REPLY
0
Entering edit mode

great tool - really

ADD REPLY

Login before adding your answer.

Traffic: 1197 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6