I have a problem with bcftools roh and hope you can guide me to what may be causing this:
I run Example-1 of the Github: http://samtools.github.io/bcftools/howtos/roh-calling.html and see the output consists of a state (0:HW, 1:AZ) , and quality for each position that was in the original vcf file (test.vcf.gz).
I ham a bam file and like to get the exact same-format output for my list of positions. Here is how I generate the original .vcf (just in case that the way I am generating this is causing the issue that I explain later).
samtools mpileup -uf $FASTA ./myfile.bam -l ./dbSNP_chr9.bed | bcftools call -m -O z > roh.test.vcf.gz
This will provide a .vcf for 309 positions that I queried, and the header/rest of the vcf looks fine:
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="List of Phred-scaled genotype likelihoods"> ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype"> ... ##INFO=<ID=MQ,Number=1,Type=Integer,Description="Average mapping quality"> ##bcftools_callVersion=1.3.1+htslib-1.3.1 ##bcftools_callCommand=call -m -O z #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT mysample 9 421981 . G . 0 . DP=470;SGB=-0.379885;RPB=1;MQB=1;MQSB=1;BQB=1;MQ0F=0;AN=2;DP4=208,196,0,1;MQ=60 GT 0/0 9 422026 . T . 0 . DP=405;MQSB=1;MQ0F=0;AN=2;DP4=148,216,0,0;MQ=60 GT 0/0 9 3206113 . A . 0 . DP=338;MQSB=1;MQ0F=0;AN=2;DP4=219,79,0,0;MQ=60 GT 0/0 ...
But when I try the command in the example, I only receive two lines (positions) in the output file as also noted in the stdout!
bcftools roh --AF-dflt 0.4 ./roh.test.vcf.gz # This file was produced by: bcftools roh(1.3.1+htslib-1.3.1) # The command line was: bcftools roh --AF-dflt 0.4 ./roh.test.vcf.gz # # Chromosome Position State (0:HW, 1:AZ) Quality 9 5009541 1 3.0 9 5125250 1 3.0 Number of lines: total/processed: 309/2
My .vcf has PL but I also tried skipping the PL scores with -G30 (as is in the example), just in case it had something to do with this but the answer was exactly the same.
I greatly appreciate if you could point me to what may be causing this.