///|
/// Ambiguous observations require an explicit cleavage policy.
pub(all) enum AmbiguityPolicy {
CertainOnly
IncludePossible
RejectPossible
} derive(Eq, Debug, ToJson, @json.FromJson)
///|
/// A fragment records the reference top strand, not a mass estimate.
pub(all) struct Fragment {
start : Int
length : Int
sequence : String
circular : Bool
} derive(Eq, Debug, ToJson, @json.FromJson)
///|
/// Reports retain recognition evidence even when a site cannot be cut.
pub(all) struct Digest {
topology : Topology
source_length : Int
sites : Array[Site]
cuts : Array[Site]
fragments : Array[Fragment]
policy : AmbiguityPolicy
} derive(Eq, Debug, ToJson, @json.FromJson)
///|
fn selected_cuts(
sites : Array[Site],
policy : AmbiguityPolicy,
) -> Array[Site] raise RestrictError {
let chosen : Array[Site] = []
for site in sites {
if site.certainty == Possible && policy == RejectPossible {
raise InvalidInput("AMBIGUOUS_SITE")
}
if site.cuttable && (site.certainty == Certain || policy == IncludePossible) {
chosen.push(site)
}
}
chosen.sort_by_key(s => s.top)
let out : Array[Site] = []
for s in chosen {
if out.length() > 0 && out[out.length() - 1].top == s.top {
if out[out.length() - 1].stagger != s.stagger {
raise CutConflict("SAME_TOP_DIFFERENT_BOTTOM")
}
} else {
out.push(s)
}
}
out
}
///|
fn partition(
dna : Dna,
cuts : Array[Site],
sites : Array[Site],
policy : AmbiguityPolicy,
) -> Digest raise RestrictError {
if dna.topology == Circular {
return circular_partition(dna, cuts, sites, policy)
}
let fragments : Array[Fragment] = []
let mut previous = 0
let mut previous_bottom = 0
let mut has_cut = false
for cut in cuts {
if (has_cut && cut.bottom <= previous_bottom) ||
cut.bottom < previous_bottom {
raise CutConflict("CROSSING_STRANDS")
}
if cut.top > previous {
fragments.push({
start: previous,
length: cut.top - previous,
sequence: dna.bases[previous:cut.top].to_owned(),
circular: false,
})
}
previous = cut.top
previous_bottom = cut.bottom
has_cut = true
}
if previous < dna.length() {
fragments.push({
start: previous,
length: dna.length() - previous,
sequence: dna.bases[previous:].to_owned(),
circular: false,
})
}
{
topology: dna.topology,
source_length: dna.length(),
sites,
cuts,
fragments,
policy,
}
}
///|
/// Complete single-model cleavage. CertainOnly never cuts merely possible sites.
pub fn digest(
dna : Dna,
enzyme : Enzyme,
policy? : AmbiguityPolicy = CertainOnly,
) -> Digest raise RestrictError {
let sites = scan(dna, enzyme)
partition(dna, selected_cuts(sites, policy), sites, policy)
}
///|
/// Reconstruct ordered top-strand fragments (rotated for cut circular DNA).
pub fn Digest::reconstruct(self : Digest) -> String {
let out = StringBuilder::new()
for fragment in self.fragments {
out.write_string(fragment.sequence)
}
out.to_string()
}
///|
/// Extract at most one revolution, always in reference top-strand direction.
fn circular_bases(dna : Dna, start : Int, length : Int) -> String {
let out = StringBuilder::new(size_hint=length)
for i in 0.. Digest raise RestrictError {
let fragments : Array[Fragment] = []
if cuts.is_empty() {
fragments.push({
start: 0,
length: dna.length(),
sequence: dna.bases,
circular: true,
})
} else {
for i in 0..