Hi everyone!

I have a query regarding variant calling from a high coverage site on the basis of the maximum likelihood variant. I have RNA-seq data mapped bam file. I called variant using the below command.

bcftools mpileup --max-depth 10000  -Oz   -f ref.fa  sample.bam  | bcftools call -mv  -Oz -o calls.vcf.gz

What I found in vcf file variant is only called at maximum DP= 597 only however I can see the variant at high coverage but it's not present in vcf file.

For example the last line of VCF file.

Gene.13009.1    12364   .   C   T   108 .   DP=464;VDB=0.00126204;SGB=-0.693147;RPB=0.480808;MQB=1;MQSB=1;BQB=7.91699e-16;MQ0F=0;AC=2;AN=2;DP4=38,19,188,66;MQ=20   GT:PL   1/1:135,177,0
Gene.13009.1    12474   .   A   G   60  .   DP=597;VDB=0.0702122;SGB=-0.693147;RPB=0.748283;MQB=1;MQSB=1;BQB=2.38978e-19;MQ0F=0;AC=2;AN=2;DP4=51,31,205,56;MQ=20    GT:PL   1/1:88,43,0

Moreover, If I see in IGV variants are present at more than "DP <=597"

I don't know why I can't able to call that variant which is high depth coverage per site.

I have attached the screenshot of the variant which I am unable to call in vcf file.

Please kindly guide me or advice me.

Thank you!