This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Proper way to call a set of BAMs with bcftools in '--gvcf' mode ?

(I know about gatk+gvcf)

i'm trying to call+merge about 1000 bams. My idea is to invoke bcftools mpileup+call with the --gvcf 1 param with 10 bams and then merge the 100 temporary bcf with bcftools merge

so for the one bed file, my command looks like;

bcftools mpileup --output-type u -o tmp.bcf --fasta-ref "ref.fasta"      --adjust-MQ 50 --min-MQ 20      --regions-file "000000043.bed" -a 'INFO/AD'  -a 'FORMAT/DP'  --gvcf "1" 1.bam 2.bam ... 10.bam


bcftools call  --gvcf "1"  --ploidy GRCh37 -f GQ --multiallelic-caller --output-type b -o "samples.bcf" tmp.bcf

I expected for each pool of 100 bcf that they would en with the same last INFO/END attribute. But it's not the case, which leads to some weird results when merging the data with bcftools merge.

Here I scan my newtflow workdir for a given bed file and print the last variant of each bcf file.

$ find /ccc/scratch/cont007/fg0156/lindenbp/20200407/work -type f -name ".command.sh" -exec grep -l "000000043.bed" '{}' ';' 2> /dev/null | head | sed 's/.command.sh/samples.bcf/' | while read F; do echo "#$F" && bcftools view -G --no-header ${F} | tail -1 ; done
#/ccc/scratch/cont007/fg0156/lindenbp/20200407/work/d2/21793ac97da462c1130263203fa150/samples.bcf
#/ccc/scratch/cont007/fg0156/lindenbp/20200407/work/d2/04e430f70b290e216dd8427740d3f9/samples.bcf
#/ccc/scratch/cont007/fg0156/lindenbp/20200407/work/d2/31fd1106cff09e31779bf587c42b35/samples.bcf
#/ccc/scratch/cont007/fg0156/lindenbp/20200407/work/82/0a14f7697d35bb30e00cd361125bc6/samples.bcf
#/ccc/scratch/cont007/fg0156/lindenbp/20200407/work/3c/6c2ea7ca7b95481433315e25393c63/samples.bcf
#/ccc/scratch/cont007/fg0156/lindenbp/20200407/work/3c/8b7b9dca964860f9d36fc2de1463d9/samples.bcf
GL000193.1      88464   .       C       .       .       .       END=88475;MinDP=20
#/ccc/scratch/cont007/fg0156/lindenbp/20200407/work/3c/7f7974666f68cabab90e145828e981/samples.bcf
#/ccc/scratch/cont007/fg0156/lindenbp/20200407/work/6d/c9dcfd74facb8b5b694068be9f4841/samples.bcf
GL000193.1      88464   .       C       .       .       .       END=88475;MinDP=1
#/ccc/scratch/cont007/fg0156/lindenbp/20200407/work/b4/8f0cf99a3240d0e1685dbad4007811/samples.bcf
#/ccc/scratch/cont007/fg0156/lindenbp/20200407/work/4d/92c85a2e5fb83a39108113f0eef70f/samples.bcf

some bcf don't have any variant at all (which was' not expected ?)

So my question is : what is the proper way to call some BAMs with gvcf and then merge the results ?

Thanks.

bcftools gvcf calling bam vcf

update: I found an empty bcf after 'bcftools mpileup' (so, it's not a problem with call) but mpileup.

my example is pretty bad because there is no read in the bed I choose. I'm deleting this question for now...

0 answers

No answers yet.

Log in to answer this question.