Hi all. I'm trying calling variant with command below
bcftools mpileup -a AD,ADF,ADR -B -q 40 -Q 18 -C 50 -O u -r [some chr names] -f my_reference input1.bam input2.bam | bcftools call -vm -f GQ,GP -O u | bcftools filter -i "INFO/MQ>=40" -O z -o output.vcf.gz
As a result, I can get out.vcf.gz, but it has only SNPs in shorter region than full chromosome length.
Specifically, size of chr.1A is about 594Mbp, but out.vcf.gz contains SNPs up to about 520Mbp.
When I see input1.bam and input2.bam on igv, there are SNPs later than 520Mbp of chr.1A.
The other chromosomes are same as well as chr.1A
I suppose the cause of this is index of input1.bam and input2.bam.
Because input1.bam and input2.bam are too large to create bai index, but bcftools seems to use bai index even if I created csi index. In addition, when I use "samtools index input1.bam" I get error message "Region 536870869..536871019 cannot be stored in a bai index". Region of this message is almost the same position with maximum position of SNPs in out.vcf.gz.
How can I solve this problem ?
Thanks.
supplementary [code of bcftools mpileup] https://github.com/samtools/bcftools/blob/develop/mpileup.c
Have you tried to create a
.csi
index file for your bam files via:?
Yes, I did. And created csi index was same directory with bam file. But bcftools didn't work.
Are you able to call variants after 520Mbp, if you specify in explicit region, where you can see a variant in igv?
Yes, I can.