Variant Calling With Bcftools Call

Question

[tayebwajb]

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

[Devon Ryan]

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.

[tayebwajb]

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.

[pd3]

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

[tayebwajb]

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!

[pd3]

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?

[tayebwajb]

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

[pd3]

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?

[tayebwajb]

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

[pd3]

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

[stephane.plaisance]

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

[stephane.plaisance]

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!

[sahargel]

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

Answer 1

[pd3]

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;

Answer 2

[Devon Ryan]

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

Answer 3

[umer.zeeshan.ijaz]

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

[stephane.plaisance]

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

[umer.zeeshan.ijaz]

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

Best Wishes,

Umer

[stephane.plaisance]

no problem and thanks for sharing your methods ++

[susan.klein]

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.

[Devon Ryan]

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.

Answer 4

[vladimir.kovacevic]

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?

Answer 5

[vladimir.kovacevic]

That did the trick! Thanks!

[vasilipankratov]

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!

[hlzi]

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"