From 216e630e450086c63d1f883ffb4e55318bd98d3d Mon Sep 17 00:00:00 2001 From: Yasset Perez-Riverol Date: Wed, 7 Oct 2026 16:24:20 +0100 Subject: [PATCH 1/2] Read large Olink CSV exports faster PAD000061's 1.7 GB Explore HT CSV took 284 s to read; now about 43 s, with identical expression, sample, feature, LOD and QC-flag tables. - the delimiter is taken from the header line and the file is parsed with pandas' C engine instead of sniffing with the Python engine; files the C engine cannot parse fall back to the Python engine - QC flags are normalised once per distinct value and placed into the sample x assay matrix by position instead of pivot_table; the first non-missing measurement still wins Peak memory on that file rises from 3.7 to 4.7 GB (whole-column parsing). Checked identical against the previous reader on all 15 local Olink files. --- src/pyprideap/io/readers/olink_csv.py | 62 +++++++++++++++++++++++---- tests/test_qc_pipeline.py | 30 +++++++++++++ 2 files changed, 84 insertions(+), 8 deletions(-) diff --git a/src/pyprideap/io/readers/olink_csv.py b/src/pyprideap/io/readers/olink_csv.py index 23dff18..e53d7f6 100644 --- a/src/pyprideap/io/readers/olink_csv.py +++ b/src/pyprideap/io/readers/olink_csv.py @@ -4,6 +4,7 @@ import warnings from pathlib import Path +import numpy as np import pandas as pd from pyprideap.core import AffinityDataset, Platform @@ -156,23 +157,68 @@ def _qc_flag_matrices(df: pd.DataFrame, sample_key: str, sample_order: object, a columns present; when two columns map to the same key the first listed wins. """ matrices: dict[str, object] = {} - for col, key in _QC_FLAG_COLUMNS.items(): - if col not in df.columns or key in matrices: + present = [(col, key) for col, key in _QC_FLAG_COLUMNS.items() if col in df.columns] + if not present: + return matrices + # Row and column position of every measurement, computed once. Placing values by + # position is much faster than pivot_table on multi-million-row Explore HT exports. + rows = pd.Index(pd.Series(sample_order)).get_indexer(df[sample_key]) + cols = pd.Index(assays).get_indexer(df["OlinkID"]) + placed = (rows >= 0) & (cols >= 0) + for col, key in present: + if key in matrices: continue - flags = df[[sample_key, "OlinkID", col]].copy() - # Olink writes WARN in current exports; accept the long form too, as filter_qc does - flags[col] = flags[col].astype("string").str.strip().str.upper().replace({"WARNING": "WARN"}) - matrix = flags.pivot_table(index=sample_key, columns="OlinkID", values=col, aggfunc="first") - matrices[key] = matrix.reindex(index=sample_order, columns=assays).reset_index(drop=True) + # Flags take a handful of distinct values: normalise each once. Olink writes + # WARN in current exports; accept the long form too, as filter_qc does. + codes, uniques = pd.factorize(df[col]) + labels = pd.Series(uniques).astype("string").str.strip().str.upper().replace({"WARNING": "WARN"}) + values = np.where(codes >= 0, labels.reindex(codes).to_numpy(dtype=object), None) + ok = placed & (codes >= 0) + matrix = np.full((len(pd.Series(sample_order)), len(assays)), None, dtype=object) + # First non-missing measurement wins (as pivot_table's "first"); dropping + # later repeats leaves each position written exactly once + r, c, v = rows[ok], cols[ok], values[ok] + first = ~pd.Series(r * len(assays) + c).duplicated(keep="first").to_numpy() + matrix[r[first], c[first]] = v[first] + matrices[key] = pd.DataFrame(matrix, columns=assays).astype("string") return matrices +_DELIMITERS = (",", ";", "\t", "|") + + +def _sniff_delimiter(path: Path) -> str | None: + """Delimiter of a delimited text file, from its header line (None if unclear).""" + with open(path, encoding="utf-8-sig", errors="replace") as f: + header = f.readline() + counts = {d: header.count(d) for d in _DELIMITERS} + best = max(counts, key=lambda d: counts[d]) + return best if counts[best] > 0 else None + + +def _read_delimited(path: Path) -> pd.DataFrame: + """Read a delimited Olink export. + + The delimiter is taken from the header line and the file is parsed with + pandas' C engine, which is about 10x faster than delimiter sniffing with the + Python engine on large exports (e.g. a 1.7 GB Explore HT CSV). Files the C + engine cannot parse (e.g. ragged rows) fall back to the Python engine. + """ + sep = _sniff_delimiter(path) + if sep is not None: + try: + return pd.read_csv(path, sep=sep, low_memory=False) + except (pd.errors.ParserError, UnicodeDecodeError) as exc: + logger.debug("C engine could not parse %s (%s); using the Python engine", path.name, exc) + return pd.read_csv(path, sep=None, engine="python") + + def read_olink_csv(path: str | Path) -> AffinityDataset: path = Path(path) if not path.exists(): raise FileNotFoundError(f"File not found: {path}") - df = _apply_olink_aliases(pd.read_csv(path, sep=None, engine="python")) + df = _apply_olink_aliases(_read_delimited(path)) missing = _REQUIRED_COLS - set(df.columns) if missing: raise ValueError(f"Missing required columns in {path.name}: {sorted(missing)}") diff --git a/tests/test_qc_pipeline.py b/tests/test_qc_pipeline.py index d61450d..90fc572 100644 --- a/tests/test_qc_pipeline.py +++ b/tests/test_qc_pipeline.py @@ -1234,3 +1234,33 @@ def test_merge_picks_the_sample_column_that_matches(self): {"source name": ["PAD000014-5530_P", "PAD000014-5531_P", "PAD000014-5532_P"], "sex": list("fmf")} ) assert merge_sdrf(ds, sdrf).samples["sex"].tolist() == ["f", "m", "f"] + + def test_csv_delimiter_detection(self, tmp_path): + from pyprideap.io.readers.olink_csv import _read_delimited, _sniff_delimiter + + for sep in (",", ";", "\t"): + path = tmp_path / "x.npx.csv" + path.write_text(sep.join(["SampleID", "OlinkID", "NPX"]) + "\n" + sep.join(["S1", "O1", "1.5"]) + "\n") + assert _sniff_delimiter(path) == sep + assert _read_delimited(path).columns.tolist() == ["SampleID", "OlinkID", "NPX"] + + def test_qc_flags_first_value_wins_and_missing_stay_empty(self, tmp_path): + from pyprideap.io.readers.olink_csv import read_olink_csv + + rows = [ + ("S1", "O1", "Pass"), + ("S1", "O1", "FAIL"), # duplicate measurement: first wins + ("S1", "O2", "Warning"), + ("S2", "O1", None), # long form normalised; missing flag stays empty + ] + df = pd.DataFrame( + [ + {"SampleID": s, "OlinkID": o, "UniProt": "P", "Assay": "A", "Panel": "X", "NPX": 1.0, "QC_Warning": f} + for s, o, f in rows + ] + ) + path = tmp_path / "flags.npx.csv" + df.to_csv(path, index=False) + m = read_olink_csv(path).metadata["sample_qc_matrix"] + assert m.loc[0, "O1"] == "PASS" and m.loc[0, "O2"] == "WARN" + assert pd.isna(m.loc[1, "O1"]) and pd.isna(m.loc[1, "O2"]) From 4d2b6226dee7fc30c4ef07f3e14b38b93ca869e4 Mon Sep 17 00:00:00 2001 From: Yasset Perez-Riverol Date: Wed, 7 Oct 2026 17:11:37 +0100 Subject: [PATCH 2/2] Olink CSV reader: pyarrow, used columns only, position-based matrices PAD000061 (1.7 GB Explore HT CSV, 6.1M rows): read in ~6.7 s with a 2.5 GB peak, against 294 s and 4.2 GB for the previous reader. - parse with pyarrow, reading only the columns the reader uses (14 of 30 here) and storing identifiers and flags as categories; files of 256 MB or more are streamed in blocks to bound peak memory; files pyarrow cannot parse fall back to pandas - NPX and LOD matrices are built by position instead of pivot_table, with the same semantics (first non-missing value, assays sorted, assays without values dropped); QC flags keep integer codes until the final table - sample IDs are always text: numeric-looking IDs keep leading zeros - NPX is parsed exactly: pandas' default fast float parser was off in the last binary digit for 87,407 of 5.9M NPX values in that file (|diff| <= 9e-16); pyarrow matches Python float() for all of them Identical output to the previous reader on all 15 local Olink files. --- src/pyprideap/io/readers/olink_csv.py | 202 ++++++++++++++++++++++---- tests/test_qc_pipeline.py | 51 +++++++ 2 files changed, 221 insertions(+), 32 deletions(-) diff --git a/src/pyprideap/io/readers/olink_csv.py b/src/pyprideap/io/readers/olink_csv.py index e53d7f6..4eb8d3a 100644 --- a/src/pyprideap/io/readers/olink_csv.py +++ b/src/pyprideap/io/readers/olink_csv.py @@ -3,6 +3,7 @@ import logging import warnings from pathlib import Path +from typing import cast import numpy as np import pandas as pd @@ -150,6 +151,57 @@ def _detect_olink_platform(olink_ids: pd.Series) -> Platform: } +def _positions(values: pd.Series, labels: pd.Index) -> np.ndarray: + """Position of each value in *labels* (-1 when absent or missing). + + Categorical columns are mapped through their categories, so millions of rows + are resolved with one lookup per distinct value. Values are compared as text + when the labels are text (identifiers read as categories are text). + """ + as_text = labels.inferred_type in ("string", "empty") + if isinstance(values.dtype, pd.CategoricalDtype): + categories = values.cat.categories + lookup = labels.get_indexer(categories.astype(str) if as_text else categories) + codes = values.cat.codes.to_numpy() + return np.where(codes >= 0, lookup[codes], -1) + if as_text: + values = values.map(lambda v: v if pd.isna(v) else str(v)) + return np.asarray(labels.get_indexer(pd.Index(values))) + + +def _first_per_cell(rows: np.ndarray, cols: np.ndarray, n_cols: int) -> np.ndarray: + """Mask keeping the first measurement of each (row, column) cell.""" + duplicated = pd.Series(rows.astype(np.int64) * n_cols + cols).duplicated(keep="first").to_numpy(dtype=bool) + return cast(np.ndarray, np.logical_not(duplicated)) + + +def _long_to_wide(df: pd.DataFrame, sample_key: str, sample_order: pd.Index, value_col: str) -> pd.DataFrame: + """Samples x assays matrix of *value_col*, equivalent to + ``pivot_table(index=sample_key, columns="OlinkID", aggfunc="first")`` reindexed to + *sample_order*: first non-missing value per cell, assays sorted, assays without + any value dropped. Built by position, which is much faster and lighter than + pivot_table on multi-million-row exports. + """ + values = pd.to_numeric(df[value_col], errors="coerce").to_numpy(dtype=float) + assay_ids = df["OlinkID"] + labels = pd.Index( + sorted(assay_ids.cat.categories.astype(str)) + if isinstance(assay_ids.dtype, pd.CategoricalDtype) + else sorted({str(v) for v in assay_ids.dropna().unique()}) + ) + rows = _positions(df[sample_key], sample_order) + cols = _positions(assay_ids, labels) + ok = (rows >= 0) & (cols >= 0) & ~np.isnan(values) + r, c, v = rows[ok], cols[ok], values[ok] + first = _first_per_cell(r, c, len(labels)) + matrix = np.full((len(sample_order), len(labels)), np.nan) + matrix[r[first], c[first]] = v[first] + keep = np.zeros(len(labels), dtype=bool) + keep[np.unique(c)] = True + wide = pd.DataFrame(matrix[:, keep], columns=pd.Index(labels[keep], name="OlinkID")) + return cast(pd.DataFrame, wide) + + def _qc_flag_matrices(df: pd.DataFrame, sample_key: str, sample_order: object, assays: pd.Index) -> dict[str, object]: """Pivot per-measurement Olink QC flags into sample x assay matrices of upper-case strings. @@ -162,8 +214,8 @@ def _qc_flag_matrices(df: pd.DataFrame, sample_key: str, sample_order: object, a return matrices # Row and column position of every measurement, computed once. Placing values by # position is much faster than pivot_table on multi-million-row Explore HT exports. - rows = pd.Index(pd.Series(sample_order)).get_indexer(df[sample_key]) - cols = pd.Index(assays).get_indexer(df["OlinkID"]) + rows = _positions(df[sample_key], pd.Index(list(sample_order))) # type: ignore[call-overload] + cols = _positions(df["OlinkID"], pd.Index(assays)) placed = (rows >= 0) & (cols >= 0) for col, key in present: if key in matrices: @@ -172,15 +224,16 @@ def _qc_flag_matrices(df: pd.DataFrame, sample_key: str, sample_order: object, a # WARN in current exports; accept the long form too, as filter_qc does. codes, uniques = pd.factorize(df[col]) labels = pd.Series(uniques).astype("string").str.strip().str.upper().replace({"WARNING": "WARN"}) - values = np.where(codes >= 0, labels.reindex(codes).to_numpy(dtype=object), None) ok = placed & (codes >= 0) - matrix = np.full((len(pd.Series(sample_order)), len(assays)), None, dtype=object) + # Work with the small integer codes; text is looked up once per cell at the end + coded = np.full((len(pd.Series(sample_order)), len(assays)), -1, dtype=np.int32) # First non-missing measurement wins (as pivot_table's "first"); dropping # later repeats leaves each position written exactly once - r, c, v = rows[ok], cols[ok], values[ok] - first = ~pd.Series(r * len(assays) + c).duplicated(keep="first").to_numpy() - matrix[r[first], c[first]] = v[first] - matrices[key] = pd.DataFrame(matrix, columns=assays).astype("string") + r, c, v = rows[ok], cols[ok], codes[ok] + first = _first_per_cell(r, c, len(assays)) + coded[r[first], c[first]] = v[first] + lookup = np.array([*labels.astype(object).where(labels.notna(), None), None], dtype=object) + matrices[key] = pd.DataFrame(lookup[coded], columns=assays).astype("string") return matrices @@ -196,21 +249,116 @@ def _sniff_delimiter(path: Path) -> str | None: return best if counts[best] > 0 else None +# Columns the reader uses; everything else in an export is skipped while parsing +# Per-measurement numbers, parsed as float. MissingFreq (per assay, sometimes +# written as "30%") is read as text and made numeric afterwards when it can be. +_NUMERIC_COLS = {"NPX", "LOD", "LODNPX", "PlateLOD"} +_USED_COLS = ( + _REQUIRED_COLS | _SAMPLE_COLS | _FEATURE_COLS | _NUMERIC_COLS | set(_QC_FLAG_COLUMNS) | set(_OLINK_COLUMN_ALIASES) +) +# pandas' default missing-value markers, so both parsers treat the same text as missing +_NA_VALUES = [ + "", + "#N/A", + "#N/A N/A", + "#NA", + "-1.#IND", + "-1.#QNAN", + "-NaN", + "-nan", + "1.#IND", + "1.#QNAN", + "", + "N/A", + "NA", + "NULL", + "NaN", + "None", + "n/a", + "nan", + "null", +] + + +# Files at least this large are streamed to bound peak memory +_STREAM_MIN_BYTES = 256 * 1024 * 1024 + + +def _read_header(path: Path, sep: str) -> list[str]: + import csv + + with open(path, encoding="utf-8-sig", errors="replace", newline="") as f: + return next(csv.reader(f, delimiter=sep), []) + + +def _read_arrow(path: Path, sep: str, columns: list[str]) -> pd.DataFrame: + """Read the file with pyarrow: only *columns*, text stored as categories. + + Repeated identifiers (sample, assay, plate, flags) are stored once per + distinct value. Large files are streamed, keeping only a few blocks of raw + text in memory at a time. + """ + import pyarrow as pa + import pyarrow.csv as pacsv + + text = pa.dictionary(pa.int32(), pa.string()) + convert = pacsv.ConvertOptions( + include_columns=columns, + column_types={c: (pa.float64() if c in _NUMERIC_COLS else text) for c in columns}, + null_values=_NA_VALUES, + strings_can_be_null=True, + ) + parse = pacsv.ParseOptions(delimiter=sep) + if path.stat().st_size < _STREAM_MIN_BYTES: + # Small files: one multithreaded read (streaming has ~1 s fixed overhead) + table = pacsv.read_csv(path, parse_options=parse, convert_options=convert) + else: + read = pacsv.ReadOptions(block_size=1 << 24) + with pacsv.open_csv(path, read_options=read, parse_options=parse, convert_options=convert) as reader: + table = pa.Table.from_batches(list(reader), schema=reader.schema) + return cast(pd.DataFrame, table.unify_dictionaries().to_pandas()) + + def _read_delimited(path: Path) -> pd.DataFrame: - """Read a delimited Olink export. + """Read a delimited Olink export, keeping only the columns the reader uses. - The delimiter is taken from the header line and the file is parsed with - pandas' C engine, which is about 10x faster than delimiter sniffing with the - Python engine on large exports (e.g. a 1.7 GB Explore HT CSV). Files the C - engine cannot parse (e.g. ragged rows) fall back to the Python engine. + The delimiter comes from the header line. The file is streamed with pyarrow + with identifiers stored as categories: on a 1.7 GB Explore HT export this + parses in ~5 s with a ~1.5 GB peak, against ~27 s and ~3.9 GB for pandas. + Files pyarrow cannot parse (e.g. ragged rows, text in a numeric column) fall + back to pandas. """ sep = _sniff_delimiter(path) if sep is not None: + header = _read_header(path, sep) + columns = [c for c in header if c in _USED_COLS] + if len(columns) == len(set(columns)) and _REQUIRED_COLS <= {_OLINK_COLUMN_ALIASES.get(c, c) for c in columns}: + try: + return _read_arrow(path, sep, columns) + except Exception as exc: # pyarrow raises several error types for malformed files + logger.debug("pyarrow could not parse %s (%s); using pandas", path.name, exc) try: - return pd.read_csv(path, sep=sep, low_memory=False) + return pd.read_csv(path, sep=sep, low_memory=False, dtype={"SampleID": str}) except (pd.errors.ParserError, UnicodeDecodeError) as exc: logger.debug("C engine could not parse %s (%s); using the Python engine", path.name, exc) - return pd.read_csv(path, sep=None, engine="python") + return pd.read_csv(path, sep=None, engine="python", dtype={"SampleID": str}) + + +def _plain(frame: pd.DataFrame) -> pd.DataFrame: + """Small sample/feature tables: categorical columns back to plain values. + + MissingFreq becomes numeric when all its values are numbers (as pandas would + infer); identifier columns stay text so leading zeros are kept. + """ + for col in frame.columns: + if isinstance(frame[col].dtype, pd.CategoricalDtype): + values = frame[col].astype(object).where(frame[col].notna(), np.nan) + if col == "MissingFreq": + numeric = pd.to_numeric(values, errors="coerce") + if not (numeric.isna() & values.notna()).any(): + values = numeric + frame[col] = values + return frame def read_olink_csv(path: str | Path) -> AffinityDataset: @@ -234,22 +382,17 @@ def read_olink_csv(path: str | Path) -> AffinityDataset: df[col] = pd.to_numeric(df[col], errors="coerce") sample_cols = [c for c in df.columns if c in _SAMPLE_COLS] - samples = df[sample_cols].drop_duplicates(subset=[sample_key]).reset_index(drop=True) + samples = _plain(df[sample_cols].drop_duplicates(subset=[sample_key]).reset_index(drop=True)) feature_cols = [c for c in df.columns if c in _FEATURE_COLS] # Drop LOD from per-assay features since it varies per plate/sample feature_cols_no_lod = [c for c in feature_cols if c != "LOD"] - features = df[feature_cols_no_lod].drop_duplicates(subset=["OlinkID"]).reset_index(drop=True) + features = _plain(df[feature_cols_no_lod].drop_duplicates(subset=["OlinkID"]).reset_index(drop=True)) - sample_order = samples[sample_key].values + # Samples whose key is missing are kept as empty rows, as pivot_table left them + sample_order = pd.Index(samples[sample_key].astype(object).map(lambda v: v if pd.isna(v) else str(v))) - expression = df.pivot_table( - index=sample_key, - columns="OlinkID", - values="NPX", - aggfunc="first", - ) - expression = expression.reindex(sample_order).reset_index(drop=True) + expression = _long_to_wide(df, sample_key, sample_order, "NPX") logger.debug("Pivot shape: %d samples x %d features", expression.shape[0], expression.shape[1]) # Align features to match expression column order (pivot_table sorts columns) @@ -264,18 +407,13 @@ def read_olink_csv(path: str | Path) -> AffinityDataset: lod_col = next((c for c in ("LODNPX", "LOD", "PlateLOD") if c in df.columns), None) if lod_col is not None: logger.debug("LOD column %s present, building LOD matrix", lod_col) - lod_matrix = df.pivot_table( - index=sample_key, - columns="OlinkID", - values=lod_col, - aggfunc="first", - ) - lod_matrix = lod_matrix.reindex(sample_order).reset_index(drop=True) + lod_matrix = _long_to_wide(df, sample_key, sample_order, lod_col) metadata["lod_matrix"] = lod_matrix else: logger.debug("No LOD column in input") metadata.update(_qc_flag_matrices(df, sample_key, sample_order, expression.columns)) + del df platform = _detect_olink_platform(features["OlinkID"]) diff --git a/tests/test_qc_pipeline.py b/tests/test_qc_pipeline.py index 90fc572..fa99607 100644 --- a/tests/test_qc_pipeline.py +++ b/tests/test_qc_pipeline.py @@ -1264,3 +1264,54 @@ def test_qc_flags_first_value_wins_and_missing_stay_empty(self, tmp_path): m = read_olink_csv(path).metadata["sample_qc_matrix"] assert m.loc[0, "O1"] == "PASS" and m.loc[0, "O2"] == "WARN" assert pd.isna(m.loc[1, "O1"]) and pd.isna(m.loc[1, "O2"]) + + @staticmethod + def _write_long(path, sep=",", npx=("1.25", "2.5", "-0.75", "3.0")): + rows = [ + { + "SampleID": sid, + "OlinkID": oid, + "UniProt": "P" + oid[-1], + "Assay": "A" + oid[-1], + "Panel": "X", + "PlateID": "P1", + "MissingFreq": "30%" if oid == "OID1" else "0.1", + "NPX": v, + "LOD": "0.5", + "QC_Warning": "PASS", + "Unused": "x" * 20, + } + for (sid, oid), v in zip([("007", "OID1"), ("007", "OID2"), ("010", "OID1"), ("010", "OID2")], npx) + ] + pd.DataFrame(rows).to_csv(path, sep=sep, index=False) + + def test_streamed_and_single_read_agree(self, tmp_path, monkeypatch): + import pyprideap.io.readers.olink_csv as oc + + path = tmp_path / "long.npx.csv" + self._write_long(path) + small = oc.read_olink_csv(path) + monkeypatch.setattr(oc, "_STREAM_MIN_BYTES", 0) + streamed = oc.read_olink_csv(path) + assert small.expression.equals(streamed.expression) + assert small.samples.equals(streamed.samples) and small.features.equals(streamed.features) + + def test_ids_keep_leading_zeros_and_unused_columns_are_skipped(self, tmp_path): + from pyprideap.io.readers.olink_csv import read_olink_csv + + path = tmp_path / "ids.npx.csv" + self._write_long(path) + ds = read_olink_csv(path) + assert ds.samples["SampleID"].tolist() == ["007", "010"] + assert ds.expression.to_numpy().tolist() == [[1.25, 2.5], [-0.75, 3.0]] + assert "Unused" not in ds.samples.columns and "Unused" not in ds.features.columns + # MissingFreq written as "30%" stays text, as before + assert ds.features["MissingFreq"].tolist() == ["30%", "0.1"] + + def test_text_in_npx_falls_back_and_becomes_missing(self, tmp_path): + from pyprideap.io.readers.olink_csv import read_olink_csv + + path = tmp_path / "shifted.npx.csv" + self._write_long(path, sep=";", npx=("1.25", "PlateID", "-0.75", "3.0")) + ds = read_olink_csv(path) + assert ds.expression.isna().to_numpy().tolist() == [[False, True], [False, False]]