The .bcf file generated by bcftools mpileup cannot be opened with BCFtools

Question

[liying]

Hello everyone! I'm encountering an issue where both "bcftools view" and "bcftools call" are unable to open the .bcf file when using the following code. Does anyone have any suggestions or advice? Thank you!

My codes:

bcftools mpileup -O u -f ~/data/GRCz11/Danio_rerio.GRCz11.dna.toplevel.fa zebrafishAligned.sortedByCoord.out.bam > zebrafish_bcftools.bcf

bcftools call -mv -O v zebrafish_bcftools.bcf -o res.vcf

The error:

Failed to open zebrafish_bcftools.bcf: unknown file type.

The bcftools version is 1.9; however, I was able to use "bcftools view" after converting the MGP mouse VCF file to BCF format.

bcftools view test.vcf -O b -o A1.bcf.gz

[Pierre Lindenbaum]

what is the output of file zebrafish_bcftools.bcf ?

[liying]

zebrafish_bcftools.bcf: data; instead of Binary Call Format (BCF) version 2.2 (just like A1.bcf)

Answer 1

[Pierre Lindenbaum]

re-try your command, do not save the output of mpileup in an uncompressed bcf file, but pipe the output into bcftools call

[liying]

Thanks a lot ! I think I may have found the reason. I used bcftools mpileup -O u -f ~/data/GRCz11/Danio_rerio.GRCz11.dna.toplevel.fa zebrafishAligned.sortedByCoord.out.bam -o zebrafish2.bcf to do the variant calling, and I found "[mpileup] 1 samples in 1 input files" in the Log file; meanwhile, I tested bcftools mpileup -O v -f ~/data/GRCz11/Danio_rerio.GRCz11.dna.toplevel.fa zebrafishAligned.sortedByCoord.out.bam > zebrafish.vcf to generate VCF directly. And found that indeed the first sentence in the VCF is "[mpileup].....". It seems that the logging information was input into the BCF file, causing this issue.

[Pierre Lindenbaum]

I found "[mpileup] 1 samples in 1 input files" in the Log file; meanwhile,

hum... did you use something like nohup ?

[liying]

Yes. I used nohup bcftools ...... &

[Pierre Lindenbaum]

nohup is evil , it combines stderr+stdout, you'd better use tmux or screen

How to run R script in linux system

How can I use nohup command when filtering SNPs?

how to use nohup with parallel

etc...

---

[liying]

Thank you!!