Question: Empty VCF file with bcftools call
0
gravatar for AP
3.5 years ago by
AP90
AP90 wrote:

Hello,

I am trying to produce a vcf file using bcftools call but it produces an empty vcf file containing only the header. In short, here is what I do:

  1. Alignment with BWA
  2. With samtools, make sorted.bam files
  3. With samtools, index the sorted.bam files
  4. Run samtools mpileup in the following way: samtools mpileup -C 50 -E -t SP -t DP -u -I –f /genome/refgenome.fa -b bam_list.txt > output.bcf
  5. Run bcftools call: bcftools call -v -c output.bcf > output.vcf

I am using versions 1.3.1 of samtools, bcftools and htslib. I tried reinstalling these programs but it did not change the issue. I also tried with versions 1.2. Same problem. As far as I know, the bcf file seems fine, it contains lots of data and is 20GB.

I tried producing a basic vcf file using bcftools view: bcftools view output.cf > output.vcf and it works. The vcf file seems completely normal.

Could anyone help me with this? Why would bcftools call produce an empty output?

Thanks

bcftools call • 2.4k views
ADD COMMENTlink modified 3.4 years ago • written 3.5 years ago by AP90

How big is your vcf file? whats is the last line of your vcf file? Have you tried to run bcftools on Galaxy with your vcf file?

ADD REPLYlink modified 3.5 years ago • written 3.5 years ago by reza.jabal330

Thanks for your answer. My BAM files varie in size from 150MB to 2GB. My bcf file is ~ 22GB. I am running my analysis on a HPC and have plenty of power. The last line of my "empty" vcf file is just the end of the header with CHROM POS ID ...etc. I did not try to run bcftools on Galaxy but I am not sure that it will help. I suspect an issue with mpileup or bcftools call.

Please let me know if you have other suggestions.

ADD REPLYlink written 3.5 years ago by AP90

Ok, I guess the bam_list.txt might be the issue. If you're running your analysis on cluster make a flat file containing your bam locations. Then include the following to your pipeline:

cat BAM_Path | while read path sample; do
samtools mpileup -C 50 -E -t SP -t DP -u -I –f /genome/refgenome.fa -b "$location" > output.bcf...

Hope this helps!

ADD REPLYlink written 3.5 years ago by reza.jabal330
2
gravatar for AP
3.4 years ago by
AP90
AP90 wrote:

For some reasons, an empty line was added at the bottom of the index genome file. Removing it solved the problem…

ADD COMMENTlink written 3.4 years ago by AP90
1

well-done on sharing this.

ADD REPLYlink written 3.4 years ago by reza.jabal330
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: 1703 users visited in the last hour