bcftools mpileup with --region option output wrong consequence
0
0
Entering edit mode
23 months ago
SSK ▴ 10

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

bcftools mpileup SNP bam • 1.1k views
ADD COMMENT
0
Entering edit mode

Have you tried to create a .csi index file for your bam files via:

$ samtools index -c input.bam

?

ADD REPLY
0
Entering edit mode

Yes, I did. And created csi index was same directory with bam file. But bcftools didn't work.

ADD REPLY
0
Entering edit mode

Are you able to call variants after 520Mbp, if you specify in explicit region, where you can see a variant in igv?

ADD REPLY
0
Entering edit mode

Yes, I can.

ADD REPLY

Login before adding your answer.

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