Question: Count mismatchs/Indel number per position in BAM file ?
0
gravatar for Chadi Saad
3.4 years ago by
Chadi Saad80
France
Chadi Saad80 wrote:

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 • 1.3k views
ADD COMMENTlink modified 3.4 years ago by Pierre Lindenbaum131k • written 3.4 years ago by Chadi Saad80

Can't you provide non-overlapping intervals?

ADD REPLYlink modified 3.4 years ago • written 3.4 years ago by h.mon31k

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

ADD REPLYlink written 3.4 years ago by Brian Bushnell17k
2
gravatar for dariober
3.4 years ago by
dariober11k
WCIP | Glasgow | UK
dariober11k wrote:

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

pysamstats -f ref.fa --type variation_strand aln.bam > aln.var.txt
ADD COMMENTlink written 3.4 years ago by dariober11k
2
gravatar for Pierre Lindenbaum
3.4 years ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum131k wrote:

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
ADD COMMENTlink written 3.4 years ago by Pierre Lindenbaum131k

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

ADD REPLYlink written 3.4 years ago by Chadi Saad80

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

ADD REPLYlink written 3.4 years ago by Pierre Lindenbaum131k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1062 users visited in the last hour