Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Cast min_mapq to str instead of bytes #912

Merged
merged 2 commits into from
Sep 22, 2024

Conversation

rach-kennedy
Copy link
Contributor

Hi, I've been evaluating cnvkit for our data and overall liking it - results are looking good! However I've been noticing that the --min-mapq parameter doesn't seem to be propagating through correctly, and believe that this one-line change would be a fix.

To replicate the issue using data already present in this repo, you can pull out this single read with MAPQ=33 into a mini bam:

> samtools view -h test/formats/na12878-chrM-Y-trunc.bam | grep -e HWI-ST1124:106:C15APACXX:1:1107:19114:135147 -e ^@ | samtools view -hb > single_read_qual_33.bam && samtools index single_read_qual_33.bam
> echo -e "chrM\t7000\t8000" > single_target.bed

And then run cnvkit.py coverage with two different thresholds, which yield identical .cnn files despite falling above and below the MAPQ of this read:

> cnvkit.py coverage -q 0 single_read_qual_33.bam single_target.bed -o cov_qual_0.cnn
> cnvkit.py coverage -q 60 single_read_qual_33.bam single_target.bed -o cov_qual_60.cnn

And poking around more in the python, here's how using str instead of bytes is fixing the behavior:

>>> import pysam
>>> bed_fname = "single_target.bed"
>>> bam_fname = "single_read_qual_33.bam"
>>> pysam.bedcov(*[bed_fname, bam_fname, "-Q", bytes(60)])
'chrM\t7000\t8000\t200\n'
>>> pysam.bedcov(*[bed_fname, bam_fname, "-Q", str(60)])
'chrM\t7000\t8000\t0\n'

Thank you, and let me know if I've misunderstood something!

@etal
Copy link
Owner

etal commented Sep 22, 2024

Thanks for the fix! This makes sense.

@etal etal merged commit b873b86 into etal:master Sep 22, 2024
8 of 10 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants