Skills › Scientific Agent Skills › depmap › depmap_data.py
depmap_data.py
201 lines · skills/depmap/scripts/depmap_data.py at 5f2bc03
1"""Validated release discovery and local DepMap analysis helpers.
3The catalogue is metadata; recent rows may deliberately omit download URLs.
4Analysis functions expect independent, model-level observations selected by the caller.
5"""
6from __future__ import annotations
8import csv
9import hashlib
10import io
11import re
12from pathlib import Path
13from urllib.request import urlopen
15import numpy as np
16import pandas as pd
17from scipy import stats
19CATALOGUE_URL = "https://depmap.org/portal/api/no-captcha/download/files"
22def fetch_catalogue() -> pd.DataFrame:
23 """Fetch the complete public CSV catalogue; reject HTML verification pages."""
24 with urlopen(CATALOGUE_URL, timeout=60) as response:
25 content_type = response.headers.get("Content-Type", "").lower()
26 text = response.read().decode("utf-8-sig")
27 if "html" in content_type or text.lstrip().startswith("<"):
28 raise ValueError("Received HTML, not DepMap catalogue CSV; use the portal downloads page")
29 frame = pd.read_csv(io.StringIO(text), dtype=str, keep_default_na=False)
30 required = {"release", "release_date", "filename", "md5_hash"}
31 if not required.issubset(frame.columns) or frame.empty:
32 raise ValueError("Unrecognized or empty DepMap catalogue")
33 return frame
36def select_file(catalogue: pd.DataFrame, release: str, filename: str) -> pd.Series:
37 """Require exactly one filename within an explicitly pinned release."""
38 rows = catalogue.loc[catalogue["release"].eq(release) & catalogue["filename"].eq(filename)]
39 if len(rows) != 1:
40 raise ValueError(f"Expected one catalogue row for {release!r}/{filename!r}; found {len(rows)}")
41 return rows.iloc[0]
44def verify_md5(path: str | Path, expected: str) -> str:
45 """Verify published file identity; MD5 is an integrity check, not authentication."""
46 if not re.fullmatch(r"[0-9a-fA-F]{32}", expected):
47 raise ValueError("A valid published MD5 is required")
48 digest = hashlib.md5(usedforsecurity=False)
49 with open(path, "rb") as stream:
50 for chunk in iter(lambda: stream.read(1024 * 1024), b""):
51 digest.update(chunk)
52 actual = digest.hexdigest()
53 if actual != expected.lower():
54 raise ValueError(f"MD5 mismatch for {Path(path).name}")
55 return actual
58def _unique_index(frame: pd.DataFrame | pd.Series) -> None:
59 if not frame.index.is_unique or frame.index.isna().any():
60 raise ValueError("Observation IDs must be unique and nonmissing")
63def _numeric(frame: pd.DataFrame) -> pd.DataFrame:
64 _unique_index(frame)
65 if frame.empty:
66 raise ValueError("Expected a nonempty measurement matrix")
67 if not frame.columns.is_unique:
68 raise ValueError("Feature identifiers must be unique")
69 numeric = frame.apply(pd.to_numeric, errors="raise")
70 if np.isinf(numeric.to_numpy(dtype=float)).any():
71 raise ValueError("Infinite measurements are invalid; retain missing values as NaN")
72 return numeric
75def load_gene_effect(path: str | Path) -> pd.DataFrame:
76 """Read model-level CRISPRGeneEffect, preserving Symbol (EntrezID) labels."""
77 with open(path, encoding="utf-8-sig", newline="") as stream:
78 header = next(csv.reader(stream))
79 # Check before pandas can silently mangle duplicate CSV column names.
80 if len(header) != len(set(header)):
81 raise ValueError("Duplicate CSV header")
82 frame = pd.read_csv(path, index_col=0)
83 if not frame.index.astype(str).str.fullmatch(r"ACH-\d+").all():
84 raise ValueError("Expected model IDs (ACH-...); map screen/condition IDs explicitly")
85 if not all(re.fullmatch(r".+ \(\d+\)", str(c)) for c in frame.columns):
86 raise ValueError("Expected complete Symbol (EntrezID) gene labels")
87 frame.index.name = "ModelID"
88 return _numeric(frame)
91def load_models(path: str | Path) -> pd.DataFrame:
92 """Read current Model.csv; legacy sample_info needs a deliberate conversion."""
93 frame = pd.read_csv(path, dtype={"ModelID": str})
94 required = {"ModelID", "CellLineName", "OncotreeLineage", "OncotreePrimaryDisease"}
95 if not required.issubset(frame.columns):
96 raise ValueError(f"Model.csv must contain {sorted(required)}")
97 frame = frame.set_index("ModelID")
98 _unique_index(frame)
99 return frame
102def default_model_rows(frame: pd.DataFrame, flag: str = "IsDefaultEntryForModel") -> pd.DataFrame:
103 """Select current omics default rows; refuse ambiguous per-model duplication.
105 Pass the actual flag column for the pinned release. A 'Yes' flag does not
106 decide which assay/datatype to analyze: subset the datatype first.
107 """
108 values = frame[flag].astype("string").str.strip().str.lower()
109 if not values.dropna().isin(["yes", "no", "true", "false", "1", "0"]).all():
110 raise ValueError("Unrecognized default-entry flag values")
111 selected = frame.loc[values.isin(["yes", "true", "1"])].copy()
112 if selected.empty:
113 raise ValueError("No default entries found")
114 selected = selected.set_index("ModelID")
115 _unique_index(selected)
116 return selected
119def gene_column(frame: pd.DataFrame, identifier: str) -> str:
120 """Resolve exact label, symbol, or Entrez ID without collapsing collisions."""
121 if not frame.columns.is_unique:
122 raise ValueError("Duplicate gene identifiers")
123 if identifier in frame.columns:
124 return identifier
125 matches = []
126 for column in frame.columns:
127 match = re.fullmatch(r"(.+) \((\d+)\)", str(column))
128 if match and identifier in match.groups():
129 matches.append(column)
130 if len(matches) != 1:
131 raise ValueError(f"Gene {identifier!r} resolves to {len(matches)} columns; use an exact label")
132 return matches[0]
135def gene_profile(effects: pd.DataFrame, models: pd.DataFrame, gene: str) -> pd.DataFrame:
136 """Return all observed effects plus a checked one-to-one Model.csv join."""
137 _unique_index(effects)
138 _unique_index(models)
139 missing = effects.index.difference(models.index)
140 if len(missing):
141 raise ValueError(f"{len(missing)} effect rows lack Model.csv metadata")
142 scores = effects[gene_column(effects, gene)].dropna().rename("gene_effect")
143 return scores.to_frame().join(models, validate="one_to_one").sort_values("gene_effect")
146def biomarker_scan(effects: pd.DataFrame, status: pd.Series, min_n: int = 5) -> pd.DataFrame:
147 """Exploratory one-sided mutant-vs-assayed-negative Mann-Whitney tests.
149 1 = prespecified biomarker present, 0 = assayed negative, NaN = unknown.
150 BH is applied to every eligible gene before any result filtering. This does
151 not adjust for lineage, batch, ploidy, or related models.
152 """
153 if min_n < 2:
154 raise ValueError("min_n must be at least 2")
155 effects = _numeric(effects)
156 _unique_index(status)
157 status = pd.to_numeric(status, errors="raise")
158 if not status.dropna().isin([0, 1]).all():
159 raise ValueError("Biomarker must be 0, 1, or missing")
160 common = effects.index.intersection(status.dropna().index)
161 statuses = status.loc[common]
162 result = []
163 for gene in effects.columns:
164 observed = effects.loc[common, gene].dropna()
165 mutated = observed.loc[statuses.reindex(observed.index).eq(1)]
166 negative = observed.loc[statuses.reindex(observed.index).eq(0)]
167 if min(len(mutated), len(negative)) < min_n:
168 continue
169 pvalue = stats.mannwhitneyu(mutated, negative, alternative="less", method="auto").pvalue
170 result.append({"gene": gene, "n_mutated": len(mutated), "n_assayed_negative": len(negative),
171 "mean_mutated": mutated.mean(), "mean_assayed_negative": negative.mean(),
172 "effect_size": negative.mean() - mutated.mean(), "pval": pvalue})
173 columns = ["gene", "n_mutated", "n_assayed_negative", "mean_mutated", "mean_assayed_negative", "effect_size", "pval", "qval"]
174 output = pd.DataFrame(result, columns=columns)
175 if not output.empty:
176 output["qval"] = stats.false_discovery_control(output["pval"].to_numpy(), method="bh")
177 output = output.sort_values(["qval", "effect_size"], ascending=[True, False])
178 return output
181def coessentiality(effects: pd.DataFrame, gene: str, min_n: int = 500) -> pd.DataFrame:
182 """Exploratory Pearson associations with pairwise n and BH-adjusted p-values."""
183 if min_n < 3:
184 raise ValueError("min_n must be at least 3")
185 effects = _numeric(effects)
186 target = gene_column(effects, gene)
187 result = []
188 for column in effects.columns:
189 if column == target:
190 continue
191 pair = effects[[target, column]].dropna()
192 if len(pair) < min_n or pair.nunique().min() < 2:
193 continue
194 correlation = stats.pearsonr(pair[target], pair[column])
195 if np.isfinite(correlation.statistic) and np.isfinite(correlation.pvalue):
196 result.append({"gene": column, "n": len(pair), "r": correlation.statistic, "pval": correlation.pvalue})
197 output = pd.DataFrame(result, columns=["gene", "n", "r", "pval", "qval"])
198 if not output.empty:
199 output["qval"] = stats.false_discovery_control(output["pval"].to_numpy(), method="bh")
200 output = output.sort_values("r", ascending=False)
201 return output