feat: port over code from worker (#1) (#2)

Author: holtgrewe

.gitattributes

@@ -0,0 +1 @@
+/tests/** filter=lfs diff=lfs merge=lfs -text

.github/workflows/release-please.yml

@@ -26,6 +26,11 @@ jobs:
           toolchain: stable
           override: true
 
+      - name: Setup protoc
+        uses: arduino/setup-protoc@v1.1.2
+        with:
+          repo-token: ${{ secrets.GITHUB_TOKEN }}
+
       - uses: Swatinem/rust-cache@v2
         if: ${{ steps.release.outputs.release_created }}
 

.github/workflows/rust.yml

@@ -22,6 +22,11 @@ jobs:
           override: true
           components: rustfmt
 
+      - name: Setup protoc
+        uses: arduino/setup-protoc@v1.1.2
+        with:
+          repo-token: ${{ secrets.GITHUB_TOKEN }}
+
       - name: Check format
         run: |
           cargo fmt -- --check
@@ -39,6 +44,11 @@ jobs:
           override: true
           components: clippy
 
+      - name: Setup protoc
+        uses: arduino/setup-protoc@v1.1.2
+        with:
+          repo-token: ${{ secrets.GITHUB_TOKEN }}
+
       - name: Lint with clippy
         uses: actions-rs/clippy-check@v1
         with:
@@ -59,6 +69,11 @@ jobs:
           toolchain: stable
           override: true
 
+      - name: Setup protoc
+        uses: arduino/setup-protoc@v1.1.2
+        with:
+          repo-token: ${{ secrets.GITHUB_TOKEN }}
+
       - uses: Swatinem/rust-cache@v2
 
       - name: Run cargo-tarpaulin

Cargo.toml

@@ -8,20 +8,76 @@ license = "Apache-2.0"
 homepage = "https://github.com/bihealth/viguno"
 readme = "README.md"
 
-[lib]
-name = "the_repo_name"
-path = "src/lib.rs"
-
 [dependencies]
+actix-files = "0.6"
+actix-web = "4.3"
+anyhow = "1.0"
+base16ct = "0.2"
+bgzip = "0.3"
+bio = "1.1"
+byte-unit = "4.0"
+chrono = { version = "0.4", features = ["serde"] }
+clap = { version = "4.1", features = ["derive", "help", "env"] }
+clap-verbosity-flag = "2.0"
+console = "0.15"
+csv = "1.1"
+derive_more = "0.99"
+enum-map = { version = "2.4", features = ["serde"] }
+env_logger = "0.10"
+fastrand = "1.9"
+flate2 = "1.0"
+hgvs = "0.8"
+hpo = "0.8"
+indexmap = { version = "1.9", features = ["serde"] }
+indicatif = { version = "0.17", features = ["rayon"] }
+itertools = "0.10"
+lazy_static = "1.4"
+linked-hash-map = { version = "0.5", features = ["serde", "serde_impl"] }
+log = "0.4"
+md-5 = "0.10"
+multimap = "0.9"
+noodles-bgzf = "0.22"
+noodles-core = "0.11"
+noodles-csi = "0.19"
+noodles-tabix = "0.22"
+noodles-vcf = "0.31"
+once_cell = "1.18"
+pretty_assertions = "1.3"
+procfs = "0.15"
+prost = "0.11"
+rayon = "1.7"
+regex = "1.7"
+result = "1.0"
+rocksdb = { version = "0.21", features = ["multi-threaded-cf"] }
+rocksdb-utils-lookup = "0.1"
+serde = { version = "1.0", features = ["serde_derive"] }
+serde-jsonlines = "0.4"
+serde_json = "1.0"
+serde_with = "3.0"
+sha2 = "0.10"
+shellexpand = "3.0"
+strum = "0.24"
+strum_macros = "0.24"
+temp_testdir = "0.2"
+tempdir = "0.3"
+test-log = "0.2"
 thiserror = "1.0"
+thousands = "0.2"
+toml = { version = "0.7", features = ["preserve_order", "parse", "display"] }
 tracing = "0.1"
+tracing-subscriber = "0.3"
+uuid = { version = "1.2", features = ["v4", "fast-rng", "serde"] }
+
+[build-dependencies]
+prost-build = "0.11"
 
 [dev-dependencies]
-anyhow = "1.0"
-clap = { version = "4.1", features = ["derive", "env"] }
-clap-verbosity-flag = {version = "2.0"}
-env_logger = "0.10"
+file_diff = "1.0"
+insta = { version = "1.29", features = ["yaml"] }
 pretty_assertions = "1.3"
+serde_test = "1.0"
 temp_testdir = "0.2"
-test-log = "0.2"
-tracing-subscriber = {version = "0.3" }
+
+[[bin]]
+name = "viguno"
+path = "src/main.rs"

build.rs

@@ -0,0 +1,13 @@
+// The custom build script, needed as we use protocolbuffers.
+
+fn main() {
+    prost_build::Config::new()
+        .protoc_arg("-Isrc/proto")
+        // Add serde serialization and deserialization to the generated code.
+        .type_attribute(".", "#[derive(serde::Serialize, serde::Deserialize)]")
+        // Skip serializing `None` values.
+        .type_attribute(".", "#[serde_with::skip_serializing_none]")
+        // Define the protobuf files to compile.
+        .compile_protos(&["viguno/v1/simulation.proto"], &["src/"])
+        .unwrap();
+}

src/algos/mod.rs

@@ -0,0 +1,3 @@
+//! Phenotype-related algorithms.
+
+pub mod phenomizer;

src/algos/phenomizer.rs

@@ -0,0 +1,91 @@
+//! Similarity computation using the Phenomizer method.
+
+use hpo::{
+    similarity::{Builtins, Similarity},
+    term::{HpoGroup, InformationContentKind},
+    Ontology,
+};
+
+// TODO: this is equivalent to `StandardCombiner::FunSimAvg` USE IT.
+
+/// Compute symmetric similarity score.
+pub fn score(q: &HpoGroup, d: &HpoGroup, o: &Ontology) -> f32 {
+    let s = Builtins::Resnik(InformationContentKind::Gene);
+    (score_dir(q, d, o, &s) + score_dir(d, q, o, &s)) / 2.0
+}
+
+/// "Directed" score part of phenomizer score.
+///
+/// # Panics
+///
+/// If there are more query terms than fit an `f32` value.
+fn score_dir(qs: &HpoGroup, ds: &HpoGroup, o: &Ontology, s: &impl Similarity) -> f32 {
+    // Handle case of empty `qs`.
+    if qs.is_empty() {
+        return 0f32;
+    }
+
+    // For each `q in qs` compute max similarity to any `d in ds`.
+    let mut tmp: Vec<f32> = Vec::new();
+    for q in qs {
+        if let Some(q) = o.hpo(q) {
+            tmp.push(
+                ds.iter()
+                    .filter_map(|d| o.hpo(d).map(|d| q.similarity_score(&d, s)))
+                    .max_by(|a, b| a.partial_cmp(b).expect("try to compare NaN"))
+                    .unwrap_or_default(),
+            );
+        }
+    }
+
+    // NB: we allow loss of precision in this function for the following statement.
+    let len: u16 = qs.len().try_into().expect("more than 2^16 query terms");
+    let len: f32 = len.try_into().expect("too many query terms for f32");
+    tmp.iter().sum::<f32>() / len
+}
+
+#[cfg(test)]
+mod test {
+    use super::*;
+    use hpo::{annotations::OmimDiseaseId, term::HpoGroup, HpoTermId, Ontology};
+
+    fn load_hpo() -> Result<Ontology, anyhow::Error> {
+        Ok(Ontology::from_standard("tests/data/hpo")?)
+    }
+
+    fn prepare(terms: &[&str]) -> HpoGroup {
+        HpoGroup::from(
+            terms
+                .iter()
+                .map(|s| HpoTermId::from((*s).to_string()))
+                .collect::<Vec<_>>(),
+        )
+    }
+
+    #[test]
+    fn phenomizer_score_gene() -> Result<(), anyhow::Error> {
+        let hpo = load_hpo()?;
+
+        let query = &[
+            // slender build
+            "HP:0001533",
+            // high, narrow palate
+            "HP:0002705",
+        ];
+        let omim_marfan = hpo
+            .omim_disease(&OmimDiseaseId::from(154_700))
+            .expect("marfan symdrome must be in HPO");
+        let hpo_marfan = omim_marfan
+            .to_hpo_set(&hpo)
+            .child_nodes()
+            .without_modifier()
+            .into_iter()
+            .collect::<HpoGroup>();
+
+        let score = score(&prepare(query), &hpo_marfan, &hpo);
+
+        assert!((score - 1.770_859_7).abs() < 0.00001, "score = {score}");
+
+        Ok(())
+    }
+}

src/common.rs

@@ -0,0 +1,69 @@
+//! Functionality shared between all commands.
+
+use clap::Parser;
+use clap_verbosity_flag::{InfoLevel, Verbosity};
+
+/// Shared command line arguments.
+#[derive(Parser, Debug)]
+pub struct Args {
+    /// Verbosity of the program
+    #[clap(flatten)]
+    pub verbose: Verbosity<InfoLevel>,
+}
+
+/// Construct the `indicatif` style for progress bars.
+///
+/// # Panics
+///
+/// In the case when writing the ETA seconds could not be written to the progress bar.
+pub fn indicatif_style() -> indicatif::ProgressStyle {
+    let tpl = "{spinner:.green} [{elapsed_precise}] [{wide_bar:.cyan/blue}] \
+    {human_pos}/{human_len} ({per_sec})";
+    indicatif::ProgressStyle::with_template(tpl)
+        .unwrap()
+        .with_key(
+            "eta",
+            |state: &indicatif::ProgressState, w: &mut dyn std::fmt::Write| {
+                write!(w, "{:.1}s", state.eta().as_secs_f64())
+                    .expect("could not write the ETA as seconds to progress bar");
+            },
+        )
+        .progress_chars("#>-")
+}
+
+/// Construct an `indicatif` progress bar with the common style.
+///
+/// Also, we will enable a steady tick every 0.1s and hide in tests.
+pub fn progress_bar(#[allow(unused_variables)] len: usize) -> indicatif::ProgressBar {
+    #[cfg(test)]
+    let pb = indicatif::ProgressBar::hidden();
+    #[cfg(not(test))]
+    let pb = indicatif::ProgressBar::new(len as u64).with_style(indicatif_style());
+    pb.enable_steady_tick(std::time::Duration::from_millis(100));
+    pb
+}
+
+/// Load HPO either from binary `$path_hpo/hpo.bin` if it exist, otherwise load as
+/// standard directory from `$path_hpo`.
+///
+/// # Errors
+///
+/// In the case of loading failure.
+pub fn load_hpo<P: AsRef<std::path::Path>>(path: P) -> Result<hpo::Ontology, anyhow::Error> {
+    if path.as_ref().join("hpo.bin").exists() {
+        tracing::info!(
+            "  attempting to load binary HPO file from {}",
+            path.as_ref().display()
+        );
+        Ok(hpo::Ontology::from_binary(path.as_ref().join("hpo.bin"))?)
+    } else {
+        tracing::info!(
+            "  attempting to load HPO from standard file {}",
+            path.as_ref().display()
+        );
+        Ok(hpo::Ontology::from_standard(&format!(
+            "{}",
+            path.as_ref().display()
+        ))?)
+    }
+}

src/convert/mod.rs

@@ -0,0 +1,58 @@
+//! Conversion of HPO text files to `hpo` binary format.
+
+use std::io::Write;
+
+use clap::Parser;
+
+/// Command line arguments for `convert` sub command.
+#[derive(Parser, Debug)]
+#[command(author, version, about = "Convert HPO text files to binary format", long_about = None)]
+pub struct Args {
+    /// Path to the directory with the HPO files.
+    #[arg(long, required = true)]
+    pub path_hpo_dir: String,
+    /// Path to the output binary file.
+    #[arg(long, required = true)]
+    pub path_out_bin: String,
+}
+
+/// Main entry point for `convert` command.
+///
+/// # Errors
+///
+/// In the case of query execution failure.
+///
+/// # Panics
+///
+/// In the case of term lookup failure.
+pub fn run(args_common: &crate::common::Args, args: &Args) -> Result<(), anyhow::Error> {
+    tracing::info!("args_common = {:?}", &args_common);
+    tracing::info!("args = {:?}", &args);
+
+    if let Some(level) = args_common.verbose.log_level() {
+        match level {
+            log::Level::Trace | log::Level::Debug => {
+                std::env::set_var("RUST_LOG", "debug");
+                env_logger::init_from_env(env_logger::Env::new().default_filter_or("info"));
+            }
+            _ => (),
+        }
+    }
+
+    tracing::info!("Loading HPO...");
+    let before_loading = std::time::Instant::now();
+    let hpo = crate::common::load_hpo(&args.path_hpo_dir)?;
+    tracing::info!("...done loading HPO in {:?}", before_loading.elapsed());
+    tracing::info!("Ontology [{}] with {} terms", hpo.hpo_version(), hpo.len());
+
+    tracing::info!("Writing binary file...");
+    let before_writing = std::time::Instant::now();
+    let filename = &args.path_out_bin;
+    let mut fh = std::fs::File::create(filename).expect("Cannot create file");
+    fh.write_all(&hpo.as_bytes())?;
+    tracing::info!("...done writing binary in {:?}", before_writing.elapsed());
+
+    tracing::info!("All done. Have a nice day!");
+
+    Ok(())
+}