uno_core/quant/directms1/
input.rs

1use super::{Config, Dataset, InputPeptide};
2use anyhow::{bail, ensure, Context, Result};
3use serde::Serialize;
4use sha2::{Digest, Sha256};
5use std::{
6    collections::{BTreeMap, BTreeSet},
7    path::{Path, PathBuf},
8};
9
10#[derive(Clone, Debug, Serialize)]
11pub struct InputHash {
12    pub path: PathBuf,
13    pub sha256: String,
14}
15#[derive(Default)]
16pub struct Reader {
17    pub hashes: Vec<InputHash>,
18}
19/// One TSV record addressed by column name. Records share their header index,
20/// so a row costs one `StringRecord` rather than a map of owned strings.
21pub struct Row {
22    columns: std::rc::Rc<BTreeMap<String, usize>>,
23    record: csv::StringRecord,
24}
25impl Row {
26    /// The value of an optional column, `None` when the file lacks it.
27    pub fn get(&self, name: &str) -> Option<&str> {
28        self.columns.get(name).map(|&i| &self.record[i])
29    }
30}
31impl std::ops::Index<&str> for Row {
32    type Output = str;
33    fn index(&self, name: &str) -> &str {
34        &self.record[self.columns[name]]
35    }
36}
37impl Reader {
38    /// Hash `path`, check its header, and hand each record to `each` in file order.
39    pub fn each(
40        &mut self,
41        path: &Path,
42        required: &[&str],
43        each: impl FnMut(&Row) -> Result<()>,
44    ) -> Result<()> {
45        let bytes = std::fs::read(path).with_context(|| format!("reading {}", path.display()))?;
46        self.hashes.push(InputHash {
47            path: path.to_path_buf(),
48            sha256: format!("{:x}", Sha256::digest(&bytes)),
49        });
50        scan(path, &bytes, required, each)
51    }
52    /// A second pass over a file this reader already hashed. The file must be
53    /// unchanged, and it is not recorded again, so the manifest's input list
54    /// is the same as for a single pass.
55    pub fn reread(
56        &self,
57        path: &Path,
58        required: &[&str],
59        each: impl FnMut(&Row) -> Result<()>,
60    ) -> Result<()> {
61        let bytes = std::fs::read(path).with_context(|| format!("reading {}", path.display()))?;
62        let sha256 = format!("{:x}", Sha256::digest(&bytes));
63        ensure!(
64            self.hashes
65                .iter()
66                .any(|h| h.path == path && h.sha256 == sha256),
67            "{} changed between passes",
68            path.display()
69        );
70        scan(path, &bytes, required, each)
71    }
72    pub fn table(&mut self, path: &Path, required: &[&str]) -> Result<Vec<Row>> {
73        self.table_where(path, required, |_| true)
74    }
75    /// Like [`Reader::table`], keeping only the rows `keep` accepts.
76    pub fn table_where(
77        &mut self,
78        path: &Path,
79        required: &[&str],
80        keep: impl Fn(&Row) -> bool,
81    ) -> Result<Vec<Row>> {
82        let mut rows = Vec::new();
83        self.each(path, required, |r| {
84            if keep(r) {
85                rows.push(Row {
86                    columns: r.columns.clone(),
87                    record: r.record.clone(),
88                });
89            }
90            Ok(())
91        })?;
92        Ok(rows)
93    }
94    pub fn lines(&mut self, path: &Path) -> Result<BTreeSet<String>> {
95        let bytes = std::fs::read(path)?;
96        self.hashes.push(InputHash {
97            path: path.to_path_buf(),
98            sha256: format!("{:x}", Sha256::digest(&bytes)),
99        });
100        Ok(std::str::from_utf8(&bytes)?
101            .lines()
102            .map(str::trim)
103            .filter(|s| !s.is_empty())
104            .map(str::to_owned)
105            .collect())
106    }
107}
108/// Check the header of the TSV `bytes` and hand each record to `each` in file order.
109fn scan(
110    path: &Path,
111    bytes: &[u8],
112    required: &[&str],
113    mut each: impl FnMut(&Row) -> Result<()>,
114) -> Result<()> {
115    let mut reader = csv::ReaderBuilder::new()
116        .delimiter(b'\t')
117        .from_reader(bytes);
118    let headers = reader.headers()?.clone();
119    ensure!(
120        headers.iter().collect::<BTreeSet<_>>().len() == headers.len(),
121        "duplicate columns in {}",
122        path.display()
123    );
124    for name in required {
125        ensure!(
126            headers.iter().any(|h| h == *name),
127            "missing column {name} in {}",
128            path.display()
129        );
130    }
131    let columns = std::rc::Rc::new(
132        headers
133            .iter()
134            .enumerate()
135            .map(|(i, h)| (h.to_owned(), i))
136            .collect::<BTreeMap<_, _>>(),
137    );
138    let mut row = Row {
139        columns,
140        record: csv::StringRecord::new(),
141    };
142    while reader.read_record(&mut row.record)? {
143        each(&row)?;
144    }
145    Ok(())
146}
147fn companion(path: &Path, suffix: &str) -> Result<PathBuf> {
148    let name = path.to_str().context("non-UTF8 reference path")?;
149    let base = name
150        .strip_suffix("_proteins_full.tsv")
151        .context("reference inputs must end in _proteins_full.tsv")?;
152    Ok(PathBuf::from(format!("{base}{suffix}")))
153}
154fn float(row: &Row, key: &str) -> Result<f64> {
155    row[key]
156        .parse()
157        .with_context(|| format!("invalid {key}: {}", &row[key]))
158}
159fn py_float(value: f64) -> String {
160    if value.is_nan() {
161        "nan".to_string()
162    } else if value.fract() == 0. {
163        format!("{value:.1}")
164    } else {
165        value.to_string()
166    }
167}
168#[derive(Default)]
169pub struct ReferenceOptions {
170    pub all_proteins: bool,
171    pub all_pfms: bool,
172    pub allowed_proteins: Option<PathBuf>,
173    pub allowed_peptides: Option<PathBuf>,
174}
175pub fn reference(
176    s1: &[PathBuf],
177    s2: &[PathBuf],
178    cfg: &Config,
179    options: &ReferenceOptions,
180) -> Result<(Dataset, Vec<InputHash>)> {
181    ensure!(
182        s1.len() >= 2 && s2.len() >= 2,
183        "at least two input files per condition required"
184    );
185    let paths: Vec<_> = s1.iter().chain(s2).collect();
186    let canonical: Result<BTreeSet<_>> = paths
187        .iter()
188        .map(|p| Ok(std::fs::canonicalize(p)?))
189        .collect();
190    ensure!(
191        canonical?.len() == paths.len(),
192        "a reference file occurs in multiple samples or conditions"
193    );
194    let mut reader = Reader::default();
195    let mut allowed_proteins = BTreeSet::new();
196    let mut allowed_peptides = BTreeSet::new();
197    for path in &paths {
198        if options.allowed_proteins.is_none() {
199            let p = if options.all_proteins {
200                (*path).clone()
201            } else {
202                companion(path, "_proteins.tsv")?
203            };
204            for row in reader.table(&p, &["dbname"])? {
205                let name = row["dbname"].to_owned();
206                if !options.all_proteins {
207                    allowed_proteins.insert(format!("{}{name}", cfg.decoy_prefix));
208                }
209                allowed_proteins.insert(name);
210            }
211        }
212        if options.allowed_peptides.is_none() {
213            for row in reader.table(
214                &companion(path, "_PFMs_ML.tsv")?,
215                &["seqs", "qpreds", "plen"],
216            )? {
217                if (options.all_pfms || float(&row, "qpreds")? <= 10.)
218                    && float(&row, "plen")? >= cfg.min_length as f64
219                {
220                    allowed_peptides.insert(row["seqs"].to_owned());
221                }
222            }
223        }
224    }
225    if let Some(p) = &options.allowed_proteins {
226        allowed_proteins = reader.lines(p)?;
227        for p in allowed_proteins.clone() {
228            allowed_proteins.insert(format!("{}{p}", cfg.decoy_prefix));
229        }
230    }
231    if let Some(p) = &options.allowed_peptides {
232        allowed_peptides = reader.lines(p)?;
233    }
234    log::info!(
235        "native input: {} allowed proteins, {} allowed peptides",
236        allowed_proteins.len(),
237        allowed_peptides.len()
238    );
239    const PFM_COLUMNS: [&str; 6] = [
240        "sequence",
241        "proteins",
242        "charge",
243        "ion_mobility",
244        "Intensity",
245        "RT",
246    ];
247    // Two passes over the PFM tables, so only one file is held at a time:
248    // the first collects the proteins of every eligible match, the second
249    // (hash-checked, not re-recorded) keeps each sample's best precursors.
250    let mut expanded = BTreeSet::new();
251    for path in &paths {
252        reader.each(&companion(path, "_PFMs.tsv")?, &PFM_COLUMNS, |row| {
253            if allowed_peptides.contains(&row["sequence"])
254                && row["proteins"]
255                    .split(';')
256                    .any(|p| allowed_proteins.contains(p))
257            {
258                expanded.extend(row["proteins"].split(';').map(str::to_owned));
259            }
260            Ok(())
261        })?;
262    }
263    log::info!(
264        "native input: {} proteins in eligible matches",
265        expanded.len()
266    );
267    let mut merged: BTreeMap<String, InputPeptide> = BTreeMap::new();
268    for (sample, path) in paths.iter().enumerate() {
269        let mut best: BTreeMap<String, (f64, InputPeptide)> = BTreeMap::new();
270        reader.reread(&companion(path, "_PFMs.tsv")?, &PFM_COLUMNS, |row| {
271            let seq = &row["sequence"];
272            if !allowed_peptides.contains(seq) {
273                return Ok(());
274            }
275            let proteins: Vec<_> = row["proteins"]
276                .split(';')
277                .filter(|p| expanded.contains(*p))
278                .map(str::to_owned)
279                .collect();
280            if proteins.is_empty() {
281                return Ok(());
282            }
283            let intensity = float(row, "Intensity")?;
284            ensure!(intensity.is_finite() && intensity>0.,"nonpositive/nonfinite reference intensity for {seq}; supply missing rows rather than zeros");
285            let charge: u8 = row["charge"].parse().context("invalid precursor charge")?;
286            ensure!(charge > 0, "charge must be positive");
287            // pandas.read_table represents absent native mobility as NaN; its
288            // string "nan" is part of the reference precursor key.
289            let im = if row["ion_mobility"].trim().is_empty() {
290                f64::NAN
291            } else {
292                float(row, "ion_mobility")?
293            };
294            ensure!(!im.is_infinite(), "invalid mobility");
295            let key = format!("{seq}{charge}{}", py_float(im));
296            let entry = InputPeptide {
297                key: key.clone(),
298                sequence: seq.to_owned(),
299                proteins,
300                values: vec![None; paths.len()],
301            };
302            if let Some((v, old)) = best.get(&key) {
303                ensure!(
304                    old.sequence == *seq,
305                    "ambiguous reference composite precursor identity"
306                );
307                if intensity <= *v {
308                    return Ok(());
309                }
310            }
311            best.insert(key, (intensity, entry));
312            Ok(())
313        })?;
314        for (key, (intensity, row)) in best {
315            let p = merged.entry(key).or_insert(row);
316            p.values[sample] = Some(intensity);
317        }
318    }
319    let samples = paths
320        .iter()
321        .map(|p| p.file_name().unwrap().to_string_lossy().to_string())
322        .collect();
323    if merged.is_empty() {
324        bail!("no eligible peptide-feature matches");
325    }
326    Ok((
327        Dataset {
328            samples,
329            n_s1: s1.len(),
330            rows: merged.into_values().collect(),
331        },
332        reader.hashes,
333    ))
334}
uno_core/quant/directms1/input.rs

uno_core/quant/directms1/
input.rs
