lfq_permutation_audit/
lfq_permutation_audit.rs

1//! Independent-permutation controls; intentionally outside the production CLI.
2#[path = "support/lfq_matrices.rs"]
3mod lfq_matrices;
4use clap::Parser;
5use koth_ms::{
6    alignment::{align_runs, RunInput},
7    config::AlignConfig,
8    input::{check_batch_mobility_scales, discover_runs, read_features, read_hills},
9    lfq::{consensus::audit_permutation_consensus, quantify_consensus},
10};
11use std::{path::PathBuf, time::Instant};
12
13#[derive(Parser)]
14struct Args {
15    batch: PathBuf,
16    config: PathBuf,
17    output: PathBuf,
18    #[arg(long)]
19    take: Option<usize>,
20    /// Extract the pooled-ten accepted set with the unchanged production pipeline.
21    #[arg(long)]
22    quantify: bool,
23}
24
25fn main() -> anyhow::Result<()> {
26    env_logger::init();
27    let args = Args::parse();
28    let cfg = AlignConfig::from_toml(&args.config)?;
29    let all = discover_runs(&args.batch)?;
30    let n = args.take.unwrap_or(all.len());
31    anyhow::ensure!(n >= 2 && n <= all.len(), "take must be 2..=number of runs");
32    let paths: Vec<_> = (0..n).map(|i| &all[i * all.len() / n]).collect();
33    let runs: Vec<_> = paths
34        .iter()
35        .map(|r| {
36            Ok(RunInput {
37                name: r.name.clone(),
38                features: read_features(&r.features_path)?,
39                hills: Vec::new(),
40                scan_times: Vec::new(),
41                rt_bounds: None,
42            })
43        })
44        .collect::<anyhow::Result<_>>()?;
45    check_batch_mobility_scales(
46        paths
47            .iter()
48            .zip(&runs)
49            .map(|(p, r)| (*p, r.features.as_slice())),
50    )?;
51    if let Some(name) = &cfg.alignment.reference_run {
52        anyhow::ensure!(
53            runs.iter().any(|r| &r.name == name),
54            "reference run is not selected"
55        );
56    }
57    std::fs::create_dir_all(&args.output)?;
58    std::fs::write(args.output.join("runs.json"),serde_json::to_vec_pretty(&paths.iter().enumerate().map(|(id,p)|
59        serde_json::json!({"run_id":id,"name":p.name,"features":p.features_path})).collect::<Vec<_>>())?)?;
60    let t = Instant::now();
61    let alignment = align_runs(&runs, &cfg.alignment);
62    let audit = audit_permutation_consensus(&runs, &alignment, &cfg.lfq);
63    let mut out = csv::WriterBuilder::new()
64        .delimiter(b'\t')
65        .from_path(args.output.join("candidates.tsv"))?;
66    out.write_record([
67        "candidate_id",
68        "massCalib",
69        "mz",
70        "charge",
71        "rtApex",
72        "im",
73        "combined_score",
74        "n_contributing_runs",
75        "weak_members",
76        "average_score",
77        "average_first",
78        "average_repeat",
79        "average_permuted",
80    ])?;
81    let mut counts = [0usize; 3];
82    let (mut intersection, mut union) = (0usize, 0usize);
83    for (i, g) in audit.groups.iter().enumerate() {
84        let keep = [audit.q_first[i], audit.q_repeat[i], audit.q_pooled[i]]
85            .map(|q| g.n_contributing_runs >= 2 && q <= cfg.lfq.consensus.max_group_qvalue);
86        for j in 0..3 {
87            counts[j] += usize::from(keep[j]);
88        }
89        intersection += usize::from(keep[0] && keep[1]);
90        union += usize::from(keep[0] || keep[1]);
91        let weak = g
92            .members
93            .iter()
94            .filter(|&&(r, f)| runs[r].features[f as usize].combined_score < 0.5)
95            .count();
96        out.write_record([
97            i.to_string(),
98            g.neutral_mass.to_string(),
99            g.ref_mz.to_string(),
100            g.charge.to_string(),
101            g.ref_rt.to_string(),
102            g.ref_im.to_string(),
103            g.seed_combined_score.to_string(),
104            g.n_contributing_runs.to_string(),
105            weak.to_string(),
106            g.group_score.to_string(),
107            audit.q_first[i].to_string(),
108            audit.q_repeat[i].to_string(),
109            audit.q_pooled[i].to_string(),
110        ])?;
111    }
112    out.flush()?;
113    let summary = serde_json::json!({"method":"independent RT permutation within run/charge/quality stratum","n_runs":n,"reference":alignment.reference_name,"elapsed_sec":t.elapsed().as_secs_f64(),"controls":audit.null_group_counts,"retained_first":counts[0],"retained_repeat":counts[1],"retained_pooled":counts[2],"retained_jaccard":if union>0 {Some(intersection as f64 / union as f64)}else{None},"gate":cfg.lfq.consensus.max_group_qvalue});
114    std::fs::write(
115        args.output.join("audit.json"),
116        serde_json::to_vec_pretty(&summary)?,
117    )?;
118    println!("{summary}");
119    if args.quantify {
120        let groups = audit
121            .groups
122            .into_iter()
123            .filter(|g| {
124                g.n_contributing_runs >= 2 && g.group_qvalue <= cfg.lfq.consensus.max_group_qvalue
125            })
126            .collect();
127        let matrix = quantify_consensus(&runs, &alignment, &cfg.lfq, groups, |i| {
128            read_hills(&paths[i].hills_path).expect("read hills")
129        });
130        let folder = args.output.join("average");
131        std::fs::create_dir_all(&folder)?;
132        lfq_matrices::matrices(&matrix, &folder)?;
133        std::fs::write(folder.join("align_config.toml"), cfg.to_toml_string()?)?;
134        std::fs::write(
135            folder.join("variant.json"),
136            serde_json::to_vec_pretty(&summary)?,
137        )?;
138    }
139    Ok(())
140}
lfq_permutation_audit/lfq_permutation_audit.rs

lfq_permutation_audit/
lfq_permutation_audit.rs
