feat: adjust "strucvars query" to read ingested VCFs (#202) (#220)

.gitattributes

@@ -2,3 +2,8 @@ tests/db/compile/** filter=lfs diff=lfs merge=lfs -text
 tests/db/hpo/** filter=lfs diff=lfs merge=lfs -text
 tests/seqvars/ingest/db/grch37/** filter=lfs diff=lfs merge=lfs -text
 tests/seqvars/ingest/db/grch38/** filter=lfs diff=lfs merge=lfs -text
+tests/strucvars/query/*.vcf filter=lfs diff=lfs merge=lfs -text
+tests/strucvars/query/grch37/** filter=lfs diff=lfs merge=lfs -text
+tests/strucvars/query/noref/** filter=lfs diff=lfs merge=lfs -text
+tests/strucvars/query/db/** filter=lfs diff=lfs merge=lfs -text
+src/strucvars/query/snapshots/*smoke_test*.snap filter=lfs diff=lfs merge=lfs -text

README.md

@@ -283,7 +283,49 @@ Convert output of [varfish-db-downloader](https://github.com/bihealth/varfish-db
 $ varfish-server-worker strucvars txt-to-bin \
     --input-type {ClinvarSv,StrucvarInhouse,...} \
     --path-input IN.txt \
-    --path-output-bin DST.bin
+    --path-output DST.bin
+```
+
+## The `strucvars query` Command
+
+Run a query on a VCF file with structural variants as created by `strucvars ingest` using a varfish worker database.
+
+```
+$ varfish-server-worker strucvars query \
+    --genome-release grch37 \
+    --path-db path/to/worker-db \
+    --path-input IN.vcf.gz \
+    --path-output OUT.jsonl
+```
+
+The worker database has the following structure:
+
+```
+$ROOT/
+    noref/
+        genes/
+            acmg.tsv              -- ACMG SF list genes
+            mim2gene.tsv          -- OMIM to NCBI mapping from clingen
+            xlink.bin             -- gene crosslinks
+    {genome_release}/             -- one per genome release
+        mehari/
+            txs.bin.zstd          -- mehari transcripts
+        features/                 -- features important for annotation
+            masked_repeat.bin     -- masked repeats
+            masked_seqdup.bin     -- masked segmental duplications
+        strucvars/                -- structural variant specific
+            bgdbs/                -- background databases
+                dbvar.bin         -- dbVar
+                dgv.bin           -- DGV
+                dgv-gs.bin        -- DGV gold standard
+                exac.bin          -- ExAC CNVs
+                g1k.bin           -- 1000 genomes CNVs
+                gnomad.bin        -- gnomAD-SVs
+            clinvar.bin           -- ClinVar SVs
+            inhouse.bin           -- inhouse SV database
+            patho-mms.bed         -- well-known pathogenic DELs/DUPs
+        tads/
+            hesc.bed              -- hESC TAD definitions
 ```
 
 # Developer Information

build.rs

@@ -1,6 +1,8 @@
 // The custom build script, needed as we use protocolbuffers.
 
 fn main() {
+    println!("cargo:rerun-if-changed=src/proto/varfish/v1/clinvar.proto");
+    println!("cargo:rerun-if-changed=src/proto/varfish/v1/sv.proto");
     prost_build::Config::new()
         .protoc_arg("-Isrc/proto")
         // Add serde serialization and deserialization to the generated code.
@@ -9,7 +11,10 @@ fn main() {
         .type_attribute(".", "#[serde_with::skip_serializing_none]")
         // Define the protobuf files to compile.
         .compile_protos(
-            &["varfish/v1/clinvar.proto", "varfish/v1/sv.proto"],
+            &[
+                "src/proto/varfish/v1/clinvar.proto",
+                "src/proto/varfish/v1/sv.proto",
+            ],
             &["src/"],
         )
         .unwrap();

misc/ok-query-trio.json

@@ -51,14 +51,13 @@
       "CNV"
     ],
     "tx_effects": [
-      "transcript_ablation",
-      "transcript_translocation",
+      "transcript_variant",
       "exon_variant",
       "splice_region_variant",
       "intron_variant",
       "upstream_variant",
       "downstream_variant",
-      "intergenic_variant",
+      "intergenic_variant"
     ],
     "gene_allowlist": null,
     "genomic_region": null,

src/common.rs

@@ -149,10 +149,10 @@ where
 )]
 pub enum GenomeRelease {
     // GRCh37 / hg19
-    #[strum(serialize = "GRCh37")]
+    #[strum(serialize = "grch37")]
     Grch37,
     /// GRCh38 / hg38
-    #[strum(serialize = "GRCh38")]
+    #[strum(serialize = "grch38")]
     Grch38,
 }
 
@@ -207,6 +207,8 @@ pub enum Genotype {
     Het,
     /// hom. alt.
     HomAlt,
+    /// other, includes no-call
+    WithNoCall,
 }
 
 impl std::str::FromStr for Genotype {
@@ -217,6 +219,7 @@ impl std::str::FromStr for Genotype {
             "0/0" | "0|0" | "0" => Genotype::HomRef,
             "0/1" | "1/0" | "0|1" | "1|0" => Genotype::Het,
             "1/1" | "1|1" | "1" => Genotype::HomAlt,
+            "./." | "./0" | "./1" | "0/." | "1/." => Genotype::WithNoCall,
             _ => anyhow::bail!("invalid genotype value: {:?}", s),
         })
     }

src/seqvars/aggregate/mod.rs

@@ -76,6 +76,7 @@ fn handle_record(
             };
 
         let carrier_genotype = match (chrom, individual.sex, genotype) {
+            (_, _, Genotype::WithNoCall) => continue,
             // on the autosomes, male/female count the same
             (Chrom::Auto, _, Genotype::HomRef) => {
                 res_counts.count_an += 2;

src/seqvars/prefilter/mod.rs

@@ -282,7 +282,6 @@ mod test {
         );
 
         let params_file = tmpdir.to_path_buf().join("params.json");
-        eprintln!("params_json = {}", &params_json);
         std::fs::write(&params_file, &params_json)?;
 
         let args = super::Args {

src/strucvars/aggregate/output.rs

@@ -154,6 +154,7 @@ impl Record {
                 };
 
             match (chrom, individual.sex, genotype) {
+                (_, _, Genotype::WithNoCall) => continue,
                 // on the autosomes, male/female count the same
                 (Chrom::Auto, _, Genotype::HomRef) => (),
                 (Chrom::Auto, _, Genotype::Het) => {

src/strucvars/ingest/mod.rs

@@ -254,7 +254,6 @@ impl mehari::annotate::seqvars::AnnotatedVcfWriter for WriterWrapper {
         }
 
         let key_callers: vcf::record::info::field::Key = "callers".parse()?;
-        eprintln!("callers = {:?}", &input_record.info().get(&key_callers));
         if let Some(Some(callers)) = input_record.info().get(&key_callers) {
             if let vcf::record::info::field::Value::Array(
                 vcf::record::info::field::value::Array::String(callers),

src/strucvars/query/bgdbs.rs

@@ -15,7 +15,7 @@ use crate::{
 };
 
 use super::{
-    records::ChromRange,
+    schema::ChromRange,
     schema::{CaseQuery, StructuralVariant, SvType},
 };
 

src/strucvars/query/clinvar.rs

@@ -9,7 +9,7 @@ use tracing::{info, warn};
 use crate::common::{reciprocal_overlap, GenomeRelease, CHROMS};
 
 use super::{
-    records::ChromRange,
+    schema::ChromRange,
     schema::{Pathogenicity, StructuralVariant, SvType},
 };
 

src/strucvars/query/genes.rs

@@ -259,7 +259,7 @@ impl OmimDb {
 }
 
 #[tracing::instrument]
-fn load_omim_db(path: &Path) -> Result<OmimDb, anyhow::Error> {
+fn load_mim2gene_db(path: &Path) -> Result<OmimDb, anyhow::Error> {
     tracing::debug!("loading OMIM TSV records from {:?}...", path);
 
     let before_loading = Instant::now();
@@ -291,7 +291,7 @@ pub struct GeneDb {
     pub ensembl: GeneRegionDb,
     pub xlink: XlinkDb,
     pub acmg: AcmgDb,
-    pub omim: OmimDb,
+    pub mim2gene: OmimDb,
 }
 
 impl GeneDb {
@@ -327,7 +327,11 @@ pub fn load_gene_db(path_db: &str, genome_release: GenomeRelease) -> Result<Gene
         )?,
         xlink: load_xlink_db(Path::new(path_db).join("noref/genes/xlink.bin").as_path())?,
         acmg: load_acmg_db(Path::new(path_db).join("noref/genes/acmg.tsv").as_path())?,
-        omim: load_omim_db(Path::new(path_db).join("noref/genes/omim.tsv").as_path())?,
+        mim2gene: load_mim2gene_db(
+            Path::new(path_db)
+                .join("noref/genes/mim2gene.tsv")
+                .as_path(),
+        )?,
     };
 
     Ok(result)