dnagent_domain/
fragment_annotations.rs

1//! Source-linked annotation projections onto each digest strand's 5′→3′ frame.
2use crate::digest::{Digest, DigestError, DigestFragment, FragmentStrand, simulate_digest};
3use crate::{DomainError, Feature, ImportedPrimer, Region, SequenceRecord, Strand, Topology};
4use serde::Serialize;
5use thiserror::Error;
6
7#[derive(Debug, Error)]
8pub enum AnnotationError {
9    #[error(transparent)]
10    Digest(#[from] DigestError),
11    #[error(transparent)]
12    Coordinates(#[from] DomainError),
13}
14
15#[derive(Debug, Clone, PartialEq, Eq, Serialize)]
16pub struct MappedPart {
17    /// Index into the original feature's ordered location parts, zero-based.
18    pub source_part: usize,
19    /// Offset within that original part, in its forward-coordinate traversal.
20    pub source_offset: usize,
21    pub source_region: Region,
22    /// Coordinates on THIS exported strand's 5′→3′ sequence, not a duplex axis.
23    pub fragment_region: Region,
24}
25
26#[derive(Debug, Clone, PartialEq, Eq, Serialize)]
27pub struct FeatureMapping {
28    pub source_feature_id: String,
29    pub fragment_strand: Strand,
30    pub source_bases: usize,
31    pub retained_bases: usize,
32    /// Coverage only, NOT a statement that a CDS remains functional or in frame.
33    pub complete: bool,
34    /// Original contiguous parts split across opposite ends of a linear product.
35    pub split_source_parts: Vec<usize>,
36    /// Original part order, then offset within each part; never genomic re-sorting.
37    pub parts: Vec<MappedPart>,
38}
39
40#[derive(Debug, Clone, PartialEq, Eq, Serialize)]
41pub struct FragmentAnnotations {
42    pub fragment_id: String,
43    pub top: Vec<FeatureMapping>,
44    pub bottom: Vec<FeatureMapping>,
45}
46
47#[derive(Debug, Clone, PartialEq, Eq, Serialize)]
48pub struct AnnotatedDigest {
49    pub source_name: String,
50    /// Unmodified originals, including qualifiers, operators and display hints.
51    pub source_features: Vec<Feature>,
52    /// Retained only as source metadata; binding sites are not inferred.
53    pub source_primers: Vec<ImportedPrimer>,
54    pub digest: Digest,
55    pub annotations: Vec<FragmentAnnotations>,
56    pub assumptions: Vec<&'static str>,
57}
58
59pub fn annotated_digest(
60    record: &SequenceRecord,
61    names: &[String],
62) -> Result<AnnotatedDigest, AnnotationError> {
63    let digest = simulate_digest(record.sequence(), record.topology(), names)?;
64    let annotations = digest
65        .fragments
66        .iter()
67        .map(|fragment| {
68            Ok(FragmentAnnotations {
69                fragment_id: fragment.id.clone(),
70                top: map_features(record, fragment, false)?,
71                bottom: map_features(record, fragment, true)?,
72            })
73        })
74        .collect::<Result<Vec<_>, DomainError>>()?;
75    Ok(AnnotatedDigest {
76        source_name: record.name().to_owned(),
77        source_features: record.features().to_vec(),
78        source_primers: record.primers().to_vec(),
79        digest,
80        annotations,
81        assumptions: vec![
82            "Annotations projected separately onto each strand's 5-prime-to-3-prime sequence; bottom projection reverses coordinates and feature strand",
83            "complete means base coverage only; split_source_parts flags original contiguous parts broken across linear fragment ends",
84            "Source feature operators, part order and qualifiers are retained without translation, frame, splicing or biological-function inference",
85            "Source qualifiers may no longer apply to clipped products; they are provenance, not new product annotations",
86            "Source primers retained without binding-site inference; opaque import metadata is not exported",
87            "FASTA is sequence-only: use the JSON report for annotations, source coordinates and duplex end geometry",
88        ],
89    })
90}
91
92fn map_features(
93    record: &SequenceRecord,
94    fragment: &DigestFragment,
95    bottom: bool,
96) -> Result<Vec<FeatureMapping>, DomainError> {
97    let strand = if bottom {
98        &fragment.bottom
99    } else {
100        &fragment.top
101    };
102    record
103        .features()
104        .iter()
105        .filter_map(|feature| {
106            match map_feature(
107                feature,
108                strand,
109                fragment.topology,
110                record.sequence().len(),
111                bottom,
112            ) {
113                Ok(mapping) if mapping.parts.is_empty() => None,
114                result => Some(result),
115            }
116        })
117        .collect()
118}
119
120pub(crate) fn map_feature(
121    feature: &Feature,
122    strand: &FragmentStrand,
123    topology: Topology,
124    n: usize,
125    bottom: bool,
126) -> Result<FeatureMapping, DomainError> {
127    let mut parts = Vec::new();
128    let mut split_source_parts = Vec::new();
129    let mut source_bases = 0;
130    let begin = strand.source_start as i128;
131    let finish = begin + strand.length as i128;
132    for (index, region) in feature.location().parts().iter().enumerate() {
133        source_bases += region.length().get();
134        let mut pieces = Vec::new();
135        // Both source and fragment spans are <= n; these three lifts cover all
136        // possible intersections, including origin-spanning source annotations.
137        for shift in [-(n as i128), 0, n as i128] {
138            let start = region.start().get() as i128 + shift;
139            let end = start + region.length().get() as i128;
140            let lo = start.max(begin);
141            let hi = end.min(finish);
142            if lo >= hi {
143                continue;
144            }
145            let source_start = lo.rem_euclid(n as i128) as usize;
146            let length = (hi - lo) as usize;
147            let source_region = if source_start + length <= n {
148                Region::linear(source_start, source_start + length, n)?
149            } else {
150                Region::circular_arc(source_start, length, n)?
151            };
152            let (a, b) = if bottom {
153                (finish - hi, finish - lo)
154            } else {
155                (lo - begin, hi - begin)
156            };
157            pieces.push(MappedPart {
158                source_part: index,
159                source_offset: (lo - start) as usize,
160                source_region,
161                fragment_region: Region::linear(a as usize, b as usize, strand.length)?,
162            });
163        }
164        pieces.sort_by_key(|p| p.source_offset);
165        if topology == Topology::Linear && pieces.len() > 1 {
166            split_source_parts.push(index);
167        }
168        parts.extend(pieces);
169    }
170    let retained_bases = parts.iter().map(|p| p.fragment_region.length().get()).sum();
171    let fragment_strand = match (feature.location().strand(), bottom) {
172        (Strand::Forward, true) => Strand::Reverse,
173        (Strand::Reverse, true) => Strand::Forward,
174        (strand, _) => strand,
175    };
176    Ok(FeatureMapping {
177        source_feature_id: feature.id().as_str().to_owned(),
178        fragment_strand,
179        source_bases,
180        retained_bases,
181        complete: retained_bases == source_bases,
182        split_source_parts,
183        parts,
184    })
185}
186
187#[cfg(test)]
188mod tests {
189    use super::*;
190    use crate::{DisplayHints, DnaSeq, FeatureId, Location, LocationOperator, Qualifier};
191    fn record(
192        seq: &str,
193        topology: Topology,
194        parts: Vec<Region>,
195        strand: Strand,
196        operator: LocationOperator,
197    ) -> SequenceRecord {
198        SequenceRecord::new(
199            "synthetic",
200            DnaSeq::new(seq).unwrap(),
201            topology,
202            vec![Feature::new(
203                FeatureId::new("f").unwrap(),
204                "CDS",
205                "synthetic feature",
206                Location::new(parts, strand, operator).unwrap(),
207                vec![Qualifier {
208                    key: "note".into(),
209                    value: None,
210                }],
211                DisplayHints::default(),
212            )],
213            vec![],
214        )
215        .unwrap()
216    }
217    #[test]
218    fn staggered_cuts_project_distinct_strands_and_preserve_provenance() {
219        let r = record(
220            "AAAGAATTCTTT",
221            Topology::Linear,
222            vec![Region::linear(2, 10, 12).unwrap()],
223            Strand::Forward,
224            LocationOperator::Contiguous,
225        );
226        let d = annotated_digest(&r, &["EcoRI".into()]).unwrap();
227        assert_eq!(d.source_features, r.features());
228        assert_eq!(d.annotations[0].top[0].retained_bases, 2);
229        assert_eq!(d.annotations[0].bottom[0].retained_bases, 6);
230        assert!(!d.annotations[0].top[0].complete);
231        assert_eq!(d.annotations[0].bottom[0].fragment_strand, Strand::Reverse);
232        assert_eq!(
233            d.annotations[0].bottom[0].parts[0].fragment_region,
234            Region::linear(0, 6, 8).unwrap()
235        );
236    }
237    #[test]
238    fn complete_coverage_does_not_hide_a_feature_split_at_linear_ends() {
239        let r = record(
240            "AAAGAATTCTTT",
241            Topology::Circular,
242            vec![Region::linear(2, 10, 12).unwrap()],
243            Strand::Reverse,
244            LocationOperator::Contiguous,
245        );
246        let d = annotated_digest(&r, &["EcoRI".into()]).unwrap();
247        let mapping = &d.annotations[0].top[0];
248        assert!(mapping.complete);
249        assert_eq!(mapping.split_source_parts, [0]);
250        assert_eq!(
251            mapping
252                .parts
253                .iter()
254                .map(|p| p.source_offset)
255                .collect::<Vec<_>>(),
256            [0, 2]
257        );
258        assert_eq!(d.annotations[0].bottom[0].fragment_strand, Strand::Forward);
259    }
260    #[test]
261    fn order_operator_and_unknown_strand_are_preserved() {
262        let r = record(
263            "AAAGAATTCTTT",
264            Topology::Linear,
265            vec![
266                Region::linear(8, 11, 12).unwrap(),
267                Region::linear(1, 3, 12).unwrap(),
268            ],
269            Strand::Unknown,
270            LocationOperator::Order,
271        );
272        let d = annotated_digest(&r, &["BamHI".into()]).unwrap();
273        assert_eq!(
274            d.source_features[0].location().operator(),
275            LocationOperator::Order
276        );
277        assert_eq!(d.annotations[0].bottom[0].fragment_strand, Strand::Unknown);
278        assert_eq!(
279            d.annotations[0].top[0]
280                .parts
281                .iter()
282                .map(|p| p.source_part)
283                .collect::<Vec<_>>(),
284            [0, 1]
285        );
286    }
287    #[test]
288    fn origin_spanning_uncut_circle_is_not_flagged_as_disrupted() {
289        let r = record(
290            "AAAGAATTCTTT",
291            Topology::Circular,
292            vec![Region::circular_arc(10, 4, 12).unwrap()],
293            Strand::Forward,
294            LocationOperator::Contiguous,
295        );
296        let d = annotated_digest(&r, &["BamHI".into()]).unwrap();
297        assert!(d.annotations[0].top[0].complete);
298        assert_eq!(d.annotations[0].top[0].split_source_parts, []);
299        assert_eq!(d.annotations[0].top[0].parts.len(), 2);
300    }
301}
dnagent_domain/fragment_annotations.rs

dnagent_domain/
fragment_annotations.rs
