twitcher 0.6.0

Find template switch mutations in genomic data
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73

use std::borrow::Borrow;

use rust_htslib::bcf::Record;

pub use crate::common::cluster_settings::ClusteringSettings;

impl ClusteringSettings {
    pub fn belongs<R: Borrow<Record>>(
        &self,
        last: Option<R>,
        candidate: R,
    ) -> anyhow::Result<bool> {
        let Some(last) = last else {
            return Ok(true);
        };
        Ok(candidate.borrow().pos().saturating_sub(last.borrow().end())
            <= self.max_gap.try_into()?)
    }

    pub fn is_cluster<R, F: Fn(&R) -> &Record, G: Fn(&R) -> Option<u32>>(
        &self,
        cluster: &[R],
        get_record: F,
        get_allele: G,
    ) -> bool {
        if cluster.is_empty() {
            return false;
        }
        let first = get_record(cluster.first().unwrap().borrow());
        let last = get_record(cluster.last().unwrap().borrow());
        // Reference extent of the cluster (inclusive).
        let l_ref = first.pos().abs_diff(last.end()) as f64 + 1.0;
        // Accumulate cluster mass and the net length change (Σ alt_len − ref_len) so we can derive
        // the query extent. Records without a usable allele pair are skipped.
        let (mass, net) = cluster
            .iter()
            .map(|r| (get_record(r), get_allele(r)))
            .filter_map(|(r, a)| self.chosen_allele_lens(r, a))
            .fold((0.0, 0.0), |(mass, net), (ref_len, alt_len)| {
                (
                    mass + self.event_mass(ref_len, alt_len),
                    net + alt_len as f64 - ref_len as f64,
                )
            });
        let span = self.span(l_ref, l_ref + net);
        self.is_valid_cluster(mass, span, cluster.len())
    }

    /// Lengths of the reference allele and the relevant alt allele for a record. When no specific
    /// allele is requested, the alt maximizing [`ClusteringSettings::event_mass`] is chosen so mass
    /// and span use the same allele. Returns `None` for records without at least two alleles.
    fn chosen_allele_lens(&self, record: &Record, allele: Option<u32>) -> Option<(usize, usize)> {
        let alls = record.alleles();
        if alls.len() < 2 {
            return None;
        }
        let ref_len = alls[0].len();
        let alt_len = if let Some(a) = allele {
            alls.get(a as usize)?.len()
        } else {
            alls[1..]
                .iter()
                .map(|a| a.len())
                .max_by(|x, y| {
                    self.event_mass(ref_len, *x)
                        .partial_cmp(&self.event_mass(ref_len, *y))
                        .unwrap()
                })
                .unwrap()
        };
        Some((ref_len, alt_len))
    }
}