pub fn read_features(path: &Path) -> Result<Vec<ScoredFeature>, KothError>Expand description
Read scored features from a file, detecting format from extension.
Examples found in repository?
examples/lfq_group_audit.rs (line 23)
11fn main() -> anyhow::Result<()> {
12 env_logger::init();
13 let args: Vec<_> = std::env::args().skip(1).map(PathBuf::from).collect();
14 anyhow::ensure!(args.len() == 3, "usage: lfq_group_audit BATCH CONFIG OUT");
15 let cfg = AlignConfig::from_toml(&args[1])?;
16 let paths = discover_runs(&args[0])?;
17 anyhow::ensure!(!paths.is_empty(), "no feature runs found");
18 let runs: Vec<_> = paths
19 .iter()
20 .map(|r| {
21 Ok(RunInput {
22 name: r.name.clone(),
23 features: read_features(&r.features_path)?,
24 hills: Vec::new(),
25 scan_times: Vec::new(),
26 rt_bounds: None,
27 })
28 })
29 .collect::<anyhow::Result<_>>()?;
30 let alignment = align_runs(&runs, &cfg.alignment);
31 let groups = build_consensus_candidates(&runs, &alignment, &cfg.lfq);
32 std::fs::create_dir_all(&args[2])?;
33 let mut out = std::io::BufWriter::new(std::fs::File::create(args[2].join("candidates.tsv"))?);
34 writeln!(out, "candidate_id\tmassCalib\tmz\tcharge\trtApex\tim\tcombined_score\tn_contributing_runs\tgroup_score\tgroup_qvalue\tweak_members\tretained")?;
35 for (i, g) in groups.iter().enumerate() {
36 let weak = g
37 .members
38 .iter()
39 .filter(|&&(r, f)| runs[r].features[f as usize].combined_score < 0.5)
40 .count();
41 writeln!(
42 out,
43 "{i}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{weak}\t{}",
44 g.neutral_mass,
45 g.ref_mz,
46 g.charge,
47 g.ref_rt,
48 g.ref_im,
49 g.seed_combined_score,
50 g.n_contributing_runs,
51 g.group_score,
52 g.group_qvalue,
53 g.n_contributing_runs >= 2 && g.group_qvalue <= cfg.lfq.consensus.max_group_qvalue
54 )?;
55 }
56 out.flush()?;
57 Ok(())
58}More examples
examples/lfq_member_geometry.rs (line 21)
8fn main() -> anyhow::Result<()> {
9 let args: Vec<_> = std::env::args().skip(1).collect();
10 anyhow::ensure!(
11 args.len() == 5,
12 "usage: lfq_member_geometry BATCH CONFIG MEMBERS_TSV OUT_TSV ID,ID"
13 );
14 let cfg = AlignConfig::from_toml(&PathBuf::from(&args[1]))?;
15 let paths = discover_runs(&PathBuf::from(&args[0]))?;
16 let runs: Vec<_> = paths
17 .iter()
18 .map(|r| {
19 Ok(RunInput {
20 name: r.name.clone(),
21 features: read_features(&r.features_path)?,
22 hills: Vec::new(),
23 scan_times: Vec::new(),
24 rt_bounds: None,
25 })
26 })
27 .collect::<anyhow::Result<_>>()?;
28 check_batch_mobility_scales(
29 paths
30 .iter()
31 .zip(&runs)
32 .map(|(p, r)| (p, r.features.as_slice())),
33 )?;
34 let alignment = align_runs(&runs, &cfg.alignment);
35 let ids: HashSet<usize> = args[4]
36 .split(',')
37 .map(str::parse)
38 .collect::<Result<_, _>>()?;
39 let mut input = csv::ReaderBuilder::new()
40 .delimiter(b'\t')
41 .from_path(&args[2])?;
42 let mut output = csv::WriterBuilder::new()
43 .delimiter(b'\t')
44 .from_path(&args[3])?;
45 output.write_record([
46 "candidate_id",
47 "run_id",
48 "feature_idx",
49 "mass",
50 "rt",
51 "im",
52 "quality",
53 "is_primary",
54 "neutron_offset",
55 ])?;
56 for row in input.records() {
57 let row = row?;
58 let id: usize = row[0].parse()?;
59 if !ids.contains(&id) {
60 continue;
61 }
62 let r: usize = row[1].parse()?;
63 let f: usize = row[2].parse()?;
64 let feat = &runs[r].features[f];
65 let mz = feat.monoisotopic_mz();
66 let rt = feat.feature.rt_apex();
67 let im = feat.feature.im_apex();
68 let (mz, rt, im) = if r == alignment.reference_idx {
69 (mz, rt, im)
70 } else {
71 let a = &alignment.alignments[&runs[r].name];
72 (a.correct_mz(mz, rt), a.warp_rt(rt), a.correct_im(im, rt))
73 };
74 let mass = (mz - 1.007276466621) * feat.feature.charge as f64;
75 output.write_record([
76 id.to_string(),
77 r.to_string(),
78 f.to_string(),
79 mass.to_string(),
80 rt.to_string(),
81 im.to_string(),
82 feat.combined_score.to_string(),
83 row[3].to_string(),
84 feat.neutron_offset.to_string(),
85 ])?;
86 }
87 output.flush()?;
88 Ok(())
89}examples/lfq_permutation_audit.rs (line 38)
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}examples/lfq_refinement.rs (line 58)
37fn main() -> anyhow::Result<()> {
38 env_logger::init();
39 let args = Args::parse();
40 anyhow::ensure!(
41 !(args.quick && args.conditioned),
42 "conditioned controls require the full audit"
43 );
44 let config = AlignConfig::from_toml(&args.config)?;
45 let all_paths = discover_runs(&args.batch)?;
46 anyhow::ensure!(!all_paths.is_empty(), "no runs");
47 let count = args.take.unwrap_or(all_paths.len());
48 anyhow::ensure!(
49 count >= 2 && count <= all_paths.len(),
50 "take must be 2..=number of runs"
51 );
52 let selected: Vec<_> = (0..count).map(|i| i * all_paths.len() / count).collect();
53 let paths: Vec<_> = selected.iter().map(|&i| &all_paths[i]).collect();
54 let mut runs = Vec::new();
55 for p in &paths {
56 runs.push(RunInput {
57 name: p.name.clone(),
58 features: read_features(&p.features_path)?,
59 hills: Vec::new(),
60 scan_times: Vec::new(),
61 rt_bounds: None,
62 });
63 }
64 check_batch_mobility_scales(
65 paths
66 .iter()
67 .zip(&runs)
68 .map(|(p, r)| (*p, r.features.as_slice())),
69 )?;
70 if let Some(name) = &config.alignment.reference_run {
71 anyhow::ensure!(
72 runs.iter().any(|r| &r.name == name),
73 "reference run is not selected"
74 );
75 }
76 std::fs::create_dir_all(&args.output)?;
77 std::fs::write(args.output.join("runs.json"),serde_json::to_vec_pretty(&paths.iter().enumerate().map(|(id,p)|
78 serde_json::json!({"run_id":id,"name":p.name,"features":p.features_path})).collect::<Vec<_>>())?)?;
79 let t = Instant::now();
80 let alignment = align_runs(&runs, &config.alignment);
81 let alignment_sec = t.elapsed().as_secs_f64();
82 let t = Instant::now();
83 let audit = audit_consensus(&runs, &alignment, &config.lfq, !args.quick);
84 let audit_sec = t.elapsed().as_secs_f64();
85 std::fs::write(
86 args.output.join("audit.json"),
87 serde_json::to_vec_pretty(&serde_json::json!({
88 "variants":VARIANT_NAMES,"controls":audit.null_group_counts,"full_controls":!args.quick,"alignment_sec":alignment_sec,"audit_sec":audit_sec,
89 "n_runs":runs.len(),"n_groups":audit.groups.len(),"reference":alignment.reference_name}))?,
90 )?;
91 let mut out = csv::WriterBuilder::new()
92 .delimiter(b'\t')
93 .from_path(args.output.join("candidates.tsv"))?;
94 let mut header: Vec<String> = [
95 "candidate_id",
96 "massCalib",
97 "mz",
98 "charge",
99 "rtApex",
100 "im",
101 "combined_score",
102 "n_contributing_runs",
103 "weak_members",
104 ]
105 .map(String::from)
106 .to_vec();
107 for v in VARIANT_NAMES {
108 for field in [
109 "score",
110 "global",
111 "repeat",
112 "density",
113 "pooled",
114 "conditioned",
115 ] {
116 header.push(format!("{v}_{field}"));
117 }
118 }
119 out.write_record(header)?;
120 let mut members = csv::WriterBuilder::new()
121 .delimiter(b'\t')
122 .from_path(args.output.join("members.tsv"))?;
123 members.write_record(["candidate_id", "run_id", "feature_idx", "is_primary"])?;
124 for (i, g) in audit.groups.iter().enumerate() {
125 let weak = g
126 .members
127 .iter()
128 .filter(|&&(r, f)| runs[r].features[f as usize].combined_score < 0.5)
129 .count();
130 let mut row = vec![
131 i.to_string(),
132 g.neutral_mass.to_string(),
133 g.ref_mz.to_string(),
134 g.charge.to_string(),
135 g.ref_rt.to_string(),
136 g.ref_im.to_string(),
137 g.seed_combined_score.to_string(),
138 g.n_contributing_runs.to_string(),
139 weak.to_string(),
140 ];
141 for v in 0..4 {
142 row.extend(
143 [
144 audit.scores[i][v],
145 audit.q_global[i][v],
146 audit.q_repeat[i][v],
147 audit.q_density[i][v],
148 audit.q_pooled[i][v],
149 audit.q_conditioned[i][v],
150 ]
151 .map(|x| x.to_string()),
152 );
153 }
154 out.write_record(row)?;
155 for &(r, f) in &g.members {
156 members.write_record([
157 i.to_string(),
158 r.to_string(),
159 f.to_string(),
160 (g.per_run_feature[r] == Some(f)).to_string(),
161 ])?;
162 }
163 }
164 out.flush()?;
165 members.flush()?;
166 let selected_q = if args.conditioned {
167 &audit.q_conditioned
168 } else {
169 &audit.q_global
170 };
171 for variant in args.quantify.split(',').filter(|s| !s.is_empty()) {
172 let v: usize = variant.parse()?;
173 anyhow::ensure!(v < 4, "variant must be 0..3");
174 let folder = args.output.join(VARIANT_NAMES[v]);
175 std::fs::create_dir_all(&folder)?;
176 let groups = audit
177 .groups
178 .iter()
179 .enumerate()
180 .filter(|(i, g)| {
181 g.n_contributing_runs >= 2
182 && selected_q[*i][v] <= config.lfq.consensus.max_group_qvalue
183 })
184 .map(|(i, g)| {
185 let mut g = g.clone();
186 g.group_score = audit.scores[i][v];
187 g.group_qvalue = selected_q[i][v];
188 g
189 })
190 .collect();
191 let t = Instant::now();
192 let matrix = quantify_consensus(&runs, &alignment, &config.lfq, groups, |i| {
193 read_hills(&paths[i].hills_path).expect("read hills")
194 });
195 let lfq_sec = t.elapsed().as_secs_f64();
196 matrices(&matrix, &folder)?;
197 std::fs::write(folder.join("align_config.toml"), config.to_toml_string()?)?;
198 let report = build_align_report(
199 &runs,
200 &alignment,
201 &matrix,
202 &config,
203 AlignTiming {
204 alignment_sec,
205 lfq_sec,
206 total_sec: alignment_sec + audit_sec + lfq_sec,
207 },
208 );
209 write_align_report(&report, &folder.join("align_report.json"))?;
210 let mut meta = std::fs::File::create(folder.join("variant.json"))?;
211 writeln!(
212 meta,
213 "{}",
214 serde_json::json!({"variant":VARIANT_NAMES[v],"group_gate":config.lfq.consensus.max_group_qvalue,"control": if args.conditioned { "ten quality-conditioned rank shifts" } else { "first five global shifts" }})
215 )?;
216 }
217 Ok(())
218}