twitcher 0.5.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

use std::borrow::Borrow;

use rust_htslib::bcf::Record;

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

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());
        let span = first.pos().abs_diff(last.end()) as f64 + 1.0;
        let complexity: f64 = cluster
            .iter()
            .map(|r| (get_record(r), get_allele(r)))
            .flat_map(|(r, a)| record_complexity(r, a))
            .sum();
        self.is_valid_cluster(complexity, span, cluster.len())
    }
}

fn record_complexity(record: &Record, allele: Option<u32>) -> anyhow::Result<f64> {
    let alls = record.alleles();
    if alls.len() < 2 {
        anyhow::bail!("Record should have at least two alleles");
    }
    if let Some(a) = allele {
        Ok(variant_complexity(alls[0].len(), alls[a as usize].len()))
    } else {
        Ok(alls[1..]
            .iter()
            .map(|a| variant_complexity(alls[0].len(), a.len()))
            .reduce(|acc, e| if e > acc { e } else { acc })
            .unwrap())
    }
}