delamerde

src/lib.rs

@@ -1222,15 +1222,17 @@ impl MyBreakPoint {
         left_mapping.target_name = target_name.clone();
         left_mapping.target_start = pos;
         left_mapping.target_end = pos;
-        
+
         let mut right_mapping = Mapping::default();
         right_mapping.strand = strand;
         right_mapping.target_name = target_name;
         right_mapping.target_start = pos;
         right_mapping.target_end = pos;
 
-        Self { id: id.to_string(), mappings: vec![left_mapping, right_mapping] }
-
+        Self {
+            id: id.to_string(),
+            mappings: vec![left_mapping, right_mapping],
+        }
     }
 
     pub fn sens(&self) -> Self {
@@ -1347,7 +1349,10 @@ impl MyBreakPoint {
 
 #[cfg(test)]
 mod tests {
-    use std::ops::{Range, RangeInclusive};
+    use std::{
+        collections::HashSet,
+        ops::{Range, RangeInclusive},
+    };
 
     use petgraph::{stable_graph::StableGraph, visit::EdgeRef, Directed};
 
@@ -1398,11 +1403,15 @@ mod tests {
         let mut graph = StableGraph::<MyBreakPoint, String, Directed>::new();
 
         // add bp to graph
+        let mut dedup_hash = HashSet::new();
         for contig in genome.contigs() {
             if let Some(mappings) = contig.contig_ref.breakpoints() {
                 let bp = MyBreakPoint::new(mappings);
-                let id = graph.add_node(bp.clone());
-                breakpoints_id.insert(bp.id, id);
+                if !dedup_hash.contains(&bp.str_sens(sens_symb, rc_symb)) {
+                    let id = graph.add_node(bp.clone());
+                    breakpoints_id.insert(bp.id.clone(), id);
+                    dedup_hash.insert(bp.str_sens(sens_symb, rc_symb));
+                }
             }
         }
 
@@ -1415,8 +1424,16 @@ mod tests {
         let g_iter = graph.clone();
         for last in g_iter.node_weights() {
             for next in g_iter.node_weights() {
-                let a = breakpoints_id.get(&last.id).unwrap();
-                let b = breakpoints_id.get(&next.id).unwrap();
+                if last == next {
+                    continue;
+                }
+                let a = *breakpoints_id.get(&last.id).unwrap();
+                let b = *breakpoints_id.get(&next.id).unwrap();
+
+                if a == start_id && b == end_id || a == end_id && b == start_id {
+                    continue;
+                }
+
 
                 let (_, (last_tail_sens, last_tail_chr, last_tail_pos)) = last.repr_sens();
                 let ((next_head_sens, next_head_chr, next_head_pos), _) = next.repr_sens();
@@ -1427,8 +1444,13 @@ mod tests {
                     } else {
                         last_tail_pos > next_head_pos
                     };
-                    if in_order {
-                        graph.add_edge(*a, *b, "".to_string());
+                    if in_order && graph.find_edge(a, b).is_none() {
+                        println!(
+                            "{}<-->{}",
+                            last.str_sens(sens_symb, rc_symb),
+                            next.str_sens(sens_symb, rc_symb)
+                        );
+                        graph.add_edge(a, b, "".to_string());
                     }
                 }
 
@@ -1440,34 +1462,105 @@ mod tests {
                     } else {
                         last_tail_pos > next_head_pos
                     };
-                    if in_order {
-                        graph.add_edge(*b, *a, "".to_string());
+                    if in_order && graph.find_edge(a, b).is_none() {
+                        println!(
+                            "{}<-->{}",
+                            last.str_sens(sens_symb, rc_symb),
+                            next.str_rc(sens_symb, rc_symb)
+                        );
+                        graph.add_edge(a, b, "".to_string());
+                    }
+                }
+
+                let (_, (last_tail_sens, last_tail_chr, last_tail_pos)) = last.repr_rc();
+                let ((next_head_sens, next_head_chr, next_head_pos), _) = next.repr_sens();
+                if last_tail_chr == next_head_chr && last_tail_sens == next_head_sens {
+                    let in_order = if last_tail_sens {
+                        last_tail_pos <= next_head_pos
+                    } else {
+                        last_tail_pos > next_head_pos
+                    };
+                    if in_order && graph.find_edge(b, a).is_none() {
+                        println!(
+                            "{}<-->{}",
+                            last.str_rc(sens_symb, rc_symb),
+                            next.str_sens(sens_symb, rc_symb)
+                        );
+                        graph.add_edge(b, a, "".to_string());
+                    }
+                }
+
+                let (_, (last_tail_sens, last_tail_chr, last_tail_pos)) = last.repr_rc();
+                let ((next_head_sens, next_head_chr, next_head_pos), _) = next.repr_rc();
+                if last_tail_chr == next_head_chr && last_tail_sens == next_head_sens {
+                    let in_order = if last_tail_sens {
+                        last_tail_pos <= next_head_pos
+                    } else {
+                        last_tail_pos > next_head_pos
+                    };
+                    if in_order  && graph.find_edge(b, a).is_none() {
+                        println!(
+                            "{}<-->{}",
+                            last.str_rc(sens_symb, rc_symb),
+                            next.str_rc(sens_symb, rc_symb)
+                        );
+                        graph.add_edge(b, a, "".to_string());
                     }
                 }
+
+
             }
         }
 
         // Remove edges in feedback arc set.
-        let fas: Vec<_> = petgraph::algo::feedback_arc_set::greedy_feedback_arc_set(&graph)
-            .map(|e| e.id())
-            .collect();
-        for edge_id in fas {
-            graph.remove_edge(edge_id);
-        }
-
+        // let fas: Vec<_> = petgraph::algo::feedback_arc_set::greedy_feedback_arc_set(&graph)
+        //     .map(|e| e.id())
+        //     .collect();
+        // for edge_id in fas {
+        //     graph.remove_edge(edge_id);
+        // }
+        //
         // remove empty nodes
-        let del: Vec<_> = graph.node_indices().filter(|id| graph.neighbors_undirected(*id).count() == 0).collect();
-        del.iter().for_each(|id| { graph.remove_node(*id); });
+        let del: Vec<_> = graph
+            .node_indices()
+            .filter(|id| graph.neighbors_undirected(*id).count() == 0)
+            .collect();
+        del.iter().for_each(|id| {
+            graph.remove_node(*id);
+        });
 
-        
         // Find all the ways from start node to end node.
-        let ways =
-            algo::all_simple_paths::<Vec<_>, _>(&graph, start_id, end_id, 0, None)
-                .collect::<Vec<_>>();
+        let ways = algo::all_simple_paths::<Vec<_>, _>(&graph, start_id, end_id, 0, None)
+            .collect::<Vec<_>>();
 
 
+
+        let mut confirmed_way = Vec::new();
+        for way in ways.iter().enumerate() {
+            let mut last_sens = true;
+            let mut confirmed = true;
+            // let mut last_sens = s.1.0.clone();
+            for window in way.windows(2) {
+                let a = graph.node_weight(window[0]).unwrap();
+                let b =  graph.node_weight(window[1]).unwrap();
+
+                let last = if last_sens { a.clone() } else { a.rc() };
+                if last.is_next(&b) {
+                    last_sens = true;
+                } else if last.is_next(&b.rc()) {
+                    last_sens = false;
+                } else {
+                    confirmed = false;
+                    break;
+                }
+            }
+            if confirmed {
+                confirmed_way.push(way.clone());
+            }
+        }
+
         // Add way number to edges.
-        for (n_way, way) in ways.iter().enumerate() {
+        for (n_way, way) in confirmed_way.iter().enumerate() {
             for window in way.windows(2) {
                 let edge_id = graph.find_edge(window[0], window[1]).unwrap();
                 let edge = graph.edge_weight_mut(edge_id).unwrap();