Question: bcftools consensus command
1
gravatar for Begonia_pavonina
5 weeks ago by
Begonia_pavonina10 wrote:

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.

A: 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?

ADD COMMENTlink modified 5 weeks ago by mhe0 • written 5 weeks ago by Begonia_pavonina10
1

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

ADD REPLYlink written 5 weeks ago by genomax78k

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 REPLYlink written 5 weeks ago by Begonia_pavonina10
1

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 REPLYlink written 5 weeks ago by genomax78k

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 REPLYlink modified 5 weeks ago by genomax78k • written 5 weeks ago by Begonia_pavonina10

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 REPLYlink modified 5 weeks ago • written 5 weeks ago by genomax78k

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 REPLYlink written 5 weeks ago by Begonia_pavonina10
1

I see. So things are working again?

ADD REPLYlink written 5 weeks ago by genomax78k

Yes it is now. Thanks for your help!

ADD REPLYlink written 4 weeks ago by Begonia_pavonina10

great tool - really

ADD REPLYlink written 5 weeks ago by mhe0
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 744 users visited in the last hour