Merge branch 'subsample_qfilt' into 'dev'

subsample_qfilt: enable filtering for proportional subsampling and change...

See merge request research/pomoxis!163

CHANGELOG.md

@@ -4,6 +4,10 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [v0.3.12] - 2023-02-09
+### Changed
+- `subsample_bam`: `--quality` filtering now uses mean error probability, not mean of quality scores as previously.
+- `subsample_bam`: enable filtering for proportional subsampling.
 ## [v0.3.11] - 2022-11-16
 ### Fixed
 - Fix crashes in `subsample_bam` with alignment filtering and `common_errors_from_bam`

Makefile

@@ -139,7 +139,7 @@ $(BINCACHEDIR)/bedtools: | $(BINCACHEDIR) $(BINBUILDDIR)
 venv: venv/bin/activate
 IN_VENV=. ./venv/bin/activate
 venv/bin/activate:
-	test -d venv || $(PYTHON) -m venv venv --prompt '(pomoxis) '
+	test -d venv || $(PYTHON) -m venv venv --prompt 'pomoxis'
 	${IN_VENV} && pip install pip --upgrade
 	${IN_VENV} && pip install -r requirements.txt
 

pomoxis/__init__.py

@@ -1,4 +1,4 @@
-__version__ = '0.3.11'
+__version__ = '0.3.12'
 
 import argparse
 import os

pomoxis/subsample_bam.py

@@ -91,9 +91,6 @@ def main():
 
 
 def subsample_region_proportionally(region, args):
-    if args.quality is not None or args.coverage is not None or args.accuracy is not None:
-        raise NotImplemented('Read filtering is not currently supported for proportion subsampling')
-
     logger = logging.getLogger(region.ref_name)
     coverage_summary = coverage_summary_of_region(region, args.bam, args.stride)
     col = 'depth_{}'.format(args.orientation) if args.orientation is not None else 'depth'
@@ -139,6 +136,8 @@ def _read_iter():
     return found_enough_depth
 
 
+QSCORES_TO_PROBS = 10 ** (-0.1 * np.array(np.arange(100)))
+
 def filter_read(r, bam, args, logger):
     """Decide whether a read should be filtered out, returning a bool"""
 
@@ -153,9 +152,9 @@ def filter_read(r, bam, args, logger):
 
     # filter quality
     if args.quality is not None:
-        mean_q = np.mean(r.query_qualities)
+        mean_q = -10 * np.log10(np.mean(QSCORES_TO_PROBS[r.query_qualities]))
         if mean_q < args.quality:
-            logger.debug("Filtering {} by quality ({:.2f}).".format(r.query_name, mean_q))
+            logger.debug(f"Filtering {r.query_name} with quality {mean_q:.2f}")
             return True
 
     # filter accuracy or alignment coverage