This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Count mismatchs/Indel number per position in BAM file ?

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 ?

alignment
• 3,325 views • link
updated by Pierre Lindenbaum Pierre Lindenbaum • written by Chadi Saad Chadi Saad

Can't you provide non-overlapping intervals?

• 0 views • link
written by h.mon h.mon

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

• 0 views • link
written by Brian Bushnell Brian Bushnell

2 answers

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

pysamstats -f ref.fa --type variation_strand aln.bam > aln.var.txt
• 0 views • link
written by dariober dariober

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
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

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

• 0 views • link
written by Chadi Saad Chadi Saad

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

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Log in to answer this question.

More posts like this