"""Cleavage Profiler.
Computes the cleavage profile of Zhou et al. (2022,
https://doi.org/10.1073/pnas.2209852119): the proportion of fragment ends at
each position over the depth at that position (as a percentage).
"""
from __future__ import annotations
import gzip
import time
import warnings
from multiprocessing import Pool
from pathlib import Path
from sys import stderr, stdin
from typing import Union
import numpy as np
import pyBigWig as pbw
from finaletoolkit.utils import (
chrom_sizes_to_dict,
chrom_sizes_to_list,
frag_array,
)
from finaletoolkit.utils.typing import FragFile
__all__ = ["cleavage_profile", "multi_cleavage_profile"]
# Structured dtype for per-position cleavage proportions.
_CLEAVAGE_DTYPE = [("contig", "U16"), ("pos", "i8"), ("proportion", "f8")]
[docs]
def cleavage_profile(
input_file: FragFile,
chrom_size: int,
contig: str,
start: int,
stop: int,
left: int = 0,
right: int = 0,
min_length: int | None = None,
max_length: int | None = None,
quality_threshold: int = 30,
verbose: Union[bool, int] = 0,
fraction_low: int | None = None,
fraction_high: int | None = None,
reference_file: str | Path | None = None,
) -> np.ndarray:
"""Compute the cleavage profile over a single interval.
Parameters
----------
input_file : str or pysam handle
BAM/CRAM/fragment input.
chrom_size : int
Length of ``contig``.
contig : str
Contig name.
start : int
0-based start coordinate.
stop : int
1-based stop coordinate.
left, right : int, optional
Amounts to expand the interval on each side (useful when only CpG
coordinates are given).
min_length, max_length : int, optional
Fragment-length filter.
quality_threshold : int, optional
Minimum mapping quality (default 30).
verbose : bool or int, optional
Print progress/timing.
fraction_low, fraction_high : int, optional
Deprecated aliases for ``min_length``/``max_length``.
reference_file : str or Path, optional
Reference genome (required for CRAM).
Returns
-------
numpy.ndarray
Structured array with fields ``('contig', 'pos', 'proportion')`` where
``proportion`` is a percentage in ``[0, 100]``.
"""
if verbose:
start_time = time.time()
stderr.write(
f"""
Calculating cleavage profile
input_file: {input_file}
contig: {contig}
start: {start}
stop: {stop}
fraction_low: {fraction_low}
fraction_high: {fraction_high}
quality_threshold: {quality_threshold}
verbose: {verbose}
"""
)
# Resolve deprecated aliases (both spellings together is an error).
if fraction_low is not None and min_length is None:
min_length = fraction_low
warnings.warn(
"fraction_low is deprecated. Use min_length instead.",
category=DeprecationWarning,
stacklevel=2,
)
elif fraction_low is not None and min_length is not None:
warnings.warn(
"fraction_low is deprecated. Use min_length instead.",
category=DeprecationWarning,
stacklevel=2,
)
raise ValueError("fraction_low and min_length cannot both be specified")
if fraction_high is not None and max_length is None:
max_length = fraction_high
warnings.warn(
"fraction_high is deprecated. Use max_length instead.",
category=DeprecationWarning,
stacklevel=2,
)
elif fraction_high is not None and max_length is not None:
warnings.warn(
"fraction_high is deprecated. Use max_length instead.",
category=DeprecationWarning,
stacklevel=2,
)
raise ValueError("fraction_high and max_length cannot both be specified")
adj_start = max(start - left, 0)
adj_stop = min(stop + right, chrom_size)
frags = frag_array(
input_file=input_file,
contig=contig,
quality_threshold=quality_threshold,
start=adj_start,
stop=adj_stop,
min_length=min_length,
max_length=max_length,
intersect_policy="any",
reference_file=reference_file,
)
positions = np.arange(adj_start, adj_stop)
# Depth: number of fragments covering each position.
fragwise_overlaps = np.logical_and(
np.greater_equal(positions[np.newaxis], frags["start"][:, np.newaxis]),
np.less(positions[np.newaxis], frags["stop"][:, np.newaxis]),
)
depth = np.sum(fragwise_overlaps, axis=0)
# Fragment ends: forward fragments end at their start (+ strand),
# reverse fragments end at their stop (- strand).
forward_ends = np.logical_and(
np.equal(positions[np.newaxis], frags["start"][:, np.newaxis]),
frags["strand"][:, np.newaxis],
)
reverse_ends = np.logical_and(
np.equal(positions[np.newaxis], frags["stop"][:, np.newaxis]),
np.logical_not(frags["strand"][:, np.newaxis]),
)
ends = np.sum(np.logical_or(forward_ends, reverse_ends), axis=0)
proportions = np.zeros_like(depth, dtype=np.float64)
non_zero_mask = depth != 0
proportions[non_zero_mask] = ends[non_zero_mask] / depth[non_zero_mask] * 100
results = np.zeros_like(proportions, dtype=_CLEAVAGE_DTYPE)
results["contig"] = contig
results["pos"] = positions
results["proportion"] = proportions
if verbose:
stderr.write(
f"cleavage_profile took {time.time() - start_time} s to complete\n"
)
return results
def _cleavage_profile_star(args):
return cleavage_profile(*args)
[docs]
def multi_cleavage_profile(
input_file: FragFile,
interval_file: Union[str, Path],
chrom_sizes: Union[str, Path],
left: int = 0,
right: int = 0,
min_length: int | None = None,
max_length: int | None = None,
quality_threshold: int = 30,
output_file: str = "-",
workers: int = 1,
verbose: Union[bool, int] = 0,
fraction_low: int | None = None,
fraction_high: int | None = None,
reference_file: str | Path | None = None,
) -> str:
"""Compute cleavage profiles over intervals in a (sorted) BED file.
Parameters
----------
input_file : str or path
BAM/CRAM/fragment input.
interval_file : str or path
Sorted BED of intervals (``"-"`` reads stdin).
chrom_sizes : str or path
``.chrom.sizes`` file (required).
left, right : int, optional
Interval expansion applied before merging overlaps.
min_length, max_length : int, optional
Fragment-length filter.
quality_threshold : int, optional
Minimum mapping quality (default 30).
output_file : str, optional
``.bw`` or ``.bed.gz``/`.bedgraph.gz` path (default ``"-"``).
workers : int, optional
Worker-process count (default 1).
verbose : bool or int, optional
Print progress/timing.
fraction_low, fraction_high : int, optional
Deprecated aliases for ``min_length``/``max_length``.
reference_file : str or Path, optional
Reference genome (required for CRAM).
Returns
-------
str
The output path.
"""
if verbose:
start_time = time.time()
stderr.write(
f"""
Calculating cleavage profile
input_file: {input_file}
interval_file: {interval_file}
chrom_sizes: {chrom_sizes}
left: {left}
right: {right}
min_length: {min_length}
max_length: {max_length}
quality_threshold: {quality_threshold}
output_file: {output_file}
workers: {workers}
verbose: {verbose}
"""
)
# Resolve deprecated aliases (both spellings together is an error).
if fraction_low is not None and min_length is None:
min_length = fraction_low
warnings.warn(
"fraction_low is deprecated. Use min_length instead.",
category=DeprecationWarning,
stacklevel=2,
)
elif fraction_low is not None and min_length is not None:
warnings.warn(
"fraction_low is deprecated. Use min_length instead.",
category=DeprecationWarning,
stacklevel=2,
)
raise ValueError("fraction_low and min_length cannot both be specified")
if fraction_high is not None and max_length is None:
max_length = fraction_high
warnings.warn(
"fraction_high is deprecated. Use max_length instead.",
category=DeprecationWarning,
stacklevel=2,
)
elif fraction_high is not None and max_length is not None:
warnings.warn(
"fraction_high is deprecated. Use max_length instead.",
category=DeprecationWarning,
stacklevel=2,
)
raise ValueError("fraction_high and max_length cannot both be specified")
if input_file == "-" and interval_file == "-":
raise ValueError("input_file and site_bed cannot both read from stdin")
if chrom_sizes is None:
raise ValueError("chrom_sizes must be specified.")
header = chrom_sizes_to_list(chrom_sizes)
chrom_dict = chrom_sizes_to_dict(chrom_sizes)
if verbose > 1:
stderr.write(f"chrom sizes {header}\n")
contigs, starts, stops = _read_intervals(
interval_file, left, right, chrom_dict
)
size_dict = dict(header)
sizes = [size_dict[contig] for contig in contigs]
count = len(contigs)
if verbose:
stderr.write("Zipping inputs\n")
interval_list = zip(
count * [input_file],
sizes,
contigs,
starts,
stops,
count * [0], # left/right precomputed to avoid double-padding
count * [0],
count * [min_length],
count * [max_length],
count * [quality_threshold],
count * [max(verbose - 1, 0)],
count * [fraction_low],
count * [fraction_high],
count * [reference_file],
)
if verbose:
stderr.write("Calculating cleavage profile...\n")
pool = Pool(workers, maxtasksperchild=500)
try:
interval_scores = pool.imap(
_cleavage_profile_star, interval_list, chunksize=100
)
if isinstance(output_file, str):
if verbose:
stderr.write(f"Output file {output_file} specified. Opening...\n")
if output_file.endswith(".bw"):
_write_bigwig(output_file, header, interval_scores)
elif (
output_file.endswith(".bed.gz")
or output_file.endswith("bedgraph.gz")
or output_file == "-"
):
_write_bedgraph_gz(output_file, interval_scores)
else:
raise ValueError(
"output_file can only have suffix .bw, .bedgraph.gz, or "
".bed.gz."
)
elif output_file is not None:
raise TypeError(
f'output_file is unsupported type "{type(input_file)}". '
"output_file should be a string specifying the path of the "
"file to output scores to."
)
finally:
pool.close()
if verbose:
end_time = time.time()
stderr.write(
f"cleavage profile took {end_time - start_time} s to complete\n"
)
return output_file
def _read_intervals(interval_file, left, right, chrom_dict):
"""Parse a sorted BED into merged, expanded intervals."""
contigs: list[str] = []
starts: list[int] = []
stops: list[int] = []
bed = stdin if interval_file == "-" else open(interval_file)
try:
prev_contig = None
prev_start = 0
prev_stop = 0
for line in bed:
contents = line.split()
contig = contents[0].strip()
start, stop = int(contents[1]), int(contents[2])
if contig not in chrom_dict:
warnings.warn(
f"Skipping interval {contig}:{start}-{stop} from "
f"interval_file ({contig} not in chrom_sizes)",
UserWarning,
)
continue
start = max(0, start - left)
stop = min(stop + right, chrom_dict[contig])
if prev_contig == contig and start < prev_stop:
prev_stop = max(prev_stop, stop)
else:
contigs.append(prev_contig)
starts.append(prev_start)
stops.append(prev_stop)
prev_contig, prev_start, prev_stop = contig, start, stop
contigs.append(prev_contig)
starts.append(prev_start)
stops.append(prev_stop)
finally:
if interval_file != "-":
bed.close()
# Drop the initial placeholder entry from the prev_* priming.
return contigs[1:], starts[1:], stops[1:]
def _write_bigwig(output_file, header, interval_scores) -> None:
"""Write per-position cleavage proportions to a bigWig file."""
with pbw.open(output_file, "w") as bigwig:
bigwig.addHeader(header)
last = "None"
for interval_score in interval_scores:
contigs = interval_score["contig"]
starts = interval_score["pos"]
scores = interval_score["proportion"]
if contigs.shape == (0,):
continue
try:
bigwig.addEntries(
contigs[0],
starts[0],
values=scores.astype(np.float64),
step=1,
span=1,
)
except RuntimeError as e:
stderr.write(f"{contigs[0]}:{starts[0]}-{starts[-1] + 1}\n")
stderr.write(
"invalid or out of order interval encountered. "
"Skipping to next.\n"
)
stderr.write(f"captured error:\n{e}\n")
stderr.write(f"current output:\n{interval_score}\n")
stderr.write(f"last output:\n{last}\n")
continue
last = interval_score
def _write_bedgraph_gz(output_file, interval_scores) -> None:
"""Write per-position cleavage proportions to a gzip-compressed bedGraph."""
with gzip.open(output_file, "wt") as bedgraph:
for interval_score in interval_scores:
contigs = interval_score["contig"]
starts = interval_score["pos"]
scores = interval_score["proportion"]
stops = starts + 1
lines = "".join(
f"{contig}\t{start}\t{stop}\t{score}\n"
for contig, start, stop, score in zip(contigs, starts, stops, scores)
)
bedgraph.write(lines)
