righor 0.2.4

Righor creates model of Ig/TCR sequences from sequencing 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
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131

use crate::shared::Dna;
use crate::shared::{DAlignment, VJAlignment};
use crate::vdj::Model;
use anyhow::{anyhow, Result};
#[cfg(all(feature = "py_binds", feature = "pyo3"))]
use pyo3::prelude::*;

#[derive(Default)]
pub struct Event<'a> {
    pub v: Option<&'a VJAlignment>,
    pub j: Option<&'a VJAlignment>,
    pub d: Option<&'a DAlignment>,
    pub delv: usize,
    pub delj: usize,
    pub deld3: usize,
    pub deld5: usize,
}

#[cfg_attr(all(feature = "py_binds", feature = "pyo3"), pyclass(get_all, set_all))]
#[derive(Default, Clone, Debug, PartialEq, Eq)]
pub struct StaticEvent {
    pub v_index: usize,
    pub v_start_gene: usize, // start of the sequence in the V gene
    pub delv: usize,
    pub j_index: usize,
    pub j_start_seq: usize, // start of the palindromic J gene (with all dels) in the sequence
    pub delj: usize,
    pub d_index: usize,
    pub d_start_seq: usize,
    pub deld3: usize,
    pub deld5: usize,
    pub insvd: Dna,
    pub insdj: Dna,
    pub errors: Vec<(usize, u8)>,
}

#[cfg(all(feature = "py_binds", feature = "pyo3"))]
#[pymethods]
impl StaticEvent {
    fn __repr__(&self) -> String {
        format!(
            "StaticEvent(\n\
		 nb. del. on V3: {},\n\
		 nb. del. on D5: {},\n\
		 nb. del. on D3: {},\n\
		 nb. del. on J5: {},\n\
		 V-D insertions: {},\n\
		 D-J insertions: {})",
            self.delv,
            self.deld5,
            self.deld3,
            self.delj,
            self.insvd.get_string(),
            self.insdj.get_string()
        )
    }
}

impl StaticEvent {
    pub fn to_sequence(&self, m: &Model) -> Dna {
        let seq_v: &Dna = m.seg_vs[self.v_index].seq_with_pal.as_ref().unwrap();
        let seq_j: &Dna = m.seg_js[self.j_index].seq_with_pal.as_ref().unwrap();
        let seq_d: &Dna = m.seg_ds[self.d_index].seq_with_pal.as_ref().unwrap();
        let mut seq: Dna = Dna::new();
        seq.extend(&seq_v.extract_subsequence(0, seq_v.len() - self.delv));
        seq.extend(&self.insvd);
        seq.extend(&seq_d.extract_subsequence(self.deld5, seq_d.len() - self.deld3));
        seq.extend(&self.insdj);
        seq.extend(&seq_j.extract_subsequence(self.delj, seq_j.len()));

        // add errors
        for (ii, nuc) in &self.errors {
            seq.seq[*ii] = *nuc
        }

        seq
    }

    pub fn extract_cdr3(&self, full_sequence: &Dna, m: &Model) -> Dna {
        let vg = &m.seg_vs[self.v_index];
        let jg = &m.seg_js[self.j_index];
        let mut end_cdr3 = full_sequence.len() - jg.seq.len() + jg.cdr3_pos.unwrap() + 3;
        if vg.cdr3_pos.unwrap() > end_cdr3 {
            end_cdr3 = vg.cdr3_pos.unwrap(); // if we cut too much we return the empty sequence
        }
        full_sequence.extract_subsequence(vg.cdr3_pos.unwrap(), end_cdr3)
    }

    pub fn to_cdr3(&self, m: &Model) -> Dna {
        let full_sequence = self.to_sequence(m);
        self.extract_cdr3(&full_sequence, m)
    }
}

impl Event<'_> {
    pub fn to_static(&self, insvd: Dna, insdj: Dna) -> Result<StaticEvent> {
        Ok(StaticEvent {
            v_index: self
                .v
                .ok_or(anyhow!("Can't move that event to static"))?
                .index,
            v_start_gene: self
                .v
                .ok_or(anyhow!("Can't move that event to static"))?
                .start_gene,
            delv: self.delv,
            j_index: self
                .j
                .ok_or(anyhow!("Can't move that event to static"))?
                .index,
            j_start_seq: self
                .j
                .ok_or(anyhow!("Can't move that event to static"))?
                .start_seq,
            delj: self.delj,
            d_index: self
                .d
                .ok_or(anyhow!("Can't move that event to static"))?
                .index,
            d_start_seq: self
                .d
                .ok_or(anyhow!("Can't move that event to static"))?
                .pos,
            deld3: self.deld3,
            deld5: self.deld5,
            insvd,
            insdj,
            ..Default::default()
        })
    }
}