diff --git a/skills/mascot-parser/SKILL.md b/skills/mascot-parser/SKILL.md index 588058c..1ba563d 100644 --- a/skills/mascot-parser/SKILL.md +++ b/skills/mascot-parser/SKILL.md @@ -1,8 +1,6 @@ --- name: mascot-parser description: Mascot Parser (msparser) Skill -version: 1.0.0 -license: Apache-2.0 --- # Mascot Parser (msparser) Skill @@ -93,7 +91,7 @@ if not resfile.isValid(): ```python params = resfile.params() # returns ms_searchparams -params.getDB() # database name +params.getDB() # PRIMARY database name (see multi-DB note below) params.getCOM() # search title params.getTOL() # peptide tolerance params.getTOLU() # tolerance units @@ -109,6 +107,26 @@ params.getCHARGE() # charge state params.getINSTRUMENT() # instrument type ``` +#### Multi-database searches: `getDB()` is not enough + +Mascot supports searching multiple databases in a single submission (`DB`, `DB2`, `DB3`, ... in the search form, common when Distiller users tick multiple databases in the search dialog). **`params.getDB()` returns only the primary (first) database — silently hiding the rest.** Reporting "the database is X" based on `getDB()` alone is a confident-sounding bug when the actual search used more than one. + +Always enumerate all databases. The indexed `getDB(i)` is **1-based** (consistent with queries, peptide ranks, fixed/variable mod indices, MS1 match iteration, etc.; only component indexing is 0-based): + +```python +n = params.getNumberOfDatabases() +dbs = [params.getDB(i) for i in range(1, n + 1)] # 1-based +``` + +| Call | Behaviour | +|------|-----------| +| `params.getDB()` | Primary DB only — equivalent to `params.getDB(1)` | +| `params.getDB(0)` | Returns `''` (no DB at index 0) | +| `params.getDB(i)` for `1 <= i <= getNumberOfDatabases()` | Database name at slot `i` | +| `params.getNumberOfDatabases()` | Total count — always check this before reporting "the database" | + +Any summary tool should print all databases when `getNumberOfDatabases() > 1`. Distiller-side note: the `` blocks inside a `.rov`'s `rover_data` stream carry empty `Database` attributes — those are GUI preference templates, not the executed search. The authoritative DB list is on the embedded `.dat`'s `params` object. + ### Creating Results: The Two-Argument Constructor (Recommended) **Use `ms_mascotresults_params` with the two-argument constructor.** The old multi-argument constructor is obsolete. @@ -314,6 +332,126 @@ while prot: prot = results.getHit(hit) ``` +**⚠️ This loop sees only family representatives — see the next section if you +need every identified protein.** + +## Reading all proteins: family members, not just `getHit()` + +`ms_peptidesummary.getHit(n)` (and `ms_proteinsummary.getHit(n)`) returns **only +the protein-family REPRESENTATIVE** of each family. When result clustering is on +(it is, by default, in the canonical reader — `setProteinFamilySwitch(1)` / +`MSRES_CLUSTER_PROTEINS`), homologous proteins that group into a family are +stored as family **MEMBERS** and are **invisible to a `getHit()`-only loop**. A +reps-only extractor silently drops them — they were identified, they pass FDR, +and they never appear in your output. + +**Verified failure (Mascot Server 3.1.241.27):** in a published *Vibrio +cholerae* 2D-TPP study, the three validated chemoreceptor hits — UniProt +**Q9KSE4** (VC1313), **Q9KQ43** (VC2161), **Q9KS54** (VC1406) — cluster into a +single 38-member family whose representative is **Q9KVD0** (VC2161 is family 251, +member 7). A `getHit()`-only loop dropped ~189 member proteins across 5 samples, +including all three biological hits, even though every one passed Percolator at +q = 0. On a routine bacterial TMT search, walking members recovered 66 extra +proteins that `getHit` alone missed (2690 vs 2624) — exactly +`getNumberOfFamilyMembers()`. + +### Walking family members + +There are two APIs. **`getNextFamilyProtein(masterHit, id)` is the simpler +per-family iterator** and is what most code should use: + +```python +def iter_all_proteins(results): + """Yield (protein, representative_accession_or_None) for EVERY protein.""" + n = results.getNumberOfHits() + for hit in range(1, n + 1): # 1-based + rep = results.getHit(hit) + if rep is None: + break + yield rep, None # the representative + rep_acc = rep.getAccession() + member_id = 1 + while True: + member = results.getNextFamilyProtein(hit, member_id) + if member is None: # no more members in this family + break + yield member, rep_acc # a homologous family member + member_id += 1 +``` + +`getNextFamilyProtein` returns `None` at `member_id = 1` when a family has no +extra members, so this degrades cleanly to a plain `getHit` walk on +non-clustered results. Family members expose the same `ms_protein` interface as +representatives (`getAccession`, `getScore`, `getNumPeptides`, +`getPeptideQuery`/`getPeptideP`, quant lookups all work on them). + +The lower-level alternative is **`getHitAndFamilyMember(prot, hitAndFamily, +rules)`**, where the 3rd argument is a `ms_mascotresults::hitAndFamily_t` cursor +(NOT a simple int) and `rules` combines `FC_PROTEIN_IGN_*` flags +(`FC_PROTEIN_IGN_SUBSETS`, `FC_PROTEIN_IGN_SAMESETS`, `FC_PROTEIN_IGN_FAMILY`, +`FC_PROTEIN_IGN_MASK`). Use it when you need fine control over which +subset/sameset/family relationships to traverse. Note +`getNumberOfFamilyMembers()` is a **GLOBAL** total across the whole result, not a +per-family count — don't use it to bound a per-family loop; iterate +`getNextFamilyProtein` until it returns `None` instead. + +### Gotcha summary + +| Call | What it returns | +|------|-----------------| +| `getHit(n)` | family REPRESENTATIVE only — members are hidden | +| `getNextFamilyProtein(hit, id)` | the `id`-th member of family `hit`; `None` past the last | +| `getNumberOfFamilyMembers()` | GLOBAL member total (all families) — not a per-family bound | +| `getNumberOfHits()` | number of representatives (= number of families) | + +## Importing Percolator results (not in the `.msr`) + +**Percolator / ML-rescored results are NOT stored in the SQLite `.msr`.** They +live in server cache files written when Percolator ran: + +``` +/data/cache////..target.pop # FDR-passing target PSMs + ..decoy.pop # decoy PSMs + ..pip # Percolator input +``` + +To get Percolator scores through msparser you must **import** them — the +canonical reader (`master_results_2.pl` → +`/perl64/site/lib/PeptideSummary/Util.pm::make_ms_mascotresults_params`) +builds the summary with clustering and Percolator on by default: +`setProteinFamilySwitch(1)`, flags `MSRES_CLUSTER_PROTEINS | MSRES_SHOW_SUBSETS`, +Percolator enabled, `MSPEPSUM_USE_CACHE`, `setIgnoreIonsScoreBelow(...)`, +`setUsePeptideSummary(1)`. msparser computes its own cache hash and often looks +under a *different* path than the server wrote, so on real installs you must +stage the server's actual `.pop`/`.pip` files where msparser expects them (glob +the server cache, prime `setPercolatorFeatures`, copy to +`getPercolatorFileNames()`, then set `MSPEPSUM_PERCOLATOR`) — see +[Common Recipes](references/common-recipes.md) recipe 20. + +### The `.target.pop` is directly usable and authoritative + +For FDR-passing PSMs you often don't need msparser's summary at all — the +`.target.pop` is a tab-separated table you can read directly. Columns: + +| Column | Meaning | +|--------|---------| +| `PSMId` | `query:N;rank:M` — **N is the msparser query number**, so the spectrum/peaks are reachable via `ms_inputquery(resfile, N)` | +| `score` | Percolator score | +| `q-value` | Percolator q-value — filter `q <= target FDR` | +| `posterior_error_prob` | PEP | +| `peptide` | `X.SEQUENCE.X` (flanking residues either side of dots) | +| `proteinIds` | space-separated accessions; a peptide mapping to exactly one accession is unique to it | + +**Recommended robust pattern — drive per-protein identification/quant off the +`.target.pop` PSMs:** read the rows, keep `q <= target FDR`, group by +`proteinIds`. This captures representatives **AND** family members automatically +and bypasses the `getHit()` blind spot entirely. For quant, parse `query:N` out +of `PSMId`, read that spectrum with `ms_inputquery(resfile, N)`, sum the reporter +ions, and add the PSM's reporter sums to every protein in `proteinIds`. A +complete working TMT extractor built this way (PSMId parsing, cache glob, +`ms_inputquery` peak access) is in the project at +`TPP2D/scripts/msr_to_quant_summary_pop.py`. + ### Getting Peptides for a Protein ```python @@ -415,14 +553,15 @@ for u in range(1, 1 + results.getNumberOfUnassigned()): 1. **`createResfile` not `createResFile`** - the 'f' is lowercase 2. **Always check `isValid()`** after creating objects, then `getLastErrorString()` for details -3. **1-based indexing** - queries, peptide ranks, peak indices all start at 1 +3. **1-based indexing** - queries, peptide ranks, peak indices, fixed/variable mod indices, MS1 match iteration, AND `params.getDB(i)` for multi-DB enumeration all start at 1. (Component indexing in quantitation is the exception — that's 0-based.) +3a. **Multi-DB searches**: `params.getDB()` returns only the primary database. Always enumerate via `params.getNumberOfDatabases()` + `params.getDB(i)` for `i in 1..N`. See [Search Parameters → Multi-database searches](#search-parameters) above. Reporting a single DB name from `getDB()` alone is a known footgun. 4. **Use `PROXY_TYPE_NO_PROXY`** for localhost connections 5. **Use the two-argument constructor** for `ms_peptidesummary` / `ms_proteinsummary` — the multi-argument constructor is obsolete -6. **Protein iteration** - `getHit(n)` returns `None` when no more hits; start at 1 +6. **Protein iteration** - `getHit(n)` returns `None` when no more hits; start at 1. **But `getHit()` returns family REPRESENTATIVES only** — homologous family members are hidden. To read every identified protein, also walk `getNextFamilyProtein(hit, id)`. See [Reading all proteins: family members, not just `getHit()`](#reading-all-proteins-family-members-not-just-gethit). 7. **Peptide duplicate check** - always check `getPeptideDuplicate(i) != DUPE_DuplicateSameQuery` before processing 8. **Error handling pattern**: check `isValid()` first, then `getLastErrorString()`, then `clearAllErrors()` 9. **ms_datfile can read from URL** - pass `ms_connection_settings` with session ID as third argument -10. **Result files** are `.dat` (text) or `.msr` (binary) format; both opened the same way +10. **Result files** are `.dat` (text) or `.msr` (binary) format; both opened the same way — **EXCEPT** the new Mascot 3.x SQLite-backed `.msr`: on **large DIA** results `ms_peptidesummary`/`ms_proteinsummary` can return 0 hits or hang, in which case read the SQLite tables directly ([SQLite-backed `.msr`](references/sqlite-msr.md)). **Do not mistake the more common "missing proteins" causes for a SQLite-format problem:** (a) family members hidden behind `getHit()` (see [Reading all proteins](#reading-all-proteins-family-members-not-just-gethit)); (b) Percolator/FDR-passing hits that live in cache `.pop` files, not the `.msr` at all (see [Importing Percolator results](#importing-percolator-results-not-in-the-msr)). Neither is fixed by switching to raw SQL. 11. **FDR target is a decimal** - `setTargetFDR(0.01)` for 1%, not `setTargetFDR(1)` 12. **Significance test** - `getIonsScore() >= getPeptideThreshold(q, 20, rank)` is equivalent to `getExpectationValue() <= 0.05` @@ -441,3 +580,4 @@ The msparser SDK ships with example scripts in `/example_python/`, - [Common Recipes](references/common-recipes.md) - copy-paste code patterns for common tasks (Python) - [Server Configuration](references/server-config.md) - directory layout, authentication flow - [Obsolete Examples](references/obsolete-examples.md) - which SDK example scripts need modernizing +- [SQLite-backed `.msr`](references/sqlite-msr.md) - Mascot 3.x writes `.msr` as a SQLite DB; `ms_peptidesummary` fails on large DIA results; full table schema + direct-SQL recipes diff --git a/skills/mascot-parser/references/common-recipes.md b/skills/mascot-parser/references/common-recipes.md index 3d9e3f7..8125d29 100644 --- a/skills/mascot-parser/references/common-recipes.md +++ b/skills/mascot-parser/references/common-recipes.md @@ -597,3 +597,202 @@ def export_mgf(resfile, output_path): f.write("END IONS\n\n") ``` + +--- + +## 18. CRITICAL: chdir to Mascot CGI dir before opening result files + +`msparser` resolves the unimod XML schema (`unimod_2.xsd`) via a *relative path* +`../html/xmlns/schema/unimod_2/unimod_2.xsd`. If your CWD is anywhere outside +`/cgi/`, the schema fails to load — `getLastErrorString()` reports +"Failed to load unimod xml file" and `pep.getPeptideStr()` returns empty +strings on every peptide (silent failure, no exception). PSM counts via +`getNumberOfHits()` come back as 0. + +**Always:** +```python +import os +os.chdir(r"C:\inetpub\mascot\cgi") # Linux: /usr/local/mascot/cgi +resfile = msparser.ms_mascotresfilebase.createResfile(path) +``` + +This is not documented in the SDK guide but is required by every script that +opens local result files. The QC dashboard (`mascot_qc_report.py`) does this +explicitly; copy that pattern. + +--- + +## 19. Get PSM/sequence counts at a target FDR (the right way) + +The naive call `getNumHitsAboveIdentity(0.05, ...)` returns counts at the *raw* +identity threshold (p ≤ 0.05). To match the CSV-export-header convention +(counts at 1 % FDR-adjusted threshold) use `getThresholdForFDRAboveIdentity`, +which returns a **list of 5 values** `[ok, achieved_FDR, sigLevel, n_target, n_decoy]`: + +```python +M = msparser.ms_mascotresults +target_fdr = 0.01 + +# 3-arg overload — db_match_type ONLY accepts DS_IDENTITY (DS_HOMOLOGY returns +# [False, -1, -1, -1, -1]). +def fdr_count(method, count_type): + r = method(target_fdr, count_type, M.DS_IDENTITY) + if not r or len(r) < 5 or not r[0]: + return 0, 0, 0.0 + return int(r[3]), int(r[4]), float(r[1]) # target, decoy, achieved_fdr + +n_psm_target, n_psm_decoy, fdr_psm = fdr_count( + results.getThresholdForFDRAboveIdentity, M.DS_COUNT_PSM) +n_seq_target, n_seq_decoy, fdr_seq = fdr_count( + results.getThresholdForFDRAboveIdentity, M.DS_COUNT_SEQUENCE) +n_psm_target_h, n_psm_decoy_h, _ = fdr_count( + results.getThresholdForFDRAboveHomology, M.DS_COUNT_PSM) +``` + +Constants: `DS_COUNT_PSM = 0`, `DS_COUNT_SEQUENCE = 1`, `DS_IDENTITY = 0`, +`DS_HOMOLOGY = 1` (the latter is used as a `db_match_type` by some methods but +**not** by `getThresholdForFDR*` — pass `DS_IDENTITY` there). + +Why this matters: `setTargetFDR(0.01)` only affects the *threshold function used +during summary creation*. The count APIs always take an explicit sig_level, +and the right one for "1 % FDR" is what `getThresholdForFDR*` returns — not 0.05. + +**There is no `pep.isDecoy()` method.** Don't try to detect decoys by walking +peptides and checking accession prefixes; use the API above, which classifies +hits via the result file's embedded decoy flag. + +--- + +## 20. Read Percolator-rescored results: stage cache files in a tmp dir + +Setting `MSPEPSUM_PERCOLATOR` on a result-params object alone is **not enough** +— `ms_peptidesummary(resfile, params)` will create with `getNumberOfHits() = 0` +and any `getPeptide()` call will raise +"Attempting to call function ... before createSummary() has completed". + +The hash msparser computes locally on `mascot.dat` may differ from the hash +the server used when it wrote the `.target.pop` / `.decoy.pop` / `.pip` files +in `/data/cache/YYYY/MM//`. The fix (mirrors what +the QC dashboard does): + +1. Open the resfile pointed at a **fresh tmp dir** as cache_dir. +2. Call `setPercolatorFeatures` with the same `ml_adapter_param` values that + the search used (e.g. `MS2Rescore.ms2pip_model=HCD2019`). +3. Read the *expected* filenames via `getPercolatorFileNames()` — returns + tuple `(pip, target.pop, decoy.pop)`. +4. Copy the real server cache files to the expected paths in the tmp dir. +5. Set `MSPEPSUM_PERCOLATOR` and create the summary. + +```python +import os, glob, shutil, tempfile +from pathlib import Path + +def open_with_percolator(res_path, adapter_params): + os.chdir(r"C:\inetpub\mascot\cgi") # see recipe 18 + + # Find server-written cache for this resfile + res_name = os.path.basename(res_path) + real_target = glob.glob(rf"C:\inetpub\mascot\data\cache\*\*\*\{res_name}.*.target.pop") + if not real_target: + return None # search wasn't Percolator-rescored + real_dir = os.path.dirname(real_target[0]) + final_hash = os.path.basename(real_target[0]).rsplit('.target.pop', 1)[0].split('.')[-1] + real_decoy = os.path.join(real_dir, f"{res_name}.{final_hash}.decoy.pop") + real_pip = os.path.join(real_dir, f"{res_name}.{final_hash}.pip") + + # Open into tmp cache dir, set features, get expected filenames + tmp = tempfile.mkdtemp(prefix="perc_cache_") + df = msparser.ms_datfile(r"C:\inetpub\mascot\config\mascot.dat") + opts = df.getMascotOptions() + rf = msparser.ms_mascotresfilebase.createResfile(res_path, 0, "", 0, tmp) + opts.setPercolatorExeFlags(opts.getPercolatorRtFlags(rf.hasRT(), opts.isPercolatorUseRT())) + vs = msparser.VectorString() + for p in adapter_params: vs.append(p) + rf.setPercolatorFeatures(opts, "", vs) + + expected_pip, expected_target, expected_decoy = rf.getPercolatorFileNames() + for src, dst in [(real_pip, expected_pip), + (real_target, expected_target), + (real_decoy, expected_decoy)]: + Path(dst).parent.mkdir(parents=True, exist_ok=True) + shutil.copyfile(src, dst) + + # Now build params and summary + rp = msparser.ms_mascotresults_params() + rf.get_ms_mascotresults_params(opts, rp) + flags2 = rp.getFlags2() | msparser.ms_peptidesummary.MSPEPSUM_PERCOLATOR + rp.setFlags2(flags2) + rp.setTargetFDR(0.01) + rp.setTargetFDRType(msparser.ms_mascotresults.DS_COUNT_SEQUENCE) + return msparser.ms_peptidesummary(rf, rp), rf +``` + +After this, `rf.getPercolatorFileNames()` returns a 3-tuple +`(pip, target.pop, decoy.pop)` (note pip is FIRST, not last). The QC +dashboard's `stage_percolator_cache` follows exactly this pattern. + +## 21. Read the `.target.pop` directly (FDR-passing PSMs, family-member-safe) + +When you only need the **FDR-passing PSMs** (and per-protein roll-ups built from +them), you can skip the `ms_peptidesummary` staging dance of recipe 20 entirely +and read the server's `.target.pop` table directly. This is the most robust path +for per-protein identification/quant because grouping by the `proteinIds` column +captures **representatives AND family members** automatically — it sidesteps the +`getHit()`-returns-representatives-only blind spot (see SKILL.md → *Reading all +proteins: family members, not just `getHit()`*). + +`.target.pop` is tab-separated with a header row. Key columns: `PSMId` +(`query:N;rank:M` — **N is the msparser query number**), `score`, `q-value`, +`posterior_error_prob`, `peptide` (`X.SEQUENCE.X`), `proteinIds` (space-separated +accessions). msparser is still used — but only for **spectrum/peak access** via +`ms_inputquery`, which is reliable on the SQLite `.msr`. + +```python +import os, glob, re +import msparser + +PSMID_RE = re.compile(r"query:(\d+);rank:(\d+)") + +def find_target_pop(cache_root, res_name): + hits = glob.glob(os.path.join(cache_root, "*", "*", "*", f"{res_name}.*.target.pop")) + return hits[0] if hits else None + +def proteins_from_target_pop(res_path, cache_root, target_fdr=0.01): + """Group FDR-passing PSMs by protein — reps AND family members included.""" + pop = find_target_pop(cache_root, os.path.basename(res_path)) + if not pop: + return None # search wasn't Percolator-rescored + resfile = msparser.ms_mascotresfilebase.createResfile(res_path) # peaks only + + by_protein = {} # acc -> {"psms": int, "seqs": set, "queries": set} + with open(pop, encoding="utf-8") as fh: + next(fh) # header + for line in fh: + f = line.rstrip("\n").split("\t") + if len(f) < 6: + continue + try: + q = float(f[2]) # q-value column + except ValueError: + continue + if q > target_fdr: + continue + m = PSMID_RE.search(f[0]) # PSMId -> query N + if not m: + continue + query = int(m.group(1)) # use with ms_inputquery(resfile, query) + pepfield = f[4] # X.SEQ.X + seq = pepfield.split(".")[1] if pepfield.count(".") >= 2 else pepfield + for acc in f[5].split(): # proteinIds, space-separated + d = by_protein.setdefault(acc, {"psms": 0, "seqs": set(), "queries": set()}) + d["psms"] += 1 + d["seqs"].add(seq) + d["queries"].add(query) + return by_protein +``` + +For TMT/iTRAQ quant, read each passing PSM's spectrum with +`ms_inputquery(resfile, query)`, sum the reporter-ion intensities (nearest peak +within tolerance per channel), and add them to every protein in `proteinIds`. A +complete working extractor is in the project at +`TPP2D/scripts/msr_to_quant_summary_pop.py`. diff --git a/skills/mascot-parser/references/sqlite-msr.md b/skills/mascot-parser/references/sqlite-msr.md new file mode 100644 index 0000000..0d41280 --- /dev/null +++ b/skills/mascot-parser/references/sqlite-msr.md @@ -0,0 +1,188 @@ +# The SQLite-backed `.msr` result file (Mascot Server 3.x) + +Mascot Server **3.1.241.x** (and later) writes its `.msr` result files as +**SQLite 3 databases** ("Mascot SQLite Results"), not the legacy flat-text +`.dat` format. The legacy `.dat` is still embedded inside, gzip-chunked, in a +`legacy__dat28_file` table — but on large DIA results that embedded copy can be +incomplete or expensive to materialise. + +> ## ⚠️ Before you reach for direct SQL: it is usually NOT the SQLite format +> ## that is hiding your proteins +> +> "Direct SQL on the `.msr`" is the right escape hatch only for the specific +> **large-DIA `ms_peptidesummary` hang / 0-hits** failure described below. It is +> **not** the fix for the two far more common "missing proteins" causes, and +> the SQLite `.msr` **cannot** supply what those need: +> +> 1. **Protein-family members are invisible to a `getHit()`-only loop.** +> `ms_peptidesummary.getHit(n)` returns only family *representatives*; +> homologous proteins clustered into a family are stored as *members* and +> are silently dropped by a reps-only walk. This is a logic gap, not a file +> format problem — switching to raw SQL does not fix it by itself, and you'd +> be reimplementing Mascot's clustering. **Walk family members instead** (see +> SKILL.md → *Reading all proteins: family members, not just `getHit()`*). +> +> 2. **Percolator / ML-rescored results are NOT in the `.msr` SQLite at all.** +> Percolator output lives in the server cache files +> `..target.pop` / `.decoy.pop` / `.pip`. No table in the schema +> below contains q-values or PEPs. If a protein "passes at q=0" but you can't +> find it, it is almost certainly a Percolator-only, FDR-passing hit you +> haven't imported — the `.msr` will never show it via SQL **or** msparser +> until you stage those `.pop` files (see SKILL.md → *Importing Percolator +> results (not in the .msr)*). For FDR-passing PSMs + per-protein quant, the +> **`.target.pop`-driven approach** in that section is the robust path and it +> sidesteps the `getHit()` family blind spot for free. +> +> Reach for the direct-SQL recipes here only after you've ruled out (1) and (2). + +## The gotcha — `ms_peptidesummary` on a big SQLite `.msr` + +On a narrow-window DIA result (hundreds of thousands of pseudo-spectra, multi-DB +search), `msparser.ms_peptidesummary(resfile, params)` and `ms_proteinsummary`: + +- **return `getNumberOfHits() == 0`** (no protein hits) — or +- **hang for many minutes** (it decompresses the ~hundreds-of-MB `legacy__dat28_file` + and rebuilds a summary over all queries × all databases) — or both, depending + on the build and the `params`. + +Meanwhile `resfile.isValid()`, `anyMSMS()`, `anyFastaMatches()`, +`hasQuantitation()` all return **True** — so the empty result *looks like* a +silent failure rather than an unsupported format. Direct SQL on the same file +shows the data is fully there (e.g. 700+ PSMs for a single biomarker protein in +one raw file). On a DIA result with Percolator/MS2Rescore FDR this overlaps with +the older `MSPEPSUM_PERCOLATOR`-makes-`getHit()`-return-`None` bug (see the main +skill's Percolator section) — clearing that flag doesn't fix the hang. + +**Practical rule:** for the new SQLite `.msr` on a large DIA result *where the +summary API genuinely hangs or returns 0 hits*, **read the SQLite tables +directly with the `sqlite3` module** instead of going through msparser's summary +API. msparser is still the right tool for `resfile.params()` (search +parameters, including multi-DB enumeration), `anyMSMS()` / `isErrorTolerant()` / +`hasQuantitation()` predicates, and `ms_inputquery` spectrum access. It's +specifically the *protein/peptide summary* on big DIA SQLite results that you +should bypass. + +**But first confirm this is actually your failure mode.** If `getNumberOfHits()` +returns a sensible non-zero count and the call returns promptly, your "missing +proteins" are almost certainly family members or Percolator-only hits, not a +SQLite-format problem — see the caveat box at the top of this file. Direct SQL +here gives you raw `psm__peptides` / `protein__data` rows but **no q-values, no +PEPs, and no family clustering** (those tables simply don't exist), so it will +not, on its own, recover Percolator-passing or family-member proteins. + +Detecting the format: open the file with `sqlite3` and read `schema_version`: + +```python +import sqlite3 +con = sqlite3.connect(f"file:{path}?mode=ro", uri=True) +maj, minr, prod, ver = con.execute( + "SELECT major_version, minor_version, product_name, product_version " + "FROM schema_version").fetchone() +# e.g. (1, 0, 'Mascot Server', '3.1.241.27') +``` + +(If `sqlite3.connect` raises `DatabaseError: file is not a database`, you have a +legacy flat `.dat`/`.msr` — use msparser normally.) + +## Table layout (schema_version 1.0) + +| Table | Rows (typical DIA file) | Purpose / key columns | +|---|---|---| +| `schema_version` | 1 | `major_version`, `minor_version`, `product_name`, `product_version` | +| `search__header` | ~17 | key/value: search title (`COM`), user, etc. | +| `search__parameters` | ~170 | key/value: every Mascot search-form field (`TOL`, `ITOL`, `CLE`, `PFA`, `MODS`, `IT_MODS`, `DECOY`, `INSTRUMENT`, …) | +| `search__databases` | one row per DB | `db_id`, `db_name`, `fasta_file`, `release`, `db_type`, `sequences`, `residues`, `decoy_type`, `et_sequences` — **this is the authoritative multi-DB list** | +| `search__fixed_mods` | one per fixed mod | `mod_num`, `mod_name`, `delta`, `neutral_loss`, `residues` | +| `search__variable_mods` | one per var mod | `mod_num`, `mod_name`, `delta` | +| `search__variable_mod_nl` | one per var-mod NL | `mod_num`, `idx`, `type`, `value` | +| `search__symbol_masses` | ~33 | `symbol` → `mass` (residue + special symbol masses) | +| `search__taxonomy` | 0+ | taxonomy filter rows | +| `search__config_files` | one per config | enzyme/mods/etc config snapshots | +| `query__data` | one per query | `query_id`, `source_index`, `title` (URL-encoded MGF title), `charge`, `rt_in_seconds`, `scans`, `raw_file`, `mass_min`/`mass_max`, `num_vals`, **`ions1` / `ions1_charge`** (the fragment peak list as `"mz:intensity,…"` / `"z,z,…"`; `ions2`/`ions3` for multi-spectrum queries), `ion_mobility`, `it_mods` | +| `query__summary` | ~2× query count | one row per (query, psm_type): `query_id`, `psm_type`, `observed_mr`, `observed_m_z`, `charge`, **`intensity`** (precursor MS1 intensity Distiller assigned during peak detection), `qmatch` (homology-threshold proxy) | +| `psm__peptides` | millions | one row per PSM: `query_id`, `rank_id`, `sequence_idx`, **`psm_type`** (0 = MS/MS-fragment-matched; 1 = secondary precursor-only, max score ~50), `missed_cleavages`, `peptide_mr`, `delta`, **`ions_score`**, **`sequence`**, **`varmods_string`** (per-position digit string, `'0'`=none), `summed_mods_string`/`local_mods_string` (for site-localisation), `ions_matched`, `ion_series_found`, `peaks_used_from_ions{1,2,3}`, `drange_start_pos`/`drange_end_pos`, **`quant_component_name`** (Distiller quant component this PSM was assigned to, when the search came from a Distiller project) | +| `psm__proteins` | ~2× PSM count | PSM → protein map: `query_id`, `rank_id`, `sequence_idx`, `psm_type`, **`protein_id`**, `frame_number`, **`start_idx`** / **`end_idx`** (1-based residue range in the protein), `multiplicity`, **`residue_before`** / **`residue_after`** (flanking residues), `sl_*` (spectral-library provenance) | +| `protein__data` | 10k–200k | **`protein_id`**, **`is_decoy`** (0/1), `protein_seq_id`, **`db_id`** (→ `search__databases`), **`accession_str`**, `mass`, `title` (description), `taxonomy_id`, `has_top_scoring_peptide` | +| `psm__substitutions` | 0+ | error-tolerant substitutions: `query_id`, `rank_id`, `psm_type`, `site`, `ambiguous_residue`, `residue` | +| `psm__et_mods`, `psm__seq_tags`, `psm__sl_mods`, `psm__linked_sites`, `psm__monolinks` | 0+ | ET mods / sequence tags / spectral-library mods / crosslink sites | +| `pmf_hit__*` | 0 in MS/MS searches | PMF (peptide-mass-fingerprint) hit tables | +| `legacy__dat28_file` | tens–hundreds of chunks | `idx`, `data` (BLOB) — the classic flat `.dat` file, gzip-compressed and split into ordered chunks. Reassemble by `ORDER BY idx`, concatenate, `gzip.decompress`. Can be hundreds of MB on a DIA result. This is what `ms_mascotresfilebase.createResfile()` actually parses — and why the summary API chokes. | + +## Working query: PSMs for one protein, one raw file + +```python +import sqlite3 +con = sqlite3.connect(f"file:{msr_path}?mode=ro", uri=True) + +SQL = """ +SELECT pp.query_id, pp.rank_id, pp.psm_type, + pp.sequence, pp.peptide_mr, pp.delta, pp.ions_score, pp.varmods_string, + qs.observed_m_z, qs.charge, qs.intensity, qs.qmatch, + qd.rt_in_seconds, qd.scans, qd.raw_file, qd.title, + pd.accession_str, pd.db_id, pd.is_decoy, + px.start_idx, px.end_idx, px.residue_before, px.residue_after +FROM psm__peptides pp +JOIN psm__proteins px + ON px.query_id = pp.query_id AND px.rank_id = pp.rank_id + AND px.psm_type = pp.psm_type AND px.sequence_idx = pp.sequence_idx +JOIN protein__data pd ON pd.protein_id = px.protein_id +JOIN query__summary qs ON qs.query_id = pp.query_id AND qs.psm_type = pp.psm_type +JOIN query__data qd ON qd.query_id = pp.query_id +WHERE pd.accession_str LIKE ? -- e.g. 'P02750%' + AND pd.is_decoy = 0 + AND pp.psm_type = 0 -- MS/MS-fragment-matched only + AND pp.ions_score >= ? -- e.g. 20.0 +ORDER BY pp.ions_score DESC +""" +for row in con.execute(SQL, ("P02750%", 20.0)): + ... +``` + +Notes: +- **Always filter `psm_type = 0`** unless you specifically want the secondary precursor-only matches. +- **`is_decoy = 0`** on `protein__data` excludes decoy hits. +- There is no `expectation_value` column — score-vs-`qmatch` is the available significance proxy. If you need true E-values, decode the embedded legacy `.dat` (`legacy__dat28_file`) and feed it to msparser — but on a big DIA file that's the slow path you're trying to avoid. +- `rt_in_seconds`, `scans`, `charge` in `query__data` are stored as **strings** (occasionally comma-lists for multi-scan queries) — cast/split as needed. +- `query__data.title` is **URL-encoded** — `urllib.parse.unquote` it to get the original MGF spectrum title (carries the raw-file path and scan number). +- `varmods_string` is one digit per peptide residue position (plus N-/C-term slots); a digit `n` means variable mod `mod_num = n` (from `search__variable_mods`) is on that residue. All-zeros = unmodified. +- The `quant_component_name` column links a PSM to a Mascot Distiller quant component, but the **integrated XIC areas themselves are not in the `.msr`** — they live in the Distiller `.rov` project's binary segment streams (see the `mascot-distiller` skill's `ROV_FILE_FORMAT.md`). For a quick label-free proxy you can trapezoidally integrate `query__summary.intensity` vs `query__data.rt_in_seconds` over the matched PSMs of a peptide form — symmetric across forms, so ratios are robust, but it underestimates the true Distiller XIC. + +## Per-PSM protein accessions: msparser methods return EMPTY on `.msr` + +Verified on Mascot 3.1.242 / msparser 3.x (2026-08): with a working +`ms_peptidesummary` over an `.msr`, the accession routes that work on `.dat` +**silently return nothing** for proteins outside the report hit list — +`summary.getProteinsWithThisPepMatch(q, r)` → empty string, +`pep.getProteins()` → empty tuple, and section access +(`getSectionValueStr(SEC_PEPTIDES, "qN_pM")`) does not exist for `.msr` at +all (`ms_mascotresfilebase` has no SEC_ constants; they live on +`ms_mascotresfile_dat`). Code ported from `.dat` workflows (e.g. +fdr_stats.pl) will run without error and classify zero matches. + +**The working hybrid pattern** — msparser for scores/expect/pretty-rank, +SQLite for accessions, joined on `(query_id, rank_id)`: + +```python +resfile = msparser.ms_mascotresfilebase.createResfile(path) +summary = msparser.ms_peptidesummary(resfile, flags, 0.05, 1, "", 0, 5, "", + msparser.ms_peptidesummary.MSPEPSUM_NONE) +expect = summary.getPeptideExpectationValue(pep.getIonsScore(), q) + +acc = {} # (query_id, rank_id) -> [accessions] +for q, r, a in sqlite3.connect(path).execute( + "SELECT p.query_id, p.rank_id, d.accession_str " + "FROM psm__proteins p JOIN protein__data d " + "ON p.protein_id = d.protein_id AND p.psm_type = 0"): + acc.setdefault((q, r), []).append(a) +``` + +(Used by `fdr_stats_compare.py` in the mascot-entrapment-tools repo; see the +`mascot-entrapment-audit` skill.) + +## When msparser *does* still work on a SQLite `.msr` + +- `resfile = ms_mascotresfilebase.createResfile(path)` — opens fine (`isValid()` True). +- `resfile.params()` — full search parameters, including `getNumberOfDatabases()` + `getDB(i)` for the multi-DB list. Use this, not the `search__databases` table, if you want the msparser-typed view. +- `resfile.anyMSMS()` / `anyPMF()` / `isErrorTolerant()` / `hasQuantitation()` — all reliable. +- `ms_inputquery(resfile, q)` — spectrum/peak access works (it can be slow to first-touch because of the legacy `.dat` materialisation). +- `ms_peptidesummary` / `ms_proteinsummary` on a **small** search (a handful of files, modest query count) — fine. It's specifically the large-DIA case that breaks.