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
39#[derive(Debug, PartialEq)]
40pub enum LoHiState {
41    Low,
42    High,
43}
44
45impl fmt::Display for RefSpecError {
46    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
47        let err_msg = match self {
48            RefSpecError::MalformedInt(mfi) => format!("Malformed integer {}", mfi),
49            RefSpecError::ZeroRef => "Ref # must be > 0".to_string(),
50            RefSpecError::RefTooLarge(max) => format!("Ref # too large (max {})", max),
51        };
52        write!(f, "{}", err_msg)
53    }
54}
55
56pub struct Alignment {
57    pub headers: Vec<String>,
58    pub sequences: Vec<String>,
59    /* The consensus sequence is now a field of Alignment, and is computed once upon creation. This
60     * contrasts with the very first implementation, in which the consensus was recomputed every
61     * time the UI was drawn... which was very inefficient but had this funny "twinkling" effect in
62     * columns with tied residue frequencies. This was due to the fact that HashMap stores its keys
63     * in an unpredictable order, and that different calls to keys() may return them indifferent
64     * orders. See best_residue().
65     */
66    /* These are properties of the whole _alignment_, or at least of whole columns. They cannot be
67     * meaningfully attributed to a sequence. */
68    pub consensus: String,
69    pub entropies: Vec<f64>,
70    pub densities: Vec<f64>,
71
72    /* By contrast, the following are properties of sequences (at least in part). Length, for
73     * example, does not depend on anything but the sequence itself, and could be a field in a
74     * struct that also contains the sequence and its header. */
75    pub id_wrt_reference: Vec<f64>, // reference is usually the consensus, but CAN be an aln seq.
76    // Recompute if ref changes.
77    // Of course the sequence length is an integer, but using an integer type like u32 would make
78    // it hard (for me, at least...) to write a function that accepts a Vec of either lengths or
79    // %IDs. Tried Box, and generics, but the extra work doesn't seem warranted.
80    pub relative_seq_len: Vec<f64>,
81    pub macromolecule_type: SeqType,
82    /* Specifies whether the reference sequence should be the consensus (see above) or one of the
83     * original sequences (identified by rank).*/
84    ref_spec: RefSpec,
85}
86
87#[derive(Debug, PartialEq)]
88struct BestResidue {
89    residue: char,
90    frequency: u64,
91}
92
93impl Alignment {
94    // Makes an Alignment from a SeqFile, which is consumed.
95    pub fn from_file(seq_file: SeqFile) -> Alignment {
96        let mut headers: Vec<String> = Vec::new();
97        let mut sequences: Vec<String> = Vec::new();
98        let mut max_len: usize = 0;
99        for record in seq_file {
100            headers.push(record.header);
101            let l = record.sequence.len();
102            sequences.push(record.sequence);
103            if l > max_len {
104                max_len = l;
105            }
106        }
107        // Pad any sequence shorter than max_len, so we are not limited to alignments with exactly
108        // identical numbers of positions (reviewer suggestion).
109        sequences
110            .iter_mut()
111            .for_each(|s| *s = format!("{:<width$}", s, width = max_len));
112        // NOTE: the 's' can also be written '&*s', which makes the automatic re-borrow explicit.
113        let first_seq = sequences.first();
114        let macromolecule_type = seq_type(first_seq.expect("No sequence found."));
115        let consensus = consensus(&sequences, macromolecule_type);
116        let entropies = entropies(&sequences);
117        let densities = densities(&sequences);
118        let id_wrt_reference = sequences
119            .iter()
120            .map(|seq| percent_identity(seq, &consensus))
121            .collect();
122        let relative_seq_len = sequences.iter().map(|seq| seq_len_nogaps(seq)).collect();
123
124        Alignment {
125            headers,
126            sequences,
127            consensus,
128            entropies,
129            densities,
130            id_wrt_reference,
131            relative_seq_len,
132            macromolecule_type,
133            ref_spec: RefSpec::Consensus,
134        }
135    }
136
137    // Makes an Alignment from a Vec of headers and a Vec of Strings, which are consumed. Mostly
138    // used for testing.
139    #[allow(dead_code)]
140    pub fn from_vecs(hdrs: Vec<String>, seqs: Vec<String>) -> Alignment {
141        assert_eq!(hdrs.len(), seqs.len());
142        let headers = hdrs;
143        let sequences = seqs;
144        let first_seq = sequences.first();
145        let macromolecule_type = seq_type(first_seq.expect("No sequence found."));
146        let consensus = consensus(&sequences, macromolecule_type);
147        let entropies = entropies(&sequences);
148        let densities = densities(&sequences);
149        let id_wrt_reference = sequences
150            .iter()
151            .map(|seq| percent_identity(seq, &consensus))
152            .collect();
153        let relative_seq_len = sequences.iter().map(|seq| seq_len_nogaps(seq)).collect();
154
155        Alignment {
156            headers,
157            sequences,
158            consensus,
159            entropies,
160            densities,
161            id_wrt_reference,
162            relative_seq_len,
163            macromolecule_type,
164            ref_spec: RefSpec::Consensus,
165        }
166    }
167
168    pub fn num_seq(&self) -> usize {
169        self.sequences.len()
170    }
171
172    // TODO: shouldn't this be aln_width?
173    pub fn aln_len(&self) -> usize {
174        self.sequences[0].len()
175    }
176
177    pub fn macromolecule_type(&self) -> SeqType {
178        self.macromolecule_type
179    }
180
181    pub fn get_ref_spec(&self) -> RefSpec {
182        self.ref_spec
183    }
184
185    pub fn set_ref_spec(&mut self, spec: RefSpec) -> Result<(), RefSpecError> {
186        match spec {
187            // Note: the rank in a RefSpec is 0-based. The conversion from user-land is done in
188            // app.rs.
189            RefSpec::Rank(rk) if rk >= self.num_seq() => {
190                return Err(RefSpecError::RefTooLarge(self.num_seq()));
191            }
192            _ => self.ref_spec = spec,
193        }
194        // Probable change of ref -> Recompute the identities WRT ref
195        let reference = self.reference();
196        self.id_wrt_reference = self
197            .sequences
198            .iter()
199            .map(|seq| percent_identity(seq, &reference))
200            .collect();
201        Ok(())
202    }
203
204    pub fn reference(&self) -> String {
205        match self.ref_spec {
206            RefSpec::Consensus => self.consensus.clone(),
207            RefSpec::Rank(rk) => self.sequences[rk].clone(),
208        }
209    }
210}
211
212fn res_count(sequences: &Vec<String>, col: usize) -> ResidueCounts {
213    let mut freqs: ResidueCounts = HashMap::new();
214    for seq in sequences {
215        let residue = seq.as_bytes()[col] as char;
216        *freqs.entry(residue).or_insert(0) += 1;
217    }
218    freqs
219}
220
221// TODO: consensus(), best_residue(), entropies(), densities(), and related functions
222// should be methods of Alignment instead of free functions. This would eliminate the need
223// to thread parameters like SeqType through signatures — they could access self.macromolecule_type
224// directly. See the B0024 fix (2026-08-25) where SeqType was threaded through consensus()
225// as motivation for this refactor.
226
227pub fn consensus(sequences: &Vec<String>, seq_type: SeqType) -> String {
228    let mut consensus = String::new();
229    for j in 0..sequences[0].len() {
230        let dist = res_count(sequences, j); // res -> count map
231        let br = best_residue(&dist, seq_type);
232        let rel_freq: f64 = (br.frequency as f64 / sequences.len() as f64) as f64;
233        if rel_freq >= UC_CONS_THRESHOLD {
234            consensus.push(br.residue.to_ascii_uppercase());
235        } else if rel_freq >= LC_CONS_THRESHOLD {
236            if br.residue.is_alphabetic() {
237                consensus.push(br.residue.to_ascii_lowercase());
238            } else {
239                //consensus.push('-');
240                consensus.push(br.residue);
241            }
242        } else {
243            consensus.push('*');
244        }
245    }
246    consensus
247}
248
249pub fn entropies(sequences: &Vec<String>) -> Vec<f64> {
250    let mut entropies: Vec<f64> = Vec::new();
251    for j in 0..sequences[0].len() {
252        let dist = res_count(sequences, j);
253        let freq = to_freq_distrib(&dist);
254        let e = entropy(&freq);
255        entropies.push(e);
256    }
257    entropies
258}
259
260pub fn col_density(sequences: &Vec<String>, col: usize) -> f64 {
261    let mut mass = 0;
262    for seq in sequences {
263        match seq.as_bytes()[col] as char {
264            'a'..='z' | 'A'..='Z' => mass += 1,
265            '-' | '.' | ' ' => {}
266            other => {
267                panic!("Character {other} unexpected in an alignment.\nThis might be due to file format, please see option -f.");
268            }
269        }
270    }
271    mass as f64 / sequences.len() as f64
272}
273
274pub fn densities(sequences: &Vec<String>) -> Vec<f64> {
275    (0..sequences[0].len())
276        .map(|col| col_density(sequences, col))
277        .collect()
278}
279
280fn iupac_ambiguity_code(amb_nt: &mut [char]) -> char {
281    let mut normalized_nt = amb_nt
282        .into_iter()
283        .map(|nt| nt.to_ascii_lowercase())
284        .collect::<Vec<char>>();
285    normalized_nt.sort();
286
287    let normalized_nt_as_string = normalized_nt.into_iter().join("");
288
289    match normalized_nt_as_string.as_str() {
290        "a" => 'a',    // Adenine
291        "ac" => 'm',   // aMino
292        "acg" => 'v',  // not-T (not-U), V follows U
293        "acgt" => 'n', // aNy
294        "act" => 'h',  // not-G, H follows G in the alphabet
295        "ag" => 'r',   // puRine
296        "agt" => 'd',  // not-C, D follows C
297        "at" => 'w',   // Weak interaction (2 H bonds)
298        "c" => 'c',    // Cytosine
299        "cg" => 's',   // Strong interaction (3 H bonds)
300        "cgt" => 'b',  // not-A, B follows A
301        "ct" => 'y',   // pYrimidine
302        "g" => 'g',    // Guanine
303        "gt" => 'k',   // Keto
304        "t" => 't',    // Thymine
305        &_ => 'n',     // Anything else folded into N - might want to bail, perhaps?
306    }
307}
308
309fn best_residue(counts: &ResidueCounts, seq_type: SeqType) -> BestResidue {
310    let max_freq = counts.values().max().unwrap();
311    let mut most_frequent_residues = counts // plural <- may be ties
312        .keys()
313        .filter(|&&k| counts.get(&k) == Some(max_freq))
314        .map(|&k| k)
315        .collect::<Vec<char>>();
316
317    let residue = if most_frequent_residues.len() == 1 {
318        most_frequent_residues[0]
319    } else if SeqType::Protein == seq_type {
320        'X'
321    } else {
322        iupac_ambiguity_code(&mut most_frequent_residues)
323    };
324
325    BestResidue {
326        residue: residue,
327        frequency: *max_freq,
328    }
329}
330
331// Convert a residue -> count map into a residue -> frequency map (relative frequency, that is).
332// While gaps are allowed (and indeed useful) in the former, they are not included in the latter
333// (in particular because they make litle sense when computing entropy).
334//
335fn to_freq_distrib(counts: &ResidueCounts) -> ResidueDistribution {
336    let total_counts: u64 = counts
337        .iter()
338        .filter(|(res, _count)| **res != '-')
339        .map(|(_res, count)| count)
340        .sum();
341    let mut distrib = ResidueDistribution::new();
342    for (residue, count) in counts.iter() {
343        if *residue == '-' {
344            continue;
345        }
346        distrib.insert(*residue, *count as f64 / total_counts as f64);
347    }
348    distrib
349}
350
351fn entropy(freqs: &ResidueDistribution) -> f64 {
352    // Discard '-'s
353    let residues: Vec<&char> = freqs.keys().filter(|&&r| r != '-').collect();
354    let sum: f64 = residues
355        .into_iter()
356        .map(|res| {
357            let p = *freqs.get(res).unwrap();
358            p * p.ln()
359        })
360        .sum();
361
362    -sum
363}
364
365fn percent_identity(s1: &str, s2: &str) -> f64 {
366    let num_identical = s1
367        .chars()
368        .zip(s2.chars())
369        .filter(|(c1, c2)| c1.eq_ignore_ascii_case(c2))
370        .count();
371    num_identical as f64 / s1.len() as f64
372}
373
374fn seq_len_nogaps(s: &str) -> f64 {
375    s.chars().filter(|c| c.is_alphabetic()).count() as f64 / s.len() as f64
376}
377
378fn seq_type(sequence: &str) -> SeqType {
379    let counts = sequence.to_lowercase().chars().counts();
380    let counts_u64: HashMap<char, u64> = counts.into_iter().map(|(k, v)| (k, v as u64)).collect();
381    let frequencies = to_freq_distrib(&counts_u64);
382    let nt_freq: f64 = *frequencies.get(&'a').unwrap_or(&0.0)
383        + *frequencies.get(&'c').unwrap_or(&0.0)
384        + *frequencies.get(&'g').unwrap_or(&0.0)
385        + *frequencies.get(&'t').unwrap_or(&0.0)
386        + *frequencies.get(&'u').unwrap_or(&0.0);
387    // A quick-and dirty heuristic, I'm afraid
388    if nt_freq > 0.75 {
389        Nucleic
390    } else {
391        Protein
392    }
393}
394
395pub fn mark_lohi(metric: &[f64], threshold: f64) -> Vec<LoHiState> {
396    assert!(!threshold.is_nan(), "threshold must not be NaN");
397    metric
398        .iter()
399        .map(|&v| {
400            assert!(!v.is_nan(), "metric value must not be NaN");
401
402            if v < threshold {
403                LoHiState::Low
404            } else {
405                LoHiState::High
406            }
407        })
408        .collect()
409}
410
411pub fn find_hi_runs(lohi_states: &[LoHiState]) -> Vec<(usize, usize)> {
412    let mut runs = Vec::new();
413    let mut run_start: Option<usize> = None;
414
415    for (i, state) in lohi_states.iter().enumerate() {
416        match (run_start, state) {
417            // start of a high-metric chunk
418            (None, LoHiState::High) => {
419                run_start = Some(i);
420            }
421            // any part of a low-metric chunk
422            (Some(start), LoHiState::Low) => {
423                runs.push((start, i - start));
424                run_start = None;
425            }
426            _ => {}
427        }
428    }
429
430    // trailing chunk
431    if let Some(start) = run_start {
432        runs.push((start, lohi_states.len() - start));
433    }
434
435    runs
436}
437
438pub fn merge_hi_runs(runs: &[(usize, usize)], threshold: usize) -> Vec<(usize, usize)> {
439    let mut merged_runs = Vec::new();
440
441    if runs.is_empty() {
442        return merged_runs;
443    }
444
445    merged_runs.push(runs[0]);
446
447    for &(cur_run_start, cur_run_len) in &runs[1..] {
448        let (prev_run_start, prev_run_len) = *merged_runs.last().unwrap();
449
450        let low_run_start = prev_run_start + prev_run_len;
451        let low_run_len = cur_run_start - low_run_start;
452
453        if low_run_len >= threshold {
454            merged_runs.push((cur_run_start, cur_run_len));
455        } else {
456            merged_runs.last_mut().unwrap().1 += low_run_len + cur_run_len;
457        }
458    }
459
460    merged_runs
461}
462
463#[cfg(test)]
464mod tests {
465    use crate::alignment::{
466        best_residue, consensus, densities, entropies, entropy, find_hi_runs, mark_lohi,
467        merge_hi_runs, percent_identity, res_count, seq_len_nogaps, seq_type, to_freq_distrib,
468        Alignment, BestResidue, LoHiState, RefSpec, ResidueCounts, ResidueDistribution, SeqType,
469        SeqType::{Nucleic, Protein},
470    };
471    use crate::seq::fasta::read_fasta_file;
472    use approx::assert_relative_eq;
473    use std::collections::HashMap;
474
475    #[test]
476    fn test_read_aln() {
477        let fasta1 = read_fasta_file("./data/test2.fas").unwrap();
478        let aln1 = Alignment::from_file(fasta1);
479        assert_eq!("seq1", aln1.headers[0]);
480        assert_eq!("seq2", aln1.headers[1]);
481        assert_eq!("seq3", aln1.headers[2]);
482        assert_eq!("TTGCCG-CGA", aln1.sequences[0]);
483        assert_eq!("TTCCCGGCGA", aln1.sequences[1]);
484        assert_eq!("TTACCG-CAA", aln1.sequences[2]);
485    }
486
487    #[test]
488    fn test_consensus() {
489        let fasta2 = read_fasta_file("data/test-cons.fas").unwrap();
490        let aln2 = Alignment::from_file(fasta2);
491        // Updated: the output changed with the IUPAC codes implementation. Position 2 and 4
492        // are now resolved through the ambiguity code matcher, which defaults to 'n' for
493        // protein residues. This is expected behavior until the function gains protein support.
494        assert_eq!("AQw-n", consensus(&aln2.sequences, SeqType::Protein));
495    }
496
497    #[test]
498    fn test_res_count() {
499        let fasta2 = read_fasta_file("data/test-cons.fas").unwrap();
500        let aln2 = Alignment::from_file(fasta2);
501        let mut d0: ResidueCounts = HashMap::new();
502        d0.insert('A', 6);
503        assert_eq!(d0, res_count(&aln2.sequences, 0));
504
505        let mut d1: ResidueCounts = HashMap::new();
506        d1.insert('Q', 5);
507        d1.insert('T', 1);
508        assert_eq!(d1, res_count(&aln2.sequences, 1));
509
510        let mut d2: ResidueCounts = HashMap::new();
511        d2.insert('W', 2);
512        d2.insert('I', 1);
513        d2.insert('S', 1);
514        d2.insert('D', 1);
515        d2.insert('F', 1);
516        assert_eq!(d2, res_count(&aln2.sequences, 2));
517
518        let mut d3: ResidueCounts = HashMap::new();
519        d3.insert('-', 3);
520        d3.insert('K', 2);
521        d3.insert('L', 1);
522        assert_eq!(d3, res_count(&aln2.sequences, 3));
523    }
524
525    #[test]
526    fn test_most_frequent_residue() {
527        let d0: ResidueCounts = HashMap::from([('A', 6)]);
528        let mut exp: BestResidue = BestResidue {
529            residue: 'A',
530            frequency: 6,
531        };
532        assert_eq!(exp, best_residue(&d0, SeqType::Nucleic));
533
534        let d1: ResidueCounts = HashMap::from([('Q', 5), ('T', 1)]);
535        exp = BestResidue {
536            residue: 'Q',
537            frequency: 5,
538        };
539        assert_eq!(exp, best_residue(&d1, SeqType::Protein));
540
541        let d2: ResidueCounts = HashMap::from([('W', 2), ('I', 1), ('S', 1), ('D', 1), ('F', 1)]);
542        exp = BestResidue {
543            residue: 'W',
544            frequency: 2,
545        };
546        assert_eq!(exp, best_residue(&d2, SeqType::Protein));
547
548        // col 3 cannot be tested <- ties
549
550        let d4: ResidueCounts = HashMap::from([('-', 3), ('K', 2), ('L', 1)]);
551        exp = BestResidue {
552            residue: '-',
553            frequency: 3,
554        };
555        assert_eq!(exp, best_residue(&d4, SeqType::Protein));
556    }
557
558    #[test]
559    fn test_to_freq_distrib() {
560        let eps = 0.001;
561        let counts: ResidueCounts = HashMap::from([('K', 3), ('L', 3), ('G', 6), ('-', 6)]);
562        let rfreqs = to_freq_distrib(&counts);
563        assert_relative_eq!(0.25, *rfreqs.get(&'K').unwrap(), epsilon = eps);
564        assert_relative_eq!(0.25, *rfreqs.get(&'L').unwrap(), epsilon = eps);
565        assert_relative_eq!(0.5, *rfreqs.get(&'G').unwrap(), epsilon = eps);
566    }
567
568    #[test]
569    fn test_entropy_1() {
570        let eps = 0.00001;
571        let distrib: ResidueDistribution = ResidueDistribution::from([('A', 1.0)]);
572        assert_relative_eq!(0.0, entropy(&distrib), epsilon = eps);
573    }
574
575    #[test]
576    fn test_entropy_2() {
577        let eps = 0.00001;
578        let distrib: ResidueDistribution = ResidueDistribution::from([('A', 0.5), ('F', 0.5)]);
579        // This should be ln(2), and as it happens Rust has a constant for this; remarkably, clippy
580        // detects the literal constant below and suggests using the (arguably more accurate)
581        // built-in definition.
582        // assert_relative_eq!(0.6931471805599453, entropy(&distrib), epsilon = eps);
583        assert_relative_eq!(std::f64::consts::LN_2, entropy(&distrib), epsilon = eps);
584    }
585
586    #[test]
587    fn test_entropy_3() {
588        let eps = 0.00001;
589        let distrib: ResidueDistribution =
590            ResidueDistribution::from([('A', 0.5), ('F', 0.25), ('T', 0.25)]);
591        assert_relative_eq!(1.0397207708399179, entropy(&distrib), epsilon = eps);
592    }
593
594    #[test]
595    fn test_entropies() {
596        let fasta2 = read_fasta_file("data/test-cons.fas").unwrap();
597        let aln2 = Alignment::from_file(fasta2);
598        let entrs = entropies(&aln2.sequences);
599        let eps = 0.001;
600        assert_relative_eq!(0.0, entrs[0], epsilon = eps);
601        assert_relative_eq!(0.4505, entrs[1], epsilon = eps);
602        assert_relative_eq!(1.5607, entrs[2], epsilon = eps);
603        assert_relative_eq!(0.6365, entrs[3], epsilon = eps);
604    }
605
606    #[test]
607    fn test_density() {
608        let fasta = read_fasta_file("data/test-density.msa").unwrap();
609        let aln = Alignment::from_file(fasta);
610        let dens = densities(&aln.sequences);
611        assert_eq!(1.0, dens[0]);
612        assert_eq!(0.8, dens[1]);
613        assert_eq!(0.6, dens[2]);
614        assert_eq!(0.4, dens[3]);
615        assert_eq!(0.2, dens[4]);
616        assert_eq!(0.0, dens[5]);
617    }
618
619    #[test]
620    fn test_order_aln() {
621        let fasta = read_fasta_file("./data/test4.aln").unwrap();
622        let aln1 = Alignment::from_file(fasta);
623        // Check original order
624        assert_eq!("Zea_001", aln1.headers[0]);
625        assert_eq!("Rana_002", aln1.headers[1]);
626        assert_eq!("Panthera_050", aln1.headers[49]);
627        assert_eq!("tgctgttcgtcaaAgtaggcc", aln1.sequences[0]);
628        assert_eq!("tgctgttAgAcaaagtaggcc", aln1.sequences[1]);
629        assert_eq!("tgctgttcgtcaaagtaggcc", aln1.sequences[49]);
630    }
631
632    #[test]
633    fn test_similarity_00() {
634        let s1 = "GAATTC";
635        assert_eq!(percent_identity(s1, s1), 1.0);
636    }
637
638    #[test]
639    fn test_similarity_05() {
640        let s1 = "GAATTC";
641        let s2 = "GAA---";
642        assert_eq!(percent_identity(s1, s2), 0.5);
643    }
644
645    #[test]
646    fn test_similarity_10() {
647        let s1 = "GAATTC";
648        let s2 = "gaattc";
649        assert_eq!(percent_identity(s1, s2), 1.0);
650    }
651
652    #[test]
653    fn test_seq_len_nogaps_00() {
654        assert_eq!(seq_len_nogaps("atgc"), 1.0);
655    }
656
657    #[test]
658    fn test_seq_len_nogaps_05() {
659        assert_eq!(seq_len_nogaps("a-gc"), 0.75);
660    }
661
662    #[test]
663    fn test_seq_len_nogaps_10() {
664        assert_eq!(seq_len_nogaps("--.-"), 0.0);
665    }
666
667    #[test]
668    fn test_seq_type_00() {
669        assert_eq!(Nucleic, seq_type("GAATTC"));
670    }
671
672    #[test]
673    fn test_seq_type_05() {
674        assert_eq!(Protein, seq_type("HGTSDA"));
675    }
676
677    #[test]
678    fn test_seq_type_10() {
679        assert_eq!(Nucleic, seq_type("cgatgcacgatgcncagtgtuucgatcga"));
680    }
681
682    #[test]
683    fn test_seq_type_15() {
684        assert_eq!(Nucleic, seq_type("UUTGAU"));
685    }
686
687    // Make sure seq files with unequal lengths get correctly padded
688    #[test]
689    fn test_unequal_seq_len() {
690        let fasta = read_fasta_file("./data/test5.aln").unwrap();
691        let _ = Alignment::from_file(fasta);
692    }
693
694    // Test the Vec constructor
695    #[test]
696    fn test_vec_ctor_00() {
697        let hdrs = vec![
698            String::from("Leo"),
699            String::from("Tigris"),
700            String::from("Pardus"),
701            String::from("Onca"),
702        ];
703        let seqs = vec![
704            String::from("catgcatatg"),
705            String::from("aatgcatatg"),
706            String::from("tatgcatatg"),
707            String::from("gatgcatatg"),
708        ];
709        let aln = Alignment::from_vecs(hdrs, seqs);
710        assert_eq!(4, aln.num_seq());
711        assert_eq!(10, aln.aln_len());
712        assert_eq!(SeqType::Nucleic, aln.macromolecule_type());
713        assert_eq!("Onca", aln.headers[3]);
714        assert_eq!("gatgcatatg", aln.sequences[3]);
715    }
716
717    // Test the reference specifier
718    #[test]
719    fn test_reference_specifier() {
720        let hdrs = vec![
721            String::from("frugilegus"),
722            String::from("monedula"),
723            String::from("corax"),
724            String::from("corone"),
725            String::from("cornix"),
726        ];
727        let seqs = vec![
728            String::from("catgcatatg"),
729            String::from("aatgcatatg"),
730            String::from("tatgcatatg"),
731            String::from("tatgcatatg"),
732            String::from("gatgcatatg"),
733        ];
734        let mut aln = Alignment::from_vecs(hdrs, seqs);
735        // By default, the reference sequence is the consensus
736        assert_eq!(RefSpec::Consensus, aln.get_ref_spec());
737        assert_eq!("tATGCATATG", aln.reference());
738        // Now set the ref to the first sequence (rank 0)
739        let _ = aln.set_ref_spec(RefSpec::Rank(0));
740        assert_eq!(RefSpec::Rank(0), aln.get_ref_spec());
741        assert_eq!("catgcatatg", aln.reference());
742        // Back to consensus
743        let _ = aln.set_ref_spec(RefSpec::Consensus);
744        assert_eq!(RefSpec::Consensus, aln.get_ref_spec());
745        assert_eq!("tATGCATATG", aln.reference());
746    }
747
748    // Tests the %id WRT ref (incl. when != consensus)
749    #[test]
750    fn test_pct_id_wrt_ref() {
751        let hdrs = vec![
752            String::from("frugilegus"),
753            String::from("monedula"),
754            String::from("corax"),
755            String::from("corone"),
756            String::from("cornix"),
757        ];
758        // consensus: ACg-
759        let seqs = vec![
760            String::from("A---"),
761            String::from("AC--"),
762            String::from("ACG-"),
763            String::from("ACGT"),
764            String::from("ACGT"),
765        ];
766        let mut aln = Alignment::from_vecs(hdrs, seqs);
767        // Check the ref, which by default is the consensus
768        assert_eq!("ACg-", aln.reference());
769        assert_eq!(vec![0.5, 0.75, 1.0, 0.75, 0.75], aln.id_wrt_reference);
770        // Now switch to seq #0 for reference
771        let _ = aln.set_ref_spec(RefSpec::Rank(0));
772        assert_eq!("A---", aln.reference());
773        assert_eq!(vec![1.0, 0.75, 0.5, 0.25, 0.25], aln.id_wrt_reference);
774        // Switch back to consensus
775        let _ = aln.set_ref_spec(RefSpec::Consensus);
776        assert_eq!("ACg-", aln.reference());
777        assert_eq!(vec![0.5, 0.75, 1.0, 0.75, 0.75], aln.id_wrt_reference);
778    }
779
780    #[test]
781    fn test_mark_lohi() {
782        let metric = vec![0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0];
783        assert_eq!(
784            mark_lohi(&metric, 0.8),
785            vec![
786                LoHiState::Low,
787                LoHiState::Low,
788                LoHiState::Low,
789                LoHiState::Low,
790                LoHiState::Low,
791                LoHiState::Low,
792                LoHiState::Low,
793                LoHiState::High,
794                LoHiState::High,
795                LoHiState::High,
796            ]
797        );
798    }
799
800    #[test]
801    fn test_find_hi_runs() {
802        let lohi_states = vec![
803            LoHiState::High,
804            LoHiState::Low,
805            LoHiState::Low,
806            LoHiState::Low,
807            LoHiState::High, // start of a hi run at pos 4 (length 4)
808            LoHiState::High,
809            LoHiState::High,
810            LoHiState::High,
811            LoHiState::Low,
812            LoHiState::Low,
813            LoHiState::Low,
814            LoHiState::Low,
815            LoHiState::High, // start of a hi run at pos 12 (length 5)
816            LoHiState::High,
817            LoHiState::High,
818            LoHiState::High,
819            LoHiState::High,
820            LoHiState::High,
821            LoHiState::Low,
822            LoHiState::Low,
823            LoHiState::Low,
824            LoHiState::High,
825            LoHiState::High,
826        ];
827        assert_eq!(
828            find_hi_runs(&lohi_states),
829            vec![(0, 1), (4, 4), (12, 6), (21, 2)]
830        );
831    }
832
833    #[test]
834    fn test_find_hi_runs_all_high() {
835        let lohi_states = vec![LoHiState::High, LoHiState::High, LoHiState::High];
836        assert_eq!(find_hi_runs(&lohi_states), vec![(0, 3)]);
837    }
838
839    #[test]
840    fn test_find_hi_runs_all_low() {
841        let lohi_states = vec![LoHiState::Low, LoHiState::Low, LoHiState::Low];
842        assert_eq!(find_hi_runs(&lohi_states), vec![]);
843    }
844
845    #[test]
846    fn test_merge_hi_runs_single_run() {
847        assert_eq!(merge_hi_runs(&[(5, 3)], 3), vec![(5, 3)]);
848    }
849
850    #[test]
851    fn test_merge_hi_runs_merges_short_gaps_only() {
852        // Runs are:
853        // 0..4, 6..8, 11..14, 16..18, 22..24
854        //
855        // Gaps are:
856        // 4..6   len 2  -> merge if threshold is 3
857        // 8..11  len 3  -> do not merge if threshold is 3
858        // 14..16 len 2  -> merge if threshold is 3
859        // 18..22 len 4  -> do not merge if threshold is 3
860
861        let runs = vec![(0, 4), (6, 2), (11, 3), (16, 2), (22, 2)];
862
863        assert_eq!(merge_hi_runs(&runs, 3), vec![(0, 8), (11, 7), (22, 2)]);
864    }
865
866    // B0024: consensus was nondeterministic when the most frequent residue at a position was tied.
867    // HashMap iteration order is randomized per process, so best_residue() would pick an arbitrary
868    // tied residue. Fixed by computing IUPAC ambiguity codes for ties.
869
870    #[test]
871    fn test_consensus_is_stable_across_calls() {
872        // Multiple calls to consensus on the same alignment should yield identical results.
873        // Without the fix, this fails ~50% of the time due to HashMap iteration order.
874        let seqs = vec![
875            "AGGCTC".to_string(),
876            "AGGCAC".to_string(),
877            "ACGTGC".to_string(),
878        ];
879        let consensus1 = consensus(&seqs, SeqType::Nucleic);
880        let consensus2 = consensus(&seqs, SeqType::Nucleic);
881        let consensus3 = consensus(&seqs, SeqType::Nucleic);
882        assert_eq!(
883            consensus1, consensus2,
884            "consensus changed between calls: '{}' vs '{}'",
885            consensus1, consensus2
886        );
887        assert_eq!(
888            consensus2, consensus3,
889            "consensus changed between calls: '{}' vs '{}'",
890            consensus2, consensus3
891        );
892    }
893
894    #[test]
895    fn test_consensus_uses_iupac_codes_for_tied_nucleotides() {
896        // Position 4 (0-indexed) has a three-way tie: T, A, G each appear once.
897        // The consensus should use the IUPAC code 'D' (not A/C/G) for this position.
898        // D = adenine, Guanine, Thymine (not C).
899        let seqs = vec![
900            "AGGCTC".to_string(),
901            "AGGCAC".to_string(),
902            "ACGTGC".to_string(),
903        ];
904        let cons = consensus(&seqs, SeqType::Nucleic);
905        // Positions: 0=A (3/3), 1=G (2/3), 2=G (3/3), 3=C (2/3), 4=d/D (tie 1/3), 5=C (3/3)
906        // The exact case depends on frequency threshold, but position 4 should be deterministic.
907        assert!(
908            cons.chars().nth(4).unwrap().to_ascii_lowercase() == 'd',
909            "position 4 should resolve to IUPAC code 'd' for A/G/T tie, got '{}'",
910            cons.chars().nth(4).unwrap()
911        );
912    }
913
914    #[test]
915    fn test_consensus_tie_breakpoint_by_frequency() {
916        // Two residues tied at the top frequency should produce an IUPAC code.
917        // (Frequency 2/4 = 0.5, above the lowercase threshold of ~20%, so renders lowercase.)
918        let seqs = vec![
919            "ACT".to_string(),
920            "ACT".to_string(),
921            "ATT".to_string(),
922            "ATT".to_string(),
923        ];
924        // Position 0: A 4/4
925        // Position 1: C 2/4, T 2/4 (tie at top frequency)
926        // Position 2: T 4/4
927        let cons = consensus(&seqs, SeqType::Nucleic);
928        // Position 1 should be a nucleotide IUPAC code (C and T -> 'y' for pYrimidine).
929        let pos1 = cons.chars().nth(1).unwrap().to_ascii_lowercase();
930        assert_eq!(
931            pos1, 'y',
932            "position 1 should resolve to IUPAC code 'y' for C/T tie, got '{}'",
933            pos1
934        );
935    }
936}
termal_alignment/alignment.rs

termal_alignment/
alignment.rs
