Coverage for src / optwps / optwps.py: 85%
233 statements
« prev ^ index » next coverage.py v7.14.0, created at 2026-05-25 22:34 +0200
1#!/usr/bin/env python
2"""
3Fast Window Protection Score (WPS) Calculator for Cell-Free DNA Analysis.
5This module provides efficient computation of Window Protection Scores from BAM files,
6designed for analyzing cell-free DNA fragmentation patterns and nucleosome positioning.
8The WPS algorithm identifies protected genomic regions by analyzing DNA fragment
9coverage patterns. It counts fragments that span versus fragments whose endpoints
10fall within a protection window around each genomic position.
12Example:
13 Basic usage::
15 from optwps import WPS
17 wps_calc = WPS(protection_size=120)
18 wps_calc.run(bamfile='input.bam', out_filepath='output.tsv')
20 Creating separate output files::
22 # Per chromosome
23 wps_calc.run(bamfile='input.bam', out_filepath='wps_{chrom}.tsv')
25 # Per target region
26 wps_calc.run(bamfile='input.bam', out_filepath='wps_{target}.tsv.gz')
28Note:
29 - Op BAM Description
30 - M 0 alignment match (can be a sequence match or mismatch)
31 - I 1 insertion to the reference
32 - D 2 deletion from the reference
33 - N 3 skipped region from the reference
34 - S 4 soft clipping (clipped sequences present in SEQ)
35 - H 5 hard clipping (clipped sequences NOT present in SEQ)
36 - P 6 padding (silent deletion from padded reference)
37 - = 7 sequence match
38 - X 8 sequence mismatch
39 - pysam always uses 0-based coordinates, converting from 1-based of the sam files!
40"""
42import os
43import sys
44import pysam
45import numpy as np
46import pandas as pd
47import joblib
49from .read_processing import (
50 collect_fragment_intervals,
51 fragment_interval,
52 iter_pysam_reads,
53 read_info,
54)
55from .read_validator import ReadValidator
56from .weighting import WeightsCalculator
57from .utils import exopen
58from tqdm.auto import tqdm
61class ROIGenerator:
62 """
63 Generate regions of interest (ROI) for processing.
65 This class generates genomic regions either from a BED file or from the entire
66 genome in the BAM file. Regions are yielded in chunks for memory-efficient processing.
68 Args:
69 bed_file (str, optional): Path to BED file containing regions to process.
70 If None, the entire genome will be processed. Default: None
71 chunk_size (int, optional): Size of chunks in base pairs for processing.
72 Default: 1e8 (100 megabases)
73 """
75 def __init__(self, bed_file=None, chunk_size=1e8, njobs=1):
76 self.bed_file = bed_file
77 self.chunk_size = chunk_size
78 self.njobs = njobs
79 if self.njobs < 0:
80 self.njobs = joblib.cpu_count() + self.njobs
81 self.njobs = max(1, self.njobs)
83 def regions(self, bam_file=None):
84 """
85 Generate regions for processing.
87 Yields genomic regions either from a BED file or from the entire genome
88 referenced in the BAM file. Regions are chunked according to chunk_size.
90 Args:
91 bam_file (str, optional): Path to BAM file. Required if no BED file
92 is provided to determine genome regions. Default: None
94 Yields:
95 tuple: (chromosome, chunk_start, chunk_end, region_id)
96 where chunk_start/chunk_end define the current chunk being processed
97 and region_id is either from the BED file or constructed from chrom/start/end.
99 Raises:
100 ValueError: If neither bed_file nor bam_file can provide regions
101 """
102 if (self.bed_file is None) or (not os.path.exists(self.bed_file)):
103 input_file = (
104 pysam.Samfile(bam_file, "rb")
105 if self.njobs == 1
106 else pysam.Samfile(bam_file, "rb", threads=self.njobs)
107 )
108 nchunks = sum(
109 (input_file.get_reference_length(chrom) - 1) // self.chunk_size + 1
110 for chrom in input_file.references
111 )
112 iterator = tqdm(total=nchunks, desc="Processing genome regions")
113 for chrom in input_file.references:
114 chrom_length = input_file.get_reference_length(chrom)
115 region_start = 0
116 while region_start < chrom_length:
117 region_end = min(region_start + self.chunk_size, chrom_length)
118 yield chrom, region_start, region_end, chrom
119 region_start = region_end
120 iterator.update(1)
121 iterator.close()
122 else:
123 # read number of lines in bed file
124 nlines = sum(1 for _ in exopen(self.bed_file, "r"))
125 with exopen(self.bed_file, "r") as bed:
126 for line in tqdm(bed, total=nlines, desc="Processing BED regions"):
127 ret = line.strip().split("\t")
128 chrom, start, end = ret[:3]
129 chrom = chrom.replace("chr", "")
130 try:
131 region_id = ret[3]
132 except IndexError:
133 region_id = f"{chrom}_{start}_{end}"
134 chunk_start = int(start)
135 chunk_end = min(chunk_start + self.chunk_size, int(end))
136 while chunk_start < int(end):
137 yield chrom, chunk_start, chunk_end, region_id
138 chunk_start = chunk_end
139 chunk_end = min(chunk_start + self.chunk_size, int(end))
142class CMWriterUnpacker:
143 """
144 Wrapper for file handles to provide consistent write and close interface.
146 This class handles both context managers and direct file handles, providing
147 a unified interface for writing data and closing files. It also prevents
148 closing stdout to avoid unexpected behavior.
150 Args:
151 cm: A context manager or file handle object
153 Attributes:
154 cm: The original context manager
155 handle: The unpacked file handle
156 """
158 def __init__(self, cm):
159 self.cm = cm
160 try:
161 self.handle = self.cm.__enter__()
162 except AttributeError:
163 self.handle = self.cm
165 def write(self, data):
166 """Write data to the file handle."""
167 self.handle.write(data)
169 def close(self):
170 """Close the file handle, avoiding closure of stdout."""
171 try:
172 if self.handle != sys.stdout:
173 self.handle.close()
174 except AttributeError:
175 self.cm.__exit__(None, None, None)
178class WPS:
179 """
180 Window Protection Score (WPS) calculator for cell-free DNA analysis.
182 This class computes Window Protection Scores from aligned sequencing reads in BAM format.
183 WPS quantifies the protection of DNA fragments around each genomic position, useful for
184 identifying nucleosome positioning and other protected regions in cell-free DNA.
186 The algorithm calculates:
187 - Outside score: fragments that completely span the protection window
188 - Inside score: fragment endpoints falling within the protection window
189 - WPS = outside - inside
191 Args:
192 bed_file (str, optional): Path to BED file with regions to process.
193 If None, processes entire genome. Default: None
194 mappability_file (str, optional): Path to BigWig file with mappability scores.
195 If None, bam will be processed without mappability filtering.
196 protection_size (int, optional): Total protection window size in base pairs.
197 This value is divided by 2 to get the window on each side of the position.
198 Default: 120
199 min_insert_size (int, optional): Minimum insert/fragment size to include.
200 If None, no minimum filter applied. Default: None
201 max_insert_size (int, optional): Maximum insert/fragment size to include.
202 If None, no maximum filter applied. Default: None
203 min_mappability (float, optional): Minimum average mappability score for fragments to include.
204 If None, no mappability filter applied. Default: 0.9
205 correct_for_bias (bool, optional): Whether to apply bias correction weights based on fragment features (length, GC content).
206 Default: False
207 bias_bins (int, optional): Number of bins used for each bias-correction feature.
208 Default: 10
209 bias_subsample (float, optional): Fraction of reads used to estimate bias-correction weights.
210 Default: 0.05
211 valid_chroms (set, optional): Set of valid chromosome names to process.
212 Default: chromosomes 1-22, X, Y
213 chunk_size (float, optional): Region chunk size for processing.
214 Default: 1e8 (100 Mb)
215 read_buffer_size (int, optional): Number of reads sent to each worker task.
216 Default: 10000
217 njobs (int, optional): Number of jobs to use for read processing and input
218 BAM decompression. If negative, uses (number of CPUs + njobs).
219 Default: -1
221 Attributes:
222 bed_file (str): Path to BED file or None
223 protection_size (int): Half of the protection window size
224 valid_chroms (set): Set of valid chromosome names
225 min_insert_size (int): Minimum fragment size filter
226 max_insert_size (int): Maximum fragment size filter
227 min_mappability (float): Minimum mappability score filter
228 chunk_size (float): Chunk size for processing
229 roi_generator (ROIGenerator): Region generator instance
230 read_validator (ReadValidator): Read validator instance
231 weights_calculator (WeightsCalculator or None): Weights calculator instance for bias correction
233 Example:
234 >>> wps = WPS(protection_size=120, min_insert_size=120, max_insert_size=180)
235 >>> wps.run(bamfile='sample.bam', out_filepath='wps_output.tsv')
236 >>> # Process specific chromosomes
237 >>> wps = WPS(valid_chroms={'1', '2', '3', 'X', 'Y'})
238 >>> wps.run(bamfile='sample.bam', out_filepath='chr_wps.tsv')
240 Note:
241 - Automatically filters duplicate, QC-failed, and unmapped reads
242 - Handles both paired-end and single-end sequencing data
243 - Supports downsampling for high-coverage samples
244 - Can write to separate files per chromosome or target region using placeholders
245 """
247 def __init__(
248 self,
249 bed_file=None,
250 mappability_file=None,
251 protection_size=120,
252 min_insert_size=None,
253 max_insert_size=None,
254 min_mappability=0.9,
255 correct_for_bias=False,
256 bias_bins=10,
257 bias_subsample=0.05,
258 valid_chroms=set(map(str, list(range(1, 23)) + ["X", "Y"])),
259 chunk_size=1e8,
260 njobs=1,
261 read_buffer_size=10000,
262 ):
263 self.bed_file = bed_file
264 self.mappability_file = mappability_file
265 self.protection_size = protection_size // 2
266 if valid_chroms is not None:
267 self.valid_chroms = [x.replace("chr", "") for x in valid_chroms]
268 else:
269 self.valid_chroms = None
270 self.min_insert_size = min_insert_size
271 self.max_insert_size = max_insert_size
272 self.min_mappability = min_mappability
273 self.chunk_size = chunk_size
274 self.njobs = njobs
275 if self.njobs < 0:
276 self.njobs = joblib.cpu_count() + self.njobs
277 self.njobs = max(1, self.njobs)
278 self.read_buffer_size = read_buffer_size
279 self.roi_generator = ROIGenerator(
280 bed_file=self.bed_file, chunk_size=self.chunk_size, njobs=self.njobs
281 )
282 self.read_validator = ReadValidator(
283 min_insert_size=self.min_insert_size,
284 max_insert_size=self.max_insert_size,
285 mappability_file=self.mappability_file,
286 min_mappability_threshold=self.min_mappability,
287 )
288 self.weights_calculator = (
289 WeightsCalculator(
290 mappability_file=self.mappability_file,
291 nbins=bias_bins,
292 subsample=bias_subsample,
293 min_insert_size=self.min_insert_size,
294 max_insert_size=self.max_insert_size,
295 min_mappability_threshold=self.min_mappability,
296 njobs=self.njobs,
297 read_buffer_size=self.read_buffer_size,
298 )
299 if correct_for_bias
300 else None
301 )
303 def __call__(self, *args, **kwargs):
304 return self.run(*args, **kwargs)
306 def run(
307 self,
308 bamfile,
309 out_filepath=None,
310 downsample_ratio=None,
311 compute_coverage=False,
312 verbose_output=False,
313 add_header=False,
314 ):
315 """
316 Calculate Window Protection Score for all regions and write to file.
318 Processes the BAM file to compute WPS values for each genomic position
319 in the specified regions (or entire genome). Results are written to a
320 tab-separated output file.
322 Args:
323 bamfile (str): Path to input BAM file (must be sorted and indexed)
324 out_filepath (str): Path to output TSV file, or stdout.
325 If it contains a formatting substring {target} or {chrom}, it will be used to create per-target or
326 per-chromosome files. Default: None (stdout)
327 downsample_ratio (float, optional): Fraction of reads to randomly keep
328 (0.0 to 1.0). Useful for high-coverage samples. Default: None (no downsampling)
329 compute_coverage (bool, optional): Whether to compute and include base coverage
330 verbose_output (bool, optional): Whether to include detailed counts
331 add_header (bool, optional): Whether to add header to the output
333 Returns:
334 None: Results are written directly to out_filepath
336 Output Format:
337 Tab-separated file with columns:
338 - chromosome: Chromosome name (without 'chr' prefix)
339 - start: Start position (0-based)
340 - end: End position (1-based, start + 1)
341 - base read coverage (if compute_coverage=True)
342 - outside: Count of fragments spanning the protection window (if verbose_output=True)
343 - inside: Count of fragment endpoints in protection window (if verbose_output=True)
344 - wps: Window Protection Score (outside - inside)
346 Raises:
347 FileNotFoundError: If bamfile does not exist
348 ValueError: If downsample_ratio is not between 0 and 1
350 Example:
351 >>> wps = WPS()
352 >>> wps.run(bamfile='input.bam', out_filepath='output.tsv')
353 >>> # With downsampling
354 >>> wps.run(bamfile='input.bam', out_filepath='output.tsv',
355 ... downsample_ratio=0.5)
356 >>> # Creating separate files per chromosome
357 >>> wps.run(bamfile='input.bam', out_filepath='wps_{chrom}.tsv')
358 >>> # Creating separate files per target region
359 >>> wps.run(bamfile='input.bam', out_filepath='wps_{target}.tsv.gz')
360 """
361 if out_filepath is None:
362 out_filepath = "stdout"
363 input_file = (
364 pysam.Samfile(bamfile, "rb")
365 if self.njobs == 1
366 else pysam.Samfile(bamfile, "rb", threads=self.njobs)
367 )
368 prefix = (
369 "chr" if any(r.startswith("chr") for r in input_file.references) else ""
370 )
371 use_partial_writer = "{target}" in out_filepath or "{chrom}" in out_filepath
372 partial_writers = dict()
373 total_outfile = None
374 header_added = set()
375 if not use_partial_writer:
376 total_outfile = CMWriterUnpacker(exopen(out_filepath, "w"))
377 try:
378 total_outfile = total_outfile.__enter__()
379 except AttributeError:
380 pass
382 if self.weights_calculator is not None:
383 self.weights_calculator.fit(input_file)
385 for chrom, start, end, region_id in self.roi_generator.regions(
386 bam_file=bamfile
387 ):
388 if "chr" in chrom:
389 chrom = chrom.replace("chr", "")
390 if self.valid_chroms is not None and chrom not in self.valid_chroms:
391 continue
392 try:
393 regionStart, regionEnd = int(start), int(end)
394 except ValueError:
395 continue
397 starts = []
398 ends = []
399 weights = []
400 reads = input_file.fetch(
401 prefix + chrom,
402 max(0, regionStart - self.protection_size - 1),
403 regionEnd + self.protection_size + 1,
404 )
405 can_parallel_reads = self.weights_calculator is None or isinstance(
406 self.weights_calculator, WeightsCalculator
407 )
408 if can_parallel_reads:
409 read_batches = []
410 batch = []
411 for read in iter_pysam_reads(reads):
412 batch.append(read_info(read))
413 if len(batch) >= self.read_buffer_size:
414 read_batches.append(batch)
415 batch = []
416 if batch:
417 read_batches.append(batch)
418 starts, ends, weights = collect_fragment_intervals(
419 read_batches,
420 min_insert_size=self.min_insert_size,
421 max_insert_size=self.max_insert_size,
422 mappability_path=self.mappability_file,
423 min_mappability_threshold=self.min_mappability,
424 upstream_limit=regionStart - self.protection_size - 1,
425 downsample_ratio=downsample_ratio,
426 bin_edges=getattr(self.weights_calculator, "bin_edges", None),
427 weight_values=getattr(self.weights_calculator, "weights", None),
428 use_weights=self.weights_calculator is not None,
429 njobs=self.njobs,
430 read_buffer_size=self.read_buffer_size,
431 )
432 else:
433 for read in iter_pysam_reads(reads):
434 if not self.read_validator.valid_read(
435 read,
436 upstream_limit=regionStart - self.protection_size - 1,
437 downsample_ratio=downsample_ratio,
438 ):
439 continue
441 rstart, rend = fragment_interval(read)
442 starts.append(rstart)
443 ends.append(rend)
444 if self.weights_calculator is not None:
445 weight = self.weights_calculator.transform(read)
446 weights.append(weight)
447 n = regionEnd - regionStart + 1
448 if len(starts) > 0:
449 starts = np.array(starts)
450 ends = np.array(ends)
451 read_weights = (
452 np.array(weights, dtype=float)
453 if self.weights_calculator is not None
454 else None
455 )
456 score_dtype = float if read_weights is not None else int
457 # Fragments fully spanning the window boundaries
458 span_start = starts + self.protection_size - regionStart
459 span_end = ends - self.protection_size - regionStart + 2
460 valid = span_end >= span_start
461 outside = np.zeros(n + 2, dtype=score_dtype)
462 np.add.at(
463 outside,
464 np.clip(span_start[valid] + 1, 0, n + 1),
465 1 if read_weights is None else read_weights[valid],
466 )
467 np.add.at(
468 outside,
469 np.clip(span_end[valid], 0, n + 1),
470 -1 if read_weights is None else -read_weights[valid],
471 )
472 np.add.at(
473 outside,
474 np.clip(span_start[~valid] + 1, 0, n + 1),
475 -1 if read_weights is None else -read_weights[~valid],
476 )
477 np.add.at(
478 outside,
479 np.clip(span_end[~valid], 0, n + 1),
480 1 if read_weights is None else read_weights[~valid],
481 )
483 outside_cum = np.cumsum(outside)[:-2]
485 # Fragments whose endpoints fall inside windows
486 all_ends = np.concatenate([starts, ends]) - regionStart
487 left = np.clip(all_ends - self.protection_size + 2, 0, n + 1)
488 right = np.clip(all_ends + self.protection_size + 1, 0, n + 1)
489 inside = np.zeros(n + 2, dtype=score_dtype)
490 endpoint_weights = (
491 1
492 if read_weights is None
493 else np.concatenate([read_weights, read_weights])
494 )
495 np.add.at(inside, left, endpoint_weights)
496 np.add.at(
497 inside,
498 right,
499 -1 if read_weights is None else -endpoint_weights,
500 )
501 inside_cum = np.cumsum(inside)[:-2]
503 wps = outside_cum - inside_cum
504 coverage = None
505 if compute_coverage:
506 coverage = np.zeros(n + 1, dtype=int)
507 np.add.at(
508 coverage,
509 np.clip(starts - regionStart, 0, n),
510 1,
511 )
512 np.add.at(
513 coverage,
514 np.clip(ends - regionStart + 1, 0, n),
515 -1,
516 )
517 coverage = np.cumsum(coverage)[:-1]
518 else:
519 score_dtype = float if self.weights_calculator is not None else int
520 outside_cum = np.zeros(n, dtype=score_dtype)
521 inside_cum = np.zeros(n, dtype=score_dtype)
522 wps = np.zeros(n, dtype=score_dtype)
523 coverage = None
524 if compute_coverage:
525 coverage = np.zeros(n, dtype=int)
526 partial_outfile = None
527 if use_partial_writer:
528 formatted_out_filepath = out_filepath.format(
529 target=region_id,
530 chrom=chrom,
531 )
532 partial_outfile = partial_writers.get(
533 formatted_out_filepath,
534 CMWriterUnpacker(exopen(formatted_out_filepath, "w")),
535 )
536 partial_writers[formatted_out_filepath] = partial_outfile
537 try:
538 partial_outfile = partial_outfile.__enter__()
539 except AttributeError:
540 pass
542 outfile = partial_outfile if use_partial_writer else total_outfile
543 st = np.arange(regionStart, regionEnd + 1)
544 en = st + 1 # add end coordinate
545 df = pd.DataFrame({"#chrom": chrom, "start": st, "end": en})
546 if compute_coverage:
547 df["coverage"] = coverage
548 if verbose_output:
549 df["outside"] = outside_cum
550 df["inside"] = inside_cum
551 df["wps"] = wps
552 header = False
553 if add_header and outfile not in header_added:
554 header_added.add(outfile)
555 header = True
556 df.to_csv(outfile, sep="\t", header=header, index=False)
558 if use_partial_writer:
559 for writer in partial_writers.values():
560 writer.close()
561 else:
562 total_outfile.close()
565def main():
566 """Main entry point for the command-line interface."""
567 from argparse import ArgumentParser
569 parser = ArgumentParser()
570 parser.add_argument(
571 "-i",
572 "--input",
573 dest="input",
574 help="Input BAM file",
575 required=True,
576 )
577 parser.add_argument(
578 "-o",
579 "--output",
580 dest="output",
581 help="The output file path for WPS results. If not provided, results will be printed to stdout.",
582 required=False,
583 )
584 parser.add_argument(
585 "-r",
586 "--regions",
587 dest="regions",
588 help="BED file with regions of interest (default: whole genome)",
589 default=None,
590 )
591 parser.add_argument(
592 "-w",
593 "--protection",
594 dest="protection",
595 help="Base pair protection window (default: 120)",
596 default=120,
597 type=int,
598 )
599 parser.add_argument(
600 "--min-insert-size",
601 dest="min_insert_size",
602 help="Minimum read length threshold to consider (Optional)",
603 default=None,
604 type=int,
605 )
606 parser.add_argument(
607 "--max-insert-size",
608 dest="max_insert_size",
609 help="Maximum read length threshold to consider (Optional)",
610 default=None,
611 type=int,
612 )
613 parser.add_argument(
614 "--mappability-file",
615 dest="mappability_file",
616 help="BigWig file with mappability scores used to filter fragments (optional)",
617 default=None,
618 )
619 parser.add_argument(
620 "--min-mappability",
621 dest="min_mappability",
622 help="Minimum average mappability score for fragments when --mappability-file is provided (default: 0.9)",
623 default=0.9,
624 type=float,
625 )
626 parser.add_argument(
627 "--correct-for-bias",
628 dest="correct_for_bias",
629 help="Apply fragment length and GC-content bias correction weights to WPS counts.",
630 action="store_true",
631 )
632 parser.add_argument(
633 "--bias-bins",
634 dest="bias_bins",
635 help="Number of bins per feature for bias-correction weights (default: 10)",
636 default=10,
637 type=int,
638 )
639 parser.add_argument(
640 "--bias-subsample",
641 dest="bias_subsample",
642 help="Fraction of reads used to estimate bias-correction weights (default: 0.05)",
643 default=0.05,
644 type=float,
645 )
646 parser.add_argument(
647 "--downsample",
648 dest="downsample",
649 help="Ratio to down sample reads (default OFF)",
650 default=None,
651 type=float,
652 )
653 parser.add_argument(
654 "--chunk-size",
655 dest="chunk_size",
656 help="Chunk size for processing in pieces, in case of low memory (default 1e8)",
657 default=1e8,
658 type=int,
659 )
660 parser.add_argument(
661 "--valid-chroms",
662 dest="valid_chroms",
663 help="Comma-separated list of valid chromosomes to include (e.g., '1,2,3,X,Y') or 'canonical' for chromosomes 1-22, X, Y. Optional",
664 default=None,
665 type=str,
666 )
667 parser.add_argument(
668 "--compute-coverage",
669 dest="compute_coverage",
670 help="If provided, output will include base read coverage as the 4th column.",
671 action="store_true",
672 )
673 parser.add_argument(
674 "--verbose-output",
675 dest="verbose_output",
676 help="If provided, output will include separate counts for 'outside' and 'inside' along with WPS.",
677 action="store_true",
678 )
679 parser.add_argument(
680 "--add-header",
681 dest="add_header",
682 help="If provided, output files will include a header line.",
683 action="store_true",
684 )
685 parser.add_argument(
686 "--njobs",
687 dest="njobs",
688 help="Number of jobs to use for read processing and input BAM decompression. If negative, uses (number of CPUs + njobs). Default: 1",
689 default=1,
690 type=int,
691 )
692 parser.add_argument(
693 "--read-buffer-size",
694 dest="read_buffer_size",
695 help="Number of reads sent to each worker task when --njobs is not 1. Default: 10000",
696 default=10000,
697 type=int,
698 )
699 args = parser.parse_args()
700 valid_chroms = None
701 if args.valid_chroms == "canonical":
702 valid_chroms = [str(i) for i in range(1, 23)] + ["X", "Y"]
703 else:
704 valid_chroms = args.valid_chroms.split(",") if args.valid_chroms else None
705 optwps = WPS(
706 bed_file=args.regions,
707 protection_size=args.protection,
708 min_insert_size=args.min_insert_size,
709 max_insert_size=args.max_insert_size,
710 mappability_file=args.mappability_file,
711 min_mappability=args.min_mappability,
712 correct_for_bias=args.correct_for_bias,
713 bias_bins=args.bias_bins,
714 bias_subsample=args.bias_subsample,
715 chunk_size=args.chunk_size,
716 valid_chroms=valid_chroms,
717 njobs=args.njobs,
718 read_buffer_size=args.read_buffer_size,
719 )
720 optwps.run(
721 bamfile=args.input,
722 out_filepath=args.output,
723 downsample_ratio=args.downsample,
724 compute_coverage=args.compute_coverage,
725 verbose_output=args.verbose_output,
726 add_header=args.add_header,
727 )
730if __name__ == "__main__":
731 main()