corems.mass_spectra.output.mgf

  1"""SIRIUS-compatible MGF export for LC-MS mass features."""
  2
  3import warnings
  4from pathlib import Path
  5
  6import numpy as np
  7
  8
  9def sirius_charge(polarity: str) -> str:
 10    """Map LCMS polarity to a SIRIUS CHARGE token.
 11
 12    Parameters
 13    ----------
 14    polarity : str
 15        ``'positive'`` or ``'negative'``.
 16
 17    Returns
 18    -------
 19    str
 20        ``'1+'`` or ``'1-'``.
 21    """
 22    value = (polarity or "").strip().lower()
 23    if value == "positive":
 24        return "1+"
 25    if value == "negative":
 26        return "1-"
 27    raise ValueError(
 28        "Polarity not set for dataset, must be a either 'positive' or 'negative'"
 29    )
 30
 31
 32def format_ion_block(headers: dict, mz, intensity) -> str:
 33    """Serialize one MGF spectrum as ``BEGIN IONS`` ... ``END IONS``.
 34
 35    Parameters
 36    ----------
 37    headers : dict
 38        MGF ``NAME=VALUE`` pairs, written in insertion order.
 39    mz, intensity : array_like
 40        Peak lists.
 41
 42    Returns
 43    -------
 44    str
 45        One ion block including a trailing newline.
 46    """
 47    lines = ["BEGIN IONS"]
 48    for key, value in headers.items():
 49        if value is None:
 50            continue
 51        lines.append(f"{key}={value}")
 52    mz = np.asarray(mz, dtype=float)
 53    intensity = np.asarray(intensity, dtype=float)
 54    for mass, abund in zip(mz, intensity):
 55        lines.append(f"{mass} {abund}")
 56    lines.append("END IONS")
 57    return "\n".join(lines) + "\n"
 58
 59
 60def ensure_mgf_path(out_file_path) -> Path:
 61    """Return a Path, appending ``.mgf`` when no suffix is given."""
 62    path = Path(out_file_path)
 63    if path.suffix == "":
 64        path = path.with_suffix(".mgf")
 65    return path
 66
 67
 68def peak_mz_abundance(spectrum):
 69    """Return (mz, abundance) from processed peaks, else raw arrays.
 70
 71    A spectrum that has been through ``process_mass_spec`` is read from
 72    ``mspeaks`` even when that list is empty. Raw ``_mz_exp`` /
 73    ``_abundance`` are used only when the spectrum has not been processed.
 74
 75    Parameters
 76    ----------
 77    spectrum : MassSpectrum or None
 78        Spectrum to read.
 79
 80    Returns
 81    -------
 82    tuple of numpy.ndarray
 83        ``(mz, abundance)``. Empty arrays if ``spectrum`` is None, the
 84        processed peak list is empty, or raw arrays are missing.
 85    """
 86    if spectrum is None:
 87        return np.array([]), np.array([])
 88    # process_mass_spec sets _dynamic_range, including 0 when no peaks remain.
 89    processed = getattr(spectrum, "_dynamic_range", None) is not None
 90    mspeaks = getattr(spectrum, "mspeaks", None)
 91    if processed or mspeaks:
 92        if not mspeaks:
 93            return np.array([]), np.array([])
 94        return np.asarray(spectrum.mz_exp, dtype=float), np.asarray(
 95            spectrum.abundance, dtype=float
 96        )
 97    mz = getattr(spectrum, "_mz_exp", None)
 98    ab = getattr(spectrum, "_abundance", None)
 99    if mz is None or ab is None:
100        return np.array([]), np.array([])
101    return np.asarray(mz, dtype=float), np.asarray(ab, dtype=float)
102
103
104def _ms1_precursor_peak(feature):
105    """Return the feature precursor as a one-peak MS1 list.
106
107    Raises
108    ------
109    ValueError
110        If the feature has no precursor m/z.
111    """
112    mz = getattr(feature, "mz", None)
113    if mz is None:
114        raise ValueError("Mass feature is missing precursor m/z; cannot export MGF")
115    intensity = getattr(feature, "intensity", None)
116    if intensity is None:
117        intensity = 0.0
118    return np.array([float(mz)]), np.array([float(intensity)])
119
120
121def _scan_of_spectrum(feature, spec):
122    scan = getattr(spec, "scan_number", None)
123    if scan is not None:
124        return scan
125    for key, value in getattr(feature, "ms2_mass_spectra", {}).items():
126        if value is spec:
127            return key
128    return None
129
130
131def _ms2_spectra(feature, ms2_mode: str):
132    """Select MS2 spectra for one feature.
133
134    ``'best'`` returns ``feature.best_ms2``: the scan with the highest
135    library entropy similarity when ``ms2_similarity_results`` is
136    non-empty, otherwise the non-chimeric scan closest to the apex.
137    ``'all'`` returns every associated scan that has peaks.
138    """
139    if ms2_mode not in ("best", "all"):
140        raise ValueError(f"ms2_mode must be 'best' or 'all', got {ms2_mode!r}")
141    if ms2_mode == "best":
142        spec = feature.best_ms2
143        if spec is None:
144            return []
145        return [(_scan_of_spectrum(feature, spec), spec)]
146    spectra = []
147    for scan, spec in feature.ms2_mass_spectra.items():
148        mz, _ = peak_mz_abundance(spec)
149        if mz.size:
150            spectra.append((scan, spec))
151    return spectra
152
153
154def _ordered_headers(headers: dict) -> dict:
155    order = [
156        "FEATURE_ID",
157        "PEPMASS",
158        "CHARGE",
159        "MSLEVEL",
160        "TITLE",
161        "RTINSECONDS",
162        "SCANS",
163    ]
164    out = {}
165    for key in order:
166        if key in headers:
167            out[key] = headers[key]
168    for key, value in headers.items():
169        if key not in out:
170            out[key] = value
171    return out
172
173
174def _shared_headers(feature, feature_id, polarity, sample_name=None, title=None):
175    if title is None:
176        if sample_name:
177            title = f"{sample_name} feature {feature_id}"
178        else:
179            title = f"feature {feature_id}"
180    rt = getattr(feature, "retention_time", None)
181    rtinseconds = None if rt is None else float(rt) * 60.0
182    return {
183        "FEATURE_ID": feature_id,
184        "PEPMASS": feature.mz,
185        "CHARGE": sirius_charge(polarity),
186        "TITLE": title,
187        "RTINSECONDS": rtinseconds,
188    }
189
190
191def feature_mgf_text(
192    feature,
193    feature_id,
194    polarity,
195    *,
196    sample_name=None,
197    ms2_mode="best",
198    title=None,
199):
200    """Return MGF text for one feature, or None if MS2 is missing.
201
202    Parameters
203    ----------
204    feature : LCMSMassFeature
205        Mass feature with associated MS2 spectra. MS1 is the precursor
206        ``mz`` / ``intensity`` on the feature, not a scanned spectrum.
207    feature_id : hashable
208        Value written as ``FEATURE_ID`` (mass-feature id or consensus cluster).
209    polarity : str
210        ``'positive'`` or ``'negative'``.
211    sample_name : str, optional
212        Used in the default TITLE.
213    ms2_mode : {'best', 'all'}, optional
214        Which associated MS2 spectra to write. ``'best'`` (default)
215        writes ``LCMSMassFeature.best_ms2``. After a spectral search
216        that is the scan with the highest library entropy similarity;
217        with no search results it is the non-chimeric scan closest to
218        the apex. ``'all'`` writes every MS2 scan that has peaks.
219    title : str, optional
220        Override TITLE.
221
222    Returns
223    -------
224    str or None
225        One or more ion blocks, or None if the feature is incomplete.
226    """
227    mz1, ab1 = _ms1_precursor_peak(feature)
228    ms2_list = _ms2_spectra(feature, ms2_mode)
229    if not ms2_list:
230        return None
231
232    shared = _shared_headers(
233        feature, feature_id, polarity, sample_name=sample_name, title=title
234    )
235    blocks = []
236    ms1_headers = dict(shared)
237    ms1_headers["MSLEVEL"] = 1
238    ms1_headers["SCANS"] = getattr(feature, "apex_scan", None)
239    blocks.append(format_ion_block(_ordered_headers(ms1_headers), mz1, ab1))
240    for scan, spec in ms2_list:
241        mz2, ab2 = peak_mz_abundance(spec)
242        if mz2.size == 0:
243            continue
244        ms2_headers = dict(shared)
245        ms2_headers["MSLEVEL"] = 2
246        ms2_headers["SCANS"] = (
247            scan if scan is not None else getattr(spec, "scan_number", None)
248        )
249        blocks.append(format_ion_block(_ordered_headers(ms2_headers), mz2, ab2))
250    if len(blocks) < 2:
251        return None
252    return "".join(blocks)
253
254
255def write_feature_records_to_mgf(
256    records,
257    polarity,
258    out_file_path,
259    *,
260    ms2_mode="best",
261    overwrite=False,
262) -> Path:
263    """Write selected feature records to an MGF file.
264
265    Parameters
266    ----------
267    records : iterable of tuple
268        ``(feature_id, feature, sample_name)`` triples. ``sample_name`` may be None.
269    polarity : str
270        ``'positive'`` or ``'negative'``.
271    out_file_path : str or Path
272        Output path. ``.mgf`` is appended when no suffix is given.
273    ms2_mode : {'best', 'all'}, optional
274        Which associated MS2 spectra to write. ``'best'`` (default)
275        writes ``LCMSMassFeature.best_ms2``. After a spectral search
276        that is the scan with the highest library entropy similarity;
277        with no search results it is the non-chimeric scan closest to
278        the apex. ``'all'`` writes every MS2 scan that has peaks.
279    overwrite : bool, optional
280        Replace an existing file. Default is False.
281
282    Returns
283    -------
284    pathlib.Path
285        Path of the written MGF file.
286
287    Raises
288    ------
289    FileExistsError
290        If the path exists and ``overwrite`` is False.
291    ValueError
292        If no complete MS1/MS2 pairs remain after skipping.
293    """
294    path = ensure_mgf_path(out_file_path)
295    if path.exists() and not overwrite:
296        raise FileExistsError(f"MGF file already exists: {path}")
297
298    chunks = []
299    n_in = 0
300    n_skip = 0
301    for feature_id, feature, sample_name in records:
302        n_in += 1
303        text = feature_mgf_text(
304            feature,
305            feature_id,
306            polarity,
307            sample_name=sample_name,
308            ms2_mode=ms2_mode,
309        )
310        if text is None:
311            n_skip += 1
312            continue
313        chunks.append(text)
314
315    if not chunks:
316        raise ValueError(
317            "No complete MS1/MS2 feature pairs to export "
318            f"(skipped {n_skip} of {n_in} selected features)"
319        )
320    if n_skip:
321        warnings.warn(
322            f"Skipped {n_skip} of {n_in} features missing a usable MS2 spectrum",
323            UserWarning,
324        )
325    path.write_text("".join(chunks), encoding="utf-8")
326    return path
327
328
329def iter_collection_mgf_records(collection, cluster_ids=None):
330    """Build ``(cluster_id, mass_feature, sample_name)`` from consensus representatives.
331
332    Parameters
333    ----------
334    collection : LCMSCollection
335        Collection with consensus clustering and loaded representatives.
336    cluster_ids : iterable, optional
337        Cluster ids to export. Default is all representative clusters.
338
339    Returns
340    -------
341    list of tuple
342        ``(cluster_id, LCMSMassFeature, sample_name)``.
343
344    Raises
345    ------
346    ValueError
347        If representatives are missing, requested clusters are unknown, or
348        representative objects are not loaded on the samples.
349    """
350    if not hasattr(collection, "get_representative_mass_features_for_all_clusters"):
351        raise ValueError("Collection cannot provide consensus representatives")
352    reps = collection.get_representative_mass_features_for_all_clusters()
353    if reps is None or len(reps) == 0:
354        raise ValueError(
355            "No consensus representatives found. Run process_consensus_features() "
356            "with load_representatives=True, add_ms2=True first."
357        )
358    if cluster_ids is None:
359        selected = reps
360    else:
361        requested = list(cluster_ids)
362        selected = reps[reps["cluster"].isin(requested)]
363        found = set(selected["cluster"].tolist())
364        missing = [c for c in requested if c not in found]
365        if missing:
366            raise ValueError(f"Consensus cluster id(s) not found: {missing}")
367
368    records = []
369    lcms_map = getattr(collection, "_lcms", {})
370    for _, row in selected.iterrows():
371        cluster = row["cluster"]
372        sample_name = None
373        if "sample_name" in selected.columns:
374            val = row["sample_name"]
375            if not (isinstance(val, (float, np.floating)) and np.isnan(val)):
376                sample_name = val
377        if sample_name is None:
378            if "sample_id" not in selected.columns:
379                raise ValueError(
380                    "Cannot resolve sample_name from consensus representatives"
381                )
382            sample_name = collection.samples[int(row["sample_id"])]
383        mf_id = row["mf_id"]
384        sample = lcms_map.get(sample_name)
385        mass_features = (
386            getattr(sample, "mass_features", {}) if sample is not None else {}
387        )
388        if mf_id not in mass_features:
389            try:
390                mf_id = int(mf_id)
391            except (TypeError, ValueError):
392                mf_id = row["mf_id"]
393        if sample is None or mf_id not in mass_features:
394            raise ValueError(
395                "Representative mass features have not been loaded into samples. "
396                "Call process_consensus_features() with load_representatives=True, "
397                "add_ms2=True before exporting MGF."
398            )
399        records.append((cluster, mass_features[mf_id], sample_name))
400    return records