m

src/lib.rs

@@ -261,27 +261,34 @@ impl SequencesGraph {
         println!("{}", String::from_utf8_lossy(&cons) );
     }
 
+    // Stack reads respecting their offsets
+    // Return the consensus sequence
     pub fn get_stretch_consensus (&self, stretch_id: usize) -> Vec<u8> {
         let v = self.get_stretch_align(stretch_id);
-        let mut ov: Vec<Vec<u8>> = Vec::with_capacity(v.len());
-        let mut res: Vec<u8> = Vec::new();
+
+        let mut reads_stack: Vec<Vec<u8>> = Vec::with_capacity(v.len());
+        let mut consensus_sequence: Vec<u8> = Vec::new();
         if v.len() > 0 {
-            let decal = -v[0].1;
+            // first read position will be used as reference
+            let first_gap = -v[0].1;
             let mut max_len = 0;
+            
+            // add each read to the stack with n offset empty space 
             for (node_index, offset) in v {
                 let seq = self.sequences.get(&node_index).unwrap();
-                let seq = format!("{}{}", " ".repeat((offset + decal) as usize).to_string(), seq).as_bytes().to_vec();
+                let seq = format!("{}{}", " ".repeat((offset + first_gap) as usize).to_string(), seq).as_bytes().to_vec();
+                // keep track of the size of reads stack
                 if max_len < seq.len() { max_len = seq.len(); }
-                ov.push(seq);
+                reads_stack.push(seq);
             }
 
+            // for each column of the reads stack
             for i in 0..max_len {
                 let mut acc: Vec<u8> = Vec::new();
                 let atcg = vec![b'A', b'T', b'C', b'G'];
                 let mut n_atcg = [0; 4];
-
-                ov.iter().for_each(|s|{
-                    if s.len() > i {
+                reads_stack.iter().for_each(|s|{
+                    if s.len() >= i {
                         acc.push(s[i].clone());
                         for nb in 0..4 {
                             if s[i] == atcg[nb] {
@@ -296,12 +303,12 @@ impl SequencesGraph {
                 if *max_base as f32 / n_base as f32 > 0.5 {
                     n_atcg.iter().enumerate().for_each(|(pos, n)| if n == max_base { base = atcg[pos] });
                 }
-                res.push(base);
+                consensus_sequence.push(base);
             }
         }
 
-        println!("RES {:?}", res);
-        res
+        println!("RES {:?}", consensus_sequence);
+        consensus_sequence
     }
 
     pub fn get_fasta (&self) -> String {