Add in combinatorial code

Cargo.toml

@@ -10,6 +10,7 @@ benchmark=[]
 cached = "0.53.1"
 chrono = "0.4.38"
 clap = { version = "4.5.8", features = ["derive"] }
+csv = "1.3.0"
 fallible-streaming-iterator = "0.1.9"
 include_dir = "0.7.4"
 intervaltree = "0.2.7"

src/exports/gfa.rs

@@ -56,7 +56,18 @@ pub fn export_gfa(
     }
 
     let mut edges = edge_set.into_iter().collect();
-    let (blocks, boundary_edges) = Edge::blocks_from_edges(conn, &edges);
+    let blocks = Edge::blocks_from_edges(conn, &edges);
+    let node_blocks_by_id: HashMap<i64, Vec<GroupBlock>> =
+        blocks
+            .clone()
+            .into_iter()
+            .fold(HashMap::new(), |mut acc, block| {
+                acc.entry(block.node_id)
+                    .and_modify(|blocks| blocks.push(block.clone()))
+                    .or_insert_with(|| vec![block.clone()]);
+                acc
+            });
+    let boundary_edges = Edge::boundary_edges_from_sequences(&node_blocks_by_id);
     edges.extend(boundary_edges.clone());
 
     let (graph, edges_by_node_pair) = Edge::build_graph(&edges, &blocks);

src/graph.rs

@@ -40,3 +40,125 @@ where
         None
     })
 }
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use petgraph::graphmap::DiGraphMap;
+    use std::collections::HashSet;
+
+    #[test]
+    fn test_path_graph() {
+        let mut graph: DiGraphMap<i64, ()> = DiGraphMap::new();
+        graph.add_node(1);
+        graph.add_node(2);
+        graph.add_node(3);
+
+        graph.add_edge(1, 2, ());
+        graph.add_edge(2, 3, ());
+
+        let paths = all_simple_paths(&graph, 1, 3).collect::<Vec<Vec<i64>>>();
+        assert_eq!(paths.len(), 1);
+        let path = paths.first().unwrap().clone();
+        assert_eq!(path, vec![1, 2, 3]);
+    }
+
+    #[test]
+    fn test_two_path_graph() {
+        let mut graph: DiGraphMap<i64, ()> = DiGraphMap::new();
+        graph.add_node(1);
+        graph.add_node(2);
+        graph.add_node(3);
+        graph.add_node(4);
+
+        graph.add_edge(1, 2, ());
+        graph.add_edge(1, 3, ());
+        graph.add_edge(2, 4, ());
+        graph.add_edge(3, 4, ());
+
+        let paths = all_simple_paths(&graph, 1, 4).collect::<Vec<Vec<i64>>>();
+        assert_eq!(paths.len(), 2);
+        assert_eq!(
+            HashSet::<Vec<i64>>::from_iter::<Vec<Vec<i64>>>(paths),
+            HashSet::from_iter(vec![vec![1, 2, 4], vec![1, 3, 4]])
+        );
+    }
+
+    #[test]
+    fn test_two_by_two_combinatorial_graph() {
+        let mut graph: DiGraphMap<i64, ()> = DiGraphMap::new();
+        graph.add_node(1);
+        graph.add_node(2);
+        graph.add_node(3);
+        graph.add_node(4);
+        graph.add_node(5);
+        graph.add_node(6);
+        graph.add_node(7);
+
+        graph.add_edge(1, 2, ());
+        graph.add_edge(1, 3, ());
+        graph.add_edge(2, 4, ());
+        graph.add_edge(3, 4, ());
+        graph.add_edge(4, 5, ());
+        graph.add_edge(4, 6, ());
+        graph.add_edge(5, 7, ());
+        graph.add_edge(6, 7, ());
+
+        let paths = all_simple_paths(&graph, 1, 7).collect::<Vec<Vec<i64>>>();
+        assert_eq!(paths.len(), 4);
+        assert_eq!(
+            HashSet::<Vec<i64>>::from_iter::<Vec<Vec<i64>>>(paths),
+            HashSet::from_iter(vec![
+                vec![1, 2, 4, 5, 7],
+                vec![1, 3, 4, 5, 7],
+                vec![1, 2, 4, 6, 7],
+                vec![1, 3, 4, 6, 7]
+            ])
+        );
+    }
+
+    #[test]
+    fn test_three_by_three_combinatorial_graph() {
+        let mut graph: DiGraphMap<i64, ()> = DiGraphMap::new();
+        graph.add_node(1);
+        graph.add_node(2);
+        graph.add_node(3);
+        graph.add_node(4);
+        graph.add_node(5);
+        graph.add_node(6);
+        graph.add_node(7);
+        graph.add_node(8);
+        graph.add_node(9);
+
+        graph.add_edge(1, 2, ());
+        graph.add_edge(1, 3, ());
+        graph.add_edge(1, 4, ());
+        graph.add_edge(2, 5, ());
+        graph.add_edge(3, 5, ());
+        graph.add_edge(4, 5, ());
+        graph.add_edge(5, 6, ());
+        graph.add_edge(5, 7, ());
+        graph.add_edge(5, 8, ());
+        graph.add_edge(6, 9, ());
+        graph.add_edge(7, 9, ());
+        graph.add_edge(8, 9, ());
+
+        let paths = all_simple_paths(&graph, 1, 9).collect::<Vec<Vec<i64>>>();
+        assert_eq!(paths.len(), 9);
+        let expected_paths = vec![
+            vec![1, 2, 5, 6, 9],
+            vec![1, 3, 5, 6, 9],
+            vec![1, 4, 5, 6, 9],
+            vec![1, 2, 5, 7, 9],
+            vec![1, 3, 5, 7, 9],
+            vec![1, 4, 5, 7, 9],
+            vec![1, 2, 5, 8, 9],
+            vec![1, 3, 5, 8, 9],
+            vec![1, 4, 5, 8, 9],
+        ];
+        assert_eq!(
+            HashSet::<Vec<i64>>::from_iter::<Vec<Vec<i64>>>(paths),
+            HashSet::from_iter(expected_paths)
+        );
+    }
+}

src/main.rs

@@ -10,6 +10,7 @@ use gen::imports::gfa::import_gfa;
 use gen::models::metadata;
 use gen::models::operations::{setup_db, Branch, OperationState};
 use gen::operation_management;
+use gen::updates::library::update_with_library;
 use gen::updates::vcf::update_with_vcf;
 use rusqlite::{types::Value, Connection};
 use std::fmt::Debug;
@@ -60,13 +61,28 @@ enum Commands {
         fasta: Option<String>,
         /// A VCF file to incorporate
         #[arg(short, long)]
-        vcf: String,
+        vcf: Option<String>,
         /// If no genotype is provided, enter the genotype to assign variants
         #[arg(short, long)]
         genotype: Option<String>,
         /// If no sample is provided, enter the sample to associate variants to
         #[arg(short, long)]
         sample: Option<String>,
+        /// A CSV with combinatorial library information
+        #[arg(short, long)]
+        library: Option<String>,
+        /// A fasta with the combinatorial library parts
+        #[arg(short, long)]
+        parts: Option<String>,
+        /// The name of the path to add the library to
+        #[arg(short, long)]
+        path_name: Option<String>,
+        /// The start coordinate for the region to add the library to
+        #[arg(short, long)]
+        start: Option<i64>,
+        /// The end coordinate for the region to add the library to
+        #[arg(short, long)]
+        end: Option<i64>,
     },
     /// Initialize a gen repository
     Init {},
@@ -214,22 +230,40 @@ fn main() {
             name,
             fasta,
             vcf,
+            library,
+            parts,
             genotype,
             sample,
+            path_name,
+            start,
+            end,
         }) => {
             conn.execute("BEGIN TRANSACTION", []).unwrap();
             let name = &name.clone().unwrap_or_else(|| {
                 get_default_collection(&operation_conn)
                     .expect("No collection specified and default not setup.")
             });
-            update_with_vcf(
-                vcf,
-                name,
-                genotype.clone().unwrap_or("".to_string()),
-                sample.clone().unwrap_or("".to_string()),
-                &conn,
-                &operation_conn,
-            );
+            if let Some(library_path) = library {
+                update_with_library(
+                    &conn,
+                    &operation_conn,
+                    name,
+                    &path_name.clone().unwrap(),
+                    start.unwrap(),
+                    end.unwrap(),
+                    &parts.clone().unwrap(),
+                    library_path,
+                );
+            } else {
+                update_with_vcf(
+                    &vcf.clone().unwrap(),
+                    name,
+                    genotype.clone().unwrap_or("".to_string()),
+                    sample.clone().unwrap_or("".to_string()),
+                    &conn,
+                    &operation_conn,
+                );
+            }
 
             conn.execute("END TRANSACTION", []).unwrap();
         }

src/models/block_group.rs

@@ -272,8 +272,20 @@ impl BlockGroup {
 
     pub fn get_all_sequences(conn: &Connection, block_group_id: i64) -> HashSet<String> {
         let mut edges = BlockGroupEdge::edges_for_block_group(conn, block_group_id);
-        let (blocks, boundary_edges) = Edge::blocks_from_edges(conn, &edges);
+        let blocks = Edge::blocks_from_edges(conn, &edges);
+        let node_blocks_by_id: HashMap<i64, Vec<GroupBlock>> =
+            blocks
+                .clone()
+                .into_iter()
+                .fold(HashMap::new(), |mut acc, block| {
+                    acc.entry(block.node_id)
+                        .and_modify(|blocks| blocks.push(block.clone()))
+                        .or_insert_with(|| vec![block.clone()]);
+                    acc
+                });
+        let boundary_edges = Edge::boundary_edges_from_sequences(&node_blocks_by_id);
         edges.extend(boundary_edges.clone());
+
         let (graph, _) = Edge::build_graph(&edges, &blocks);
 
         let mut start_nodes = vec![];
@@ -288,7 +300,6 @@ impl BlockGroup {
                 end_nodes.push(node);
             }
         }
-
         let blocks_by_id = blocks
             .clone()
             .into_iter()

src/models/edge.rs

@@ -320,7 +320,7 @@ impl Edge {
             .collect::<Vec<i64>>()
     }
 
-    pub fn blocks_from_edges(conn: &Connection, edges: &Vec<Edge>) -> (Vec<GroupBlock>, Vec<Edge>) {
+    pub fn blocks_from_edges(conn: &Connection, edges: &Vec<Edge>) -> Vec<GroupBlock> {
         let mut node_ids = HashSet::new();
         let mut edges_by_source_node_id: HashMap<i64, Vec<&Edge>> = HashMap::new();
         let mut edges_by_target_node_id: HashMap<i64, Vec<&Edge>> = HashMap::new();
@@ -346,7 +346,6 @@ impl Edge {
 
         let mut blocks = vec![];
         let mut block_index = 0;
-        let mut boundary_edges = vec![];
         // we sort by keys to exploit the external sequence cache which keeps the most recently used
         // external sequence in memory.
         for (node_id, sequence) in sequences_by_node_id
@@ -358,21 +357,6 @@ impl Edge {
                 edges_by_target_node_id.get(node_id),
                 sequence.length,
             );
-            for block_boundary in &block_boundaries {
-                // NOTE: Most of this data is bogus, the Edge struct is just a convenient wrapper
-                // for the data we need to set up boundary edges in the block group graph
-                boundary_edges.push(Edge {
-                    id: -1,
-                    source_node_id: *node_id,
-                    source_coordinate: *block_boundary,
-                    source_strand: Strand::Forward,
-                    target_node_id: *node_id,
-                    target_coordinate: *block_boundary,
-                    target_strand: Strand::Forward,
-                    chromosome_index: 0,
-                    phased: 0,
-                });
-            }
 
             if !block_boundaries.is_empty() {
                 let start = 0;
@@ -422,7 +406,7 @@ impl Edge {
             0,
         );
         blocks.push(end_block);
-        (blocks, boundary_edges)
+        blocks
     }
 
     pub fn build_graph(
@@ -483,6 +467,30 @@ impl Edge {
 
         (graph, edges_by_node_pair)
     }
+
+    pub fn boundary_edges_from_sequences(
+        node_blocks_by_id: &HashMap<i64, Vec<GroupBlock>>,
+    ) -> Vec<Edge> {
+        let mut boundary_edges = vec![];
+        for node_blocks in node_blocks_by_id.values() {
+            for (previous_block, next_block) in node_blocks.iter().tuple_windows() {
+                // NOTE: Most of this data is bogus, the Edge struct is just a convenient wrapper
+                // for the data we need to set up boundary edges in the block group graph
+                boundary_edges.push(Edge {
+                    id: -1,
+                    source_node_id: previous_block.node_id,
+                    source_coordinate: previous_block.end,
+                    source_strand: Strand::Forward,
+                    target_node_id: next_block.node_id,
+                    target_coordinate: next_block.start,
+                    target_strand: Strand::Forward,
+                    chromosome_index: 0,
+                    phased: 0,
+                });
+            }
+        }
+        boundary_edges
+    }
 }
 
 #[cfg(test)]