diff --git a/.github/pull_request_template.md b/.github/pull_request_template.md new file mode 100644 index 00000000..8a886273 --- /dev/null +++ b/.github/pull_request_template.md @@ -0,0 +1,7 @@ +## Summary + + + +## Checklist + +- [ ] Checked whether this PR needs a version bump and updated it if needed. diff --git a/Cargo.lock b/Cargo.lock index 4fcc85ea..27e6a319 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -204,7 +204,7 @@ dependencies = [ [[package]] name = "alloc_track" -version = "0.34.1" +version = "0.35.0" [[package]] name = "allocator-api2" @@ -320,7 +320,7 @@ checksum = "330a5ed07fa54e4702c9d6c4174f74427fc0ef6e214bbd677ae50a5099946470" [[package]] name = "apex_sim" -version = "0.34.1" +version = "0.35.0" dependencies = [ "eframe", "egui", @@ -371,7 +371,7 @@ dependencies = [ [[package]] name = "array2d" -version = "0.34.1" +version = "0.35.0" dependencies = [ "serde", ] @@ -1746,7 +1746,7 @@ dependencies = [ [[package]] name = "calib_dash" -version = "0.34.1" +version = "0.35.0" dependencies = [ "calibrt", "insta", @@ -1758,7 +1758,7 @@ dependencies = [ [[package]] name = "calibrt" -version = "0.34.1" +version = "0.35.0" dependencies = [ "array2d", "rand 0.9.3", @@ -4693,7 +4693,7 @@ dependencies = [ [[package]] name = "micromzpaf" -version = "0.34.1" +version = "0.35.0" dependencies = [ "serde", "serde_json", @@ -4791,7 +4791,7 @@ dependencies = [ [[package]] name = "msspeculator-core" version = "0.1.0" -source = "git+https://github.com/TalusBio/msspeculator?rev=14634ad065ddcb9750d3c1bb0c817a2c3e332e07#14634ad065ddcb9750d3c1bb0c817a2c3e332e07" +source = "git+https://github.com/TalusBio/msspeculator?rev=3fb9dddbd848f4e7ccf3ccfc049690d2f10431a9#3fb9dddbd848f4e7ccf3ccfc049690d2f10431a9" dependencies = [ "anyhow", "blake2", @@ -4807,7 +4807,7 @@ dependencies = [ [[package]] name = "msspeculator-inference" version = "0.1.0" -source = "git+https://github.com/TalusBio/msspeculator?rev=14634ad065ddcb9750d3c1bb0c817a2c3e332e07#14634ad065ddcb9750d3c1bb0c817a2c3e332e07" +source = "git+https://github.com/TalusBio/msspeculator?rev=3fb9dddbd848f4e7ccf3ccfc049690d2f10431a9#3fb9dddbd848f4e7ccf3ccfc049690d2f10431a9" dependencies = [ "anyhow", "blake2", @@ -6666,7 +6666,7 @@ dependencies = [ [[package]] name = "rescore_dash" -version = "0.34.1" +version = "0.35.0" dependencies = [ "ratatui", "rayon", @@ -7845,7 +7845,7 @@ dependencies = [ [[package]] name = "tims_stage" -version = "0.34.1" +version = "0.35.0" dependencies = [ "bytes", "filetime", @@ -7861,7 +7861,7 @@ dependencies = [ [[package]] name = "timscentroid" -version = "0.34.1" +version = "0.35.0" dependencies = [ "arrow 59.2.0", "async-trait", @@ -7896,7 +7896,7 @@ dependencies = [ [[package]] name = "timsquery" -version = "0.34.1" +version = "0.35.0" dependencies = [ "array2d", "arrow 59.2.0", @@ -7925,7 +7925,7 @@ dependencies = [ [[package]] name = "timsquery_cli" -version = "0.34.1" +version = "0.35.0" dependencies = [ "base64 0.22.1", "clap", @@ -7945,7 +7945,7 @@ dependencies = [ [[package]] name = "timsquery_pyo3" -version = "0.34.1" +version = "0.35.0" dependencies = [ "numpy", "pyo3", @@ -7956,7 +7956,7 @@ dependencies = [ [[package]] name = "timsquery_viewer" -version = "0.34.1" +version = "0.35.0" dependencies = [ "calibrt", "clap", @@ -8091,7 +8091,7 @@ dependencies = [ [[package]] name = "timsseek" -version = "0.34.1" +version = "0.35.0" dependencies = [ "array2d", "arrow 59.2.0", @@ -8119,7 +8119,7 @@ dependencies = [ [[package]] name = "timsseek_cli" -version = "0.34.1" +version = "0.35.0" dependencies = [ "alloc_track", "anyhow", @@ -8149,7 +8149,7 @@ dependencies = [ [[package]] name = "timsseek_macros" -version = "0.34.1" +version = "0.35.0" dependencies = [ "proc-macro2", "quote", diff --git a/Cargo.toml b/Cargo.toml index 769b583b..c28c314f 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -37,7 +37,7 @@ default-members = [ ] [workspace.package] -version = "0.34.1" +version = "0.35.0" edition = "2024" authors = ["Sebastian Paez"] license = "Apache-2.0" diff --git a/README.md b/README.md index efee2d65..8c6df8d6 100644 --- a/README.md +++ b/README.md @@ -69,6 +69,19 @@ cargo run --release --bin timsseek -- build-library \ --max-fragments 10 \ -o $SPECLIB_NAME +# Or predict exact modified peptides and charge states from a TSV. +# protein_ids separates multiple IDs with semicolons; decoy_group is an integer. +cat > peptides.tsv <<'EOF' +proforma protein_ids decoy decoy_group +PEC[UNIMOD:4]TIDEK/2 P1 false 1 +PECTIDEK/3 P1 false 2 +PEK[UNIMOD:259]TIDEK/2 PRTC false 3 +PECTDIEK/3 P1 true 2 +EOF +timsseek build-library --peptides peptides.tsv \ + --decoy-method shuffle --decoy-seed 42 \ + -o peptides.mzspeclib.txt.gz + # Run timsseek. Config is optional; defaults work for most runs. # To tweak tolerances: `timsseek --write-default-config config.toml`, edit, pass with `-c`. # TOML and JSON both accepted (sniffed by extension). @@ -138,11 +151,10 @@ A proteome build takes minutes and reports as it goes: a progress bar per phase Not supported by `build-library` right now, having been dropped along with the separate `speclib_build_cli`: -- **Remote paths.** `--fasta` and `--out` are filesystem paths; a remote URI is rejected by name rather than staged. Build locally and copy. +- **Remote paths.** `--fasta`, `--peptides` and `--out` are filesystem paths; a remote URI is rejected by name rather than staged. Build locally and copy. - **Acquisition and chromatography context.** A build uses the model artifact's own defaults; picking a different one is a decision about the model. -- **A peptide list instead of a FASTA.** Digestion is the only input path. - **Fragment and precursor filters** beyond `--min-intensity` and `--max-fragments`: no minimum transition count, and no precursor or fragment *m/z* bounds. -- **A choice of decoy method.** Decoys are pseudo-reversed; the old `reverse` / `edge_mutate` selection is gone. They are predicted by default, because a search derives ±CH2 mass-shift decoys for a library that ships none and a mass shift is a far weaker null than a reversed sequence through the same model. `--no-decoys` opts out. +- **Decoy generation.** FASTA decoys are pseudo-reversed. Peptide TSV builds support `--decoy-method pseudo-reverse|shuffle` and `--decoy-seed`. Generated decoys are predicted by default; `--no-decoys` disables generation but retains supplied decoys. Pseudo-reverse retries with successively shorter interior spans; shuffle retries with its seed. When no collision-free decoy exists, the builder warns and keeps the target unpaired. Supplied target and decoy rows must share an integer `decoy_group` and charge; each group has one target and at most one supplied decoy. Targets without a group get unique numeric IDs. Peptide TSV builds require an mzSpecLib output suffix so the groups survive serialization. Auth via AWS default chain. MinIO/R2: set `AWS_ENDPOINT_URL`. See `docs/development.md` for `[staging]` config + env var list. @@ -181,4 +193,3 @@ Some of the forms of contributing to the current state of the project could be: - Code - We welcome pull requests! We would really appreciate if an issue is open to discuss potential changes before they are merged. - diff --git a/pyproject.toml b/pyproject.toml index 36faf61f..e1a76c9a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "timsseek-workspace" -version = "0.34.1" +version = "0.35.0" requires-python = ">=3.11,<3.13" dependencies = [] @@ -31,7 +31,7 @@ preview = true select = ["E", "F", "T20", "I"] [tool.bumpver] -current_version = "0.34.1" +current_version = "0.35.0" version_pattern = "MAJOR.MINOR.PATCH[-PYTAGNUM]" tag_message = "v{new_version}" commit_message = "chore: bump version to {new_version}" diff --git a/python/timsquery_pyo3/pyproject.toml b/python/timsquery_pyo3/pyproject.toml index 56e98c94..2f615406 100644 --- a/python/timsquery_pyo3/pyproject.toml +++ b/python/timsquery_pyo3/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "maturin" [project] name = "timsquery_pyo3" -version = "0.34.1" +version = "0.35.0" requires-python = ">=3.11,<3.14" classifiers = [ "Programming Language :: Rust", diff --git a/rust/timsquery/src/models/target_columns.rs b/rust/timsquery/src/models/target_columns.rs index bb868d53..76a26dd7 100644 --- a/rust/timsquery/src/models/target_columns.rs +++ b/rust/timsquery/src/models/target_columns.rs @@ -136,12 +136,13 @@ pub mod test_handles { } } -/// Interned competition groups: `codes[row]` points at `labels`, so rows that -/// compete share one label rather than each storing a copy. +/// Interned declared groups. Codes below `singleton_limit` identify ungrouped +/// rows directly; larger codes index `labels` after subtracting that limit. #[derive(Debug, Clone)] pub struct DecoyGroups { codes: Vec, labels: SourceIds, + singleton_limit: usize, } /// One row's worth of input to [`TargetColumnsBuilder::push_row`]. @@ -187,12 +188,9 @@ pub struct Row<'a, L: KeyLike> { /// is true of every row, and when the groups it did get are all singletons /// -- both cases are what deriving the group from the row already does. /// - /// One reader supplies this, and that is a fact about the formats rather - /// than about this field: mzSpecLib's `related spectrum keys` is the only - /// place a library states that one entry is the counterpart of another. A - /// DIA-NN/Spectronaut/Skyline `transition_group_id` names the row it is on, - /// not a partner, so those readers have nothing to pair by -- a decoy there - /// is a row flagged as a decoy and nothing more. + /// mzSpecLib's `related spectrum keys` and msspeculator's `decoy_group` + /// declare competition groups. A DIA-NN/Spectronaut/Skyline + /// `transition_group_id` names the row it is on, not a partner. pub decoy_group: Option, } @@ -247,14 +245,11 @@ pub struct TargetColumns { /// Built at [`Self::seal`], minted there for formats that carry no ids, so /// a sealed arena always has one per row. pub(crate) source_ids: SourceIds, - /// The competition groups the input declared, interned: a code per row, - /// plus the labels those codes point at. Rows that compete share a code, so - /// the label is stored once per group rather than once per row. + /// The competition groups the input declared: a code per row, plus one + /// label per declared group. Ungrouped rows use their row code and output id. /// - /// `None` when the input declared none, which is every format today. A - /// group is then a singleton -- the row competes with its own decoy - /// variants and nothing else -- so both the code and the label are - /// derivable from the row and nothing is stored. + /// `None` when every group is a singleton. A row then competes only with + /// its own decoy variants, so both code and label are derived from the row. pub(crate) decoy_groups: Option, // CSR prefix offsets (n+1) pub(crate) frag_off: Vec, @@ -432,25 +427,22 @@ impl TargetColumns { let mut claimed: std::collections::HashMap<(GroupCode, u8), usize> = std::collections::HashMap::new(); for (row, group) in self.pending_groups.drain(..).enumerate() { - let group = match group { - Some(g) => g, - // No id column yet -- `seal` builds it after this -- so mirror - // what `output_id` will report for an unnamed row. - None => match self.pending_ids.get(row) { - Some(Some(id)) => id.clone(), - _ => OwnedSourceId::Numeric(row as u64), - }, + let Some(group) = group else { + codes.push(GroupCode::new(row as u32)); + continue; }; let code = *seen.entry(group.clone()).or_insert_with(|| { labels.push(group); - GroupCode::new((labels.len() - 1) as u32) + GroupCode::new( + u32::try_from(n + labels.len() - 1).expect("decoy group codes exceed u32"), + ) }); if !self.is_decoy[row] { let charge = self.charge[row]; match claimed.entry((code, charge)) { std::collections::hash_map::Entry::Occupied(first) => { return Err(SourceIdError::GroupHasTwoTargets { - group: labels[code.get()].to_string(), + group: labels[code.get() - n].to_string(), first_row: *first.get(), row, charge, @@ -465,7 +457,7 @@ impl TargetColumns { } self.pending_groups.shrink_to_fit(); - if labels.len() == n { + if seen.len() + codes.iter().filter(|code| code.get() < n).count() == n { // Every row alone in its own group, which is what `decoy_group_code` // derives for a groupless arena. return Ok(()); @@ -475,6 +467,7 @@ impl TargetColumns { self.decoy_groups = Some(DecoyGroups { codes, labels: SourceIds::owned(labels, n_labels)?, + singleton_limit: n, }); Ok(()) } @@ -496,11 +489,12 @@ impl TargetColumns { /// every result, which is the whole point of the code above. pub fn decoy_group(&self, tgt: RowIdx) -> SourceId<'_> { match &self.decoy_groups { - Some(g) => g + Some(g) if g.codes[tgt.get()].get() >= g.singleton_limit => g .labels - .get(g.codes[tgt.get()].get()) - .expect("every code indexes a label"), + .get(g.codes[tgt.get()].get() - g.singleton_limit) + .expect("every declared code indexes a label"), None => self.output_id(tgt), + Some(_) => self.output_id(tgt), } } @@ -599,8 +593,8 @@ impl TargetColumns { } .into()); } - // Groups first: an undeclared row falls back to its own id, which is - // still in `pending_ids` at this point. + // Groups first. Ungrouped rows retain their row code and use their + // output id once source ids have been built. self.build_decoy_groups()?; self.build_source_ids()?; let ships_decoys = self.is_decoy.iter().any(|&d| d); @@ -947,9 +941,8 @@ mod tests { assert_ne!(c.decoy_group_code(a), c.decoy_group_code(b)); } - /// A row the file left out of a group has to get a code anyway once any row - /// declared one, since a code indexes the label column directly. It competes - /// alone, under its own id. + /// A row the file left out of a group keeps its row code and output id. + /// It does not need an allocated group label even when other rows link up. #[test] fn a_row_outside_every_declared_group_competes_alone() { let mut c = TargetColumnsBuilder::with_capabilities(TargetCapabilities::default_diann()); @@ -979,6 +972,9 @@ mod tests { assert_eq!(c.decoy_group_code(rows[0]), c.decoy_group_code(rows[1])); assert_ne!(c.decoy_group_code(rows[0]), c.decoy_group_code(rows[2])); assert_eq!(c.decoy_group(rows[2]), SourceId::Text("c")); + let labels = &c.decoy_groups.as_ref().unwrap().labels; + assert_eq!(labels.get(0), Some(SourceId::Text("pair"))); + assert_eq!(labels.get(1), None); } /// Without declared groups a row is its own group, derived rather than diff --git a/rust/timsquery/src/serde/mzspeclib_io.rs b/rust/timsquery/src/serde/mzspeclib_io.rs index 5a9091fb..a69da747 100644 --- a/rust/timsquery/src/serde/mzspeclib_io.rs +++ b/rust/timsquery/src/serde/mzspeclib_io.rs @@ -51,6 +51,7 @@ use mzannotate::fragment::{ use mzannotate::mzdata::params::{ ControlledVocabulary, ParamValue, + Value, }; use mzannotate::mzdata::prelude::PeakCollection; use mzannotate::mzspeclib::record::{ @@ -87,6 +88,7 @@ use crate::models::capabilities::{ UnannotatedPeaks, }; use crate::models::{ + OwnedSourceId, Row, TargetCapabilities, TargetColumnsBuilder, @@ -126,9 +128,9 @@ const SPECTRUM_AGGREGATION_TYPE: u32 = 1_003_065; const INVERSE_REDUCED_MOBILITY: u32 = 1_002_815; /// `MS:1002476|ion mobility drift time`, the drift-tube spelling. const ION_MOBILITY_DRIFT_TIME: u32 = 1_002_476; -/// msspeculator's decoy-to-target link, carried as a project-defined name/value -/// pair because no PSI-MS term identifies the target a decoy was derived from. -const MSSPECULATOR_PAIR_ID: &str = "msspeculator:decoy_pair_id"; +/// msspeculator's numeric competition group, carried as a project-defined +/// name/value pair because no PSI-MS term identifies it. +const MSSPECULATOR_DECOY_GROUP: &str = "msspeculator:decoy_group"; /// Whether `accession` names a decoy spectrum. /// @@ -233,6 +235,8 @@ fn read_counting_degradation( let mut geom = TargetColumnsBuilder::with_capabilities(TargetCapabilities::default_diann()); let mut frag_intens: Vec = Vec::new(); + let mut spectrum_keys = Vec::new(); + let mut related_keys = HashSet::new(); let mut degradation = Degradation::default(); // Read off the header before the iterator takes the parser, which is also // the only place it can be read: it is what an entry declaring nothing about @@ -277,16 +281,17 @@ fn read_counting_degradation( degradation.rows_without_fragments += 1; } - // Both halves of a pair have to land in one namespace. A decoy names its - // target by `` key, so every row names its own key and a - // decoy overrides that with its target's -- not its source id, which is - // a different string entirely. A library that declares no pair therefore - // ends up with one key per row, all singletons, and `seal` drops the - // column rather than storing a group per row. - let group = row - .declared_group - .unwrap_or_else(|| row.key.to_string()) - .into(); + // SpectraST declares a link only on the decoy. Record the target key + // so its row can receive the same group after the full file is read. + let group = match row.declared_group { + Some(DeclaredGroup::Direct(group)) => Some(OwnedSourceId::Numeric(group)), + Some(DeclaredGroup::RelatedSpectrum(key)) => { + related_keys.insert(key); + Some(OwnedSourceId::Numeric(key as u64)) + } + None => None, + }; + spectrum_keys.push(row.key); frag_intens.extend_from_slice(&row.intensities); geom.push_row(Row { @@ -299,10 +304,18 @@ fn read_counting_degradation( entry_name: row.entry_name.as_deref(), is_decoy: row.is_decoy, id: Some(row.id.as_str().into()), - decoy_group: Some(group), + decoy_group: group, }); } + if !related_keys.is_empty() { + for (row, key) in spectrum_keys.into_iter().enumerate() { + if related_keys.contains(&key) && geom.inner.pending_groups[row].is_none() { + geom.inner.pending_groups[row] = Some(OwnedSourceId::Numeric(key as u64)); + } + } + } + let geom = geom.seal(policy.decoys)?; degradation.report(geom.n_rows()); @@ -434,7 +447,7 @@ struct SpectrumRow { charge: u8, mobility: f32, is_decoy: bool, - declared_group: Option, + declared_group: Option, analyte: crate::chemistry::analyte::Analyte, /// Parallel to `intensities`, and pushed straight into the arena. frags: Vec<(IonAnnot, f64)>, @@ -562,7 +575,7 @@ impl SpectrumRow { charge: charge(spectrum, peptidoform), mobility: mobility(spectrum), is_decoy: is_decoy(spectrum)?, - declared_group: declared_group(spectrum), + declared_group: declared_group(spectrum)?, analyte, frags, intensities, @@ -895,29 +908,53 @@ fn is_ms_attribute(attribute: &Attribute, accession: u32) -> bool { } } -/// Which entries compete. msspeculator's project-defined pair id first, then -/// SpectraST's `related spectrum keys`. -/// -/// Both name the target a decoy was derived from, so the target and its decoy -/// land in one group. `None` leaves the row as its own group. -fn declared_group(spectrum: &Spectrum) -> Option { +enum DeclaredGroup { + Direct(u64), + RelatedSpectrum(u32), +} + +/// Which entries compete. SpectraST names a spectrum key on the decoy, while +/// msspeculator writes the group on every member. Both values are integers. +fn declared_group(spectrum: &Spectrum) -> Result, TargetReadingError> { if let Some(param) = spectrum .description .params .iter() - .find(|p| p.name == MSSPECULATOR_PAIR_ID) + .find(|p| p.name == MSSPECULATOR_DECOY_GROUP) { - return Some(param.value.to_string()); + let group = match ¶m.value { + Value::Int(n) => u64::try_from(*n).ok(), + Value::String(s) => s.parse::().ok(), + _ => None, + } + .ok_or_else(|| { + TargetReadingError::SpeclibParse(format!( + "msspeculator:decoy_group must be an integer, got {}", + param.value + )) + })?; + return Ok(Some(DeclaredGroup::Direct(group))); } - spectrum + let related = spectrum .attributes .iter() .flatten() - .find(|a| a.name.accession == mzcv::curie!(MS:1003259)) - .map(|a| match &a.value { - AttributeValue::Scalar(v) => v.to_string(), - other => other.to_string(), - }) + .find(|a| a.name.accession == mzcv::curie!(MS:1003259)); + let Some(related) = related else { + return Ok(None); + }; + let key = match &related.value { + AttributeValue::Scalar(Value::Int(n)) => u32::try_from(*n).ok(), + AttributeValue::Scalar(Value::String(s)) => s.parse::().ok(), + _ => None, + } + .ok_or_else(|| { + TargetReadingError::SpeclibParse(format!( + "MS:1003259|related spectrum keys must name one numeric spectrum key, got {}", + related.value + )) + })?; + Ok(Some(DeclaredGroup::RelatedSpectrum(key))) } /// The theoretical m/z for an annotated peak. @@ -1538,6 +1575,16 @@ mod tests { let groups: std::collections::HashSet<_> = geom.rows().map(|r| geom.decoy_group_code(r)).collect(); assert_eq!(groups.len(), 5, "each decoy competes with its own target"); + let rows: Vec<_> = geom.rows().collect(); + assert_eq!( + geom.decoy_group(rows[6]), + crate::models::SourceId::Numeric(7) + ); + assert_eq!( + geom.decoy_group(rows[9]), + crate::models::SourceId::Numeric(7) + ); + assert_eq!(geom.decoy_group(rows[8]), geom.output_id(rows[8])); } /// The two exports disagree on nearly every term, so reading both is what @@ -1547,6 +1594,10 @@ mod tests { // DIA-NN writes `MS:1000744|selected ion m/z` and no retention time. let diann = arena("diann_export.mzspeclib.txt"); assert!(diann.n_rows() > 0); + assert!( + diann.decoy_groups.is_none(), + "unlinked rows store no group column" + ); let first = diann.rows().next().unwrap(); assert!((diann.precursor_mz(first) - 778.412_96).abs() < 1e-4); assert_eq!(diann.charge(first), 2); @@ -1796,7 +1847,7 @@ mod tests { /// it. It exercises every term msspeculator actually emits rather than the /// terms a vendor happens to use: `MS:1002815|inverse reduced ion mobility` /// for the timsTOF axis, `MS:1000896` for retention, and the - /// project-defined `msspeculator:decoy_pair_id` for the target/decoy + /// project-defined `msspeculator:decoy_group` for the target/decoy /// pairing. /// /// Stored as plain text, not gzipped, so it can be read in a review; the @@ -1811,13 +1862,17 @@ mod tests { let geom = arena("msspeculator_built.mzspeclib.txt"); assert_eq!(geom.n_rows(), 4, "two targets and their two decoys"); assert_eq!(geom.rows().filter(|r| geom.is_decoy(*r)).count(), 2); + assert!( + geom.rows() + .all(|r| matches!(geom.decoy_group(r), crate::models::SourceId::Numeric(_))) + ); let groups: std::collections::HashSet<_> = geom.rows().map(|r| geom.decoy_group_code(r)).collect(); assert_eq!( groups.len(), 2, - "`decoy_pair_id` puts each decoy in its target's group" + "`decoy_group` puts each decoy in its target's group" ); // Four peaks per entry, decoys included: the writer's `Decoy` set diff --git a/rust/timsquery/tests/mzspeclib_files/msspeculator_built.mzspeclib.txt b/rust/timsquery/tests/mzspeclib_files/msspeculator_built.mzspeclib.txt index 38187c3e..cb93eddd 100644 --- a/rust/timsquery/tests/mzspeclib_files/msspeculator_built.mzspeclib.txt +++ b/rust/timsquery/tests/mzspeclib_files/msspeculator_built.mzspeclib.txt @@ -70,7 +70,7 @@ MS:1003065|spectrum aggregation type=MS:1003074|predicted spectrum MS:1003072|spectrum origin type=MS:1003195|unnatural peptidoform decoy spectrum MS:1003061|library spectrum name=VLSAAKPEDR/2 -[3]MS:1003275|other attribute name=msspeculator:decoy_pair_id +[3]MS:1003275|other attribute name=msspeculator:decoy_group [3]MS:1003276|other attribute value=2 MS:1000041|charge state=2 MS:1000744|selected ion m/z=543.30111279 @@ -93,7 +93,7 @@ MS:1000885|protein accession=sp|P00001|TEST1_HUMAN MS:1003061|library spectrum name=VDEPKAASLR/2 MS:1003212|library attribute set name=Decoy -[3]MS:1003275|other attribute name=msspeculator:decoy_pair_id +[3]MS:1003275|other attribute name=msspeculator:decoy_group [3]MS:1003276|other attribute value=2 MS:1000041|charge state=2 MS:1000744|selected ion m/z=543.30111279 @@ -115,7 +115,7 @@ MS:1000885|protein accession=DECOY_sp|P00001|TEST1_HUMAN MS:1003061|library spectrum name=TAWQFYEEKPAGLLSDR/2 -[3]MS:1003275|other attribute name=msspeculator:decoy_pair_id +[3]MS:1003275|other attribute name=msspeculator:decoy_group [3]MS:1003276|other attribute value=1 MS:1000041|charge state=2 MS:1000744|selected ion m/z=1005.99962197 @@ -138,7 +138,7 @@ MS:1000885|protein accession=sp|P00001|TEST1_HUMAN MS:1003061|library spectrum name=TDSLLGAPKEEYFQWAR/2 MS:1003212|library attribute set name=Decoy -[3]MS:1003275|other attribute name=msspeculator:decoy_pair_id +[3]MS:1003275|other attribute name=msspeculator:decoy_group [3]MS:1003276|other attribute value=1 MS:1000041|charge state=2 MS:1000744|selected ion m/z=1005.99962197 diff --git a/rust/timsseek_cli/Cargo.toml b/rust/timsseek_cli/Cargo.toml index f4ced5f3..2fab6c49 100644 --- a/rust/timsseek_cli/Cargo.toml +++ b/rust/timsseek_cli/Cargo.toml @@ -15,10 +15,10 @@ tempfile = { workspace = true } # Library prediction. The digestion, modification and prediction logic lives # there; this crate owns only the flag surface. Pinned to a revision rather than # a branch, and unpublished, so it is a git reference until it reaches crates.io. -msspeculator-inference = { git = "https://github.com/TalusBio/msspeculator", rev = "14634ad065ddcb9750d3c1bb0c817a2c3e332e07" } +msspeculator-inference = { git = "https://github.com/TalusBio/msspeculator", rev = "3fb9dddbd848f4e7ccf3ccfc049690d2f10431a9" } # For `builtin::BUILTIN_PREFIX` and `builtin::names()`, which the inference # crate does not re-export. Same revision, already in the graph. -msspeculator-core = { git = "https://github.com/TalusBio/msspeculator", rev = "14634ad065ddcb9750d3c1bb0c817a2c3e332e07" } +msspeculator-core = { git = "https://github.com/TalusBio/msspeculator", rev = "3fb9dddbd848f4e7ccf3ccfc049690d2f10431a9" } calib_dash = { path = "../calib_dash", optional = true } # The error type of msspeculator's `LibrarySink`, which this crate implements. diff --git a/rust/timsseek_cli/src/build_library.rs b/rust/timsseek_cli/src/build_library.rs index be7f5fe3..45c3518b 100644 --- a/rust/timsseek_cli/src/build_library.rs +++ b/rust/timsseek_cli/src/build_library.rs @@ -15,23 +15,29 @@ use std::path::{ use msspeculator_inference::{ BuiltinModel, + DecoyMethod, Difference, LibraryCheck, LibraryOptions, LibraryProvenance, ModelSource, + PeptideLibraryOptions, ProgressFn, StreamOptions, check_against, check_library, stream_library, write_library, + write_peptide_library, }; use timsseek::DecoyPolicy; use tracing::info; use crate::build_progress::BuildProgress; -use crate::cli::BuildLibraryArgs; +use crate::cli::{ + BuildLibraryArgs, + CliDecoyMethod, +}; use crate::config::{ BuildConfig, LibraryConfig, @@ -116,6 +122,9 @@ impl ResolvedPrediction { #[derive(Debug, Clone, PartialEq)] pub struct ResolvedBuild { pub prediction: ResolvedPrediction, + pub peptides: Option, + pub decoy_method: DecoyMethod, + pub decoy_seed: u64, pub out: PathBuf, pub config_out: Option, pub overwrite: bool, @@ -124,8 +133,19 @@ pub struct ResolvedBuild { /// Fold the command line over the configuration file's `[library]` section. pub fn resolve_build(args: &BuildLibraryArgs, config: &BuildConfig) -> ResolvedBuild { let library = library_overrides(args, config.library.clone().unwrap_or_default()); + let input = args + .fasta + .clone() + .or_else(|| args.peptides.clone()) + .expect("clap requires FASTA or peptide TSV"); ResolvedBuild { - prediction: resolve_prediction(args.fasta.clone(), &library), + prediction: resolve_prediction(input, &library), + peptides: args.peptides.clone(), + decoy_method: match args.decoy_method.unwrap_or(CliDecoyMethod::PseudoReverse) { + CliDecoyMethod::PseudoReverse => DecoyMethod::PseudoReverse, + CliDecoyMethod::Shuffle => DecoyMethod::Shuffle, + }, + decoy_seed: args.decoy_seed.unwrap_or(0), out: args.out.clone(), config_out: sidecar_path(args), // No configuration counterpart: replacing a file is a decision about @@ -219,6 +239,9 @@ pub(crate) fn resolve_search_build( ) -> ResolvedBuild { ResolvedBuild { prediction: resolve_search_prediction(fasta, library), + peptides: None, + decoy_method: DecoyMethod::PseudoReverse, + decoy_seed: 0, config_out: Some(default_sidecar(&out)), out, overwrite, @@ -234,7 +257,17 @@ pub(crate) fn resolve_search_build( /// limitation is stated instead of discovered. fn reject_remote_paths(resolved: &ResolvedBuild) -> Result<(), CliError> { for (flag, path) in [ - ("--fasta", &resolved.prediction.fasta), + ( + if resolved.peptides.is_some() { + "--peptides" + } else { + "--fasta" + }, + resolved + .peptides + .as_ref() + .unwrap_or(&resolved.prediction.fasta), + ), ("--out", &resolved.out), ] { let uri = path.to_string_lossy(); @@ -439,9 +472,13 @@ pub fn run(resolved: &ResolvedBuild) -> Result<(), CliError> { reject_existing_output(resolved)?; let model = parse_model_source(&resolved.prediction.model)?; + let input = resolved + .peptides + .as_ref() + .unwrap_or(&resolved.prediction.fasta); info!( "Predicting a library from {} with {}", - resolved.prediction.fasta.display(), + input.display(), resolved.prediction.model ); // Held in a binding: the callback borrows it, so a temporary would be @@ -465,12 +502,28 @@ pub fn run(resolved: &ResolvedBuild) -> Result<(), CliError> { report_rebuild_check(&resolved.out, resolved.config_out.as_deref(), provenance); } }; - let stats = write_library(&LibraryOptions { - out: &pending_out, - config_out: pending_config.as_deref(), - stream: stream_options(&resolved.prediction, model, &report), - before_writing: Some(&report_rebuild), - }) + let stats = if let Some(peptides) = resolved.peptides.as_deref() { + write_peptide_library(&PeptideLibraryOptions { + model, + peptides, + out: &pending_out, + config_out: pending_config.as_deref(), + min_intensity: resolved.prediction.min_intensity, + max_fragments: resolved.prediction.max_fragments, + generate_decoys: resolved.prediction.decoys, + decoy_method: resolved.decoy_method, + decoy_seed: resolved.decoy_seed, + progress: Some(&report), + before_writing: Some(&report_rebuild), + }) + } else { + write_library(&LibraryOptions { + out: &pending_out, + config_out: pending_config.as_deref(), + stream: stream_options(&resolved.prediction, model, &report), + before_writing: Some(&report_rebuild), + }) + } .map_err(|e| CliError::LibraryBuild { source: format!("{}: {e:#}", resolved.out.display()), })?; diff --git a/rust/timsseek_cli/src/cli.rs b/rust/timsseek_cli/src/cli.rs index 4336b539..34a7bce4 100644 --- a/rust/timsseek_cli/src/cli.rs +++ b/rust/timsseek_cli/src/cli.rs @@ -70,7 +70,7 @@ pub enum Command { /// Deprecated spelling of a bare invocation; use `timsseek `. /// Removed after 2026-12-31. Search(Box), - /// Predict a spectral library from a FASTA and write it. + /// Predict a spectral library from a FASTA or peptide TSV and write it. BuildLibrary(Box), } @@ -294,13 +294,28 @@ pub struct SearchArgs { #[derive(clap::Args, Debug)] pub struct BuildLibraryArgs { /// Sequence database to digest. - #[arg(long, value_name = "PATH")] - pub fasta: PathBuf, + #[arg( + long, + value_name = "PATH", + required_unless_present = "peptides", + conflicts_with = "peptides" + )] + pub fasta: Option, + + /// TSV with proforma and protein_ids columns; optional decoy and decoy_group. + /// ProForma must include charge, e.g. PEPC[UNIMOD:4]IDEK/2. + #[arg( + long, + value_name = "PATH", + required_unless_present = "fasta", + conflicts_with = "fasta" + )] + pub peptides: Option, /// Where to write the library. The suffix picks the format: /// `.mzspeclib.txt` writes mzSpecLib, anything else writes DIA-NN TSV, and /// a trailing `.gz` compresses either. `.mzspeclib.txt.gz` is the - /// recommended spelling. + /// recommended spelling. Peptide TSV input requires mzSpecLib output. #[arg(long, short = 'o', value_name = "PATH")] pub out: PathBuf, @@ -313,31 +328,39 @@ pub struct BuildLibraryArgs { #[arg(long, value_name = "MODEL")] pub model: Option, - #[arg(long, value_name = "N")] + #[arg(long, value_name = "N", conflicts_with = "peptides")] pub missed_cleavages: Option, - #[arg(long, value_name = "N")] + #[arg(long, value_name = "N", conflicts_with = "peptides")] pub min_length: Option, - #[arg(long, value_name = "N")] + #[arg(long, value_name = "N", conflicts_with = "peptides")] pub max_length: Option, - #[arg(long, value_name = "Z")] + #[arg(long, value_name = "Z", conflicts_with = "peptides")] pub min_charge: Option, - #[arg(long, value_name = "Z")] + #[arg(long, value_name = "Z", conflicts_with = "peptides")] pub max_charge: Option, /// Fixed modification rule, repeatable. - #[arg(long = "fixed-mod", value_name = "TARGETS[MOD]")] + #[arg( + long = "fixed-mod", + value_name = "TARGETS[MOD]", + conflicts_with = "peptides" + )] fixed_mods: Option>, /// Predict with no fixed modifications at all. /// /// Needed because the default is carbamidomethyl and an empty /// `--fixed-mod` list cannot be spelled: repeating a flag zero times is /// indistinguishable from not passing it. - #[arg(long, conflicts_with = "fixed_mods")] + #[arg(long, conflicts_with_all = ["fixed_mods", "peptides"])] pub no_fixed_mods: bool, /// Variable modification rule, repeatable. - #[arg(long = "variable-mod", value_name = "TARGETS[MOD]")] + #[arg( + long = "variable-mod", + value_name = "TARGETS[MOD]", + conflicts_with = "peptides" + )] pub variable_mods: Option>, - #[arg(long, value_name = "N")] + #[arg(long, value_name = "N", conflicts_with = "peptides")] pub max_variable_mods: Option, // No acquisition or chromatography context. The model artifact has its own @@ -363,6 +386,14 @@ pub struct BuildLibraryArgs { #[arg(long, overrides_with = "decoys")] no_decoys: bool, + /// Decoy method for peptide TSV input. + #[arg(long, value_enum, requires = "peptides")] + pub decoy_method: Option, + + /// Reproducible seed for peptide TSV decoy generation. + #[arg(long, value_name = "N", requires = "peptides")] + pub decoy_seed: Option, + /// Replace the library if it already exists. /// /// Predicting a proteome takes minutes, so overwriting one by accident is @@ -380,6 +411,12 @@ pub struct BuildLibraryArgs { pub no_config_out: bool, } +#[derive(Debug, Clone, Copy, PartialEq, Eq, ValueEnum)] +pub enum CliDecoyMethod { + PseudoReverse, + Shuffle, +} + impl BuildLibraryArgs { /// `Some(true)`/`Some(false)` when either decoy flag was given, `None` when /// neither was -- which is what lets the configuration file supply the value diff --git a/rust/timsseek_cli/src/predicted_library.rs b/rust/timsseek_cli/src/predicted_library.rs index def2e99b..1496ad1b 100644 --- a/rust/timsseek_cli/src/predicted_library.rs +++ b/rust/timsseek_cli/src/predicted_library.rs @@ -35,6 +35,7 @@ use timsquery::ion::{ NeutralLoss, }; use timsquery::models::{ + OwnedSourceId, Row, TargetCapabilities, TargetColumnsBuilder, @@ -223,10 +224,9 @@ struct PredictedRow { /// same precursor under `MS:1003061|library spectrum name`, so a row keeps /// one name across both routes. Stored independently of chemistry. id: String, - /// The source peptide's pair id, `None` when decoys are off. All modified - /// forms and all charges of one peptide share it, so a group is a peptide - /// rather than a target/decoy couple. - group: Option, + /// Numeric group ID shared by a target and its decoy. Each modified form + /// and charge has its own group; absent when decoys are off. + group: Option, analyte: timsquery::chemistry::analyte::Analyte, precursor_mz: f64, charge: u8, @@ -262,7 +262,9 @@ impl PredictedRow { // this is the decoy's own sequence rather than its target's. analyte: timsquery::chemistry::analyte::Analyte::from_sequence(row.proforma), id, - group: row.decoy_pair_id.map(|pair| pair.to_string()), + group: row + .decoy_group + .map(|group| OwnedSourceId::Numeric(group as u64)), precursor_mz: row.precursor_mz, charge, // `rt` is the quantity the mzSpecLib writer puts under the term the @@ -337,7 +339,10 @@ impl ArenaRows { if !self.decoys.accepts(row.is_decoy) { return; } - let group = row.group.clone().unwrap_or_else(|| row.id.clone()); + let group = row + .group + .clone() + .unwrap_or_else(|| OwnedSourceId::Text(row.id.clone())); self.frag_intens.extend_from_slice(&row.intensities); self.geom.push_row(Row { precursor_mz: row.precursor_mz, @@ -356,7 +361,7 @@ impl ArenaRows { entry_name: Some(&row.id), is_decoy: row.is_decoy, id: Some(row.id.clone().into()), - decoy_group: Some(group.into()), + decoy_group: Some(group), }); } @@ -490,7 +495,7 @@ mod tests { peptide: &self.peptide, proforma: &self.proforma, decoy, - decoy_pair_id: pair, + decoy_group: pair, charge, precursor_mz: 500.25 + charge as f64, neutral_mass: 998.5, @@ -618,7 +623,7 @@ mod tests { } #[test] - fn a_target_and_its_decoy_that_share_a_pair_id_land_in_one_competition_group() { + fn a_target_and_its_decoy_share_one_numeric_competition_group() { let target = Fixture::new("PEPTIDEK", "PEPTIDEK"); let decoy = Fixture::new("PDITPEEK", "PDITPEEK"); let lib = build( @@ -1015,7 +1020,7 @@ mod tests { /// /// Every such quantity rather than a sample of them, because the two routes /// spell each one independently and any of them could drift. Two charge - /// states of one peptide, which msspeculator gives one `decoy_pair_id`, so + /// states of one peptide, which msspeculator gives distinct `decoy_group` IDs, so /// the group column is stored rather than derived and the comparison of it /// means something. /// diff --git a/rust/timsseek_cli/tests/build_library.rs b/rust/timsseek_cli/tests/build_library.rs index 35efe339..c6228aef 100644 --- a/rust/timsseek_cli/tests/build_library.rs +++ b/rust/timsseek_cli/tests/build_library.rs @@ -59,3 +59,214 @@ fn build_library_writes_a_readable_library_and_sidecar() { assert!(stderr.contains("built with different settings"), "{stderr}"); assert!(stderr.contains("fragments.max_fragments"), "{stderr}"); } + +#[test] +fn peptide_tsv_preserves_modified_charges_and_supplied_decoy_groups() { + let dir = tempfile::tempdir().unwrap(); + let input = dir.path().join("peptides.tsv"); + std::fs::write( + &input, + concat!( + "proforma\tprotein_ids\tdecoy\tdecoy_group\n", + "PEC[UNIMOD:4]TIDEK/2\tP1;P2\tfalse\t1\n", + "PECTIDEK/3\tP1\tfalse\t2\n", + "PEK[UNIMOD:259]TIDEK/2\tPRTC\tfalse\t3\n", + "PECTDIEK/3\tP1\ttrue\t2\n", + ), + ) + .unwrap(); + let out = dir.path().join("peptides.mzspeclib.txt"); + let result = Command::new(env!("CARGO_BIN_EXE_timsseek")) + .args([ + "build-library", + "--peptides", + input.to_str().unwrap(), + "--out", + out.to_str().unwrap(), + "--no-decoys", + "--max-fragments", + "4", + ]) + .output() + .unwrap(); + assert!( + result.status.success(), + "{}", + String::from_utf8_lossy(&result.stderr) + ); + let text = std::fs::read_to_string(&out).unwrap(); + for peptide in [ + "PEC[UNIMOD:4]TIDEK/2", + "PECTIDEK/3", + "PEK[UNIMOD:259]TIDEK/2", + ] { + assert!(text.contains(peptide), "missing {peptide}"); + } + assert_eq!( + text.matches("[3]MS:1003276|other attribute value=2\n") + .count(), + 2 + ); + let table = timsquery::serde::read_targets_with( + &out, + timsquery::models::capabilities::LoadPolicy::default(), + ) + .unwrap(); + let timsquery::serde::TargetTable::Mzpaf { geom, .. } = table else { + panic!("predicted fragments should have mzPAF annotations"); + }; + assert_eq!(geom.n_rows(), 4); + assert_eq!(geom.rows().filter(|r| geom.is_decoy(*r)).count(), 1); + assert!( + geom.rows() + .all(|r| matches!(geom.decoy_group(r), timsquery::models::SourceId::Numeric(_))) + ); + assert_eq!( + geom.rows() + .map(|r| geom.decoy_group_code(r)) + .collect::>() + .len(), + 3 + ); + assert!(geom.rows().any(|r| geom.charge(r) == 3)); +} + +#[test] +fn peptide_tsv_generates_seeded_decoys_and_assigns_groups() { + let dir = tempfile::tempdir().unwrap(); + let input = dir.path().join("peptides.tsv"); + std::fs::write( + &input, + "proforma\tprotein_ids\nPEC[UNIMOD:4]TIDEK/2\tP1\nPECTIDEK/2\tP1\n", + ) + .unwrap(); + let out = dir.path().join("peptides.mzspeclib.txt"); + let result = Command::new(env!("CARGO_BIN_EXE_timsseek")) + .args([ + "build-library", + "--peptides", + input.to_str().unwrap(), + "--out", + out.to_str().unwrap(), + "--decoy-method", + "shuffle", + "--decoy-seed", + "42", + "--max-fragments", + "4", + ]) + .output() + .unwrap(); + assert!( + result.status.success(), + "{}", + String::from_utf8_lossy(&result.stderr) + ); + let table = timsquery::serde::read_targets_with( + &out, + timsquery::models::capabilities::LoadPolicy::default(), + ) + .unwrap(); + let timsquery::serde::TargetTable::Mzpaf { geom, .. } = table else { + panic!("predicted fragments should have mzPAF annotations"); + }; + assert_eq!(geom.n_rows(), 4); + assert_eq!(geom.rows().filter(|r| geom.is_decoy(*r)).count(), 2); + assert!( + geom.rows() + .all(|r| matches!(geom.decoy_group(r), timsquery::models::SourceId::Numeric(_))) + ); + assert_eq!( + geom.rows() + .map(|r| geom.decoy_group_code(r)) + .collect::>() + .len(), + 2 + ); + let sidecar: serde_json::Value = + serde_json::from_slice(&std::fs::read(format!("{}.config.json", out.display())).unwrap()) + .unwrap(); + assert_eq!(sidecar["decoys"]["method"], "shuffle"); + assert_eq!(sidecar["decoys"]["seed"], 42); +} + +#[test] +fn peptide_tsv_retries_pseudo_reverse_after_a_target_collision() { + let dir = tempfile::tempdir().unwrap(); + let input = dir.path().join("peptides.tsv"); + std::fs::write( + &input, + "proforma\tprotein_ids\nPEPTIDEK/2\tP1\nPEDITPEK/3\tP2\n", + ) + .unwrap(); + let out = dir.path().join("peptides.mzspeclib.txt"); + let result = Command::new(env!("CARGO_BIN_EXE_timsseek")) + .args([ + "build-library", + "--peptides", + input.to_str().unwrap(), + "--out", + out.to_str().unwrap(), + "--decoy-method", + "pseudo-reverse", + ]) + .output() + .unwrap(); + assert!( + result.status.success(), + "{}", + String::from_utf8_lossy(&result.stderr) + ); + assert!(!String::from_utf8_lossy(&result.stderr).contains("skipped")); + let table = timsquery::serde::read_targets_with( + &out, + timsquery::models::capabilities::LoadPolicy::default(), + ) + .unwrap(); + let timsquery::serde::TargetTable::Mzpaf { geom, .. } = table else { + panic!("predicted fragments should have mzPAF annotations"); + }; + assert_eq!(geom.n_rows(), 4); + assert_eq!(geom.rows().filter(|r| geom.is_decoy(*r)).count(), 2); + assert!( + std::fs::read_to_string(&out) + .unwrap() + .contains("PEPITDEK/2") + ); +} + +#[test] +fn peptide_tsv_keeps_unpaired_targets_when_every_reversal_collides() { + let dir = tempfile::tempdir().unwrap(); + let input = dir.path().join("peptides.tsv"); + std::fs::write(&input, "proforma\tprotein_ids\nPEEEEEEK/2\tP1\n").unwrap(); + let out = dir.path().join("peptides.mzspeclib.txt"); + let result = Command::new(env!("CARGO_BIN_EXE_timsseek")) + .args([ + "build-library", + "--peptides", + input.to_str().unwrap(), + "--out", + out.to_str().unwrap(), + "--decoy-method", + "pseudo-reverse", + ]) + .output() + .unwrap(); + assert!( + result.status.success(), + "{}", + String::from_utf8_lossy(&result.stderr) + ); + assert!(String::from_utf8_lossy(&result.stderr).contains("skipped 1 generated decoys")); + let table = timsquery::serde::read_targets_with( + &out, + timsquery::models::capabilities::LoadPolicy::default(), + ) + .unwrap(); + let timsquery::serde::TargetTable::Mzpaf { geom, .. } = table else { + panic!("predicted fragments should have mzPAF annotations"); + }; + assert_eq!(geom.n_rows(), 1); + assert_eq!(geom.rows().filter(|r| geom.is_decoy(*r)).count(), 0); +} diff --git a/uv.lock b/uv.lock index 47d2cf6c..2acf9f41 100644 --- a/uv.lock +++ b/uv.lock @@ -676,7 +676,7 @@ wheels = [ [[package]] name = "timsseek-workspace" -version = "0.34.1" +version = "0.35.0" source = { virtual = "." } [package.dev-dependencies]