neofold 0.1.0

NeoFold Program, Program for Maximum-expected-accuracy Estimations of RNA Secondary Structures with Incorporating Homologous-RNA sequences
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
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

extern crate rnafamprob;

pub use rnafamprob::*;

pub type Mea = Prob;
#[derive(Clone)]
pub struct MeaSs {
  pub bp_pos_pair_seqs_inside_pos_pairs: PosPairSeqsWithPosPairs,
  pub ea: Mea,
}
pub type Meas = Vec<Mea>;
pub type MeaMat = HashMap<PosPair, Mea, Hasher>;
pub type Poss = Vec<Pos>;
pub type PosSeqsWithPoss = HashMap<Pos, Poss, Hasher>;
pub type PosPairs = Vec<PosPair>;
pub type PosPairSeqsWithPosPairs = HashMap<PosPair, PosPairs, Hasher>;
pub type MeaSsChar = u8;
pub type MeaSsStr = Vec<MeaSsChar>;

impl MeaSs {
  pub fn new() -> MeaSs {
    MeaSs {
      bp_pos_pair_seqs_inside_pos_pairs: PosPairSeqsWithPosPairs::default(),
      ea: 0.,
    }
  }
}

pub const UNPAIRING_BASE: MeaSsChar = '.' as MeaSsChar;
pub const BASE_PAIRING_LEFT_BASE: MeaSsChar = '(' as MeaSsChar;
pub const BASE_PAIRING_RIGHT_BASE: MeaSsChar = ')' as MeaSsChar;

pub fn neofold(bpp_mat: &SparseProbMat, upp_mat: &Probs, seq_len: usize, gamma: Prob) -> MeaSs {
  let mut mea_mat_4_bp_pos_pairs = MeaMat::default();
  let mut pos_seqs_with_poss_4_forward_bps = PosSeqsWithPoss::default();
  for sub_seq_len in 2 .. seq_len + 3 {
    for i in 0 .. seq_len + 3 - sub_seq_len {
      let j = i + sub_seq_len - 1;
      let pos_pair = (i, j);
      match bpp_mat.get(&pos_pair) {
        Some(&bpp) => {
          let meas_4_bp_pos_pair = get_meas_4_bp_pos_pair(&pos_pair, &mea_mat_4_bp_pos_pairs, &pos_seqs_with_poss_4_forward_bps, upp_mat);
          mea_mat_4_bp_pos_pairs.insert(pos_pair, meas_4_bp_pos_pair[j - i - 1] + gamma * bpp);
          let poss_exist = match pos_seqs_with_poss_4_forward_bps.get(&j) {
            Some(_) => {true},
            None => {false},
          };
          if poss_exist {
            pos_seqs_with_poss_4_forward_bps.get_mut(&j).expect("Failed to get an element of a hash map.").push(i);
          } else {
            pos_seqs_with_poss_4_forward_bps.insert(j, vec![i]);
          }
        },
        None => {},
      }
    }
  }
  let mut mea_ss = MeaSs::new();
  let pseudo_pos_pair = (0, seq_len + 1);
  let mut pos_pair_stack = vec![pseudo_pos_pair];
  while pos_pair_stack.len() > 0 {
    let pos_pair_1 = pos_pair_stack.pop().expect("Failed to pop an element of a vector.");
    let meas_4_bp_pos_pair = get_meas_4_bp_pos_pair(&pos_pair_1, &mea_mat_4_bp_pos_pairs, &pos_seqs_with_poss_4_forward_bps, upp_mat);
    let (i, j) = pos_pair_1;
    let mea = meas_4_bp_pos_pair[j - i - 1];
    if mea == 0. {continue;}
    let mut n = j - 1;
    while meas_4_bp_pos_pair[n - i] > 0. {
      let mea = meas_4_bp_pos_pair[n - i];
      if mea == meas_4_bp_pos_pair[n - i - 1] + upp_mat[n] {
        n = n - 1;
      } else {
        match pos_seqs_with_poss_4_forward_bps.get(&n) {
          Some(poss) => {
            for &m in poss {
              if m <= i {continue;}
              let pos_pair_2 = (m, n);
              if mea == meas_4_bp_pos_pair[m - i - 1] + mea_mat_4_bp_pos_pairs[&pos_pair_2] {
                let bp_pos_pairs_exist = match mea_ss.bp_pos_pair_seqs_inside_pos_pairs.get(&pos_pair_1) {
                  Some(_) => {true},
                  None => {false},
                };
                if bp_pos_pairs_exist {
                  mea_ss.bp_pos_pair_seqs_inside_pos_pairs.get_mut(&pos_pair_1).expect("Failed to get an element of a hash map.").push(pos_pair_2);
                } else {
                  mea_ss.bp_pos_pair_seqs_inside_pos_pairs.insert(pos_pair_1, vec![pos_pair_2]);
                }
                pos_pair_stack.push(pos_pair_2);
                n = m - 1;
                break;
              }
            }
          },
          None => {},
        }
      }
    }
  }
  mea_ss.ea = mea_mat_4_bp_pos_pairs[&pseudo_pos_pair];
  mea_ss
}

fn get_meas_4_bp_pos_pair(pos_pair: &PosPair, mea_mat_4_bp_pos_pairs: &MeaMat, pos_seqs_with_poss_4_forward_bps: &PosSeqsWithPoss, upp_mat: &Probs) -> Meas {
  let (i, j) = *pos_pair;
  let sub_seq_len = j - i + 1;
  let mut meas_4_bp_pos_pair = vec![0.; sub_seq_len - 1];
  for n in i + 2 .. j {
    let mut mea = 0.;
    match pos_seqs_with_poss_4_forward_bps.get(&n) {
      Some(poss) => {
        for &m in poss {
          if m <= i {continue;}
          let ea = meas_4_bp_pos_pair[m - i - 1] + mea_mat_4_bp_pos_pairs[&(m, n)];
          if ea > mea {
            mea = ea;
          }
        }
      },
      None => {},
    };
    let ea = meas_4_bp_pos_pair[n - i - 1] + upp_mat[n];
    if ea > mea {
      mea = ea;
    }
    meas_4_bp_pos_pair[n - i] = mea;
  }
  meas_4_bp_pos_pair
}