Skip to content

Commit e277da5

Browse files
authored
Merge pull request #725 from bigbio/fix/dia-box-plot-stats
dia: the remaining #717 fixes that did not reach dev (late #722 commits + #724)
2 parents c540dca + e08f531 commit e277da5

7 files changed

Lines changed: 326 additions & 53 deletions

File tree

pmultiqc/modules/common/dia_utils.py

Lines changed: 36 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
import itertools
22
import numpy as np
3+
from pmultiqc.modules.common.plots.general import run_to_sample_codes
34
import pandas as pd
45
import re
56
from collections import OrderedDict
@@ -346,8 +347,7 @@ def _sample_identification_counts(report_data: pd.DataFrame, file_df: pd.DataFra
346347
"""Counts per sample, de-duplicated across the sample's runs (what the per-run sets used to feed)."""
347348
if file_df is None or file_df.empty or not {"Sample", "Run"} <= set(file_df.columns):
348349
return dict()
349-
run_to_sample = file_df[["Run", "Sample"]].drop_duplicates().set_index("Run")["Sample"].astype(int)
350-
sample = report_data["Run"].astype(str).map(run_to_sample)
350+
sample = run_to_sample_codes(report_data["Run"], file_df)
351351
keep = sample.notna()
352352
if not keep.any():
353353
return dict()
@@ -402,24 +402,37 @@ def _handle_files_without_psm(ms_paths, ms_with_psm, cal_num_table_data):
402402

403403

404404
def _get_peptide_length(df):
405-
405+
"""{run: {length: count}} without materialising the sequences.
406+
407+
``.str.len()`` on the categorical column converts every row back to a
408+
Python string (231 M on PXD030304) and the per-run value_counts loop runs
409+
5,798 times - a transient that OOM-killed the summary after mod_plot_dict
410+
with no log line of its own (bigbio/pmultiqc#717). Length is computed once
411+
per distinct sequence and broadcast through the codes; the histogram is a
412+
single grouped size.
413+
"""
406414
if "Stripped.Sequence" not in df.columns:
407415
return None
408-
409-
df_sub = df[["Run", "Stripped.Sequence"]].copy()
410-
df_sub["length"] = df_sub["Stripped.Sequence"].str.len()
411-
416+
seqs = df["Stripped.Sequence"]
417+
if isinstance(seqs.dtype, pd.CategoricalDtype):
418+
per_category = seqs.cat.categories.astype(str).str.len().to_numpy()
419+
codes = seqs.cat.codes.to_numpy()
420+
lengths = np.where(codes >= 0, per_category[np.clip(codes, 0, None)], -1)
421+
else:
422+
lengths = seqs.astype(str).str.len().to_numpy()
423+
hist = (
424+
pd.DataFrame({"Run": df["Run"].to_numpy() if not hasattr(df["Run"], "cat") else df["Run"].cat.codes.to_numpy(), "length": lengths})
425+
.query("length >= 0")
426+
.groupby(["Run", "length"], sort=True)
427+
.size()
428+
)
429+
run_labels = df["Run"].cat.categories if hasattr(df["Run"], "cat") else None
412430
plot_data = {}
413-
for run, group in df_sub.groupby("Run", observed=True):
414-
stats_dict = group["length"].value_counts().sort_index().to_dict()
415-
plot_data[run] = stats_dict
416-
431+
for (run, length), count in hist.items():
432+
key = run_labels[run] if run_labels is not None else run
433+
plot_data.setdefault(key, {})[int(length)] = int(count)
417434
return plot_data
418435

419-
420-
## Removed draw_dia_heatmap wrapper; call cal_dia_heatmap and dia_plots.draw_heatmap directly.
421-
422-
423436
def draw_dia_intensitys(sub_section, report_df, sdrf_file_df):
424437
"""Draw the precursor intensity distribution and standard-deviation plots."""
425438
# Both consumers below only ever read these columns, so narrow the frame before
@@ -1070,14 +1083,14 @@ def dia_sample_level_modifications(df, sdrf_file_df):
10701083
if sdrf_file_df is None or sdrf_file_df.empty:
10711084
return {}
10721085

1073-
report_data = df.copy()
1074-
1075-
report_data = report_data.merge(
1076-
right=sdrf_file_df[["Sample", "Run"]].drop_duplicates(),
1077-
on="Run"
1078-
)
1079-
1080-
report_data["Sample"] = report_data["Sample"].astype(int)
1086+
# No merge: it upcast the categorical Run key to object for every row and
1087+
# was the transient that OOM-killed the summary in this stage (#717).
1088+
sample = run_to_sample_codes(df["Run"], sdrf_file_df)
1089+
keep = sample.notna()
1090+
# Keep Run: drop_duplicates below counts a peptidoform once per run within a
1091+
# sample, as the merge-based version did.
1092+
report_data = df.loc[keep, ["Run", "Modified.Sequence", "Modifications", "Protein.Group"]].copy()
1093+
report_data["Sample"] = sample[keep].astype(int).to_numpy()
10811094

10821095
mod_plot = dict()
10831096
for sample, group in report_data.groupby("Sample", sort=True, observed=True):

pmultiqc/modules/common/plots/dia.py

Lines changed: 53 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,9 @@
55

66
from pmultiqc.modules.common.plots.general import (
77
summarise_box_data,
8+
box_stats_by_group,
9+
run_to_sample_codes,
10+
FLAT_THRESHOLD,
811
plot_html_check,
912
plot_data_check
1013
)
@@ -75,38 +78,50 @@ def draw_heatmap(sub_section, hm_colors, heatmap_data):
7578
# Intensity Distribution
7679
def draw_dia_intensity_dis(sub_section, df, sdrf_file_df):
7780

78-
df_sub = df[["Run", "Modified.Sequence", "Protein.Group", "log_intensity"]].copy()
81+
large = len(df) >= FLAT_THRESHOLD
7982

8083
if not sdrf_file_df.empty:
81-
82-
df_sub = df_sub.merge(
83-
sdrf_file_df[["Sample", "Run"]].drop_duplicates(),
84-
on="Run"
84+
run_to_sample = (
85+
sdrf_file_df[["Sample", "Run"]].drop_duplicates().set_index("Run")["Sample"].astype(int)
8586
)
86-
87-
df_sub["Sample"] = df_sub["Sample"].astype(int)
88-
89-
box_data = [
90-
{
91-
(
92-
f"Sample {str(run)}"
93-
if data_type == "Sample"
94-
else str(run)
95-
): group["log_intensity"].dropna().tolist()
96-
for run, group in df_sub.groupby(data_type, sort=True, observed=True)
97-
}
98-
for data_type in ["Run", "Sample"]
99-
]
87+
if large:
88+
# No merge and no per-group lists. On PXD030304 the merge upcast the
89+
# categorical Run key to object for 231 M rows and the float lists
90+
# were built twice (by run, by sample), both alive at once: the
91+
# spike that OOM-killed the summary at 72 GB (bigbio/pmultiqc#717).
92+
# Map Run->Sample on the key and aggregate.
93+
sample = run_to_sample_codes(df["Run"], sdrf_file_df)
94+
keep = sample.notna()
95+
by_sample = df.loc[keep, ["log_intensity"]].assign(Sample=sample[keep].astype(int).to_numpy())
96+
box_data = [
97+
box_stats_by_group(df, "log_intensity", "Run"),
98+
box_stats_by_group(by_sample, "log_intensity", "Sample", label=lambda k: f"Sample {int(k)}"),
99+
]
100+
else:
101+
df_sub = df[["Run", "log_intensity"]].copy()
102+
df_sub["Sample"] = df_sub["Run"].astype(object).map(run_to_sample)
103+
df_sub = df_sub[df_sub["Sample"].notna()]
104+
df_sub["Sample"] = df_sub["Sample"].astype(int)
105+
box_data = [
106+
{
107+
(f"Sample {str(run)}" if data_type == "Sample" else str(run)): group["log_intensity"].dropna().tolist()
108+
for run, group in df_sub.groupby(data_type, sort=True, observed=True)
109+
}
110+
for data_type in ["Run", "Sample"]
111+
]
100112

101113
plot_label = ["by Run", "by Sample"]
102114

103115
else:
104-
box_data = [
105-
{
106-
str(run): group["log_intensity"].dropna().tolist()
107-
for run, group in df.groupby("Run", observed=True)
108-
}
109-
]
116+
if large:
117+
box_data = [box_stats_by_group(df, "log_intensity", "Run")]
118+
else:
119+
box_data = [
120+
{
121+
str(run): group["log_intensity"].dropna().tolist()
122+
for run, group in df.groupby("Run", observed=True)
123+
}
124+
]
110125
plot_label = ["by Run"]
111126

112127
draw_config = {
@@ -149,10 +164,14 @@ def draw_dia_intensity_dis(sub_section, df, sdrf_file_df):
149164
# Ms1.Area non-normalised MS1 peak area
150165
def draw_dia_ms1_area(sub_section, df):
151166

152-
box_data = {
153-
str(run): group["log_ms1_area"].dropna().tolist()
154-
for run, group in df.groupby("Run", observed=True)
155-
}
167+
if len(df) >= FLAT_THRESHOLD:
168+
# No per-run Python lists on large reports (bigbio/pmultiqc#717).
169+
box_data = box_stats_by_group(df, "log_ms1_area", "Run")
170+
else:
171+
box_data = {
172+
str(run): group["log_ms1_area"].dropna().tolist()
173+
for run, group in df.groupby("Run", observed=True)
174+
}
156175

157176
draw_config = {
158177
"id": "ms1_area_distribution_box",
@@ -164,6 +183,9 @@ def draw_dia_ms1_area(sub_section, df):
164183
"save_data_file": False,
165184
}
166185

186+
# 5,798 runs x raw MS1 areas is >4 GiB of serialised points and panics polars
187+
# (bigbio/pmultiqc#717). Above the flat threshold hand MultiQC the box statistics.
188+
box_data = summarise_box_data(box_data)
167189
box_html = box.plot(list_of_data_by_sample=box_data, pconfig=draw_config)
168190

169191
# box_html.flat
@@ -312,6 +334,8 @@ def draw_dia_intensity_std(sub_section, df, sdrf_file_df):
312334
"save_data_file": False,
313335
}
314336

337+
# Same failure mode as draw_dia_ms1_area: 42.7 M points on the big DIA run.
338+
box_data = summarise_box_data(box_data)
315339
box_html = box.plot(
316340
list_of_data_by_sample=box_data,
317341
pconfig=draw_box_config,

pmultiqc/modules/common/plots/general.py

Lines changed: 62 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -505,3 +505,65 @@ def summarise_box_data(plot_data, whisker_iqr=1.5, max_points=None):
505505
summarised.append(out)
506506

507507
return summarised if was_list else summarised[0]
508+
509+
510+
def run_to_sample_codes(runs: pd.Series, file_df: pd.DataFrame) -> pd.Series:
511+
"""Sample id per row from the Run column, without materialising strings.
512+
513+
Merging the report with the SDRF sample table on ``Run`` was measured to be
514+
the largest transient in the summary: it upcasts the categorical key to
515+
object for every row (231 M on PXD030304) and adds an object Sample column.
516+
Mapping through the category codes touches one small array instead.
517+
Rows whose run is not in ``file_df`` get NaN.
518+
"""
519+
run_to_sample = file_df[["Run", "Sample"]].drop_duplicates().set_index("Run")["Sample"]
520+
run_to_sample.index = run_to_sample.index.astype(str)
521+
if isinstance(runs.dtype, pd.CategoricalDtype):
522+
per_category = run_to_sample.reindex(runs.cat.categories.astype(str)).to_numpy(dtype="float64", na_value=np.nan)
523+
codes = runs.cat.codes.to_numpy()
524+
values = np.where(codes >= 0, per_category[np.clip(codes, 0, None)], np.nan)
525+
return pd.Series(values, index=runs.index)
526+
return runs.astype(str).map(run_to_sample).astype("float64")
527+
528+
529+
def box_stats_by_group(df, value_col, key, whisker_iqr=1.5, label=str):
530+
"""Box statistics per group without materialising per-group Python lists.
531+
532+
Same numbers as ``summarise_box_data`` applied to ``{group: values.tolist()}``
533+
-- percentile q1/median/q3, Tukey whiskers (min/max within q1-1.5*IQR ..
534+
q3+1.5*IQR, or the plain min/max when IQR is 0), and the mean -- computed
535+
with vectorised groupbys. On PXD030304 the list form was 231 M Python
536+
floats built twice (by run and by sample) on top of a merge that upcast the
537+
categorical run key to object: the spike that OOM-killed the summary at
538+
72 GB (bigbio/pmultiqc#717).
539+
540+
Returns ``{label(group): {min, q1, median, q3, max, mean}}`` for groups
541+
with at least one finite value.
542+
"""
543+
values = pd.to_numeric(df[value_col], errors="coerce")
544+
finite = np.isfinite(values.to_numpy(dtype="float64", na_value=np.nan))
545+
sub = pd.DataFrame({"k": df[key].to_numpy() if not hasattr(df[key], "cat") else df[key].astype(object).to_numpy(),
546+
"v": values.to_numpy(dtype="float64", na_value=np.nan)})[finite]
547+
if sub.empty:
548+
return {}
549+
g = sub.groupby("k", sort=True)
550+
q = g["v"].quantile([0.25, 0.5, 0.75]).unstack()
551+
q.columns = ["q1", "median", "q3"]
552+
stats = q.assign(mean=g["v"].mean(), lo=g["v"].min(), hi=g["v"].max())
553+
iqr = stats["q3"] - stats["q1"]
554+
low_bound = (stats["q1"] - whisker_iqr * iqr)
555+
high_bound = (stats["q3"] + whisker_iqr * iqr)
556+
lb = sub["k"].map(low_bound).to_numpy(); hb = sub["k"].map(high_bound).to_numpy()
557+
v = sub["v"].to_numpy()
558+
fenced_min = sub[v >= lb].groupby("k")["v"].min()
559+
fenced_max = sub[v <= hb].groupby("k")["v"].max()
560+
use_fence = iqr > 0
561+
stats["min"] = np.where(use_fence, fenced_min.reindex(stats.index), stats["lo"])
562+
stats["max"] = np.where(use_fence, fenced_max.reindex(stats.index), stats["hi"])
563+
out = {}
564+
for k, row in stats.iterrows():
565+
out[label(k)] = {
566+
"min": float(row["min"]), "q1": float(row["q1"]), "median": float(row["median"]),
567+
"q3": float(row["q3"]), "max": float(row["max"]), "mean": float(row["mean"]),
568+
}
569+
return out

pmultiqc/modules/quantms/quantms.py

Lines changed: 8 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -664,8 +664,15 @@ def draw_plots(self):
664664
name="draw_quantms_time_section"
665665
)
666666

667-
if self.qpx_source is None and self.msstats_input_valid:
667+
# parse_msstats_input only builds the peptide/protein quantification
668+
# tables. Do not read the MSstats input when tables are disabled: on
669+
# PXD030304 it is a 20.6 GiB, 231 M-row CSV that pd.read_csv turned
670+
# into 93 GB of Python strings for output that was then discarded
671+
# (bigbio/pmultiqc#717).
672+
if self.qpx_source is None and self.msstats_input_valid and not config.kwargs.get("disable_table", False):
668673
self.parse_msstats_input()
674+
elif self.msstats_input_valid and config.kwargs.get("disable_table", False):
675+
log.info("Tables disabled; skipping the MSstats input parse.")
669676
quant_method = config.kwargs.get("quantification_method", None)
670677

671678
if (

tests/test_dia_box_plots.py

Lines changed: 90 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,90 @@
1+
"""DIA box plots must hand MultiQC summary statistics, not raw points, on large runs (#717).
2+
3+
On PXD030304 (5,798 runs) draw_dia_ms1_area serialised every raw MS1 area into the
4+
box plot and polars panicked on a >4 GiB buffer; draw_dia_intensity_std carried
5+
42.7 M points. Above the flat threshold both must pass {min,q1,median,q3,max,mean}.
6+
"""
7+
8+
import numpy as np
9+
import pandas as pd
10+
import pytest
11+
12+
from pmultiqc.modules.common.plots import dia as dia_plots
13+
from pmultiqc.modules.common.plots import general
14+
15+
16+
@pytest.fixture
17+
def capture_box(monkeypatch):
18+
seen = []
19+
monkeypatch.setattr(dia_plots.box, "plot", lambda list_of_data_by_sample, pconfig=None: seen.append(list_of_data_by_sample) or "html")
20+
monkeypatch.setattr(dia_plots, "add_sub_section", lambda **kw: None)
21+
monkeypatch.setattr(dia_plots, "plot_html_check", lambda h: h, raising=False)
22+
return seen
23+
24+
25+
def _big_ms1(n_runs=20, per_run=None):
26+
per_run = per_run or (general.FLAT_THRESHOLD // n_runs + 1)
27+
rng = np.random.default_rng(0)
28+
return pd.DataFrame({
29+
"Run": np.repeat([f"run{i}" for i in range(n_runs)], per_run),
30+
"log_ms1_area": rng.normal(20, 2, n_runs * per_run),
31+
})
32+
33+
34+
def test_ms1_area_uses_summary_stats_above_threshold(capture_box):
35+
dia_plots.draw_dia_ms1_area(None, _big_ms1())
36+
(data,) = capture_box
37+
assert data, "box.plot received no data"
38+
for run, stats in data.items():
39+
assert isinstance(stats, dict), f"{run}: raw list reached box.plot"
40+
assert {"min", "q1", "median", "q3", "max", "mean"} <= set(stats)
41+
42+
43+
def test_ms1_area_keeps_raw_points_below_threshold(capture_box):
44+
dia_plots.draw_dia_ms1_area(None, _big_ms1(n_runs=2, per_run=10))
45+
(data,) = capture_box
46+
assert all(isinstance(v, list) for v in data.values()), "small reports should keep raw points"
47+
48+
49+
def test_intensity_std_uses_summary_stats_above_threshold(capture_box, monkeypatch):
50+
n = general.FLAT_THRESHOLD + 10
51+
fake = [{"Sample 1": list(np.linspace(0, 1, n))}]
52+
monkeypatch.setattr(dia_plots, "calculate_dia_intensity_std", lambda df, sdrf: fake)
53+
dia_plots.draw_dia_intensity_std(None, pd.DataFrame(), pd.DataFrame())
54+
(data,) = capture_box
55+
ds = data[0] if isinstance(data, list) else data
56+
assert isinstance(ds["Sample 1"], dict)
57+
assert {"min", "q1", "median", "q3", "max", "mean"} <= set(ds["Sample 1"])
58+
59+
60+
def test_box_stats_by_group_matches_summarise_box_data():
61+
"""Vectorised per-group statistics must equal the list-based summary (#717)."""
62+
rng = np.random.default_rng(7)
63+
n = general.FLAT_THRESHOLD + 500
64+
keys = rng.choice([f"r{i}" for i in range(12)], n)
65+
vals = np.concatenate([rng.normal(20, 2, n - 40), np.full(20, 20.0), rng.normal(20, 40, 20)]) # ties + outliers
66+
vals[:5] = np.nan
67+
df = pd.DataFrame({"Run": pd.Categorical(keys), "x": vals})
68+
lists = {str(k): g["x"].dropna().tolist() for k, g in df.groupby("Run", observed=True)}
69+
expected = general.summarise_box_data(lists)
70+
got = general.box_stats_by_group(df, "x", "Run")
71+
assert got.keys() == expected.keys()
72+
for k in expected:
73+
for stat in ("min", "q1", "median", "q3", "max", "mean"):
74+
assert np.isclose(got[k][stat], expected[k][stat]), (k, stat, got[k][stat], expected[k][stat])
75+
76+
77+
def test_intensity_dis_by_sample_without_merge(capture_box):
78+
n = general.FLAT_THRESHOLD + 10
79+
df = pd.DataFrame({
80+
"Run": pd.Categorical(np.repeat(["a", "b"], n // 2)),
81+
"Modified.Sequence": pd.Categorical(["P"] * n),
82+
"Protein.Group": pd.Categorical(["G"] * n),
83+
"log_intensity": np.linspace(10, 30, n),
84+
})
85+
sdrf = pd.DataFrame({"Run": ["a", "b"], "Sample": [1, 1]})
86+
dia_plots.draw_dia_intensity_dis(None, df, sdrf)
87+
(data,) = capture_box
88+
assert isinstance(data, list) and len(data) == 2
89+
assert set(data[0]) == {"a", "b"} and set(data[1]) == {"Sample 1"}
90+
assert all(isinstance(v, dict) for ds in data for v in ds.values())

0 commit comments

Comments
 (0)