genome / bam-readcount

Count bases in BAM/CRAM files
MIT License
298 stars 95 forks source link

Inconsistent output #87

Open mnshgl0110 opened 2 years ago

mnshgl0110 commented 2 years ago

Hi,

This issue has been reported earlier https://github.com/genome/bam-readcount/issues/55, but there is no answer so I submit a new issue.

In short, the output is different when using bed file to select a position vs when all positions are used. Example:

18:17 goel@pc-t7-130 dep_mercier:pheno_mut$ bam-readcount -b 30 -q 10 -w 0 -f $refcur -l tmp.bed tmp.bam 
Minimum mapping quality is set to 10
CUR2G   9659721 T   528 =:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00  A:2:23.00:37.00:11.50:0:2:0.03:0.06:251.50:0:0.00:150.00:0.99   C:1:23.00:37.00:23.00:0:1:0.03:0.06:259.00:0:0.00:150.00:0.99   G:208:23.76:37.00:18.68:140:68:0.07:0.07:265.44:140:0.03:147.78:0.34    T:317:39.59:37.00:8.48:162:155:0.47:0.01:27.19:162:0.51:139.31:0.52 N:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00

18:17 goel@pc-t7-130 dep_mercier:pheno_mut$ bam-readcount -b 30 -q 10 -w 0 -f $refcur tmp.bam | grep 9659721
Minimum mapping quality is set to 10
CUR2G   9659721 T   410 =:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00  A:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00  C:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00  G:119:24.10:37.00:22.36:104:15:0.00:0.07:0.00:0:0.00:0.00:0.00  T:291:39.56:37.00:8.65:150:141:0.00:0.01:0.00:0:0.00:0.00:0.00  N:0:0.00:0.00:0.00:0:0:0.00:0.00:0.00:0:0.00:0.00:0.00

It will be really helpful if this can be resolved.

apldx commented 2 years ago

Which version of bam-readcount are you using? I'm mostly curious if it's the new 1.0.0 release, which switches to the 1.x series of htslib, and if the error occurs under old and new htslib.

Can you attach the BED file and a subset of the BAM? (say 500 flanking bases around the position of interest), or provide them for download somewhere? Then I can look into reproducing and troubleshooting this.

mnshgl0110 commented 2 years ago

The bam-readcount version is 0.8.0. Link to bam file: https://websafe.mpipz.mpg.de/d/1SDShnRT3r/tmp.bam

Link to BED file: https://websafe.mpipz.mpg.de/d/LSrYPjQLBh/TMP.bed