Alternative to blocks

Cargo.toml

@@ -7,6 +7,8 @@ edition = "2021"
 bio = "2.0.0"
 clap = { version = "4.5.8", features = ["derive"] }
 include_dir = "0.7.4"
+intervaltree = "0.2.7"
+itertools = "0.13.0"
 rusqlite = { version = "0.31.0", features = ["bundled", "array"] }
 rusqlite_migration = { version = "1.2.0" , features = ["from-directory"]}
 sha2 = "0.10.8"

migrations/01-initial/up.sql

@@ -9,7 +9,7 @@ CREATE TABLE sample (
 CREATE TABLE sequence (
   hash TEXT PRIMARY KEY NOT NULL,
   sequence_type TEXT NOT NULL,
-  sequence TEXT,
+  sequence TEXT NOT NULL,
   "length" INTEGER NOT NULL
 );
 
@@ -79,4 +79,41 @@ CREATE TABLE change_log (
   FOREIGN KEY(path_id) REFERENCES path(id),
   FOREIGN KEY(sequence_hash) REFERENCES sequence(hash)
 );
-CREATE UNIQUE INDEX change_log_uidx ON change_log(hash);
+CREATE UNIQUE INDEX change_log_uidx ON change_log(hash);
+
+CREATE TABLE new_edges (
+  id INTEGER PRIMARY KEY NOT NULL,
+  source_hash TEXT NOT NULL,
+  source_coordinate INTEGER NOT NULL,
+  source_strand TEXT NOT NULL,
+  target_hash TEXT NOT NULL,
+  target_coordinate INTEGER NOT NULL,
+  target_strand TEXT NOT NULL,
+  chromosome_index INTEGER NOT NULL,
+  phased INTEGER NOT NULL,
+  FOREIGN KEY(source_hash) REFERENCES sequence(hash),
+  FOREIGN KEY(target_hash) REFERENCES sequence(hash),
+  constraint chk_phased check (phased in (0, 1))
+);
+CREATE UNIQUE INDEX new_edge_uidx ON new_edges(source_hash, source_coordinate, source_strand, target_hash, target_coordinate, target_strand, chromosome_index, phased);
+
+CREATE TABLE path_edges (
+  id INTEGER PRIMARY KEY NOT NULL,
+  path_id INTEGER NOT NULL,
+  index_in_path INTEGER NOT NULL,
+  edge_id INTEGER NOT NULL,
+  FOREIGN KEY(edge_id) REFERENCES new_edges(id),
+  FOREIGN KEY(path_id) REFERENCES path(id)
+);
+CREATE UNIQUE INDEX path_edges_uidx ON path_edges(path_id, edge_id);
+
+CREATE TABLE block_group_edges (
+  id INTEGER PRIMARY KEY NOT NULL,
+  block_group_id INTEGER NOT NULL,
+  edge_id INTEGER NOT NULL,
+  FOREIGN KEY(block_group_id) REFERENCES block_group(id),
+  FOREIGN KEY(edge_id) REFERENCES new_edges(id)
+);
+CREATE UNIQUE INDEX block_group_edges_uidx ON block_group_edges(block_group_id, edge_id);
+
+INSERT INTO sequence (hash, sequence_type, sequence, "length") values ("start-node-yyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyy", "OTHER", "start-node-yyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyy", 64), ("end-node-zzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz", "OTHER", "end-node-zzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz", 64);

src/main.rs

@@ -6,7 +6,17 @@ use std::path::PathBuf;
 
 use bio::io::fasta;
 use gen::migrations::run_migrations;
-use gen::models::{self, block::Block, edge::Edge, path::Path, sequence::Sequence, BlockGroup};
+use gen::models::{
+    self,
+    block::Block,
+    block_group_edge::BlockGroupEdge,
+    edge::Edge,
+    new_edge::NewEdge,
+    path::{NewBlock, Path},
+    path_edge::PathEdge,
+    sequence::Sequence,
+    BlockGroup,
+};
 use gen::{get_connection, parse_genotype};
 use noodles::vcf;
 use noodles::vcf::variant::record::samples::series::value::genotype::Phasing;
@@ -94,6 +104,54 @@ fn import_fasta(fasta: &String, name: &str, shallow: bool, conn: &mut Connection
     }
 }
 
+fn new_import_fasta(fasta: &String, name: &str, shallow: bool, conn: &mut Connection) {
+    // TODO: support gz
+    let mut reader = fasta::Reader::from_file(fasta).unwrap();
+
+    if !models::Collection::exists(conn, name) {
+        let collection = models::Collection::create(conn, name);
+
+        for result in reader.records() {
+            let record = result.expect("Error during fasta record parsing");
+            let sequence = String::from_utf8(record.seq().to_vec()).unwrap();
+            let seq_hash = Sequence::create(conn, "DNA", &sequence, !shallow);
+            let block_group = BlockGroup::create(conn, &collection.name, None, record.id());
+            let edge_into = NewEdge::create(
+                conn,
+                NewEdge::PATH_START_HASH.to_string(),
+                0,
+                "+".to_string(),
+                seq_hash.to_string(),
+                0,
+                "+".to_string(),
+                0,
+                0,
+            );
+            let edge_out_of = NewEdge::create(
+                conn,
+                seq_hash.to_string(),
+                sequence.len() as i32,
+                "+".to_string(),
+                NewEdge::PATH_END_HASH.to_string(),
+                0,
+                "+".to_string(),
+                0,
+                0,
+            );
+            BlockGroupEdge::bulk_create(conn, block_group.id, vec![edge_into.id, edge_out_of.id]);
+            Path::new_create(
+                conn,
+                record.id(),
+                block_group.id,
+                vec![edge_into.id, edge_out_of.id],
+            );
+        }
+        println!("Created it");
+    } else {
+        println!("Collection {:1} already exists", name);
+    }
+}
+
 fn update_with_vcf(
     vcf_path: &String,
     collection_name: &String,
@@ -233,6 +291,156 @@ fn update_with_vcf(
     }
 }
 
+fn new_update_with_vcf(
+    vcf_path: &String,
+    collection_name: &String,
+    fixed_genotype: String,
+    fixed_sample: String,
+    conn: &mut Connection,
+) {
+    run_migrations(conn);
+
+    let mut reader = vcf::io::reader::Builder::default()
+        .build_from_path(vcf_path)
+        .expect("Unable to parse");
+    let header = reader.read_header().unwrap();
+    let sample_names = header.sample_names();
+    for name in sample_names {
+        models::Sample::create(conn, name);
+    }
+    if !fixed_sample.is_empty() {
+        models::Sample::create(conn, &fixed_sample);
+    }
+    let mut genotype = vec![];
+    if !fixed_genotype.is_empty() {
+        genotype = parse_genotype(&fixed_genotype);
+    }
+
+    for result in reader.records() {
+        let record = result.unwrap();
+        let seq_name = record.reference_sequence_name().to_string();
+        let ref_allele = record.reference_bases();
+        // this converts the coordinates to be zero based, start inclusive, end exclusive
+        let ref_start = record.variant_start().unwrap().unwrap().get() - 1;
+        let ref_end = record.variant_end(&header).unwrap().get();
+        let alt_bases = record.alternate_bases();
+        let alt_alleles: Vec<_> = alt_bases.iter().collect::<io::Result<_>>().unwrap();
+        // TODO: fix this duplication of handling an insert
+        if !fixed_sample.is_empty() && !genotype.is_empty() {
+            for (chromosome_index, genotype) in genotype.iter().enumerate() {
+                if let Some(gt) = genotype {
+                    if gt.allele != 0 {
+                        let alt_seq = alt_alleles[chromosome_index - 1];
+                        let phased = match gt.phasing {
+                            Phasing::Phased => 1,
+                            Phasing::Unphased => 0,
+                        };
+                        // TODO: new sequence may not be real and be <DEL> or some sort. Handle these.
+                        let new_sequence_hash = Sequence::create(conn, "DNA", alt_seq, true);
+                        let sequence =
+                            Sequence::sequence_from_hash(conn, &new_sequence_hash).unwrap();
+                        let sample_bg_id = BlockGroup::get_or_create_sample_block_group(
+                            conn,
+                            collection_name,
+                            &fixed_sample,
+                            &seq_name,
+                        );
+                        let sample_paths = Path::get_paths(
+                            conn,
+                            "select * from path where block_group_id = ?1 AND name = ?2",
+                            vec![
+                                SQLValue::from(sample_bg_id),
+                                SQLValue::from(seq_name.clone()),
+                            ],
+                        );
+                        let new_block = NewBlock {
+                            id: 0,
+                            sequence: sequence.clone(),
+                            block_sequence: alt_seq.to_string(),
+                            sequence_start: 0,
+                            sequence_end: alt_seq.len() as i32,
+                            path_start: ref_start as i32,
+                            path_end: ref_end as i32,
+                            strand: "+".to_string(),
+                        };
+                        BlockGroup::new_insert_change(
+                            conn,
+                            sample_bg_id,
+                            &sample_paths[0],
+                            ref_start as i32,
+                            ref_end as i32,
+                            &new_block,
+                            chromosome_index as i32,
+                            phased,
+                        );
+                    }
+                }
+            }
+        } else {
+            for (sample_index, sample) in record.samples().iter().enumerate() {
+                let genotype = sample.get(&header, "GT");
+                if genotype.is_some() {
+                    if let Value::Genotype(genotypes) = genotype.unwrap().unwrap().unwrap() {
+                        for (chromosome_index, gt) in genotypes.iter().enumerate() {
+                            if gt.is_ok() {
+                                let (allele, phasing) = gt.unwrap();
+                                let phased = match phasing {
+                                    Phasing::Phased => 1,
+                                    Phasing::Unphased => 0,
+                                };
+                                let allele = allele.unwrap();
+                                if allele != 0 {
+                                    let alt_seq = alt_alleles[allele - 1];
+                                    // TODO: new sequence may not be real and be <DEL> or some sort. Handle these.
+                                    let new_sequence_hash =
+                                        Sequence::create(conn, "DNA", alt_seq, true);
+                                    let sequence =
+                                        Sequence::sequence_from_hash(conn, &new_sequence_hash)
+                                            .unwrap();
+                                    let sample_bg_id = BlockGroup::get_or_create_sample_block_group(
+                                        conn,
+                                        collection_name,
+                                        &sample_names[sample_index],
+                                        &seq_name,
+                                    );
+                                    let sample_paths = Path::get_paths(
+                                        conn,
+                                        "select * from path where block_group_id = ?1 AND name = ?2",
+                                        vec![
+                                            SQLValue::from(sample_bg_id),
+                                            SQLValue::from(seq_name.clone()),
+                                        ],
+                                    );
+                                    let new_block = NewBlock {
+                                        id: 0,
+                                        sequence: sequence.clone(),
+                                        block_sequence: alt_seq.to_string(),
+                                        sequence_start: 0,
+                                        sequence_end: alt_seq.len() as i32,
+                                        path_start: ref_start as i32,
+                                        path_end: ref_end as i32,
+                                        strand: "+".to_string(),
+                                    };
+                                    BlockGroup::new_insert_change(
+                                        conn,
+                                        sample_bg_id,
+                                        &sample_paths[0],
+                                        ref_start as i32,
+                                        ref_end as i32,
+                                        &new_block,
+                                        chromosome_index as i32,
+                                        phased,
+                                    );
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
 fn main() {
     let cli = Cli::parse();
 
@@ -242,15 +450,15 @@ fn main() {
             name,
             db,
             shallow,
-        }) => import_fasta(fasta, name, *shallow, &mut get_connection(db)),
+        }) => new_import_fasta(fasta, name, *shallow, &mut get_connection(db)),
         Some(Commands::Update {
             name,
             db,
             fasta,
             vcf,
             genotype,
             sample,
-        }) => update_with_vcf(
+        }) => new_update_with_vcf(
             vcf,
             name,
             genotype.clone().unwrap_or("".to_string()),