Count mismatchs/Indel number per position in BAM file ?

Question

[Chadi Saad]

how to count mismatchs number at each position in BAM file ? I found BAM-readcound, but it seems that the bug isn't solved yet : https://github.com/genome/bam-readcount/issues/19 any alternative method ?

[h.mon]

Can't you provide non-overlapping intervals?

[Brian Bushnell]

To clarify, are you talking about positions in the read, or in the reference?

Answer 1

[dariober]

I think pysamstats will do the job here with something like:

pysamstats -f ref.fa --type variation_strand aln.bam > aln.var.txt

Answer 2

[Pierre Lindenbaum]

my tool FindAllCoverageAtPosition do this job:

http://lindenb.github.io/jvarkit/FindAllCoverageAtPosition.html

$ find ./testdata/ -type f -name "*.bam" | \
 java -jar dist/findallcoverageatposition.jar -p rotavirus:100


#File              CHROM      POS  SAMPLE  DEPTH  M    I  D  N  S   H  P  EQ  X  Base(A)  Base(C)  Base(G)  Base(T)  Base(N)  Base(^)  Base(-)
./testdata/S4.bam  rotavirus  100  S4      126    126  0  0  0  29  0  0  0   0  5        0        0        121      0        0        0
./testdata/S1.bam  rotavirus  100  S1      317    317  1  0  0  50  0  0  0   0  27       0        1        289      0        1        0
./testdata/S2.bam  rotavirus  100  S2      311    311  0  1  0  60  0  0  0   0  29       1        0        281      0        0        1
./testdata/S3.bam  rotavirus  100  S3      446    446  1  0  0  86  0  0  0   0  39       0        1        406      0        1        0

or you can use sam2tsv: http://lindenb.github.io/jvarkit/Sam2Tsv.html

samtools view your.bam "chr1:1234-1234" | java -jar dist/sam2tsv.jar    -r  ref.fa

[Chadi Saad]

thank you ! but can I do it for all positions without -p option ?

[Pierre Lindenbaum]

at this point you should just parse the output of samtools mpileup....