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
105
106
107
108
109
110
111
112
113
114
115
116
117
118

use crate::common::Args as CommonArgs;
use crate::db::keys::Var;
use crate::db::{ContigIdMap, PipelineConfig, write_contig_dictionary};
use crate::pbs::seqvars::{SpliceAiPrediction, SpliceAiRecord};
use anyhow::Error;
use clap::Parser;
use prost::Message;
use std::collections::{HashMap, HashSet};
use std::sync::Arc;

/// Arguments for the SpliceAI database construction command.
#[derive(Parser, Debug, Clone)]
#[command(about = "Construct SpliceAI score 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: "spliceai",
        schema_version: "1.0",
        extra_meta: HashMap::new(),
    };

    let chrom_to_id = ContigIdMap::default();
    let chrom_to_id_closure = Arc::clone(&chrom_to_id);

    crate::db::run_vcf_pipeline(
        config,
        None::<fn(&mut noodles::vcf::Header)>,
        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 spliceai_val = match record.info().get("SpliceAI").flatten() {
                Some(val) => val,
                None => return Ok((kvs, HashSet::new())),
            };

            let spliceai_str = match crate::db::get_info_string(spliceai_val) {
                Some(s) => s,
                None => 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();
            let mut predictions_by_allele: HashMap<String, Vec<SpliceAiPrediction>> =
                HashMap::new();

            for pred_str in spliceai_str.split(',') {
                let fields: Vec<&str> = pred_str.split('|').collect();
                if fields.len() < 10 {
                    anyhow::bail!(
                        "Malformed SpliceAI prediction block: expected at least 10 pipe-separated fields, got {} in token '{}'",
                        fields.len(),
                        pred_str
                    );
                }

                let allele = fields[0].to_string();
                let prediction = SpliceAiPrediction {
                    allele: allele.clone(),
                    symbol: fields[1].to_string(),
                    ds_ag: fields[2].parse()?,
                    ds_al: fields[3].parse()?,
                    ds_dg: fields[4].parse()?,
                    ds_dl: fields[5].parse()?,
                    dp_ag: fields[6].parse()?,
                    dp_al: fields[7].parse()?,
                    dp_dg: fields[8].parse()?,
                    dp_dl: fields[9].parse()?,
                };
                predictions_by_allele
                    .entry(allele)
                    .or_default()
                    .push(prediction);
            }

            for (allele, predictions) in predictions_by_allele {
                let var = Var::new(chrom_std.clone(), pos, reference.to_string(), allele);
                let key = var.encode_with_id(chrom_id);
                let record_pb = SpliceAiRecord { predictions };

                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 SpliceAI contig index metadata mapping into the meta CF...");
    write_contig_dictionary(&args.common.output, "spliceai", &chrom_to_id)?;

    Ok(())
}