SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
bcftools view generates an empty output manore Bioinformatics 14 07-08-2016 02:19 AM
Varscan2 multi-variant samples VCF output with empty sample call FORMAT column? BhariD Bioinformatics 0 04-30-2015 01:56 PM
samtools/bcftools: How to inspect SNPs not made it into vcf file muenalan Bioinformatics 2 09-06-2012 05:13 AM
mpileup | bcftools 1000 genomes - empty vcf cdias Bioinformatics 2 09-02-2012 06:10 AM
Conflict between mpileup/bcftools and GATK in VCF file ericarcher Bioinformatics 0 09-25-2011 04:33 PM

Reply
 
Thread Tools
Old 05-04-2016, 11:24 AM   #1
AP38
Member
 
Location: Canada

Join Date: May 2016
Posts: 14
Default Empty VCF file with bcftools call

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
AP38 is offline   Reply With Quote
Old 05-04-2016, 11:27 AM   #2
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,783
Default

For reference cross-posted: https://www.biostars.org/p/189996
GenoMax is offline   Reply With Quote
Old 05-04-2016, 11:31 AM   #3
AP38
Member
 
Location: Canada

Join Date: May 2016
Posts: 14
Default

Thanks, yes I also asked the question on biostars as it may hit more people. If it's not appropriate to cross post, I will remove it from seqanswers.
AP38 is offline   Reply With Quote
Old 05-04-2016, 11:40 AM   #4
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,783
Default

Quote:
Originally Posted by AP38 View Post
Thanks, yes I also asked the question on biostars as it may hit more people. If it's not appropriate to cross post, I will remove it from seqanswers.
It is ok to cross-post on SeqAnswers. I included a link to your post on Biostars for reference.

If you get an answer over @Biostars then please come back and indicate that here.
GenoMax is offline   Reply With Quote
Old 05-04-2016, 08:23 PM   #5
SNPsaurus
Registered Vendor
 
Location: Eugene, OR

Join Date: May 2013
Posts: 450
Default

Shouldn't you have the -g option with mpileup to compute genotype likelihoods?
__________________
Providing nextRAD genotyping and PacBio sequencing services. http://snpsaurus.com
SNPsaurus is offline   Reply With Quote
Old 05-05-2016, 05:10 AM   #6
AP38
Member
 
Location: Canada

Join Date: May 2016
Posts: 14
Default

Yes indeed, if one wants to compute genotype likelihoods, the -g option is required. However, this should not help solve the problem and make a difference.
AP38 is offline   Reply With Quote
Old 05-05-2016, 10:25 AM   #7
SNPsaurus
Registered Vendor
 
Location: Eugene, OR

Join Date: May 2013
Posts: 450
Default

Sorry... don't know what I was thinking!

I tried your commands one a recent project and it also gave a header-only vcf file.

This minimal pipeline worked fine and produced a normal vcf from the same data:
mpileup -gu -Q 10 -t DP,DPR -f ref.fasta -b samples.txt | bcftools call -cv - > test.vcf
(it also worked without the -g!)

So from that I would conclude it is not a problem with your versions, file list or reference.

When I generate a bcf file using the minimal pipeline, it reports the reference allele. Your mpileup does not for some reason.
__________________
Providing nextRAD genotyping and PacBio sequencing services. http://snpsaurus.com

Last edited by SNPsaurus; 05-05-2016 at 10:35 AM.
SNPsaurus is offline   Reply With Quote
Old 05-05-2016, 10:50 AM   #8
AP38
Member
 
Location: Canada

Join Date: May 2016
Posts: 14
Default

Thanks! This is very good to know. It is possible that the -C 50 option is causing the issue because it downgrades mapping quality for excessive mismatches. I am working here with libraries of fairly small coverage so I might want to remove that option. I'll try that and stay in touch about the results.
AP38 is offline   Reply With Quote
Old 05-10-2016, 09:24 AM   #9
AP38
Member
 
Location: Canada

Join Date: May 2016
Posts: 14
Default

For some reasons, an empty line was added at the bottom of the index genome file. Removing it solved the problem…
AP38 is offline   Reply With Quote
Reply

Tags
bcftools, samtools

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 04:53 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2018, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO