How to call variant by --max-depth for RNAseq
0
0
Entering edit mode
3.1 years ago
Info.shi ▴ 30

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.

enter image description here

Please kindly guide me or advice me.

Thank you!

bcftools Samtools • 705 views
ADD COMMENT

Login before adding your answer.

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