mehari 0.45.0

Variant effect prediction all in Rust
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104

use crate::common::Args as CommonArgs;
use crate::db::keys::Var;
use crate::db::{PipelineConfig, write_contig_dictionary};
use crate::pbs::seqvars::DbsnpRecord;
use anyhow::Error;
use clap::Parser;
use itertools::Itertools;
use noodles::vcf::variant::record::Ids;
use prost::Message;
use rustc_hash::FxHashMap;
use std::collections::{HashMap, HashSet};
use std::sync::{Arc, RwLock};

/// Arguments for the dbSNP database construction command.
#[derive(Parser, Debug, Clone)]
#[command(about = "Construct dbSNP RocksDB database", long_about = None)]
pub struct Args {
    #[command(flatten)]
    pub common: crate::db::CommonPipelineArgs,
}

pub mod cli {
    pub use super::Args;
}

pub fn run(_common: &CommonArgs, args: &Args) -> Result<(), Error> {
    let config = PipelineConfig {
        assembly: &args.common.assembly,
        input: &args.common.input,
        output: &args.common.output,
        batch_size: args.common.batch_size,
        no_progress: args.common.no_progress,
        threads: args.common.threads,
        db_type: "dbsnp",
        schema_version: "3.0",
        extra_meta: HashMap::new(),
    };

    let header_modifier = |header: &mut noodles::vcf::Header| {
        if let Some(info) = header.infos_mut().get_mut("RS") {
            *info.type_mut() = noodles::vcf::header::record::value::map::info::Type::String;
        }
    };

    let chrom_to_id = Arc::new(RwLock::new(FxHashMap::default()));
    let chrom_to_id_closure = Arc::clone(&chrom_to_id);

    crate::db::run_vcf_pipeline(
        config,
        Some(header_modifier),
        move |record, contig_manager| {
            let mut kvs = Vec::new();
            let chrom = record.reference_sequence_name();
            let pos = match record.variant_start() {
                Some(start) => start.get() as i32,
                None => return Ok((kvs, HashSet::new())),
            };

            let rs_id = if record.ids().is_empty() {
                "".to_string()
            } else {
                record.ids().iter().join(";")
            };
            if rs_id.is_empty() {
                return Ok((kvs, HashSet::new()));
            }

            let (chrom_std, chrom_id) =
                crate::db::get_or_intern_contig(chrom, contig_manager, &chrom_to_id_closure);

            let reference = record.reference_bases();

            for alt in record.alternate_bases().as_ref() {
                let alt_str = alt.to_string();
                let var = Var::new(
                    chrom_std.clone(),
                    pos,
                    reference.to_string(),
                    alt_str.clone(),
                );
                let key = var.encode_with_id(chrom_id);
                let record_pb = DbsnpRecord {
                    allele: alt_str,
                    rs_id: rs_id.clone(),
                };

                let mut value = Vec::new();
                record_pb.encode(&mut value)?;

                let var_label = format!(
                    "{}:{}{}>{}",
                    var.chrom, var.pos, var.reference, var.alternative
                );
                kvs.push((key, value, var_label));
            }
            Ok((kvs, HashSet::new()))
        },
    )?;

    tracing::info!("Writing dbSNP contig index metadata mapping into the meta CF...");
    write_contig_dictionary(&args.common.output, "dbsnp", &chrom_to_id)?;

    Ok(())
}