termal_alignment/
alignment.rs

1// SPDX-License-Identifier: MIT
2// Copyright (c) 2025-2026 Thomas Junier
3
4mod permutation;
5
6use std::{collections::HashMap, fmt};
7
8use itertools::Itertools;
9
10use crate::seq::file::SeqFile;
11
12use crate::alignment::SeqType::{Nucleic, Protein};
13
14// Whether to show the most frequent residue as LC or UC
15const UC_CONS_THRESHOLD: f64 = 0.8; // uppercase if at least this
16const LC_CONS_THRESHOLD: f64 = 0.2; // lowercase if at least this (else '*')
17
18type ResidueDistribution = HashMap<char, f64>;
19type ResidueCounts = HashMap<char, u64>;
20
21#[derive(PartialEq, Clone, Copy, Debug)]
22pub enum SeqType {
23    Nucleic,
24    Protein,
25}
26
27#[derive(Clone, Copy, Debug, PartialEq)]
28pub enum RefSpec {
29    Consensus,
30    Rank(usize),
31}
32
33pub enum RefSpecError {
34    MalformedInt(String),
35    ZeroRef,
36    RefTooLarge(usize),
37}
38
39impl fmt::Display for RefSpecError {
40    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
41        let err_msg = match self {
42            RefSpecError::MalformedInt(mfi) => format!("Malformed integer {}", mfi),
43            RefSpecError::ZeroRef => "Ref # must be > 0".to_string(),
44            RefSpecError::RefTooLarge(max) => format!("Ref # too large (max {})", max),
45        };
46        write!(f, "{}", err_msg)
47    }
48}
49
50pub struct Alignment {
51    pub headers: Vec<String>,
52    pub sequences: Vec<String>,
53    /* The consensus sequence is now a field of Alignment, and is computed once upon creation. This
54     * contrasts with the very first implementation, in which the consensus was recomputed every
55     * time the UI was drawn... which was very inefficient but had this funny "twinkling" effect in
56     * columns with tied residue frequencies. This was due to the fact that HashMap stores its keys
57     * in an unpredictable order, and that different calls to keys() may return them indifferent
58     * orders. See best_residue().
59     */
60    /* These are properties of the whole _alignment_, or at least of whole columns. They cannot be
61     * meaningfully attributed to a sequence. */
62    pub consensus: String,
63    pub entropies: Vec<f64>,
64    pub densities: Vec<f64>,
65
66    /* By contrast, the following are properties of sequences (at least in part). Length, for
67     * example, does not depend on anything but the sequence itself, and could be a field in a
68     * struct that also contains the sequence and its header. */
69    pub id_wrt_reference: Vec<f64>, // reference is usually the consensus, but CAN be an aln seq.
70    // Recompute if ref changes.
71    // Of course the sequence length is an integer, but using an integer type like u32 would make
72    // it hard (for me, at least...) to write a function that accepts a Vec of either  lengths or
73    // %IDs. Tried Box, and generics, but the extra work doesn't seem warranted.
74    pub relative_seq_len: Vec<f64>,
75    pub macromolecule_type: SeqType,
76    /* Specifies whether the reference sequence should be the consensus (see above) or one of the
77     * original sequences (identified by rank).*/
78    ref_spec: RefSpec,
79}
80
81#[derive(Debug, PartialEq)]
82struct BestResidue {
83    residue: char,
84    frequency: u64,
85}
86
87impl Alignment {
88    // Makes an Alignment from a SeqFile, which is consumed.
89    pub fn from_file(seq_file: SeqFile) -> Alignment {
90        let mut headers: Vec<String> = Vec::new();
91        let mut sequences: Vec<String> = Vec::new();
92        let mut max_len: usize = 0;
93        for record in seq_file {
94            headers.push(record.header);
95            let l = record.sequence.len();
96            sequences.push(record.sequence);
97            if l > max_len {
98                max_len = l;
99            }
100        }
101        // Pad any sequence shorter than max_len, so we are not limited to alignments with exactly
102        // identical numbers of positions (reviewer suggestion).
103        sequences
104            .iter_mut()
105            .for_each(|s| *s = format!("{:<width$}", s, width = max_len));
106        // NOTE: the 's' can also be written '&*s', which makes the automatic re-borrow explicit.
107        let consensus = consensus(&sequences);
108        let entropies = entropies(&sequences);
109        let densities = densities(&sequences);
110        let id_wrt_reference = sequences
111            .iter()
112            .map(|seq| percent_identity(seq, &consensus))
113            .collect();
114        let relative_seq_len = sequences.iter().map(|seq| seq_len_nogaps(seq)).collect();
115        let first_seq = sequences.first();
116        let macromolecule_type = seq_type(first_seq.expect("No sequence found."));
117
118        Alignment {
119            headers,
120            sequences,
121            consensus,
122            entropies,
123            densities,
124            id_wrt_reference,
125            relative_seq_len,
126            macromolecule_type,
127            ref_spec: RefSpec::Consensus,
128        }
129    }
130
131    // Makes an Alignment from a Vec of headers and a Vec of Strings, which are consumed. Mostly
132    // used for testing.
133    #[allow(dead_code)]
134    pub fn from_vecs(hdrs: Vec<String>, seqs: Vec<String>) -> Alignment {
135        assert_eq!(hdrs.len(), seqs.len());
136        let headers = hdrs;
137        let sequences = seqs;
138        let consensus = consensus(&sequences);
139        let entropies = entropies(&sequences);
140        let densities = densities(&sequences);
141        let id_wrt_reference = sequences
142            .iter()
143            .map(|seq| percent_identity(seq, &consensus))
144            .collect();
145        let relative_seq_len = sequences.iter().map(|seq| seq_len_nogaps(seq)).collect();
146        let first_seq = sequences.first();
147        let macromolecule_type = seq_type(first_seq.expect("No sequence found."));
148
149        Alignment {
150            headers,
151            sequences,
152            consensus,
153            entropies,
154            densities,
155            id_wrt_reference,
156            relative_seq_len,
157            macromolecule_type,
158            ref_spec: RefSpec::Consensus,
159        }
160    }
161
162    pub fn num_seq(&self) -> usize {
163        self.sequences.len()
164    }
165
166    // TODO: shouldn't this be aln_width?
167    pub fn aln_len(&self) -> usize {
168        self.sequences[0].len()
169    }
170
171    pub fn macromolecule_type(&self) -> SeqType {
172        self.macromolecule_type
173    }
174
175    pub fn get_ref_spec(&self) -> RefSpec {
176        self.ref_spec
177    }
178
179    pub fn set_ref_spec(&mut self, spec: RefSpec) -> Result<(), RefSpecError> {
180        match spec {
181            // Note: the rank in a RefSpec is 0-based. The conversion from user-land is done in
182            // app.rs.
183            RefSpec::Rank(rk) if rk >= self.num_seq() => {
184                return Err(RefSpecError::RefTooLarge(self.num_seq()));
185            }
186            _ => self.ref_spec = spec,
187        }
188        // Probable change of ref -> Recompute the identities WRT ref
189        let reference = self.reference();
190        self.id_wrt_reference = self
191            .sequences
192            .iter()
193            .map(|seq| percent_identity(seq, &reference))
194            .collect();
195        Ok(())
196    }
197
198    pub fn reference(&self) -> String {
199        match self.ref_spec {
200            RefSpec::Consensus => self.consensus.clone(),
201            RefSpec::Rank(rk) => self.sequences[rk].clone(),
202        }
203    }
204}
205
206// TODO should these be methods of Alignment?
207
208fn res_count(sequences: &Vec<String>, col: usize) -> ResidueCounts {
209    let mut freqs: ResidueCounts = HashMap::new();
210    for seq in sequences {
211        let residue = seq.as_bytes()[col] as char;
212        *freqs.entry(residue).or_insert(0) += 1;
213    }
214    freqs
215}
216
217pub fn consensus(sequences: &Vec<String>) -> String {
218    let mut consensus = String::new();
219    for j in 0..sequences[0].len() {
220        let dist = res_count(sequences, j); // res -> count map
221        let br = best_residue(&dist);
222        let rel_freq: f64 = (br.frequency as f64 / sequences.len() as f64) as f64;
223        if rel_freq >= UC_CONS_THRESHOLD {
224            consensus.push(br.residue.to_ascii_uppercase());
225        } else if rel_freq >= LC_CONS_THRESHOLD {
226            if br.residue.is_alphabetic() {
227                consensus.push(br.residue.to_ascii_lowercase());
228            } else {
229                //consensus.push('-');
230                consensus.push(br.residue);
231            }
232        } else {
233            consensus.push('*');
234        }
235    }
236    consensus
237}
238
239pub fn entropies(sequences: &Vec<String>) -> Vec<f64> {
240    let mut entropies: Vec<f64> = Vec::new();
241    for j in 0..sequences[0].len() {
242        let dist = res_count(sequences, j);
243        let freq = to_freq_distrib(&dist);
244        let e = entropy(&freq);
245        entropies.push(e);
246    }
247    entropies
248}
249
250pub fn col_density(sequences: &Vec<String>, col: usize) -> f64 {
251    let mut mass = 0;
252    for seq in sequences {
253        match seq.as_bytes()[col] as char {
254            'a'..='z' | 'A'..='Z' => mass += 1,
255            '-' | '.' | ' ' => {}
256            other => {
257                panic!("Character {other} unexpected in an alignment.\nThis might be due to file format, please see option -f.");
258            }
259        }
260    }
261    mass as f64 / sequences.len() as f64
262}
263
264pub fn densities(sequences: &Vec<String>) -> Vec<f64> {
265    (0..sequences[0].len())
266        .map(|col| col_density(sequences, col))
267        .collect()
268}
269
270fn best_residue(dist: &ResidueCounts) -> BestResidue {
271    let max_freq = dist.values().max().unwrap();
272    let most_frequent_residue = dist
273        .keys()
274        .find(|&&k| dist.get(&k) == Some(max_freq))
275        .unwrap();
276
277    BestResidue {
278        residue: *most_frequent_residue,
279        frequency: *max_freq,
280    }
281}
282
283// Convert a residue -> count map into a residue -> frequency map (relative frequency, that is).
284// While gaps are allowed (and indeed useful) in the former, they are not included in the latter
285// (in particular because they make litle sense when computing entropy).
286//
287fn to_freq_distrib(counts: &ResidueCounts) -> ResidueDistribution {
288    let total_counts: u64 = counts
289        .iter()
290        .filter(|(res, _count)| **res != '-')
291        .map(|(_res, count)| count)
292        .sum();
293    let mut distrib = ResidueDistribution::new();
294    for (residue, count) in counts.iter() {
295        if *residue == '-' {
296            continue;
297        }
298        distrib.insert(*residue, *count as f64 / total_counts as f64);
299    }
300    distrib
301}
302
303fn entropy(freqs: &ResidueDistribution) -> f64 {
304    // Discard '-'s
305    let residues: Vec<&char> = freqs.keys().filter(|&&r| r != '-').collect();
306    let sum: f64 = residues
307        .into_iter()
308        .map(|res| {
309            let p = *freqs.get(res).unwrap();
310            p * p.ln()
311        })
312        .sum();
313
314    -sum
315}
316
317fn percent_identity(s1: &str, s2: &str) -> f64 {
318    let num_identical = s1
319        .chars()
320        .zip(s2.chars())
321        .filter(|(c1, c2)| c1.eq_ignore_ascii_case(c2))
322        .count();
323    num_identical as f64 / s1.len() as f64
324}
325
326fn seq_len_nogaps(s: &str) -> f64 {
327    s.chars().filter(|c| c.is_alphabetic()).count() as f64 / s.len() as f64
328}
329
330fn seq_type(sequence: &str) -> SeqType {
331    let counts = sequence.to_lowercase().chars().counts();
332    let counts_u64: HashMap<char, u64> = counts.into_iter().map(|(k, v)| (k, v as u64)).collect();
333    let frequencies = to_freq_distrib(&counts_u64);
334    let nt_freq: f64 = *frequencies.get(&'a').unwrap_or(&0.0)
335        + *frequencies.get(&'c').unwrap_or(&0.0)
336        + *frequencies.get(&'g').unwrap_or(&0.0)
337        + *frequencies.get(&'t').unwrap_or(&0.0)
338        + *frequencies.get(&'u').unwrap_or(&0.0);
339    // A quick-and dirty heuristic, I'm afraid
340    if nt_freq > 0.75 {
341        Nucleic
342    } else {
343        Protein
344    }
345}
346
347#[cfg(test)]
348mod tests {
349    use crate::alignment::{
350        best_residue, consensus, densities, entropies, entropy, percent_identity, res_count,
351        seq_len_nogaps, seq_type, to_freq_distrib, Alignment, BestResidue, RefSpec, ResidueCounts,
352        ResidueDistribution, SeqType,
353        SeqType::{Nucleic, Protein},
354    };
355    use crate::seq::fasta::read_fasta_file;
356    use approx::assert_relative_eq;
357    use std::collections::HashMap;
358
359    #[test]
360    fn test_read_aln() {
361        let fasta1 = read_fasta_file("./data/test2.fas").unwrap();
362        let aln1 = Alignment::from_file(fasta1);
363        assert_eq!("seq1", aln1.headers[0]);
364        assert_eq!("seq2", aln1.headers[1]);
365        assert_eq!("seq3", aln1.headers[2]);
366        assert_eq!("TTGCCG-CGA", aln1.sequences[0]);
367        assert_eq!("TTCCCGGCGA", aln1.sequences[1]);
368        assert_eq!("TTACCG-CAA", aln1.sequences[2]);
369    }
370
371    #[test]
372    fn test_consensus() {
373        let fasta2 = read_fasta_file("data/test-cons.fas").unwrap();
374        let aln2 = Alignment::from_file(fasta2);
375        assert_eq!("AQw-n", consensus(&aln2.sequences));
376    }
377
378    #[test]
379    fn test_res_count() {
380        let fasta2 = read_fasta_file("data/test-cons.fas").unwrap();
381        let aln2 = Alignment::from_file(fasta2);
382        let mut d0: ResidueCounts = HashMap::new();
383        d0.insert('A', 6);
384        assert_eq!(d0, res_count(&aln2.sequences, 0));
385
386        let mut d1: ResidueCounts = HashMap::new();
387        d1.insert('Q', 5);
388        d1.insert('T', 1);
389        assert_eq!(d1, res_count(&aln2.sequences, 1));
390
391        let mut d2: ResidueCounts = HashMap::new();
392        d2.insert('W', 2);
393        d2.insert('I', 1);
394        d2.insert('S', 1);
395        d2.insert('D', 1);
396        d2.insert('F', 1);
397        assert_eq!(d2, res_count(&aln2.sequences, 2));
398
399        let mut d3: ResidueCounts = HashMap::new();
400        d3.insert('-', 3);
401        d3.insert('K', 2);
402        d3.insert('L', 1);
403        assert_eq!(d3, res_count(&aln2.sequences, 3));
404    }
405
406    #[test]
407    fn test_most_frequent_residue() {
408        let d0: ResidueCounts = HashMap::from([('A', 6)]);
409        let mut exp: BestResidue = BestResidue {
410            residue: 'A',
411            frequency: 6,
412        };
413        assert_eq!(exp, best_residue(&d0));
414
415        let d1: ResidueCounts = HashMap::from([('Q', 5), ('T', 1)]);
416        exp = BestResidue {
417            residue: 'Q',
418            frequency: 5,
419        };
420        assert_eq!(exp, best_residue(&d1));
421
422        let d2: ResidueCounts = HashMap::from([('W', 2), ('I', 1), ('S', 1), ('D', 1), ('F', 1)]);
423        exp = BestResidue {
424            residue: 'W',
425            frequency: 2,
426        };
427        assert_eq!(exp, best_residue(&d2));
428
429        // col 3 cannot be tested <- ties
430
431        let d4: ResidueCounts = HashMap::from([('-', 3), ('K', 2), ('L', 1)]);
432        exp = BestResidue {
433            residue: '-',
434            frequency: 3,
435        };
436        assert_eq!(exp, best_residue(&d4));
437    }
438
439    #[test]
440    fn test_to_freq_distrib() {
441        let eps = 0.001;
442        let counts: ResidueCounts = HashMap::from([('K', 3), ('L', 3), ('G', 6), ('-', 6)]);
443        let rfreqs = to_freq_distrib(&counts);
444        assert_relative_eq!(0.25, *rfreqs.get(&'K').unwrap(), epsilon = eps);
445        assert_relative_eq!(0.25, *rfreqs.get(&'L').unwrap(), epsilon = eps);
446        assert_relative_eq!(0.5, *rfreqs.get(&'G').unwrap(), epsilon = eps);
447    }
448
449    #[test]
450    fn test_entropy_1() {
451        let eps = 0.00001;
452        let distrib: ResidueDistribution = ResidueDistribution::from([('A', 1.0)]);
453        assert_relative_eq!(0.0, entropy(&distrib), epsilon = eps);
454    }
455
456    #[test]
457    fn test_entropy_2() {
458        let eps = 0.00001;
459        let distrib: ResidueDistribution = ResidueDistribution::from([('A', 0.5), ('F', 0.5)]);
460        // This should be ln(2), and as it happens Rust has a constant for this; remarkably, clippy
461        // detects the literal constant below and suggests using the (arguably more accurate)
462        // built-in definition.
463        // assert_relative_eq!(0.6931471805599453, entropy(&distrib), epsilon = eps);
464        assert_relative_eq!(std::f64::consts::LN_2, entropy(&distrib), epsilon = eps);
465    }
466
467    #[test]
468    fn test_entropy_3() {
469        let eps = 0.00001;
470        let distrib: ResidueDistribution =
471            ResidueDistribution::from([('A', 0.5), ('F', 0.25), ('T', 0.25)]);
472        assert_relative_eq!(1.0397207708399179, entropy(&distrib), epsilon = eps);
473    }
474
475    #[test]
476    fn test_entropies() {
477        let fasta2 = read_fasta_file("data/test-cons.fas").unwrap();
478        let aln2 = Alignment::from_file(fasta2);
479        let entrs = entropies(&aln2.sequences);
480        let eps = 0.001;
481        assert_relative_eq!(0.0, entrs[0], epsilon = eps);
482        assert_relative_eq!(0.4505, entrs[1], epsilon = eps);
483        assert_relative_eq!(1.5607, entrs[2], epsilon = eps);
484        assert_relative_eq!(0.6365, entrs[3], epsilon = eps);
485    }
486
487    #[test]
488    fn test_density() {
489        let fasta = read_fasta_file("data/test-density.msa").unwrap();
490        let aln = Alignment::from_file(fasta);
491        let dens = densities(&aln.sequences);
492        assert_eq!(1.0, dens[0]);
493        assert_eq!(0.8, dens[1]);
494        assert_eq!(0.6, dens[2]);
495        assert_eq!(0.4, dens[3]);
496        assert_eq!(0.2, dens[4]);
497        assert_eq!(0.0, dens[5]);
498    }
499
500    #[test]
501    fn test_order_aln() {
502        let fasta = read_fasta_file("./data/test4.aln").unwrap();
503        let aln1 = Alignment::from_file(fasta);
504        // Check original order
505        assert_eq!("Zea_001", aln1.headers[0]);
506        assert_eq!("Rana_002", aln1.headers[1]);
507        assert_eq!("Panthera_050", aln1.headers[49]);
508        assert_eq!("tgctgttcgtcaaAgtaggcc", aln1.sequences[0]);
509        assert_eq!("tgctgttAgAcaaagtaggcc", aln1.sequences[1]);
510        assert_eq!("tgctgttcgtcaaagtaggcc", aln1.sequences[49]);
511    }
512
513    #[test]
514    fn test_similarity_00() {
515        let s1 = "GAATTC";
516        assert_eq!(percent_identity(s1, s1), 1.0);
517    }
518
519    #[test]
520    fn test_similarity_05() {
521        let s1 = "GAATTC";
522        let s2 = "GAA---";
523        assert_eq!(percent_identity(s1, s2), 0.5);
524    }
525
526    #[test]
527    fn test_similarity_10() {
528        let s1 = "GAATTC";
529        let s2 = "gaattc";
530        assert_eq!(percent_identity(s1, s2), 1.0);
531    }
532
533    #[test]
534    fn test_seq_len_nogaps_00() {
535        assert_eq!(seq_len_nogaps("atgc"), 1.0);
536    }
537
538    #[test]
539    fn test_seq_len_nogaps_05() {
540        assert_eq!(seq_len_nogaps("a-gc"), 0.75);
541    }
542
543    #[test]
544    fn test_seq_len_nogaps_10() {
545        assert_eq!(seq_len_nogaps("--.-"), 0.0);
546    }
547
548    #[test]
549    fn test_seq_type_00() {
550        assert_eq!(Nucleic, seq_type("GAATTC"));
551    }
552
553    #[test]
554    fn test_seq_type_05() {
555        assert_eq!(Protein, seq_type("HGTSDA"));
556    }
557
558    #[test]
559    fn test_seq_type_10() {
560        assert_eq!(Nucleic, seq_type("cgatgcacgatgcncagtgtuucgatcga"));
561    }
562
563    #[test]
564    fn test_seq_type_15() {
565        assert_eq!(Nucleic, seq_type("UUTGAU"));
566    }
567
568    // Make sure seq files with unequal lengths get correctly padded
569    #[test]
570    fn test_unequal_seq_len() {
571        let fasta = read_fasta_file("./data/test5.aln").unwrap();
572        let _ = Alignment::from_file(fasta);
573    }
574
575    // Test the Vec constructor
576    #[test]
577    fn test_vec_ctor_00() {
578        let hdrs = vec![
579            String::from("Leo"),
580            String::from("Tigris"),
581            String::from("Pardus"),
582            String::from("Onca"),
583        ];
584        let seqs = vec![
585            String::from("catgcatatg"),
586            String::from("aatgcatatg"),
587            String::from("tatgcatatg"),
588            String::from("gatgcatatg"),
589        ];
590        let aln = Alignment::from_vecs(hdrs, seqs);
591        assert_eq!(4, aln.num_seq());
592        assert_eq!(10, aln.aln_len());
593        assert_eq!(SeqType::Nucleic, aln.macromolecule_type());
594        assert_eq!("Onca", aln.headers[3]);
595        assert_eq!("gatgcatatg", aln.sequences[3]);
596    }
597
598    // Test the reference specifier
599    #[test]
600    fn test_reference_specifier() {
601        let hdrs = vec![
602            String::from("frugilegus"),
603            String::from("monedula"),
604            String::from("corax"),
605            String::from("corone"),
606            String::from("cornix"),
607        ];
608        let seqs = vec![
609            String::from("catgcatatg"),
610            String::from("aatgcatatg"),
611            String::from("tatgcatatg"),
612            String::from("tatgcatatg"),
613            String::from("gatgcatatg"),
614        ];
615        let mut aln = Alignment::from_vecs(hdrs, seqs);
616        // By default, the reference sequence is the consensus
617        assert_eq!(RefSpec::Consensus, aln.get_ref_spec());
618        assert_eq!("tATGCATATG", aln.reference());
619        // Now set the ref to the first sequence (rank 0)
620        let _ = aln.set_ref_spec(RefSpec::Rank(0));
621        assert_eq!(RefSpec::Rank(0), aln.get_ref_spec());
622        assert_eq!("catgcatatg", aln.reference());
623        // Back to consensus
624        let _ = aln.set_ref_spec(RefSpec::Consensus);
625        assert_eq!(RefSpec::Consensus, aln.get_ref_spec());
626        assert_eq!("tATGCATATG", aln.reference());
627    }
628
629    // Tests the %id WRT ref (incl. when != consensus)
630    #[test]
631    fn test_pct_id_wrt_ref() {
632        let hdrs = vec![
633            String::from("frugilegus"),
634            String::from("monedula"),
635            String::from("corax"),
636            String::from("corone"),
637            String::from("cornix"),
638        ];
639        // consensus: ACg-
640        let seqs = vec![
641            String::from("A---"),
642            String::from("AC--"),
643            String::from("ACG-"),
644            String::from("ACGT"),
645            String::from("ACGT"),
646        ];
647        let mut aln = Alignment::from_vecs(hdrs, seqs);
648        // Check the ref, which by default is the consensus
649        assert_eq!("ACg-", aln.reference());
650        assert_eq!(vec![0.5, 0.75, 1.0, 0.75, 0.75], aln.id_wrt_reference);
651        // Now switch to seq #0 for reference
652        let _ = aln.set_ref_spec(RefSpec::Rank(0));
653        assert_eq!("A---", aln.reference());
654        assert_eq!(vec![1.0, 0.75, 0.5, 0.25, 0.25], aln.id_wrt_reference);
655        // Switch back to consensus
656        let _ = aln.set_ref_spec(RefSpec::Consensus);
657        assert_eq!("ACg-", aln.reference());
658        assert_eq!(vec![0.5, 0.75, 1.0, 0.75, 0.75], aln.id_wrt_reference);
659    }
660}
termal_alignment/alignment.rs

termal_alignment/
alignment.rs
