from copy import deepcopy
from dataclasses import dataclass
from typing import Sequence, Literal
from .interval import Interval
from .genome_annotation import Transcript
from .genome import Genome
@dataclass(frozen=True)
class _CoordinateMetadata:
name: str | None
reference_genome: str
chromosome: str
transcript_strand: Literal["+", "-"]
@dataclass(frozen=True)
class _SegmentMetadata:
name: str | None
on_coordinate_strand: bool
[docs]
class DisjointIntervalSequence:
"""A flattened coordinate system over a sequence of disjoint genomic Intervals.
A DIS represents two layers:
- A **coordinate space** defined by a sequence of non-overlapping genomic
:py:class:`~genome_kit.Interval` objects (e.g. the exons of a transcript),
which are flattened into a contiguous 0-based index space. Index 0 is
always at the coordinate space's 5' end, regardless of genomic strand.
- A **segment** within that coordinate space, defined by a 5' and 3' index.
The segment may lie on the same, or opposite, strand as the coordinate space.
Use :py:meth:`from_transcript` or :py:meth:`from_intervals` to construct
instances rather than calling the constructor directly.
"""
[docs]
def __init__(
self,
coordinate_intervals: Sequence[Interval],
*,
coord_name: str | None = None,
segment_name: str | None = None,
on_coordinate_strand: bool = True,
start: int | None = None,
end: int | None = None,
):
"""Low-level constructor.
Prefer :py:meth:`from_transcript` or :py:meth:`from_intervals` for
public construction.
Parameters
----------
coordinate_intervals
Non-empty sequence of non-overlapping Intervals on the same
chromosome, strand, and reference genome.
coord_name
Optional name for the coordinate space.
segment_name
Optional name for the segment.
on_coordinate_strand
Whether the segment is on the same strand as the coordinate
intervals. Defaults to True. Can be used to represent a sequence that binds
to the transcript if set to False.
start
start index of the segment in the coordinate space. Defaults to 0
end
end index of the segment in the coordinate space. Defaults to the length
of the coordinate space.
Raises
------
ValueError
If coordinate intervals are empty, inconsistent, overlapping, or if start
is greater than end.
TypeError
If any element is not an Interval.
"""
if len(coordinate_intervals) == 0:
raise ValueError("coordinate_intervals must be non-empty")
for i, iv in enumerate(coordinate_intervals):
if not isinstance(iv, Interval):
raise TypeError(
f"coordinate_intervals[{i}] is {type(iv).__name__}, expected Interval"
)
if iv.anchor is not None:
raise ValueError(
f"coordinate_intervals[{i}] has an anchor set; "
f"anchored Intervals are not supported"
)
# Consistent chromosome, strand, reference_genome
iv0 = coordinate_intervals[0]
for iv in coordinate_intervals[1:]:
if iv.chromosome != iv0.chromosome:
raise ValueError(
f"All intervals must share the same chromosome, "
f"got {iv0.chromosome!r} and {iv.chromosome!r}"
)
if iv.strand != iv0.strand:
raise ValueError(
f"All intervals must share the same strand, "
f"got {iv0.strand!r} and {iv.strand!r}"
)
if iv.reference_genome != iv0.reference_genome:
raise ValueError(
f"All intervals must share the same reference genome, "
f"got {iv0.reference_genome!r} and {iv.reference_genome!r}"
)
# Sort 5'->3'
if iv0.strand == "+":
sorted_intervals = sorted(coordinate_intervals, key=lambda iv: iv.start)
else:
# On negative strand end is the 5' end since start < end.
# Sort by -end to get 5'->3' order.
sorted_intervals = sorted(coordinate_intervals, key=lambda iv: -iv.end)
# No overlaps (adjacent/touching OK)
for i in range(len(sorted_intervals) - 1):
cur_iv, next_iv = sorted_intervals[i], sorted_intervals[i + 1]
plus_strand_overlap = iv0.strand == "+" and cur_iv.end > next_iv.start
minus_strand_overlap = iv0.strand == "-" and cur_iv.start < next_iv.end
if plus_strand_overlap or minus_strand_overlap:
raise ValueError(
f"Intervals must not overlap: [{cur_iv.start}, {cur_iv.end}) and [{next_iv.start}, {next_iv.end})"
)
# Merge adjacent/touching intervals (e.g. [10, 20) and [20, 30) -> [10, 30)).
# Touching is detected on genomic coordinates (cur_iv.end == next_iv.start),
# which holds regardless of strand since start < end for every Interval.
merged_intervals: list[Interval] = [sorted_intervals[0]]
for next_iv in sorted_intervals[1:]:
last_iv = merged_intervals[-1]
if last_iv.end == next_iv.start or next_iv.end == last_iv.start:
merged_intervals[-1] = Interval(
iv0.chromosome,
iv0.strand,
min(last_iv.start, next_iv.start),
max(last_iv.end, next_iv.end),
iv0.reference_genome,
)
else:
merged_intervals.append(next_iv)
self._coordinate_intervals: tuple[Interval, ...] = tuple(merged_intervals)
self._coord_metadata = _CoordinateMetadata(
name=coord_name,
reference_genome=iv0.reference_genome,
chromosome=iv0.chromosome,
transcript_strand=iv0.strand,
)
self._segment_metadata = _SegmentMetadata(
name=segment_name,
on_coordinate_strand=on_coordinate_strand,
)
# Default segment start/end to span the full coordinate
if start is None:
start = 0
if end is None:
end = self.coordinate_length
# Validate that start is less than or equal to end
if start > end:
raise ValueError(
f"start index {start} cannot be greater than end index {end}"
)
self._start: int = start
self._end: int = end
[docs]
@classmethod
def from_intervals(
cls,
intervals: Sequence[Interval],
*,
coord_name: str | None = None,
segment_name: str | None = None,
) -> "DisjointIntervalSequence":
"""Construct a DIS from a sequence of Intervals
(or :py:class:`~genome_kit.Exon`/:py:class:`~genome_kit.Cds`/:py:class:`~genome_kit.Utr` objects).
Parameters
----------
intervals
Sequence of Interval or annotation objects.
coord_name
Optional name for the coordinate space.
segment_name
Optional name for the segment.
"""
coord_intervals = [
iv.interval if hasattr(iv, "interval") else iv for iv in intervals
]
return cls(coord_intervals, coord_name=coord_name, segment_name=segment_name)
[docs]
@classmethod
def from_transcript(
cls,
transcript: Transcript,
*,
region: Literal["exons", "cds", "utr5", "utr3"] = "exons",
coord_name: str | None = None,
segment_name: str | None = None,
) -> "DisjointIntervalSequence":
"""Construct a DIS from a transcript's exons, CDS, or UTR regions.
Parameters
----------
transcript
The source Transcript object.
region
Which region to extract — ``"exons"``, ``"cds"``,
``"utr5"``, or ``"utr3"``.
coord_name
Optional name for the coordinate space. Defaults to ``transcript.id``.
segment_name
Optional name for the segment. Defaults to ``transcript.id``.
Raises
------
ValueError
If region is not one of the allowed values.
"""
match region:
case "exons":
region_elements = transcript.exons
case "cds":
region_elements = transcript.cdss
case "utr5":
region_elements = transcript.utr5s
case "utr3":
region_elements = transcript.utr3s
case _:
raise ValueError(f"Invalid region: {region!r}")
coord_intervals = [element.interval for element in region_elements]
if coord_name is None:
coord_name = transcript.id
if segment_name is None:
segment_name = transcript.id
return cls(coord_intervals, coord_name=coord_name, segment_name=segment_name)
@property
def coord_name(self) -> str | None:
"""Name of the coordinate space, or None."""
return self._coord_metadata.name
@property
def reference_genome(self) -> str:
"""Reference genome of the coordinate intervals."""
return self._coord_metadata.reference_genome
@property
def chromosome(self) -> str:
"""Chromosome of the coordinate intervals."""
return self._coord_metadata.chromosome
@property
def coord_strand(self) -> Literal["+", "-"]:
"""Strand of the coordinate intervals (the transcript strand)."""
return self._coord_metadata.transcript_strand
@property
def name(self) -> str | None:
"""Name of the segment, or None."""
return self._segment_metadata.name
@property
def on_coordinate_strand(self) -> bool:
"""True if the segment is on the same strand as the coordinate intervals."""
return self._segment_metadata.on_coordinate_strand
@property
def strand(self) -> Literal["+", "-"]:
"""Effective strand of the segment, accounting for on_coordinate_strand."""
if self.on_coordinate_strand:
return self.coord_strand
# Segment is on opposite strand
if self.coord_strand == "+":
return "-"
return "+"
@property
def end5_index(self) -> int:
"""5' index of the segment."""
if self._upstream_index_step() == -1:
return self._start
return self._end
@property
def end3_index(self) -> int:
"""3' index of the segment."""
if self._upstream_index_step() == -1:
return self._end
return self._start
@property
def start(self) -> int:
"""Start index of the segment in the coordinate space."""
return self._start
@property
def end(self) -> int:
"""End index of the segment in the coordinate space."""
return self._end
def _at_index(
self, idx: int, on_coordinate_strand: bool
) -> "DisjointIntervalSequence":
"""Return a 0-length DIS at the given index position."""
return DisjointIntervalSequence(
self._coordinate_intervals,
coord_name=self._coord_metadata.name,
on_coordinate_strand=on_coordinate_strand,
start=idx,
end=idx,
)
@property
def end5(self) -> "DisjointIntervalSequence":
"""0-length DIS at the segment's 5' end."""
return self._at_index(
self.end5_index, on_coordinate_strand=self.on_coordinate_strand
)
@property
def end3(self) -> "DisjointIntervalSequence":
"""0-length DIS at the segment's 3' end."""
return self._at_index(
self.end3_index, on_coordinate_strand=self.on_coordinate_strand
)
@property
def coord_end5(self) -> "DisjointIntervalSequence":
"""0-length DIS at the coordinate space's 5' end."""
return self._at_index(0, on_coordinate_strand=True)
@property
def coord_end3(self) -> "DisjointIntervalSequence":
"""0-length DIS at the coordinate space's 3' end."""
return self._at_index(self.coordinate_length, on_coordinate_strand=True)
@property
def coordinate_intervals(self) -> tuple[Interval, ...]:
"""The underlying genomic intervals of the coordinate-space, sorted 5'->3'."""
# Deepcopy to preserve imutability of this DIS
return deepcopy(self._coordinate_intervals)
@property
def coordinate_length(self) -> int:
"""Total length of the coordinate space in bases."""
return sum(len(iv) for iv in self._coordinate_intervals)
@property
def length(self) -> int:
"""Length of the segment on the coordinate space."""
return self.end - self.start
def _lower_coord(self, coord: int) -> list[int]:
"""Convert a DIS coordinate index to the genomic coordinate(s).
Returns a 1-element list for indices interior to a single coord
interval, at the outer 5' (``coord == 0``) and 3'
(``coord == coordinate_length``) edges of the DIS, and for any
extrapolated index (``coord < 0`` or ``coord > coordinate_length``).
Extrapolation is linear from the nearest interval boundary and may
yield negative values or values exceeding the chromosome length.
Returns a 2-element list ``[iv_upstream_coord, iv_downstream_coord]`` sorted
5' -> 3' when ``coord`` falls exactly on an internal boundary between two
adjacent coord intervals — where ``iv_upstream_coord`` is the upstream-most
coordinate boundary, and ``iv_downstream_coord`` is the downstream-most
coordinate boundary (regardless of strand). Touching intervals produce
two equal values.
Parameters
----------
coord
Index in the DIS coordinate space.
"""
ivs = self._coordinate_intervals
on_plus = self.coord_strand == "+"
coord_len = self.coordinate_length
# Upstream of first interval — extrapolate (no upstream neighbor)
if coord <= 0:
if on_plus:
return [ivs[0].start - abs(coord)]
return [ivs[0].end + abs(coord)]
# Downstream of last interval — extrapolate (no downstream neighbor)
if coord >= coord_len:
overshoot = coord - coord_len
if on_plus:
return [ivs[-1].end + overshoot]
return [ivs[-1].start - overshoot]
# Walk intervals, consuming coord with a decrementing delta
delta = coord
interval_index = -1
while delta >= 0:
interval_index += 1
delta -= len(ivs[interval_index])
iv = ivs[interval_index]
if delta == -len(iv):
# coord landed on the cumulative end of ivs[interval_index - 1],
# i.e. the internal boundary between two coord intervals.
upstream = ivs[interval_index - 1]
if on_plus:
return [upstream.end, iv.start] # indices ordered 5' -> 3' w.r.t segment strand
else:
return [upstream.start, iv.end] # indices ordered 5' -> 3' w.r.t segment strand
if on_plus:
return [iv.end - abs(delta)]
return [iv.start + abs(delta)]
def _upstream_index_step(self, on_coordinate_strand: bool | None = None) -> int:
"""Return +1 or -1 indicating the upstream direction in index space.
Parameters
----------
on_coordinate_strand
Override for which strand to compute the step for. Defaults to
this segment's ``on_coordinate_strand``.
"""
if on_coordinate_strand is None:
on_coordinate_strand = self.on_coordinate_strand
return -1 if on_coordinate_strand else 1
def _validate_same_coordinate_space(
self, other: "DisjointIntervalSequence"
) -> None:
"""Raise if other does not share the same coordinate space."""
if not isinstance(other, DisjointIntervalSequence):
raise TypeError(
f"Expected DisjointIntervalSequence, got {type(other).__name__}"
)
if self._coordinate_intervals != other._coordinate_intervals:
raise ValueError("DIS objects must share the same coordinate intervals")
def _from_end_indices(self, end5: int, end3: int) -> "DisjointIntervalSequence":
"""Return a new DIS with the same coordinate space but different segment indices."""
# Validate end5 is upstream of or equal to end3
if self._upstream_index_step() == -1:
if end5 > end3:
raise ValueError(
f"Invalid indices: end5 index {end5} is downstream of end3 index {end3}"
)
if self._upstream_index_step() == 1:
if end5 < end3:
raise ValueError(
f"Invalid indices: end5 index {end5} is downstream of end3 index {end3}"
)
return DisjointIntervalSequence(
self._coordinate_intervals,
coord_name=self._coord_metadata.name,
segment_name=self._segment_metadata.name,
on_coordinate_strand=self.on_coordinate_strand,
start=min(end5, end3),
end=max(end5, end3),
)
[docs]
def shift(self, amount: int) -> "DisjointIntervalSequence":
"""Shift the segment downstream by amount (negative shifts upstream).
The coordinate space is unchanged. Only the segment indices move.
Parameters
----------
amount
Bases to shift the segment downstream. Negative values shift upstream.
"""
downstream_step = -self._upstream_index_step()
delta = amount * downstream_step
return self._from_end_indices(
self.end5_index + delta,
self.end3_index + delta,
)
[docs]
def expand(
self, upstream: int, dnstream: int | None = None
) -> "DisjointIntervalSequence":
"""Expand the segment upstream and/or downstream.
Negative values contract the segment.
Parameters
----------
upstream
Bases to expand (or contract if negative) toward the 5' end.
dnstream
Bases to expand (or contract if negative) toward the 3' end.
Defaults to upstream (symmetric).
Raises
------
ValueError
If contraction would result in end5 being downstream of end3.
"""
if dnstream is None:
dnstream = upstream
up_step = self._upstream_index_step()
down_step = -up_step
new_end5 = self.end5_index + (upstream * up_step)
new_end3 = self.end3_index + (dnstream * down_step)
# Validate end5 is still upstream of or equal to end3
if (new_end5 - new_end3) * up_step < 0:
raise ValueError(
"Invalid expansion: end5 would be downstream of end3 "
f"(end5={new_end5}, end3={new_end3})"
)
return self._from_end_indices(new_end5, new_end3)
[docs]
def expand_coord(
self, upstream: int, dnstream: int | None = None
) -> "DisjointIntervalSequence":
"""Expand the coordinate space and possibly segment at its 5' and/or 3' ends.
The outer 5' edge of the first coord interval is extended ``upstream``
bases and the outer 3' edge of the last coord interval is extended
``dnstream`` bases. If the segment spans the coordinate intervals exactly,
it is expanded an equal amount in the upstream/downstream direction as
the coordinate intervals, otherwise it is left unexpanded.
Parameters
----------
upstream
Bases to add at the coordinate space's 5' end.
dnstream
Bases to add at the coordinate space's 3' end.
Defaults to upstream (symmetric).
Raises
------
ValueError
If either argument is negative.
"""
if dnstream is None:
dnstream = upstream
if upstream < 0 or dnstream < 0:
raise ValueError(
"expand_coord requires non-negative values; "
f"got upstream={upstream}, dnstream={dnstream}"
)
# Segment does not span coord intervals exactly: do not expand segment
if self.start != 0 or self.end != self.coordinate_length:
# Re-indexing of coord space only occurs on upstream expansion, so only
# adjust indices based on this change. Downstream expansion leaves existing
# indices intact, so no adjustment needed
new_start = self._start + upstream
new_end = self._end + upstream
else:
# The interval's lower index stays at self._start in the new coord space due
# to re-indexing of the coord space implicitly 'expanding' start upstream.
# end index is adjusted to account for extra bases introduced by upstream and
# downstream expansion
new_start = self._start
new_end = self._end + upstream + dnstream
coord_ivs = list(self._coordinate_intervals)
iv0, ivn = coord_ivs[0], coord_ivs[-1]
if len(coord_ivs) == 1:
coord_ivs = [
iv0.expand(upstream, dnstream) # iv0 is ivn
]
else:
coord_ivs[0] = iv0.expand(upstream, 0)
coord_ivs[-1] = ivn.expand(0, dnstream)
return DisjointIntervalSequence(
coord_ivs,
coord_name=self._coord_metadata.name,
segment_name=self._segment_metadata.name,
on_coordinate_strand=self.on_coordinate_strand,
start=new_start,
end=new_end
)
[docs]
def upstream_of(self, other: "DisjointIntervalSequence") -> bool:
"""True if self is strictly upstream of other (no overlap).
Requires the same coordinate space and same on_coordinate_strand.
"""
self._validate_same_coordinate_space(other)
if self.on_coordinate_strand != other.on_coordinate_strand:
raise ValueError(
f"Cannot compare: self is on "
f"{'same' if self.on_coordinate_strand else 'opposite'} "
f"strand but other is on "
f"{'same' if other.on_coordinate_strand else 'opposite'} strand"
)
if self.length == 0 and other.length == 0 and self.start == other.start:
return False
if self._upstream_index_step() == -1:
return self._end <= other.start
return self._start >= other.end
[docs]
def dnstream_of(self, other: "DisjointIntervalSequence") -> bool:
"""True if self is strictly downstream of other (no overlap).
Requires the same coordinate space and same on_coordinate_strand.
"""
self._validate_same_coordinate_space(other)
if self.on_coordinate_strand != other.on_coordinate_strand:
raise ValueError(
f"Cannot compare: self is on "
f"{'same' if self.on_coordinate_strand else 'opposite'} "
f"strand but other is on "
f"{'same' if other.on_coordinate_strand else 'opposite'} strand"
)
if self.length == 0 and other.length == 0 and self.start == other.start:
return False
if self._upstream_index_step() == -1:
return self._start >= other.end
return self._end <= other.start
[docs]
def within(self, other: "DisjointIntervalSequence") -> bool:
"""True if self's segment is contained within other's segment.
Requires the same coordinate space and same on_coordinate_strand.
"""
self._validate_same_coordinate_space(other)
if self.on_coordinate_strand != other.on_coordinate_strand:
raise ValueError(
f"Cannot compare: self is on "
f"{'same' if self.on_coordinate_strand else 'opposite'} "
f"strand but other is on "
f"{'same' if other.on_coordinate_strand else 'opposite'} strand"
)
return self._start >= other.start and self._end <= other.end
[docs]
def is_same_strand(self) -> bool:
"""True if the segment is on the same strand as the coordinate intervals.
"""
return self.on_coordinate_strand
[docs]
def is_positive_strand(self) -> bool:
"""True if the segment is on the positive strand.
"""
if self.strand == "+":
return True
return False
[docs]
def as_positive_strand(self) -> "DisjointIntervalSequence":
"""Return a DIS with the segment on the positive strand.
Returns ``self`` if already on the positive strand. The coordinate
intervals are unchanged; only the segment strand is affected.
"""
if self.is_positive_strand():
return self
return self.flip_strand()
[docs]
def as_negative_strand(self) -> "DisjointIntervalSequence":
"""Return a DIS with the segment on the negative strand.
Returns ``self`` if already on the negative strand. The coordinate
intervals are unchanged; only the segment strand is affected.
"""
if not self.is_positive_strand():
return self
return self.flip_strand()
[docs]
def as_opposite_strand(self) -> "DisjointIntervalSequence":
"""Return a DIS with the segment on the opposite strand.
Returns ``self`` if already on the opposite strand. The coordinate
intervals are unchanged; only the segment strand is affected.
"""
if not self.on_coordinate_strand:
return self
return self.flip_strand()
[docs]
def as_same_strand(self) -> "DisjointIntervalSequence":
"""Return a DIS with the segment on the coordinate strand.
Returns ``self`` if already on the coordinate strand. The coordinate
intervals are unchanged; only the segment strand is affected.
"""
if self.on_coordinate_strand:
return self
return self.flip_strand()
[docs]
def flip_strand(self) -> "DisjointIntervalSequence":
"""Return a new DIS with ``on_coordinate_strand`` toggled.
The coordinate intervals are unchanged. The segment's
``on_coordinate_strand`` is flipped.
"""
return DisjointIntervalSequence(
self._coordinate_intervals,
coord_name=self._coord_metadata.name,
segment_name=self._segment_metadata.name,
on_coordinate_strand=not self.on_coordinate_strand,
start=self._start,
end=self._end,
)
[docs]
def lower(self) -> list[Interval]:
"""Project the interval back to genomic :py:class:`~genome_kit.Interval` objects.
Returns the genomic representation of the DIS interval as one or more
Intervals, in 5'->3' order. Multiple intervals are returned when the
DIS interval spans boundaries between coordinate intervals. All
returned Intervals use the effective strand of the DIS interval.
Indices outside the coordinate space are extrapolated linearly from
the nearest boundary, so returned Intervals may extend beyond the
chromosome (their start may be negative or their end may exceed the
chromosome length).
When lowering a 0-length interval on an internal boundary, the returned
value corresponds to either an upstream, or downstream interval boundary
(more specifically, whichever value is an Interval.start,
since Interval.end is exclusive).
"""
chrom = self.chromosome
strand = self.strand
refg = self.reference_genome
coord_strand = self.coord_strand
on_plus = coord_strand == "+"
start, end = self._start, self._end
if start == end:
positions = self._lower_coord(start)
# On an internal boundary _lower_coord returns two flanking values;
# collapse to a single 0-length Interval corresponding to the value of
# the index that is an Interval.start (we do this because Interval.end
# is exclusive).
pos = positions[0] if len(positions) == 1 else max(positions)
return [Interval(chrom, strand, pos, pos, refg)]
# _lower_coord returns 1 value interior/edge/extrapolated, or 2 values
# in 5'->3' order at an internal boundary. Cumulative len(_lower_coord(start)) +
# len(_lower_coord(end)) is 2 (neither on boundary), 3 (one on boundary), or 4
# (both on boundary). In all cases, the start_pos should be the downstream-most
# value and the end_pos should be the upstream-most value.
start_positions = self._lower_coord(start)
end_positions = self._lower_coord(end)
start_pos = start_positions[-1] # Since sorted 5' -> 3', -1 is downstream
end_pos = end_positions[0] # 0 is the upstream value
start_iv = Interval(chrom, coord_strand, start_pos, start_pos, refg)
end_iv = Interval(chrom, coord_strand, end_pos, end_pos, refg)
coord_ivs = self._coordinate_intervals
last_interval_index = len(coord_ivs) - 1
first_idx = 0
while first_idx < last_interval_index and start_iv.dnstream_of(
coord_ivs[first_idx]
):
first_idx += 1
last_idx = last_interval_index
while last_idx > 0 and end_iv.upstream_of(coord_ivs[last_idx]):
last_idx -= 1
result: list[Interval] = []
for i in range(first_idx, last_idx + 1):
iv = coord_ivs[i]
is_first = i == first_idx
is_last = i == last_idx
if on_plus:
genomic_start = start_pos if is_first else iv.start
genomic_end = end_pos if is_last else iv.end
else:
# swap start and end positions since going from coordinate-space that
# runs 5' -> 3' (DIS) to genomic space (3' -> 5' on - strand)
genomic_start = end_pos if is_last else iv.start
genomic_end = start_pos if is_first else iv.end
result.append(Interval(chrom, strand, genomic_start, genomic_end, refg))
# The loop emits intervals in coord 5'->3' order. When the segment is on
# the opposite strand, that is the reverse of segment 5'->3'.
if not self.on_coordinate_strand:
result.reverse()
return result
[docs]
def genomic_span(self) -> Interval:
"""Return a single :py:class:`~genome_kit.Interval` spanning the
segment's genomic extent (from its 5'-most to 3'-most positions),
ignoring gaps between coordinate intervals.
"""
start_positions = self._lower_coord(self.start)
end_positions = self._lower_coord(self.end)
start_pos = start_positions[-1] # Since sorted 5' -> 3', -1 is downstream
end_pos = end_positions[0] # 0 is the upstream value
return Interval(
self.chromosome,
self.strand,
min(start_pos, end_pos),
max(start_pos, end_pos),
self.reference_genome,
)
def _lift_position(self, pos: int) -> int:
"""Map a genomic position to a DIS index in this coordinate space.
Positions outside the coord intervals are linearly extrapolated from
the nearest outer edge. Positions in a gap between coord intervals
are clipped to the cumulative end of the previous interval (i.e. the
boundary index).
"""
ivs = self._coordinate_intervals
coord_len = self.coordinate_length
if self.coord_strand == "+":
if pos < ivs[0].start:
return pos - ivs[0].start
if pos > ivs[-1].end:
return coord_len + (pos - ivs[-1].end)
cumulative = 0
for iv in ivs:
if iv.start <= pos <= iv.end:
return cumulative + (pos - iv.start)
if pos < iv.start:
return cumulative
cumulative += len(iv)
assert False, "Position not found in any interval"
# minus
if pos > ivs[0].end:
return ivs[0].end - pos
if pos < ivs[-1].start:
return coord_len + (ivs[-1].start - pos)
cumulative = 0
for iv in ivs:
if iv.start <= pos <= iv.end:
return cumulative + (iv.end - pos)
if pos > iv.end:
return cumulative
cumulative += len(iv)
assert False, "Position not found in any interval"
[docs]
def lift_interval(self, other: Interval) -> "DisjointIntervalSequence | None":
"""Lift a genomic :py:class:`~genome_kit.Interval` onto this DIS's segment.
The interval's genomic span is mapped into DIS coordinate-space indices
and intersected with this DIS's segment. The returned DIS represents
that intersection as a segment in this coordinate space.
``other`` may lie on either strand. The returned DIS's
``on_coordinate_strand`` reflects whether ``other`` is on the
coordinate strand, so an interval on the strand opposite this DIS's
segment lifts to an opposite-strand segment.
Raises
------
ValueError
If ``other`` is not on the same chromosome and reference genome.
"""
if other.chromosome != self.chromosome:
raise ValueError(
f"Interval chromosome {other.chromosome!r} does not match DIS "
f"chromosome {self.chromosome!r}"
)
if other.reference_genome != self.reference_genome:
raise ValueError(
f"Interval reference_genome {other.reference_genome!r} does not "
f"match DIS reference_genome {self.reference_genome!r}"
)
if self.coord_strand == "+":
seg_start = self._lift_position(other.start)
seg_end = self._lift_position(other.end)
else:
# On minus, lower genomic position maps to higher DIS index.
seg_start = self._lift_position(other.end)
seg_end = self._lift_position(other.start)
# Clip to self's segment via half-open intersection.
intersected_start = max(seg_start, self._start)
intersected_end = min(seg_end, self._end)
if intersected_start >= intersected_end:
return None
return DisjointIntervalSequence(
self._coordinate_intervals,
coord_name=self._coord_metadata.name,
on_coordinate_strand=(other.strand == self.coord_strand),
start=intersected_start,
end=intersected_end,
)
[docs]
def intersect(
self, other: "DisjointIntervalSequence"
) -> "DisjointIntervalSequence | None":
"""Segment-wise intersection with another DIS in the same coord space.
Both DIS objects must share the same coordinate intervals and the same
``on_coordinate_strand``. 0-length intersections return None.
Raises
------
ValueError
If the DIS objects do not share the same coordinate space or
``on_coordinate_strand`` differs.
"""
self._validate_same_coordinate_space(other)
if self.on_coordinate_strand != other.on_coordinate_strand:
raise ValueError(
f"Cannot compare: self is on "
f"{'same' if self.on_coordinate_strand else 'opposite'} "
f"strand but other is on "
f"{'same' if other.on_coordinate_strand else 'opposite'} strand"
)
intersected_start = max(self._start, other._start)
intersected_end = min(self._end, other._end)
if intersected_start >= intersected_end:
return None
return DisjointIntervalSequence(
self._coordinate_intervals,
coord_name=self._coord_metadata.name,
on_coordinate_strand=self.on_coordinate_strand,
start=intersected_start,
end=intersected_end,
)
[docs]
def dna(self, allow_outside_coord: bool = True) -> str:
"""Return the DNA sequence corresponding to the segment.
Bases outside the coordinate intervals (segment indices ``< 0`` or
``>= coordinate_length``) are returned as ``N`` when
``allow_outside_coord`` is True; otherwise a ValueError is raised.
Returns an empty string for a 0-length segment. DNA is returned in
5' -> 3' order.
Parameters
----------
allow_outside_coord
When True (default), pad out-of-coord regions with ``N``. When
False, raise if the segment extends past the coord intervals.
Raises
------
ValueError
If the segment extends past the coord intervals and
``allow_outside_coord`` is False.
"""
if self._start == self._end:
return ""
coord_len = self.coordinate_length
if self._upstream_index_step() == -1:
upstream_pad = 0 if self._start >= 0 else abs(self._start)
downstream_pad = max(self._end - coord_len, 0)
else:
upstream_pad = max(self._end - coord_len, 0)
downstream_pad = 0 if self._start >= 0 else abs(self._start)
if (upstream_pad or downstream_pad) and not allow_outside_coord:
raise ValueError(
f"DIS segment [{self._start}, {self._end}) extends outside "
f"coord range [0, {coord_len})"
)
# Clip start to being at least start of coord space
clipped_start = max(self._start, 0)
# Clip end to being at most end of coord space
clipped_end = min(self._end, coord_len)
# If the clipped segment has a positive length, get the DNA
if clipped_start < clipped_end:
clipped = DisjointIntervalSequence(
self._coordinate_intervals,
on_coordinate_strand=self.on_coordinate_strand,
start=clipped_start,
end=clipped_end,
)
genome = Genome(self.reference_genome)
clipped_dna = "".join(genome.dna(iv) for iv in clipped.lower())
else:
clipped_dna = ""
return "N" * upstream_pad + clipped_dna + "N" * downstream_pad
[docs]
def __len__(self) -> int:
"""Return the length of the segment."""
return self.length
[docs]
def __repr__(self) -> str:
"""Return a human-readable representation."""
return (
f"DisjointIntervalSequence("
f"coord_name={self._coord_metadata.name!r}, "
f"name={self._segment_metadata.name!r}, "
f"{self.chromosome}:{self.coord_strand}, "
f"len={self.length}, "
f"coord_intervals={self._coordinate_intervals}, "
f"start={self._start}, "
f"end={self._end}, "
f"end5={self.end5_index}, "
f"end3={self.end3_index})"
)
[docs]
def __eq__(self, other: object) -> bool:
"""Equality based on coordinate intervals, metadata, and index values."""
if not isinstance(other, DisjointIntervalSequence):
return NotImplemented
return (
self._coord_metadata == other._coord_metadata
and self._segment_metadata == other._segment_metadata
and self._start == other._start
and self._end == other._end
and self._coordinate_intervals == other._coordinate_intervals
)