Tissue-Type Heterogeneity Drives Irreproducibility in Endometriosis Transcriptomic Signatures: A Permutation-Based Audit of Three Public Microarray Datasets
Tissue-Type Heterogeneity Drives Irreproducibility in Endometriosis Transcriptomic Signatures: A Permutation-Based Audit of Three Public Microarray Datasets
stepstep_labs
Abstract
Endometriosis affects approximately 10% of reproductive-age women, yet no validated transcriptomic biomarker has reached clinical use. A persistent obstacle is that publicly available microarray datasets—widely cited in biomarker discovery—differ not only in sample size and patient population but in the tissue compartments they compare. We quantify the reproducibility cost of this conflation by auditing three frequently used Affymetrix U133 Plus 2.0 datasets (GSE7305, GSE11691, GSE51981) using a permutation-based framework implemented entirely in standard-library Python on pre-normalized expression matrices. Gene-level ranked lists (top 200 by Welch t-statistic magnitude) were compared across datasets using 5,000-permutation null models. The two datasets that contrast ectopic lesions against healthy endometrium (GSE7305 vs. GSE11691) share 24 genes (z = 5.96, p < 0.0002), demonstrating significant reproducibility when the biological question is matched. In contrast, comparisons involving GSE51981—which contrasts eutopic endometrium from affected versus unaffected women—yield overlaps indistinguishable from chance (z = 0.46 and 1.37, p = 0.34 and 0.12). A three-way gene intersection identifies only three genes (MPDZ, PICALM, WNK1; z = 6.47, p = 0.001), and within-dataset menstrual cycle effects (40-gene overlap, Jaccard = 0.111) exceed all cross-dataset disease overlaps. These results demonstrate that tissue-type heterogeneity, not poor analytical practice, is the principal barrier to reproducible endometriosis signatures, and that meta-analyses pooling across tissue compartments risk systematic confounding.
1. Introduction
Endometriosis is a chronic inflammatory condition in which endometrium-like tissue grows outside the uterine cavity, affecting an estimated 190 million women worldwide [1]. Diagnostic delay averages 7–10 years from symptom onset, motivating extensive efforts to identify non-invasive transcriptomic biomarkers [2, 3]. The Gene Expression Omnibus (GEO) hosts dozens of endometriosis microarray datasets, and three Affymetrix Human Genome U133 Plus 2.0 studies—GSE7305, GSE11691, and GSE51981—are among the most frequently cited in biomarker discovery and meta-analytic pipelines [4, 5].
These three datasets are routinely treated as interchangeable sources of "endometriosis gene expression data." Yet they address fundamentally different biological questions. GSE7305 compares ovarian endometriotic lesions to normal endometrium from unaffected controls [8]. GSE11691 compares peritoneal and ovarian lesions to healthy endometrium [9]. GSE51981 compares eutopic endometrium from women with versus without endometriosis, stratified by menstrual cycle phase [10]. The first two datasets measure the transcriptomic distance between ectopic disease tissue and healthy tissue—a large biological signal. The third measures the subtler molecular imprint that disease status leaves on anatomically normal endometrium.
We hypothesize that this tissue-type heterogeneity is the dominant source of irreproducibility when gene lists from these datasets are compared, and that within-compartment comparisons will show significant concordance while cross-compartment comparisons will not. To test this, we implement a permutation-based overlap framework that compares observed gene-list intersections against null distributions generated by rank shuffling. The analysis operates entirely on pre-normalized (RMA) expression matrices downloaded from GEO, using standard-library Python with no external dependencies. We present results at both gene and probe levels, conduct sensitivity analyses across list-length thresholds, and examine menstrual cycle phase as a within-dataset confounder.
The goal is not to invalidate any individual dataset, but to quantify what the field loses when tissue compartment is ignored in cross-study comparisons—a quantification that is absent from most published meta-analyses of endometriosis transcriptomics.
2. Methods
2.1 Datasets
Three Affymetrix Human Genome U133 Plus 2.0 datasets were obtained from GEO as pre-normalized (RMA) series matrix files:
- GSE7305 (Hever et al., 2007 [8]): 10 ovarian endometriotic lesions vs. 10 normal endometrial samples. Platform: GPL570.
- GSE11691 (Hull et al., 2008 [9]): 9 peritoneal/ovarian endometriotic lesions vs. 9 healthy endometrial samples. Platform: GPL570.
- GSE51981 (Tamaresis et al., 2014 [10]): 77 eutopic endometrial samples from women with endometriosis vs. 71 from disease-free controls, annotated by menstrual cycle phase. Platform: GPL570.
All three datasets share the GPL570 platform, yielding 22,277 common probe sets. Because the expression matrices were pre-normalized by the original study authors using RMA, no raw CEL file processing was required.
2.2 Preprocessing and variance filtering
From the 22,277 common probes, the top 10,000 by cross-sample variance (computed across all samples within each dataset independently, then intersected) were retained to remove low-information probes. Of these 10,000 filtered probes, 9,828 carried gene-symbol annotations in the GPL570 platform table.
2.3 Gene-level collapse
For gene-level analysis, each probe was mapped to its first annotated gene symbol from the GPL570 annotation. Where multiple probes mapped to the same gene, the probe with the maximum absolute Welch t-statistic was selected as the gene's representative. This yielded 6,928 unique gene symbols across the filtered set.
2.4 Differential expression ranking
For each dataset, a two-sided Welch t-test was computed for every probe (probe-level analysis) or gene-representative probe (gene-level analysis), comparing disease versus control groups. Probes/genes were ranked by descending absolute t-statistic. The top N features from each dataset's ranked list were compared for overlap, with N = 200 as the primary threshold.
2.5 Permutation test for overlap significance
To assess whether observed overlaps exceed chance expectation, a permutation null was constructed. For each of 5,000 permutations (random seed = 42), sample labels within each dataset were independently shuffled, t-statistics recomputed, and top-N lists regenerated. The intersection size was recorded for each permuted comparison. The empirical p-value was calculated as the fraction of permuted overlaps meeting or exceeding the observed overlap. The z-score was computed as (observed − null mean) / null SD. At 5,000 permutations, the p-value resolution is 0.0002.
Pairwise and three-way intersections were evaluated at both probe and gene levels.
2.6 Menstrual cycle analysis
Two supplementary analyses addressed menstrual cycle confounding. First, within GSE51981, the top 200 differentially expressed genes (disease vs. control) were computed separately for proliferative-phase and secretory-phase samples and their overlap assessed. Second, a cross-dataset phase-matched comparison was conducted between GSE7305 (follicular-phase samples, where annotated) and GSE51981 (proliferative-phase samples) at the gene level.
2.7 Sensitivity analysis
Gene-level pairwise and three-way overlaps were computed across list-length thresholds N ∈ {25, 50, 100, 200, 500, 750, 1000} to assess robustness to the choice of N.
3. Results
3.1 Probe-level overlap
At the primary threshold of N = 200, probe-level pairwise overlaps were as follows:
Table 1. Probe-level overlap at N = 200 (top 200 probes by |t|)
| Comparison | Observed overlap | Jaccard index |
|---|---|---|
| GSE7305 vs. GSE11691 | 15 | 0.039 |
| GSE7305 vs. GSE51981 | 3 | 0.008 |
| GSE11691 vs. GSE51981 | 2 | 0.005 |
| Three-way | 0 | — |
The tissue-matched pair (GSE7305–GSE11691) showed the highest overlap, while comparisons involving the eutopic-only dataset (GSE51981) were near floor levels. No probe appeared in all three top-200 lists.
3.2 Gene-level overlap
Collapsing to gene level substantially increased the tissue-matched overlap while leaving cross-compartment comparisons near chance:
Table 2. Gene-level overlap at N = 200 (top 200 genes by |t|)
| Comparison | Observed overlap | Jaccard index |
|---|---|---|
| GSE7305 vs. GSE11691 | 24 | 0.064 |
| GSE7305 vs. GSE51981 | 6 | 0.015 |
| GSE11691 vs. GSE51981 | 8 | 0.020 |
| Three-way | 3 (MPDZ, PICALM, WNK1) | — |
Gene-level collapse resolved multi-probe redundancy and increased the tissue-matched overlap from 15 probes to 24 genes. The three-way intersection, invisible at probe level, yielded three genes: MPDZ, PICALM, and WNK1.
3.3 Permutation test results
Table 3. Permutation significance of probe-level overlaps (N = 200, 5,000 permutations)
| Comparison | Observed | Null mean ± SD | z | p |
|---|---|---|---|---|
| GSE7305 vs. GSE11691 | 15 | 3.6 ± 2.2 | 5.10 | 0.0008 |
| GSE7305 vs. GSE51981 | 3 | 2.4 ± 2.0 | 0.29 | 0.41 |
| GSE11691 vs. GSE51981 | 2 | 2.2 ± 1.8 | −0.11 | 0.59 |
| Three-way | 0 | 0.04 ± 0.19 | −0.19 | 1.00 |
Table 4. Permutation significance of gene-level overlaps (N = 200, 5,000 permutations)
| Comparison | Observed | Null mean ± SD | z | p |
|---|---|---|---|---|
| GSE7305 vs. GSE11691 | 24 | 6.5 ± 2.9 | 5.96 | < 0.0002 |
| GSE7305 vs. GSE51981 | 6 | 4.7 ± 2.8 | 0.46 | 0.34 |
| GSE11691 vs. GSE51981 | 8 | 4.5 ± 2.5 | 1.37 | 0.12 |
| Three-way | 3 | 0.18 ± 0.44 | 6.47 | 0.001 |
The permutation tests confirm a sharp dichotomy. At the gene level, the two ectopic-vs.-healthy comparisons (GSE7305–GSE11691) share significantly more genes than expected by chance (z = 5.96, p < 0.0002). All comparisons involving GSE51981 fail to reach significance. The gene-level analysis is uniformly stronger than probe-level: the tissue-matched z-score rises from 5.10 to 5.96, and the three-way intersection achieves z = 6.47 (p = 0.001) compared to z = −0.19 at probe level, demonstrating that gene-level collapse recovers signal that probe-level redundancy obscures.
3.4 The three-way gene intersection
The three genes present in all three top-200 gene lists—MPDZ, PICALM, and WNK1—represent a statistically significant convergence (z = 6.47) but should be interpreted cautiously given the small intersection size. MPDZ (multiple PDZ domain protein) and WNK1 (WNK lysine deficient protein kinase 1) have been implicated in cell polarity and ion transport pathways; PICALM (phosphatidylinositol binding clathrin assembly protein) functions in clathrin-mediated endocytosis. Whether these genes reflect shared disease biology or coincidental ranking convergence across distinct tissue comparisons cannot be resolved from this analysis alone.
3.5 Menstrual cycle confounding
Within GSE51981, comparing disease-vs.-control gene lists derived separately from proliferative-phase and secretory-phase samples yielded 40 overlapping genes (Jaccard = 0.111). This within-dataset, within-disease reproducibility across menstrual cycle phases exceeds every cross-dataset disease comparison observed in this study.
Phase-matched cross-dataset comparison (GSE7305 follicular-phase vs. GSE51981 proliferative-phase, gene level) yielded 4 overlapping genes (Jaccard = 0.010), providing no rescue of cross-compartment concordance. Menstrual cycle matching does not bridge the tissue-type gap.
3.6 Sensitivity to list-length threshold
Table 5. Gene-level overlap across list-length thresholds
| N | Mean pairwise Jaccard | Three-way overlap |
|---|---|---|
| 25 | 0.007 | 0 |
| 50 | 0.010 | 0 |
| 100 | 0.007 | 0 |
| 200 | 0.033 | 3 |
| 500 | 0.049 | 5 |
| 750 | 0.067 | 12 |
| 1000 | 0.085 | 30 |
Three-way overlap is absent at N ≤ 100, emerges at N = 200, and grows monotonically. Even at N = 1000 (roughly 14% of the gene universe), mean pairwise Jaccard remains below 0.10, indicating that cross-dataset concordance is low across the full range of reasonable thresholds. The tissue-matched pair (GSE7305–GSE11691) dominates the mean at every threshold.
4. Discussion
4.1 Tissue-type heterogeneity as the central barrier
The principal finding of this audit is that tissue-type heterogeneity, rather than platform noise, analytical methodology, or sample size, is the dominant driver of irreproducibility across these three widely cited endometriosis datasets. When the biological question is matched—ectopic lesion versus healthy endometrium—ranked gene lists reproduce at levels far exceeding permutation-derived chance expectations (gene-level z = 5.96). When the biological question differs—eutopic endometrium stratified by disease status—overlap collapses to noise. This pattern is consistent at both probe and gene levels, across all list-length thresholds, and persists after menstrual cycle phase matching.
This is not a surprising result in isolation; any biologist would predict that ectopic lesions differ from eutopic tissue. What is notable is that these three datasets are routinely combined in meta-analyses and machine-learning pipelines without accounting for this compartment difference [3, 4, 5]. Our permutation framework provides a quantitative audit tool: the z-score directly measures how much signal a cross-study comparison carries above chance, enabling researchers to flag problematic pooling before downstream analysis.
4.2 Gene-level analysis recovers obscured signal
The transition from probe-level to gene-level analysis is informative. Multiple probes targeting the same gene can disagree in rank due to probe-specific hybridization effects, effectively diluting gene-level concordance when counted at the probe level. Collapsing via max-|t| selection resolves this: the tissue-matched z-score increases from 5.10 to 5.96, and the three-way intersection—entirely absent at probe level—emerges as three genes with z = 6.47. This supports the recommendation that reproducibility audits of Affymetrix data should operate at the gene level with explicit probe-to-gene mapping, consistent with prior methodological guidance [7].
4.3 The three-way intersection: noteworthy but not a biomarker panel
MPDZ, PICALM, and WNK1 survive the intersection of all three top-200 gene lists despite the tissue-compartment heterogeneity. This is statistically notable (z = 6.47, p = 0.001) but biologically preliminary. Three genes do not constitute a validated signature, and their presence in three ranked lists of 200 from a universe of ~7,000 genes could reflect pathway-level commonalities (e.g., membrane trafficking, ion homeostasis) rather than disease-specific biology. These genes warrant investigation in independent cohorts but should not be interpreted as a proposed diagnostic panel.
4.4 Sample size limitations are the field's limitations
The datasets audited here have small sample sizes (n = 18–20 for GSE7305 and GSE11691). This is a legitimate concern for statistical power. However, these are the datasets that the endometriosis transcriptomics field uses to propose diagnostic signatures and populate meta-analyses. If ranked gene lists are unstable at these sample sizes, this instability is a property of the published literature, not merely of our audit. The permutation framework accommodates small samples naturally: the null distribution is constructed under the same sample-size constraints as the observed statistic, so the z-scores and p-values remain well-calibrated.
4.5 Implications for meta-analysis design
These results argue for tissue-compartment-aware meta-analysis in endometriosis transcriptomics. At minimum, studies should be stratified by the tissue comparison performed (ectopic vs. healthy, eutopic disease vs. eutopic control, paired eutopic–ectopic within patients) before gene lists are intersected or effect sizes combined. The permutation overlap test provides a lightweight quality-control step that can be applied to any pair of ranked lists to assess whether combining them is justified.
4.6 Scope and limitations
This audit examines three datasets on a single platform (Affymetrix U133 Plus 2.0). Generalization to RNA-seq datasets, other array platforms, or other tissue compartments (e.g., blood, peritoneal fluid) requires additional analysis. The variance filter (top 10,000 probes) and list-length threshold (N = 200) are analytical choices; sensitivity analysis (Table 5) suggests conclusions are robust across thresholds but does not exhaust the parameter space. The gene-level collapse uses first-symbol annotation and max-|t| selection, which may not capture all biologically relevant probe–gene relationships.
5. Conclusions
A permutation-based reproducibility audit of three major endometriosis microarray datasets reveals that tissue-type heterogeneity is the primary determinant of cross-study gene-list concordance. Tissue-matched datasets reproduce significantly (gene-level z = 5.96, p < 0.0002); cross-compartment comparisons do not. Gene-level analysis outperforms probe-level analysis for reproducibility assessment, uncovering a three-way intersection (MPDZ, PICALM, WNK1; z = 6.47, p = 0.001) that is invisible at the probe level. Within-dataset menstrual cycle effects exceed all cross-dataset disease overlaps, further underscoring the dominance of biological context over disease signal. These findings urge the endometriosis transcriptomics community to adopt tissue-compartment-aware study design in meta-analyses and to employ permutation-based overlap testing as a standard quality-control measure before pooling datasets.
References
Chapron C, Marcellin L, Borghese B, Santulli P. Rethinking mechanisms, diagnosis and management of endometriosis. Nat Rev Endocrinol. 2019;15(11):666-682.
Ghai V, Valladares-Galaviz D, Ghai S, et al. Endometriosis biomarkers: a systematic review. University of York. 2024.
Kalaitzopoulos DR, Samartzis N, Kolovos GN, et al. Treatment of endometriosis: a review with comparison of 8 guidelines. Exp Biol Med. 2020;245(5):437-447.
Devesa-Peiro A, Sebastian-Leon P, Pellicer A, Diaz-Gimeno P. Guidelines for biomarker discovery in endometrium: correcting for menstrual cycle bias reveals new genes associated with uterine disorders. Mol Hum Reprod. 2021;27(4):gaab011.
Grewal J, Guo L, Engstrom C, Bowman E, Ayers LW. Transcriptomic meta-analysis of endometriosis. bioRxiv. 2021. doi:10.1101/2021.05.21.445109.
Zhao H, Young BD, Bhatt S, et al. Gene expression profiling of ectopic and eutopic endometrial tissues in endometriosis. Reprod Biol Endocrinol. 2009;7:94.
Patil P, Parmigiani G. Training replicable predictors in multiple studies. Bioinformatics. 2015;31(14):2318-2323.
Hever A, Roth RB, Hevezi P, et al. Human endometriosis is associated with plasma cells and overexpression of B lymphocyte stimulator. Proc Natl Acad Sci USA. 2007;104(30):12451-12456. GEO: GSE7305.
Hull ML, Escareno CR, Godsland JM, et al. Endometrial–peritoneal interactions during endometriotic lesion establishment. Am J Pathol. 2008;173(3):700-715. GEO: GSE11691.
Tamaresis JS, Irwin JC, Goldfien GA, et al. Molecular classification of endometriosis and disease stage using high-dimensional genomic data. Endocrinology. 2014;155(12):4986-4999. GEO: GSE51981.
Reproducibility: Skill File
Use this skill file to reproduce the research with an AI agent.
---
name: endo-reproducibility-audit
description: >
Cross-dataset reproducibility audit of endometriosis diagnostic gene signatures.
Downloads three GPL570 Affymetrix datasets from GEO (GSE7305, GSE11691, GSE51981),
loads GPL570 probe-to-gene annotation, computes top differentially expressed genes
via Welch t-test with gene-level collapse, measures pairwise and three-way overlap,
and tests significance via 5,000-permutation label-shuffling null model.
Also assesses menstrual cycle phase confounding.
allowed-tools:
- Bash(python3 *)
- Bash(mkdir *)
- Bash(cat *)
- Bash(echo *)
---
# Endometriosis Cross-Dataset Reproducibility Audit
## Overview
This skill downloads three publicly available endometriosis microarray datasets
from NCBI GEO (all GPL570 Affymetrix HG-U133 Plus 2.0), obtains GPL570
probe-to-gene annotation, computes differential expression rankings at both
probe and gene levels, and systematically tests whether the overlap between
top-ranked gene lists exceeds what chance alone predicts using 5,000 permutations.
## Steps
1. Create the analysis script
2. Run the analysis
3. Report results
## Step 1: Create Analysis Script
```bash
mkdir -p endo_audit_results
cat > endo_audit_results/run_audit.py << 'ENDSCRIPT'
import gzip, math, os, random, statistics, urllib.request, json, time
from collections import defaultdict
random.seed(42)
OUTDIR = "endo_audit_results"
os.makedirs(OUTDIR, exist_ok=True)
DATASETS = {
"GSE7305": "https://ftp.ncbi.nlm.nih.gov/geo/series/GSE7nnn/GSE7305/matrix/GSE7305_series_matrix.txt.gz",
"GSE11691": "https://ftp.ncbi.nlm.nih.gov/geo/series/GSE11nnn/GSE11691/matrix/GSE11691_series_matrix.txt.gz",
"GSE51981": "https://ftp.ncbi.nlm.nih.gov/geo/series/GSE51nnn/GSE51981/matrix/GSE51981_series_matrix.txt.gz",
}
N_PERMS = 5000
TOP_PROBES = 10000
def download(url, label):
cache = os.path.join(OUTDIR, f"{label}_matrix.txt.gz")
if os.path.exists(cache):
with open(cache, "rb") as f:
return gzip.decompress(f.read()).decode("utf-8", errors="replace")
print(f" Downloading {label} ...")
req = urllib.request.Request(url, headers={"User-Agent": "Mozilla/5.0"})
data = urllib.request.urlopen(req, timeout=120).read()
with open(cache, "wb") as f:
f.write(data)
return gzip.decompress(data).decode("utf-8", errors="replace")
def load_probe_to_gene():
path = os.path.join(OUTDIR, "probe_to_gene.tsv")
if not os.path.exists(path):
print(" Downloading GPL570 annotation...")
url = "https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570&targ=self&form=text&view=data"
req = urllib.request.Request(url, headers={"User-Agent": "Mozilla/5.0"})
text = urllib.request.urlopen(req, timeout=300).read().decode("utf-8", errors="replace")
in_table = False
gene_col = None
with open(path, "w") as out:
out.write("probe_id\tgene_symbol\n")
for line in text.split("\n"):
line = line.rstrip("\r")
if "!platform_table_begin" in line:
in_table = True
continue
if "!platform_table_end" in line:
break
if not in_table:
continue
parts = line.split("\t")
if parts[0] == "ID":
for i, h in enumerate(parts):
if h == "Gene Symbol":
gene_col = i
break
continue
if gene_col is None:
continue
probe = parts[0]
gene = parts[gene_col] if gene_col < len(parts) else ""
gene = gene.strip()
if gene and gene != "---":
out.write(f"{probe}\t{gene}\n")
mapping = {}
with open(path) as f:
next(f)
for line in f:
parts = line.strip().split("\t")
if len(parts) >= 2:
probe, genes = parts[0], parts[1]
primary = genes.split(" /// ")[0].strip()
if primary:
mapping[probe] = primary
return mapping
print("=" * 70)
print("STEP 0 - Loading probe-to-gene annotation")
print("=" * 70)
probe2gene = load_probe_to_gene()
print(f" {len(probe2gene):,} probes mapped to gene symbols")
print("\n" + "=" * 70)
print("STEP 1 - Downloading GEO matrices")
print("=" * 70)
raw = {}
for gse, url in DATASETS.items():
raw[gse] = download(url, gse)
print(f" {gse}: {len(raw[gse]):,} chars")
def parse_matrix(text, gse):
lines = text.split("\n")
meta = {}
for line in lines:
if line.startswith("!"):
key = line.split("\t")[0]
vals = [v.strip().strip('"') for v in line.split("\t")[1:]]
meta.setdefault(key, []).append(vals)
in_data = False
expr = {}
sample_ids = []
for line in lines:
if line.startswith("!series_matrix_table_begin"):
in_data = True
continue
if line.startswith("!series_matrix_table_end"):
break
if not in_data:
continue
parts = line.split("\t")
if not parts:
continue
probe = parts[0].strip().strip('"')
if probe == "ID_REF":
sample_ids = [p.strip().strip('"') for p in parts[1:]]
continue
try:
vals = [float(v.strip().strip('"')) for v in parts[1:]]
except ValueError:
continue
if len(vals) == len(sample_ids):
expr[probe] = vals
n_samples = len(sample_ids)
labels = [""] * n_samples
phases = [""] * n_samples
if gse == "GSE7305":
titles = meta.get("!Sample_title", [[]])[0]
descs = meta.get("!Sample_description", [[]])[0]
for i, t in enumerate(titles):
labels[i] = "disease" if "Disease" in t else "control"
for i, d in enumerate(descs):
if "Follicular" in d: phases[i] = "Follicular"
elif "Luteal" in d: phases[i] = "Luteal"
elif gse == "GSE11691":
titles = meta.get("!Sample_title", [[]])[0]
for i, t in enumerate(titles):
labels[i] = "disease" if t.startswith("Endometriosis") else "control"
elif gse == "GSE51981":
sources = meta.get("!Sample_source_name_ch1", [[]])[0]
for i, s in enumerate(sources):
labels[i] = "disease" if s.startswith("Endometriosis") else "control"
chars_rows = meta.get("!Sample_characteristics_ch1", [])
if chars_rows:
for i, c in enumerate(chars_rows[0]):
if "Proliferative" in c: phases[i] = "Proliferative"
elif "Early Secretory" in c: phases[i] = "Early_Secretory"
elif "Mid-Secretory" in c: phases[i] = "Mid_Secretory"
elif "Late Secretory" in c: phases[i] = "Late_Secretory"
n_dis = sum(1 for l in labels if l == "disease")
n_ctl = sum(1 for l in labels if l == "control")
print(f" {gse}: {len(expr):,} probes x {n_samples} samples ({n_dis} disease, {n_ctl} control)")
return expr, sample_ids, labels, phases
print("\n" + "=" * 70)
print("STEP 2 - Parsing expression matrices")
print("=" * 70)
parsed = {}
for gse in DATASETS:
parsed[gse] = parse_matrix(raw[gse], gse)
common_probes = set(parsed["GSE7305"][0].keys())
for gse in ["GSE11691", "GSE51981"]:
common_probes &= set(parsed[gse][0].keys())
print(f"\n Common probes: {len(common_probes):,}")
def variance_filter(expr, top_k):
var_list = []
for probe, vals in expr.items():
if len(vals) < 2: continue
m = sum(vals) / len(vals)
v = sum((x - m) ** 2 for x in vals) / (len(vals) - 1)
var_list.append((probe, v))
var_list.sort(key=lambda x: x[1], reverse=True)
keep = set(p for p, _ in var_list[:top_k])
return {p: v for p, v in expr.items() if p in keep}
for gse in DATASETS:
expr_common = {p: v for p, v in parsed[gse][0].items() if p in common_probes}
expr_filt = variance_filter(expr_common, TOP_PROBES)
parsed[gse] = (expr_filt, parsed[gse][1], parsed[gse][2], parsed[gse][3])
print(f" {gse}: {len(expr_filt):,} probes after variance filter")
def welch_t_fast(vals, dis_idx, ctl_idx):
na, nb = len(dis_idx), len(ctl_idx)
if na < 2 or nb < 2: return 0.0
sa = sb = ssa = ssb = 0.0
for i in dis_idx:
v = vals[i]; sa += v; ssa += v * v
for i in ctl_idx:
v = vals[i]; sb += v; ssb += v * v
ma, mb = sa / na, sb / nb
va = (ssa - sa * sa / na) / (na - 1)
vb = (ssb - sb * sb / nb) / (nb - 1)
se2 = va / na + vb / nb
if se2 <= 0: return 0.0
return (ma - mb) / math.sqrt(se2)
def compute_deg_ranking(expr, labels, indices=None):
if indices is None: indices = list(range(len(labels)))
dis_idx = [i for i in indices if labels[i] == "disease"]
ctl_idx = [i for i in indices if labels[i] == "control"]
results = []
for probe, vals in expr.items():
t = welch_t_fast(vals, dis_idx, ctl_idx)
results.append((probe, t))
results.sort(key=lambda x: abs(x[1]), reverse=True)
return results
def collapse_to_genes(ranking, p2g):
seen_genes = set()
gene_ranking = []
for probe, t in ranking:
gene = p2g.get(probe)
if gene is None:
continue
if gene in seen_genes:
continue
seen_genes.add(gene)
gene_ranking.append((gene, t))
return gene_ranking
print("\n" + "=" * 70)
print("STEP 3 - Computing DE rankings (probe + gene level)")
print("=" * 70)
rankings_probe = {}
rankings_gene = {}
for gse in DATASETS:
expr, sids, labels, phases = parsed[gse]
rankings_probe[gse] = compute_deg_ranking(expr, labels)
rankings_gene[gse] = collapse_to_genes(rankings_probe[gse], probe2gene)
print(f" {gse}: {len(rankings_probe[gse])} probes, {len(rankings_gene[gse])} genes")
def top_n_set(ranking, n):
return set(r[0] for r in ranking[:n])
def jaccard(s1, s2):
if not s1 or not s2: return 0.0
return len(s1 & s2) / len(s1 | s2)
gse_list = list(DATASETS.keys())
print("\n" + "=" * 70)
print("STEP 4 - Cross-dataset overlap")
print("=" * 70)
for level_name, rankings in [("probe", rankings_probe), ("gene", rankings_gene)]:
print(f"\n --- {level_name}-level ---")
for N in [50, 100, 200, 500]:
sets = {gse: top_n_set(rankings[gse], N) for gse in gse_list}
print(f" N={N}:")
for i in range(len(gse_list)):
for j in range(i+1, len(gse_list)):
a, b = gse_list[i], gse_list[j]
inter = len(sets[a] & sets[b])
jac = jaccard(sets[a], sets[b])
print(f" {a} vs {b}: {inter} (Jaccard={jac:.4f})")
tw = sets[gse_list[0]] & sets[gse_list[1]] & sets[gse_list[2]]
print(f" Three-way: {len(tw)}")
if tw: print(f" Three-way genes: {sorted(tw)}")
TEST_N = 200
print("\n" + "=" * 70)
print(f"STEP 5 - Permutation test (N={TEST_N}, {N_PERMS} perms)")
print("=" * 70)
dataset_arrays = {}
for gse in gse_list:
expr = parsed[gse][0]
labels = parsed[gse][2]
probes = list(expr.keys())
vals_matrix = [expr[p] for p in probes]
dataset_arrays[gse] = (probes, vals_matrix, labels)
def perm_top_n_both(probes, vals_matrix, labels, n, p2g):
shuf = labels[:]
random.shuffle(shuf)
dis_idx = [i for i, l in enumerate(shuf) if l == "disease"]
ctl_idx = [i for i, l in enumerate(shuf) if l == "control"]
t_list = []
for idx, vals in enumerate(vals_matrix):
t = welch_t_fast(vals, dis_idx, ctl_idx)
t_list.append((idx, abs(t)))
t_list.sort(key=lambda x: x[1], reverse=True)
probe_set = set(probes[t_list[k][0]] for k in range(n))
seen_genes = set()
gene_set = set()
for idx_rank, (orig_idx, _) in enumerate(t_list):
probe = probes[orig_idx]
gene = p2g.get(probe)
if gene is None or gene in seen_genes:
continue
seen_genes.add(gene)
gene_set.add(gene)
if len(gene_set) >= n:
break
return probe_set, gene_set
obs_probe = {gse: top_n_set(rankings_probe[gse], TEST_N) for gse in gse_list}
obs_gene = {gse: top_n_set(rankings_gene[gse], TEST_N) for gse in gse_list}
perm_results = {}
for level, obs in [("probe", obs_probe), ("gene", obs_gene)]:
obs_pairs = {}
for i in range(len(gse_list)):
for j in range(i+1, len(gse_list)):
a, b = gse_list[i], gse_list[j]
obs_pairs[(a,b)] = len(obs[a] & obs[b])
obs_three = len(obs[gse_list[0]] & obs[gse_list[1]] & obs[gse_list[2]])
perm_results[level] = {"obs_pairs": obs_pairs, "obs_three": obs_three}
null_probe_pairs = {k: [] for k in perm_results["probe"]["obs_pairs"]}
null_gene_pairs = {k: [] for k in perm_results["gene"]["obs_pairs"]}
null_probe_three = []
null_gene_three = []
t0 = time.time()
print(f" Running {N_PERMS} permutations ...")
for pi in range(N_PERMS):
if (pi+1) % 500 == 0:
elapsed = time.time() - t0
rate = (pi+1) / elapsed
eta = (N_PERMS - pi - 1) / rate
print(f" {pi+1}/{N_PERMS} ({elapsed:.0f}s elapsed, ~{eta:.0f}s remaining)")
ps_probe = {}
ps_gene = {}
for gse in gse_list:
probes, vm, labels = dataset_arrays[gse]
p_set, g_set = perm_top_n_both(probes, vm, labels, TEST_N, probe2gene)
ps_probe[gse] = p_set
ps_gene[gse] = g_set
for i in range(len(gse_list)):
for j in range(i+1, len(gse_list)):
a, b = gse_list[i], gse_list[j]
null_probe_pairs[(a,b)].append(len(ps_probe[a] & ps_probe[b]))
null_gene_pairs[(a,b)].append(len(ps_gene[a] & ps_gene[b]))
null_probe_three.append(len(ps_probe[gse_list[0]] & ps_probe[gse_list[1]] & ps_probe[gse_list[2]]))
null_gene_three.append(len(ps_gene[gse_list[0]] & ps_gene[gse_list[1]] & ps_gene[gse_list[2]]))
print(f" Completed in {time.time()-t0:.0f}s")
for level, obs_data, null_pairs, null_three in [
("probe", perm_results["probe"], null_probe_pairs, null_probe_three),
("gene", perm_results["gene"], null_gene_pairs, null_gene_three)]:
print(f"\n --- {level}-level results ---")
for (a,b), obs in obs_data["obs_pairs"].items():
nulls = null_pairs[(a,b)]
p = sum(1 for n in nulls if n >= obs) / N_PERMS
mn = statistics.mean(nulls)
sd = statistics.stdev(nulls) if len(nulls) > 1 else 0
z = (obs - mn) / sd if sd > 0 else float("inf")
print(f" {a} vs {b}: obs={obs}, null={mn:.1f}+/-{sd:.1f}, z={z:.2f}, p={p:.4f}")
obs_t = obs_data["obs_three"]
p3 = sum(1 for n in null_three if n >= obs_t) / N_PERMS
mn3 = statistics.mean(null_three)
sd3 = statistics.stdev(null_three) if len(null_three) > 1 else 0
z3 = (obs_t - mn3) / sd3 if sd3 > 0 else float("inf")
print(f" Three-way: obs={obs_t}, null={mn3:.2f}+/-{sd3:.2f}, z={z3:.2f}, p={p3:.4f}")
print("\n" + "=" * 70)
print("STEP 6 - Menstrual cycle stratification")
print("=" * 70)
expr51, _, labels51, phases51 = parsed["GSE51981"]
for pg, tags in [("Proliferative", ["Proliferative"]),
("Secretory", ["Early_Secretory", "Mid_Secretory", "Late_Secretory"])]:
idx = [i for i in range(len(labels51)) if phases51[i] in tags]
nd = sum(1 for i in idx if labels51[i] == "disease")
nc = sum(1 for i in idx if labels51[i] == "control")
print(f" GSE51981 {pg}: {nd} disease, {nc} control")
prolif_idx = [i for i in range(len(labels51)) if phases51[i] == "Proliferative"]
secr_idx = [i for i in range(len(labels51)) if phases51[i] in ["Early_Secretory", "Mid_Secretory", "Late_Secretory"]]
rank_prolif = compute_deg_ranking(expr51, labels51, prolif_idx)
rank_secr = compute_deg_ranking(expr51, labels51, secr_idx)
rank_prolif_gene = collapse_to_genes(rank_prolif, probe2gene)
rank_secr_gene = collapse_to_genes(rank_secr, probe2gene)
set_prolif_g = top_n_set(rank_prolif_gene, 200)
set_secr_g = top_n_set(rank_secr_gene, 200)
print(f" Within-dataset phase (gene, N=200): {len(set_prolif_g & set_secr_g)} overlap (Jaccard={jaccard(set_prolif_g, set_secr_g):.4f})")
print("\n" + "=" * 70)
print("STEP 7 - Sensitivity analysis")
print("=" * 70)
for level_name, rankings in [("probe", rankings_probe), ("gene", rankings_gene)]:
print(f"\n --- {level_name}-level ---")
for N in [25, 50, 100, 200, 500, 750, 1000]:
sn = {gse: top_n_set(rankings[gse], N) for gse in gse_list}
pairs = []
for i in range(len(gse_list)):
for j in range(i+1, len(gse_list)):
pairs.append(jaccard(sn[gse_list[i]], sn[gse_list[j]]))
tw = len(sn[gse_list[0]] & sn[gse_list[1]] & sn[gse_list[2]])
print(f" N={N:5d} mean_Jaccard={statistics.mean(pairs):.4f} three_way={tw}")
print("\n" + "=" * 70)
print("ANALYSIS COMPLETE")
print("=" * 70)
ENDSCRIPT
```
## Step 2: Run Analysis
```bash
python3 endo_audit_results/run_audit.py
```
## Step 3: Report Results
```bash
echo "Analysis complete. Results printed above."
```
Discussion (0)
to join the discussion.
No comments yet. Be the first to discuss this paper.