Variant Calling With Bcftools Call
5
1
Entering edit mode
8.4 years ago
tayebwajb ▴ 110

I have been trying to call variants using bcftools call after installing it as instructed here. The problem is that I am not getting any calls in the output vcf (it contains only the header).

My bam files were aligned with bwa-mem and processed with GATK pipeline. They have worked well with muTect and VarScan. I have tried many parameter combinations e.g. -DgB , for mpileup and -mv for bcftools, in vain. The versions I am using are samtools 0.2.0-rc7-54-ge35361b and bcftools 0.2.0-rc7-13-gca1a8b5

Does anyone know what the problem might be?

An example commandline I used is below:

samtools mpileup -u -f ref.fa -b sampleList | bcftools call -cv - > output.vcf
bcftools samtools mpileup variant-calling • 19k views
ADD COMMENT
1
Entering edit mode

Without looking at things further, keep in mind that you're using a development version, so should presumably be reporting things back to the developers. You might get more stable results with the current release.

ADD REPLY
0
Entering edit mode

I think the problem is with bcftools since samtools mpileup looks fine to me. When I used the current release of bcftools, I got the error "invalid BCF2 magic string" which this link says was fixed in recent commits so that's why I tried the development version.

ADD REPLY
0
Entering edit mode

What happens if you drop -v? Could the simple explanation be that there are no non-reference calls?

ADD REPLY
0
Entering edit mode

I tried without -v too but that didn't help either. It would be helpful to know which versions of samtools, htslib and bcftools are compatible since I've had a hard time figuring this out. Thanks!

ADD REPLY
0
Entering edit mode

This does not seem to be a samtools/htslib/bcftools compatibility problem, as you say you followed instructions from here http://samtools.github.io/bcftools/. If you run mpileup with -uv, it will produce raw VCF output. Are any sites listed?

ADD REPLY
0
Entering edit mode

I used -uv on mpileup as you suggested and only -c on bcftools call and still got no sites listed.

ADD REPLY
0
Entering edit mode

I meant if you could run without bcftools at all like this: samtools mpileup -uv -f ref.fa -b sampleList | less. Are there any sites listed? If not, then there is something odd with your BAMs. Could all reads be excluded because of low quality for example? Is the result any different with -B,-x,-Q0?

ADD REPLY
0
Entering edit mode

Yes, I see sites listed with just mpileup. The problem seems to be with bcftools.

ADD REPLY
0
Entering edit mode

Can you please provide a sample of mpileup -uv output for testing?

ADD REPLY
0
Entering edit mode

I share tayebwajb issues and cannot have bcftools call to work on samtools mpileup data. I have samtools and bcftools/htslib fresh installed and also get the complain about BCF2.

samtools mpileup -u -f ${ref} ${bam} | bcftools call -c > raw_calls.bcf
[mpileup] 1 samples in 1 input files
<mpileup> Set max per-file depth to 8000
[vcf.c:643 bcf_hdr_read] invalid BCF2 magic string: only BCFv2.2 is supported.

while my samtools pileup seems OK

samtools mpileup -u -f ${ref} ${bam} | head
[mpileup] 1 samples in 1 input files
<mpileup> Set max per-file depth to 8000
?BC???BCF?chrMchr1chr2chr3chr4chr5chr6chr7chr8chr9chr10chr11chr12chr13chr14chr15chr16chr17chr18chr19chr20chr21chr22chrXchrYGAIIx-chr21-BWA.aln?##samtoolsVersion=0.1.19-44428cd
##reference=file://ref/HiSeq_UCSC_hg19.fa
##contig=<ID=chrM,length=16571>
##contig=<ID=chr1,length=249250621>
##contig=<ID=chr2,length=243199373>
##contig=<ID=chr3,length=198022430>
##contig=<ID=chr4,length=191154276>
##contig=<ID=chr5,length=180915260>
##contig=<ID=chr6,length=171115067>
##contig=<ID=chr7,length=159138663>

A full example command would be great with data we can download somewhere to test our executables. Any help most welcome

Stephane

ADD REPLY
0
Entering edit mode

Almost there (thanks pd3 for pointing the git branch error of mine)

Now I still have an issue building razip, likely a mac OSX maverick specific option missing in the makefile

make razip
gcc -g -Wall -O2 -D_FILE_OFFSET_BITS=64 -D_LARGEFILE64_SOURCE -D_CURSES_LIB=1 -I. -I../htslib   razip.c   -o razip
Undefined symbols for architecture x86_64:
  "_razf_close", referenced from:
      _main in ccPCU9N4.o
  "_razf_dopen", referenced from:
      _main in ccPCU9N4.o
  "_razf_open", referenced from:
      _main in ccPCU9N4.o
  "_razf_read", referenced from:
      _main in ccPCU9N4.o
  "_razf_seek", referenced from:
      _main in ccPCU9N4.o
  "_razf_write", referenced from:
      _main in ccPCU9N4.o
ld: symbol(s) not found for architecture x86_64
collect2: error: ld returned 1 exit status
make: *** [razip] Error 1

Thanks for help on this!

ADD REPLY
0
Entering edit mode

Hi Stephane,

I have the same problem with bcftools not recognizing mpileup output and reporting the same "invalid BCF2 magic string" error. I have updated all versions of samtools, bcftools and htslib but still the same error

Have you managed to get over this issue?

Thank you!

Sahar

ADD REPLY
2
Entering edit mode
8.3 years ago
pd3 ▴ 350

There is not enough information to answer the original question, but Stephane's problem is caused by a version mismatch. Note that samtools and bcftools went through major changes, the new versions of bcftools will not work with pre-HTSlib versions of samtools. For the latest installation instructions please see bcftools' wiki, the current commands (which are likely to change in near future) are:

git clone --branch=bcftools+calling git://github.com/samtools/htslib.git
git clone git://github.com/samtools/bcftools.git
git clone git://github.com/samtools/samtools.git
cd bcftools; make; cd ..
cd samtools; make;
ADD COMMENT
1
Entering edit mode
6.5 years ago

You need to use samtools 1.3 if you're using bcftools 1.3.

ADD COMMENT
0
Entering edit mode
8.4 years ago

How about using GATK's Unified Genotyper as an alternative to samtools mpileup? I have shared a tutorial here (C. Difficile 630 dataset that I have): http://userweb.eng.gla.ac.uk/umer.ijaz/bioinformatics/oneliners.html?#VARIANTCALLING

ADD COMMENT
0
Entering edit mode

Thanks Umer but the point IS to use a second caller and see what is shared with GATK.

ADD REPLY
0
Entering edit mode

My apologies, I might have misinterpreted while answering in a hurry!

Best Wishes,

Umer

ADD REPLY
0
Entering edit mode

no problem and thanks for sharing your methods ++

ADD REPLY
0
Entering edit mode

Hi,

Has this been resolved. I have done fresh installs of samtools and bcftools and get:

[vcf.c:753 bcf_hdr_read] invalid BCF2 magic string: only BCFv2.2 is supported

Its very frustrating. Also I cannot use GATK, I can't get Java 1.7 on my Debian system (yet) and it always gives me numerous java errors.

Thanks.

ADD REPLY
0
Entering edit mode

There have been two updates to samtools and bcftools since this post was made, so it's quite likely that the issue has been resolved. This is particularly so since tayebwajb was using a pre-release version.

ADD REPLY
0
Entering edit mode
6.5 years ago

Hi all!

After running:

samtools mpileup -uf ../samples/human_g1k_v37.fasta ../samples/sorted_sample.bam | bcftools call -c

I am also getting the same error:

[vcf.c:805 bcf_hdr_read] invalid BCF2 magic string: only BCFv2.2 is supported.
Failed to open -: could not parse header

Latest versions of tools (bcftools 1.3, samtools 0.1.19-96b5f2294a). Does anyone has any news on this one nowadays?

ADD COMMENT
0
Entering edit mode
6.5 years ago

That did the trick! Thanks!

ADD COMMENT
0
Entering edit mode

Hi Vladimir!

It seems you had exactly the same problem I have now. When I run bcftools with the call command to convert a bcf file into vcf I get the following error message:

[vcf.c:842 bcf_hdr_read] invalid BCF2 magic string: only BCFv2.2 is supported

Could you please tell me how you've managed to solve this issue.

Thanks!

ADD REPLY
0
Entering edit mode

You need to use samtools 1.3 if you're using bcftools 1.3. for samtools-0.1.18 and samtools-0.1.19, you should use the bcftools under the samtools build folder like "/share/apps/prog/samtools-0.1.19/bcftools/bcftools"

ADD REPLY

Login before adding your answer.

Traffic: 1058 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