Coverage for src / optwps / optwps.py: 96%
227 statements
« prev ^ index » next coverage.py v7.14.0, created at 2026-05-26 00:25 +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_batches,
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 bias_prior_count (float, optional): Pseudocount used to shrink rare-bin
212 bias weights.
213 Default: 20.0
214 bias_min_weight (float, optional): Minimum fragment bias weight. Use None
215 to disable.
216 Default: 0.2
217 bias_max_weight (float, optional): Maximum fragment bias weight. Use None
218 to disable.
219 Default: 5.0
220 valid_chroms (set, optional): Set of valid chromosome names to process.
221 Default: chromosomes 1-22, X, Y
222 chunk_size (float, optional): Region chunk size for processing.
223 Default: 1e8 (100 Mb)
224 read_buffer_size (int, optional): Number of reads sent to each worker task.
225 Default: 10000
226 njobs (int, optional): Number of jobs to use for read processing and input
227 BAM decompression. If negative, uses (number of CPUs + njobs).
228 Default: -1
230 Attributes:
231 bed_file (str): Path to BED file or None
232 protection_size (int): Half of the protection window size
233 valid_chroms (set): Set of valid chromosome names
234 min_insert_size (int): Minimum fragment size filter
235 max_insert_size (int): Maximum fragment size filter
236 min_mappability (float): Minimum mappability score filter
237 chunk_size (float): Chunk size for processing
238 roi_generator (ROIGenerator): Region generator instance
239 read_validator (ReadValidator): Read validator instance
240 weights_calculator (WeightsCalculator or None): Weights calculator instance for bias correction
242 Example:
243 >>> wps = WPS(protection_size=120, min_insert_size=120, max_insert_size=180)
244 >>> wps.run(bamfile='sample.bam', out_filepath='wps_output.tsv')
245 >>> # Process specific chromosomes
246 >>> wps = WPS(valid_chroms={'1', '2', '3', 'X', 'Y'})
247 >>> wps.run(bamfile='sample.bam', out_filepath='chr_wps.tsv')
249 Note:
250 - Automatically filters duplicate, QC-failed, and unmapped reads
251 - Handles both paired-end and single-end sequencing data
252 - Supports downsampling for high-coverage samples
253 - Can write to separate files per chromosome or target region using placeholders
254 """
256 def __init__(
257 self,
258 bed_file=None,
259 mappability_file=None,
260 protection_size=120,
261 min_insert_size=None,
262 max_insert_size=None,
263 min_mappability=0.9,
264 correct_for_bias=False,
265 bias_bins=10,
266 bias_subsample=0.05,
267 bias_prior_count=20.0,
268 bias_min_weight=0.2,
269 bias_max_weight=5.0,
270 valid_chroms=set(map(str, list(range(1, 23)) + ["X", "Y"])),
271 chunk_size=1e8,
272 njobs=1,
273 read_buffer_size=10000,
274 ):
275 self.bed_file = bed_file
276 self.mappability_file = mappability_file
277 self.protection_size = protection_size // 2
278 if valid_chroms is not None:
279 self.valid_chroms = [x.replace("chr", "") for x in valid_chroms]
280 else:
281 self.valid_chroms = None
282 self.min_insert_size = min_insert_size
283 self.max_insert_size = max_insert_size
284 self.min_mappability = min_mappability
285 self.chunk_size = chunk_size
286 self.njobs = njobs
287 if self.njobs < 0:
288 self.njobs = joblib.cpu_count() + self.njobs
289 self.njobs = max(1, self.njobs)
290 self.read_buffer_size = read_buffer_size
291 self.roi_generator = ROIGenerator(
292 bed_file=self.bed_file, chunk_size=self.chunk_size, njobs=self.njobs
293 )
294 self.read_validator = ReadValidator(
295 min_insert_size=self.min_insert_size,
296 max_insert_size=self.max_insert_size,
297 mappability_file=self.mappability_file,
298 min_mappability_threshold=self.min_mappability,
299 )
300 self.weights_calculator = (
301 WeightsCalculator(
302 mappability_file=self.mappability_file,
303 nbins=bias_bins,
304 subsample=bias_subsample,
305 min_insert_size=self.min_insert_size,
306 max_insert_size=self.max_insert_size,
307 min_mappability_threshold=self.min_mappability,
308 prior_count=bias_prior_count,
309 min_weight=bias_min_weight,
310 max_weight=bias_max_weight,
311 njobs=self.njobs,
312 read_buffer_size=self.read_buffer_size,
313 )
314 if correct_for_bias
315 else None
316 )
318 def __call__(self, *args, **kwargs):
319 return self.run(*args, **kwargs)
321 def run(
322 self,
323 bamfile,
324 out_filepath=None,
325 downsample_ratio=None,
326 compute_coverage=False,
327 verbose_output=False,
328 add_header=False,
329 ):
330 """
331 Calculate Window Protection Score for all regions and write to file.
333 Processes the BAM file to compute WPS values for each genomic position
334 in the specified regions (or entire genome). Results are written to a
335 tab-separated output file.
337 Args:
338 bamfile (str): Path to input BAM file (must be sorted and indexed)
339 out_filepath (str): Path to output TSV file, or stdout.
340 If it contains a formatting substring {target} or {chrom}, it will be used to create per-target or
341 per-chromosome files. Default: None (stdout)
342 downsample_ratio (float, optional): Fraction of reads to randomly keep
343 (0.0 to 1.0). Useful for high-coverage samples. Default: None (no downsampling)
344 compute_coverage (bool, optional): Whether to compute and include base coverage
345 verbose_output (bool, optional): Whether to include detailed counts
346 add_header (bool, optional): Whether to add header to the output
348 Returns:
349 None: Results are written directly to out_filepath
351 Output Format:
352 Tab-separated file with columns:
353 - chromosome: Chromosome name (without 'chr' prefix)
354 - start: Start position (0-based)
355 - end: End position (1-based, start + 1)
356 - base read coverage (if compute_coverage=True)
357 - outside: Count of fragments spanning the protection window (if verbose_output=True)
358 - inside: Count of fragment endpoints in protection window (if verbose_output=True)
359 - wps: Window Protection Score (outside - inside)
361 Raises:
362 FileNotFoundError: If bamfile does not exist
363 ValueError: If downsample_ratio is not between 0 and 1
365 Example:
366 >>> wps = WPS()
367 >>> wps.run(bamfile='input.bam', out_filepath='output.tsv')
368 >>> # With downsampling
369 >>> wps.run(bamfile='input.bam', out_filepath='output.tsv',
370 ... downsample_ratio=0.5)
371 >>> # Creating separate files per chromosome
372 >>> wps.run(bamfile='input.bam', out_filepath='wps_{chrom}.tsv')
373 >>> # Creating separate files per target region
374 >>> wps.run(bamfile='input.bam', out_filepath='wps_{target}.tsv.gz')
375 """
376 if out_filepath is None:
377 out_filepath = "stdout"
378 input_file = (
379 pysam.Samfile(bamfile, "rb")
380 if self.njobs == 1
381 else pysam.Samfile(bamfile, "rb", threads=self.njobs)
382 )
383 prefix = (
384 "chr" if any(r.startswith("chr") for r in input_file.references) else ""
385 )
386 use_partial_writer = "{target}" in out_filepath or "{chrom}" in out_filepath
387 partial_writers = dict()
388 total_outfile = None
389 header_added = set()
390 if not use_partial_writer:
391 total_outfile = CMWriterUnpacker(exopen(out_filepath, "w"))
392 try:
393 total_outfile = total_outfile.__enter__()
394 except AttributeError:
395 pass
397 if self.weights_calculator is not None:
398 self.weights_calculator.fit(input_file)
400 for chrom, start, end, region_id in self.roi_generator.regions(
401 bam_file=bamfile
402 ):
403 if "chr" in chrom:
404 chrom = chrom.replace("chr", "")
405 if self.valid_chroms is not None and chrom not in self.valid_chroms:
406 continue
407 try:
408 regionStart, regionEnd = int(start), int(end)
409 except ValueError:
410 continue
412 starts = []
413 ends = []
414 weights = []
415 reads = input_file.fetch(
416 prefix + chrom,
417 max(0, regionStart - self.protection_size - 1),
418 regionEnd + self.protection_size + 1,
419 )
420 can_parallel_reads = self.weights_calculator is None or isinstance(
421 self.weights_calculator, WeightsCalculator
422 )
423 if can_parallel_reads:
424 starts, ends, weights = collect_fragment_intervals(
425 read_info_batches(iter_pysam_reads(reads), self.read_buffer_size),
426 min_insert_size=self.min_insert_size,
427 max_insert_size=self.max_insert_size,
428 mappability_path=self.mappability_file,
429 min_mappability_threshold=self.min_mappability,
430 upstream_limit=regionStart - self.protection_size - 1,
431 downsample_ratio=downsample_ratio,
432 bin_edges=getattr(self.weights_calculator, "bin_edges", None),
433 weight_values=getattr(self.weights_calculator, "weights", None),
434 use_weights=self.weights_calculator is not None,
435 njobs=self.njobs,
436 )
437 else:
438 for read in iter_pysam_reads(reads):
439 if not self.read_validator.valid_read(
440 read,
441 upstream_limit=regionStart - self.protection_size - 1,
442 downsample_ratio=downsample_ratio,
443 ):
444 continue
446 rstart, rend = fragment_interval(read)
447 starts.append(rstart)
448 ends.append(rend)
449 if self.weights_calculator is not None:
450 weight = self.weights_calculator.transform(read)
451 weights.append(weight)
452 n = regionEnd - regionStart + 1
453 if len(starts) > 0:
454 starts = np.array(starts)
455 ends = np.array(ends)
456 read_weights = (
457 np.array(weights, dtype=float)
458 if self.weights_calculator is not None
459 else None
460 )
461 score_dtype = float if read_weights is not None else int
462 # Fragments fully spanning the window boundaries
463 span_start = starts + self.protection_size - regionStart
464 span_end = ends - self.protection_size - regionStart + 2
465 valid = span_end >= span_start
466 outside = np.zeros(n + 2, dtype=score_dtype)
467 np.add.at(
468 outside,
469 np.clip(span_start[valid] + 1, 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_end[valid], 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_start[~valid] + 1, 0, n + 1),
480 -1 if read_weights is None else -read_weights[~valid],
481 )
482 np.add.at(
483 outside,
484 np.clip(span_end[~valid], 0, n + 1),
485 1 if read_weights is None else read_weights[~valid],
486 )
488 outside_cum = np.cumsum(outside)[:-2]
490 # Fragments whose endpoints fall inside windows
491 all_ends = np.concatenate([starts, ends]) - regionStart
492 left = np.clip(all_ends - self.protection_size + 2, 0, n + 1)
493 right = np.clip(all_ends + self.protection_size + 1, 0, n + 1)
494 inside = np.zeros(n + 2, dtype=score_dtype)
495 endpoint_weights = (
496 1
497 if read_weights is None
498 else np.concatenate([read_weights, read_weights])
499 )
500 np.add.at(inside, left, endpoint_weights)
501 np.add.at(
502 inside,
503 right,
504 -1 if read_weights is None else -endpoint_weights,
505 )
506 inside_cum = np.cumsum(inside)[:-2]
508 wps = outside_cum - inside_cum
509 coverage = None
510 if compute_coverage:
511 coverage = np.zeros(n + 1, dtype=int)
512 np.add.at(
513 coverage,
514 np.clip(starts - regionStart, 0, n),
515 1,
516 )
517 np.add.at(
518 coverage,
519 np.clip(ends - regionStart + 1, 0, n),
520 -1,
521 )
522 coverage = np.cumsum(coverage)[:-1]
523 else:
524 score_dtype = float if self.weights_calculator is not None else int
525 outside_cum = np.zeros(n, dtype=score_dtype)
526 inside_cum = np.zeros(n, dtype=score_dtype)
527 wps = np.zeros(n, dtype=score_dtype)
528 coverage = None
529 if compute_coverage:
530 coverage = np.zeros(n, dtype=int)
531 partial_outfile = None
532 if use_partial_writer:
533 formatted_out_filepath = out_filepath.format(
534 target=region_id,
535 chrom=chrom,
536 )
537 partial_outfile = partial_writers.get(
538 formatted_out_filepath,
539 CMWriterUnpacker(exopen(formatted_out_filepath, "w")),
540 )
541 partial_writers[formatted_out_filepath] = partial_outfile
542 try:
543 partial_outfile = partial_outfile.__enter__()
544 except AttributeError:
545 pass
547 outfile = partial_outfile if use_partial_writer else total_outfile
548 st = np.arange(regionStart, regionEnd + 1)
549 en = st + 1 # add end coordinate
550 df = pd.DataFrame({"#chrom": chrom, "start": st, "end": en})
551 if compute_coverage:
552 df["coverage"] = coverage
553 if verbose_output:
554 df["outside"] = outside_cum
555 df["inside"] = inside_cum
556 df["wps"] = wps
557 header = False
558 if add_header and outfile not in header_added:
559 header_added.add(outfile)
560 header = True
561 df.to_csv(outfile, sep="\t", header=header, index=False)
563 if use_partial_writer:
564 for writer in partial_writers.values():
565 writer.close()
566 else:
567 total_outfile.close()
570def main():
571 """Main entry point for the command-line interface."""
572 from argparse import ArgumentParser
574 parser = ArgumentParser()
575 parser.add_argument(
576 "-i",
577 "--input",
578 dest="input",
579 help="Input BAM file",
580 required=True,
581 )
582 parser.add_argument(
583 "-o",
584 "--output",
585 dest="output",
586 help="The output file path for WPS results. If not provided, results will be printed to stdout.",
587 required=False,
588 )
589 parser.add_argument(
590 "-r",
591 "--regions",
592 dest="regions",
593 help="BED file with regions of interest (default: whole genome)",
594 default=None,
595 )
596 parser.add_argument(
597 "-w",
598 "--protection",
599 dest="protection",
600 help="Base pair protection window (default: 120)",
601 default=120,
602 type=int,
603 )
604 parser.add_argument(
605 "--min-insert-size",
606 dest="min_insert_size",
607 help="Minimum read length threshold to consider (Optional)",
608 default=None,
609 type=int,
610 )
611 parser.add_argument(
612 "--max-insert-size",
613 dest="max_insert_size",
614 help="Maximum read length threshold to consider (Optional)",
615 default=None,
616 type=int,
617 )
618 parser.add_argument(
619 "--mappability-file",
620 dest="mappability_file",
621 help="BigWig file with mappability scores used to filter fragments (optional)",
622 default=None,
623 )
624 parser.add_argument(
625 "--min-mappability",
626 dest="min_mappability",
627 help="Minimum average mappability score for fragments when --mappability-file is provided (default: 0.9)",
628 default=0.9,
629 type=float,
630 )
631 parser.add_argument(
632 "--correct-for-bias",
633 dest="correct_for_bias",
634 help="Apply fragment length and GC-content bias correction weights to WPS counts.",
635 action="store_true",
636 )
637 parser.add_argument(
638 "--bias-bins",
639 dest="bias_bins",
640 help="Number of bins per feature for bias-correction weights (default: 10)",
641 default=10,
642 type=int,
643 )
644 parser.add_argument(
645 "--bias-subsample",
646 dest="bias_subsample",
647 help="Fraction of reads used to estimate bias-correction weights (default: 0.05)",
648 default=0.05,
649 type=float,
650 )
651 parser.add_argument(
652 "--bias-prior-count",
653 dest="bias_prior_count",
654 help="Pseudocount used to shrink rare-bin bias weights (default: 20.0)",
655 default=20.0,
656 type=float,
657 )
658 parser.add_argument(
659 "--bias-min-weight",
660 dest="bias_min_weight",
661 help=(
662 "Minimum fragment bias weight. Use a negative value to disable "
663 "(default: 0.2)"
664 ),
665 default=0.2,
666 type=float,
667 )
668 parser.add_argument(
669 "--bias-max-weight",
670 dest="bias_max_weight",
671 help=(
672 "Maximum fragment bias weight. Use a negative value to disable "
673 "(default: 5.0)"
674 ),
675 default=5.0,
676 type=float,
677 )
678 parser.add_argument(
679 "--downsample",
680 dest="downsample",
681 help="Ratio to down sample reads (default OFF)",
682 default=None,
683 type=float,
684 )
685 parser.add_argument(
686 "--chunk-size",
687 dest="chunk_size",
688 help="Chunk size for processing in pieces, in case of low memory (default 1e8)",
689 default=1e8,
690 type=int,
691 )
692 parser.add_argument(
693 "--valid-chroms",
694 dest="valid_chroms",
695 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",
696 default=None,
697 type=str,
698 )
699 parser.add_argument(
700 "--compute-coverage",
701 dest="compute_coverage",
702 help="If provided, output will include base read coverage as the 4th column.",
703 action="store_true",
704 )
705 parser.add_argument(
706 "--verbose-output",
707 dest="verbose_output",
708 help="If provided, output will include separate counts for 'outside' and 'inside' along with WPS.",
709 action="store_true",
710 )
711 parser.add_argument(
712 "--add-header",
713 dest="add_header",
714 help="If provided, output files will include a header line.",
715 action="store_true",
716 )
717 parser.add_argument(
718 "--njobs",
719 dest="njobs",
720 help="Number of jobs to use for read processing and input BAM decompression. If negative, uses (number of CPUs + njobs). Default: 1",
721 default=1,
722 type=int,
723 )
724 parser.add_argument(
725 "--read-buffer-size",
726 dest="read_buffer_size",
727 help="Number of reads sent to each worker task when --njobs is not 1. Default: 10000",
728 default=10000,
729 type=int,
730 )
731 args = parser.parse_args()
732 valid_chroms = None
733 if args.valid_chroms == "canonical":
734 valid_chroms = [str(i) for i in range(1, 23)] + ["X", "Y"]
735 else:
736 valid_chroms = args.valid_chroms.split(",") if args.valid_chroms else None
737 optwps = WPS(
738 bed_file=args.regions,
739 protection_size=args.protection,
740 min_insert_size=args.min_insert_size,
741 max_insert_size=args.max_insert_size,
742 mappability_file=args.mappability_file,
743 min_mappability=args.min_mappability,
744 correct_for_bias=args.correct_for_bias,
745 bias_bins=args.bias_bins,
746 bias_subsample=args.bias_subsample,
747 bias_prior_count=args.bias_prior_count,
748 bias_min_weight=None if args.bias_min_weight < 0 else args.bias_min_weight,
749 bias_max_weight=None if args.bias_max_weight < 0 else args.bias_max_weight,
750 chunk_size=args.chunk_size,
751 valid_chroms=valid_chroms,
752 njobs=args.njobs,
753 read_buffer_size=args.read_buffer_size,
754 )
755 optwps.run(
756 bamfile=args.input,
757 out_filepath=args.output,
758 downsample_ratio=args.downsample,
759 compute_coverage=args.compute_coverage,
760 verbose_output=args.verbose_output,
761 add_header=args.add_header,
762 )
765if __name__ == "__main__":
766 main()