cealign 0.1.0

Implementation of the Combinatorial Extension (cealign) algorithm to align protein structures
Documentation 
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
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90

use plotters::prelude::*;

/// Plot the alignment path, mirroring Figure 2 from Shindyalov & Bourne (1998).
///
/// Each AFP is drawn as a thick diagonal segment from (i_a, i_b) to
/// (i_a+win_size-1, i_b+win_size-1). Consecutive AFPs are joined by a thin
/// grey connector. Axes are residue indices into protein A (x) and B (y).
/// Output is written to `plot.png`.
pub fn plot(
    path: &[(usize, usize)],
    win_size: usize,
    len_a: usize,
    len_b: usize,
) -> Result<(), Box<dyn std::error::Error>> {
    let root = BitMapBackend::new("plot.png", (640, 640)).into_drawing_area();
    root.fill(&WHITE)?;
    let root = root.margin(30, 20, 20, 30);

    let mut chart = ChartBuilder::on(&root)
        .caption("Alignment path", ("sans-serif", 18))
        .x_label_area_size(30)
        .y_label_area_size(35)
        .build_cartesian_2d(0u32..len_a as u32, 0u32..len_b as u32)?;

    chart
        .configure_mesh()
        .x_desc("Protein A residue")
        .y_desc("Protein B residue")
        .x_labels(10)
        .y_labels(10)
        .x_label_formatter(&|x| format!("{}", x))
        .y_label_formatter(&|y| format!("{}", y))
        .draw()?;

    // Thin grey connectors between the end of one AFP and the start of the next
    if path.len() > 1 {
        let connectors: Vec<(u32, u32)> = path
            .windows(2)
            .flat_map(|w| {
                let (i_a, i_b) = w[0];
                let (j_a, j_b) = w[1];
                vec![
                    ((i_a + win_size - 1) as u32, (i_b + win_size - 1) as u32),
                    (j_a as u32, j_b as u32),
                ]
            })
            .collect();
        chart.draw_series(LineSeries::new(connectors, &BLACK.mix(0.35)))?;
    }

    // Thick diagonal segment for each AFP
    for &(i_a, i_b) in path {
        let segment: Vec<(u32, u32)> = (0..win_size)
            .map(|k| ((i_a + k) as u32, (i_b + k) as u32))
            .collect();
        chart.draw_series(LineSeries::new(
            segment,
            ShapeStyle::from(&BLACK).stroke_width(2),
        ))?;
    }

    root.present()?;
    Ok(())
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_plot_single_afp() {
        let path = vec![(0usize, 0usize)];
        assert!(plot(&path, 8, 46, 46).is_ok());
        std::fs::remove_file("plot.png").ok();
    }

    #[test]
    fn test_plot_multiple_afps() {
        // Multiple AFPs exercise both the diagonal segments and the connectors
        let path = vec![(0, 0), (10, 10), (20, 20)];
        assert!(plot(&path, 8, 46, 46).is_ok());
        std::fs::remove_file("plot.png").ok();
    }

    #[test]
    fn test_plot_empty_path() {
        assert!(plot(&[], 8, 46, 46).is_ok());
        std::fs::remove_file("plot.png").ok();
    }
}