ccRCC Network Rewiring

Statistical methodology and research implementation for detecting changes in gene co-expression structure in clear-cell renal cell carcinoma (ccRCC).

The repository develops a two-condition network-rewiring test against synthetic planted truth, applies it to TCGA KIRC bulk RNA-seq, explores malignant-versus-proximal-tubule structure in GSE159115 single-cell RNA-seq, and follows selected signals in CPTAC proteomic data and TCGA survival data.

Status as of June 26, 2026: remediation pass complete. The ten known violations below have been addressed (see METHODS_DECISIONS.md, ANALYSIS_PLAN.md, PROVENANCE.md). Single-cell inference is now patient-level (scrna_patient_level.py), TCGA samples are de-duplicated at source, the TOST margin is per-assay standardized, equivalence is multiplicity-reported, an external CPTAC replication and proportional-hazards diagnostics were added, and downstream artifacts were regenerated from one method commit. Key scientific update: at the patient level the glycolysis-coordination transition is robust (7/7 patients, exact p=0.008) but PT-identity coordination is not (1/7, p=0.19); PT-identity loss remains strongly prognostic as a bulk expression module (HR 0.64) and survives adjustment for FBP1 + HIF/glycolysis, while a per-patient network-coordination score adds no prognostic value over mean expression.

Sign audit update, June 30, 2026: module_axes.py found no scoring or module-definition inversion bug. PT-identity is scored in the identity-positive direction, and the positive HIF/glycolysis/angiogenesis correlations persist when FBP1 is removed from the PT score. The old shorthand "PT-low = HIF/glycolytic" is therefore rejected for bulk module abundance. The defensible claim is that PT-high tumors are differentiated/metabolically expressive and angiogenic in TCGA bulk data, glycolysis is mostly a global metabolic-amplitude coupling, and the robust PT-low pole is myeloid/stromal/EMT.

Scientific Questions

The project asks:

1. Can differential co-expression be detected with calibrated error rates in the presence of differential expression and cell-composition confounding?
2. Which genes change network neighborhoods across ccRCC grade, stage, tumor/normal status, or driver-mutation status?
3. Which rewired genes also have evidence that their mean expression is practically equivalent between conditions?
4. Do candidate network states reproduce across contrasts, patients, omic layers, or clinical outcomes?

The term silent shifter is used here for a gene that satisfies both a rewiring criterion and a positive expression-equivalence criterion. It does not mean merely "not statistically significant for differential expression."

Statistical Target

Let:

• X be an n x G expression matrix;
• y_i in {0,1} be the condition label for sample i;
• C^(0) and C^(1) be condition-specific gene-correlation matrices; and
• A be a shared binary network skeleton.

For gene g, the primary estimand is the total change in its incident co-expression relationships over the shared skeleton:

$$R_g = \sum_{h=1}^{G} A_{gh} \left|\operatorname{atanh}(C^{(0)}_{gh}) - \operatorname{atanh}(C^{(1)}_{gh})\right|.$$

This is a node-level summary of differential correlation. It measures the magnitude of neighborhood change, not its direction and not a causal regulatory effect.

The implemented decision rule is:

$$\text{silent-shifter}_g = I(q^{RW}_g \le 0.10)\,I(p^{EQ}_g \le 0.05),$$

where q_RW is the BH-adjusted rewiring p-value and p_EQ is the TOST equivalence p-value for the gene's mean-expression difference.

Analysis Overview

expression matrix + two-group labels
        |
        v
composition / latent-factor adjustment
        |
        v
pooled top-k correlation skeleton
        |
        v
condition-specific correlations and Fisher-z edge changes
        |
        v
per-gene edge-sum rewiring statistic
        |
        +--> label-permutation null --> pooled standardized p-values --> BH q
        |
        +--> Welch estimates on input expression --> TOST equivalence p-value
        |
        v
rewired + expression-equivalent = silent-shifter call

The implementation is in methodology_benchmark/methodology.py.

Input Data and Preprocessing

TCGA KIRC bulk RNA-seq

The bulk pipeline uses UCSC Xena TCGA KIRC STAR TPM and clinical tables.

Processing in prep_ccrcc.py:

1. Ensembl version suffixes are removed.
2. Duplicate Ensembl identifiers retain the first row.
3. Samples are classified as adjacent normal when the TCGA sample code is 11; all other retained samples are classified as tumor.
4. AJCC stage is read from the GDC clinical table.
5. Histologic grade is read from the legacy Xena clinical matrix.
6. Genes with mean log2-TPM below 1.0 are removed.
7. The 1,200 remaining genes with greatest pooled variance are selected.

The current artifact contains 538 tumor and 72 adjacent-normal samples. Five tumor rows are repeat samples from patients already represented. Current bulk contrast scripts treat rows as independent and do not model patient clustering.

GSE159115 single-cell RNA-seq

Processing in prep_scrna.py:

1. Tumor-annotated ccRCC cells define the malignant population.
2. PT-A cells from normal kidney samples define the proximal-tubule reference.
3. Counts are library-size normalized to 10,000 counts per cell and transformed with log(1+x).
4. Genes must be detected in at least 10% of cells in both populations.
5. The 1,500 eligible genes with greatest pooled cell-level variance are saved.
6. Downstream rewiring scripts use the first 800 genes from this ordered panel.

The current artifact contains 9,139 malignant cells from seven tumor patients and 1,941 normal proximal-tubule cells. Three normal donors contribute almost all normal cells; a fourth contributes three cells.

Before network analysis, the combined cell matrix is residualized on:

• an intercept;
• standardized log(1 + total_UMI);
• standardized log(1 + detected_gene_count);
• standardized mitochondrial percentage; and
• standardized mean ribosomal-protein expression.

If T is the technical design matrix, the residualized expression is:

$$X_{res} = X - T(T^T T)^{-1}T^T X,$$

implemented numerically with least squares.

Core Rewiring Method

1. Composition and latent-factor adjustment

The current bulk configuration uses pooled principal-component removal with k = 3.

For centered expression X_c, let its singular value decomposition be:

$$X_c = U S V^T.$$

The adjusted matrix is:

$$X_{adj} = X_c - U_k U_k^T X_c,$$

where U_k contains the first three left singular vectors. The projection is estimated without using condition labels.

This is a generic latent-factor adjustment. It does not identify those components specifically as tumor purity or cell composition, and it can remove biological signal as well as confounding signal.

The code also contains experimental alternatives:

• none: pooled mean centering only;
• proportions: estimate nonnegative cell-type proportions from a reference, then regress centered proportions from expression;
• deconvolve: subtract estimated non-target mixtures and divide by estimated target proportion; and
• bmind: a regularized target-compartment approximation that shrinks the recovered residual toward a reference profile.

These alternatives are available for experiments but are not the current bulk analysis configuration. Single-cell analyses set adjustment to none after technical residualization.

2. Correlation estimation

Each gene is centered and divided by its population standard deviation (ddof=0). For a group matrix with n_c rows:

$$C^{(c)} = \frac{Z_c^T Z_c}{n_c}.$$

This is numerically the Pearson correlation matrix under the same population standardization. Correlations are clipped to [-1,1].

3. Shared network skeleton

A pooled correlation matrix is calculated from all adjusted samples without using group labels. For each gene, the 120 largest absolute correlations are retained. The directed neighbor selections are symmetrized by union:

$$A_{gh}=1 \quad\text{if gene }h\text{ selects }g \text{ or gene }g\text{ selects }h.$$

The diagonal is excluded. Because of union symmetrization, a gene may have more than 120 final neighbors.

The same skeleton is held fixed for the observed statistic and every label permutation. This avoids rebuilding a label-dependent feature-selection step inside only the observed contrast.

4. Edge-sum rewiring statistic

Condition-specific correlations are Fisher transformed:

$$Z^{(c)}_{gh} = \operatorname{atanh} \left(\operatorname{clip}(C^{(c)}_{gh},-0.999999,0.999999)\right).$$

The edge change is:

$$D_{gh}=|Z^{(0)}_{gh}-Z^{(1)}_{gh}|,$$

and the gene-level statistic is:

$$R_g=\sum_h A_{gh}D_{gh}.$$

The absolute value makes the test two-sided with respect to correlation change: correlation gains, losses, and sign reversals all increase the statistic.

5. Label-permutation null

The current method performs B = 150 random permutations of the condition labels. Group sizes are preserved. The adjusted matrix and pooled skeleton are held fixed.

For each gene, the null mean and standard deviation are calculated across permutations:

$$\mu_g = \frac{1}{B}\sum_b R_{gb}^{perm}, \qquad s_g = \operatorname{sd}_b(R_{gb}^{perm}).$$

The observed and permuted statistics are standardized gene by gene:

$$Z_g^{obs}=\frac{R_g^{obs}-\mu_g}{s_g}, \qquad Z_{gb}^{perm}=\frac{R_{gb}^{perm}-\mu_g}{s_g}.$$

All B x G standardized permutation statistics are pooled into one empirical null distribution. The rewiring p-value is:

$$p_g^{RW}= \frac{1 + \#\{(h,b):Z_{hb}^{perm}\ge Z_g^{obs}\}} {1 + BG}.$$

Pooling gives much finer numerical p-value resolution than 150 gene-specific permutations. It also assumes that standardized null statistics are sufficiently exchangeable across genes. That assumption has been evaluated only in the repository's simulator, not established for real transcriptomic data.

6. Multiple-testing correction

Rewiring p-values are corrected across all genes in the analyzed panel using the Benjamini-Hochberg procedure. A gene is called rewired when:

$$q_g^{RW} \le 0.10.$$

The reported threshold is an FDR target within the selected panel, conditional on the validity of the pooled empirical p-values. It is not a genome-wide error guarantee.

7. Differential-expression summary

For context, the method also calculates a two-sided Welch test on the expression matrix passed to run_methodology.

For groups 0 and 1:

$$d_g=\bar X_{0g}-\bar X_{1g}, \qquad SE_g=\sqrt{\frac{s_{0g}^2}{n_0}+\frac{s_{1g}^2}{n_1}},$$

with Welch-Satterthwaite degrees of freedom. These p-values receive BH correction and are returned as de_qvalue, but the current silent-shifter call does not use de_qvalue.

8. Positive expression-equivalence test

The method uses a two-one-sided test (TOST) with fixed margin Delta = 0.30 on the input expression scale.

The hypotheses are:

$$H_{01}: d_g \le -\Delta, \qquad H_{02}: d_g \ge +\Delta.$$

The implemented statistics are:

$$t_{low}=\frac{d_g+\Delta}{SE_g}, \qquad t_{high}=\frac{d_g-\Delta}{SE_g},$$

and:

$$p_g^{EQ}=\max\left[ P(T_{df}\ge t_{low}), P(T_{df}\le t_{high}) \right].$$

Expression equivalence is declared at p_EQ <= 0.05. These equivalence p-values are not multiplicity adjusted.

The numeric margin is applied to different scales in different analyses: log2-TPM for TCGA and technically residualized log-normalized expression for single-cell data. It has not yet been biologically calibrated separately for those assays.

9. Silent-shifter decision

A gene is called a silent shifter only when both conditions hold:

rewiring BH q-value <= 0.10
TOST equivalence p-value <= 0.05

The method therefore requires positive evidence of practical mean equivalence. Failure to detect differential expression is not sufficient.

Alternative Rewiring Statistics

The code supports two non-default alternatives.

LIONESS-style node strength

For each sample i, a leave-one-out correlation matrix is used to estimate:

$$Q_i = |nC_{all}-(n-1)C_{-i}|.$$

Values outside the shared skeleton are set to zero and incident edge weights are summed per gene. The rewiring statistic is the absolute difference in mean node strength between groups.

Spectral displacement

For each condition, the absolute skeleton-weighted correlation network is degree normalized. The leading eigenvectors, scaled by square roots of absolute eigenvalues, define a node embedding. One condition is orthogonally aligned to the other by Procrustes rotation, and per-gene Euclidean displacement is used as the statistic.

Neither alternative is used for the current tracked bulk result tables.

Synthetic Ground-Truth Benchmark

The benchmark is intended to prevent method development from optimizing real ccRCC discovery counts or p-values.

Simulator

simulate.py generates pseudobulk samples from four compartments, with compartment 0 representing malignant cells:

$$X_i=\sum_{c=1}^{4}p_{ic}X_{ic}.$$

Cell fractions are drawn from condition-specific Dirichlet distributions. The mean malignant fraction is 0.62 in condition A and 0.45 in condition B, creating a composition shift.

The default simulation contains:

Component	Default
Genes	600
Cell types	4
Background malignant modules	20
Rewired modules	2 x 30 genes
DE-confounder genes	40
Rewiring loading	0.65
Background loading	0.60
Cell-type marker fraction	0.55
Marker mean boost	2.0

One rewired module is correlated in condition A and the other in condition B. DE-confounder genes receive a malignant-compartment mean shift in condition B. Non-malignant marker genes serve as composition confounders. Null simulations remove planted rewiring while retaining the other structures.

Evaluation panel

benchmark.py evaluates signal and null simulations at 50, 100, and 150 samples per condition, using three fixed seeds at each sample size. A noisy cell-type reference is supplied to methods that use one.

The headline score is:

$$\text{benchmark score}=0.5\times AUROC+0.5\times sensitivity.$$

The run is accepted only if all gates pass:

• null-background call rate <= 0.10;
• DE-confounder call rate <= 0.10; and
• composition-confounder call rate <= 0.10.

Current audited execution:

Metric	Value
Benchmark score	0.5684
AUROC	0.8035
Sensitivity at q <= 0.10	0.3333
Null-background call rate	0.0230
DE-confounder call rate	0.0056
Composition-confounder call rate	0.0273

All three gates pass. This is performance on the fixed development simulator and fixed seeds. It is not an untouched external simulation validation.

Real-Data Study Designs

TCGA aggressiveness contrasts

analyze_contrast.py uses tumor samples for three aggressiveness comparisons:

Contrast	Group 0	Group 1
Grade	G1/G2	G3/G4
Stage	Stage I/II	Stage III/IV
Concordant extremes	low grade and early stage	high grade and late stage

Adjacent normal versus tumor is analyzed as a positive-control contrast. The minimum permitted group size is 10. Labels are permuted without covariate or patient stratification.

Current outputs:

Contrast	Rewired	Silent shifters
Grade	200 / 1,200	36
Stage	323 / 1,200	40
Concordant extremes	250 / 1,200	13
Tumor versus normal	1,099 / 1,200	8

Only FPR3, SLC16A3, and SNORD17 satisfy the current silent-shifter rule in all three aggressiveness contrasts. Fifteen genes occur in at least two.

Driver-mutation contrasts

analyze_mutation_contrast.py defines a patient as mutant for a driver if at least one mutation in that gene has an effect outside a hard-coded silent/noncoding set. Wild-type eligibility is restricted to patients represented in the mutation table.

Current outputs:

Driver	Mutant	Wild-type	Rewired	Silent shifters
VHL	174	203	0	0
PBRM1	155	222	123	2
BAP1	41	336	0	0
SETD2	46	331	0	0

The matched-driver cache and mutation_specificity.csv were produced under an older method state and are not synchronized with these result tables.

Single-cell malignant-transformation contrast

analyze_scrna.py compares normal proximal-tubule cells with malignant cells after technical residualization. The pooled malignant group is capped at 3,000 randomly selected cells. Separate runs compare each tumor patient's malignant cells against the same pooled normal reference. Patients with fewer than 100 malignant cells are skipped.

The current pooled run calls 363 of 800 genes silent shifters. Fifty-six are called in all seven tumor-patient comparisons, and the neighborhood analysis classifies 798 of 800 genes as rewired.

These numbers are not confirmatory gene-level discoveries because:

• condition labels are permuted across cells, not patients;
• cells from one patient are correlated;
• tumor/normal status is confounded with donor and sample source;
• each patient run reuses the same normal-cell pool; and
• patient recurrence counts therefore do not represent independent replication.

The single-cell outputs should be interpreted as descriptive network prioritization and pathway-level hypothesis generation.

Neighborhood and Module Analyses

Single-cell neighborhoods

scrna_neighborhood.py computes the 25 largest absolute-correlation partners of selected PT-identity, HIF/hypoxia, and glycolysis genes in normal and malignant networks. Gained and lost neighbor sets are annotated against manually curated programs.

Glycolysis coordination is summarized as mean absolute pairwise correlation among glycolysis genes:

$$M(C,S)=\frac{1}{|E_S|}\sum_{g<h,\;g,h\in S}|C_{gh}|.$$

The same quantity is calculated for each tumor patient's malignant cells and compared descriptively with the pooled normal baseline. No patient-level test or confidence interval is calculated.

TCGA module scores

bulk_module_test.py converts curated hypoxia/glycolysis and PT-identity gene sets into abundance scores. For module S, each gene is standardized across tumors and the score is:

$$score_i(S)=\frac{1}{|S|}\sum_{g\in S} \frac{X_{ig}-\bar X_g}{s_g}.$$

A purity proxy is the average of standardized immune and stromal module scores. Grade and stage are tested with Welch tests before and after linear residualization of the module score on this proxy. Driver associations are also tested with Welch tests.

These analyses test module abundance, not patient-specific network rewiring. The collection of module, endpoint, and driver tests does not currently have a family-wide multiplicity correction.

module_axes.py applies the shared modules.py definitions to the pre-declared Aim-1 module-axis family. The sign audit in that report rejects the original expectation that PT-identity should anti-correlate with HIF/glycolysis/angiogenesis in bulk. The adjusted TCGA correlations are positive for HIF, glycolysis, and angiogenesis; removing FBP1 from the PT score leaves the signs positive. After adding a global tumor-metabolic-amplitude covariate, glycolysis is near zero/slightly negative, whereas HIF and angiogenesis remain PT-high. The project should therefore state PT-low as the myeloid/stromal/EMT pole, not as a HIF-glycolytic pole.

Survival Methodology

survival_analysis.py evaluates a legacy 60-gene aggressiveness set that was frozen before survival outcomes were used. This set predates the current PCA-adjusted TOST results.

For tumor expression standardized across samples, the implemented per-sample node score is:

$$S_{ig}=|Z_{ig}|\sum_h A_{gh}|Z_{ih}|.$$

Each selected gene's node score is standardized across patients, then averaged to form the patient network score.

Overall survival analyses include:

• Kaplan-Meier curves and log-rank test at the median score;
• univariate Cox regression per score standard deviation;
• multivariable Cox regression with age, sex, binary grade, and binary stage;
• likelihood-ratio comparison of clinical and clinical-plus-score models;
• apparent concordance indices; and
• 300 bootstrap resamples comparing clinical and clinical-plus-score models on out-of-bag patients.

The legacy score is null univariately and has an out-of-bag mean C-index gain of +0.0009, with percentile interval [-0.0160,+0.0106]. The nominal adjusted Cox coefficient should not be interpreted as validated incremental prognosis. Proportional-hazards assumptions are not currently tested.

CPTAC Cross-Omic Methodology

Legacy silent-shifter validation

cptac_validate.py tests low-grade versus high-grade CPTAC tumors separately for protein features and phosphosites:

1. A Welch test is calculated for each feature with at least five observations per grade group.
2. BH correction is applied at the feature level.
3. Multiple features mapping to one gene are collapsed to the minimum q-value.
4. Genes with minimum q <= 0.10 are called grade differential.
5. A one-sided Fisher exact test asks whether legacy TCGA silent shifters are enriched among grade-differential genes within the measured panel.

Taking the minimum adjusted feature value per gene is a pragmatic aggregation, not a formal site-to-gene hierarchical testing procedure.

The tracked analysis uses the older 60-gene set and found no enrichment for protein or phosphosite changes.

PT-identity protein module

cptac_pt_validate.py evaluates 11 PT-identity genes.

• Protein abundance is z-scored feature by feature.
• Module abundance is the mean protein z-score.
• Protein coherence is mean pairwise correlation among module proteins.
• Grade association uses a Welch test of the module score.
• Individual protein and phosphosite tests use nominal Welch p-values.

The current report gives mean protein coherence 0.395 and a nonsignificant protein module grade association (p = 0.337). Missing protein z-scores are filled with zero for the coherence calculation. No uncertainty interval is reported for coherence, and individual feature tests are not multiplicity corrected.

Statistical Assumptions and Limitations

Core-method assumptions

• Samples are exchangeable under label permutation within each contrast.
• Rows are independent experimental units.
• The pooled skeleton is an acceptable fixed feature-selection step.
• Standardized null statistics are exchangeable enough to pool across genes.
• BH behavior is adequate under dependence among gene-level p-values.
• Three pooled principal components capture unwanted bulk variation without removing the network signal of interest.
• A margin of 0.30 represents practical expression equivalence on the relevant assay scale.

Known violations — resolution log (remediation pass, 2026-06-26)

1. Single-cell row independence / exchangeability — RESOLVED. scrna_patient_level.py uses the patient/donor as the unit (exact Mann-Whitney on per-unit module coordination; within-donor paired deltas for SS_2017/SS_2022). Rationale: Squair et al., Nat Commun 2021.
2. Five duplicate TCGA tumor rows — RESOLVED. prep_ccrcc.py collapses to one primary-tumor sample per patient (538→533); headline Cox unchanged (HR 0.644).
3. Shared equivalence margin across assays — RESOLVED. EQUIV_MARGIN_MODE=standardized: margin = 0.30 x per-assay pooled SD, computed separately for bulk and single-cell.
4. Equivalence p-values uncorrected — RESOLVED. equiv_qvalue (BH) computed and reported in every result; primary call uses raw TOST (BH-on-TOST is margin-fragile), corrected counts reported alongside.
5. Panel-selection conditioning — ADDRESSED. Confirmatory tests use a-priori gene sets (modules.py); scrna_panel_split.py shows the glycolysis result and the genome-wide rewired fraction are stable under a selection/inference cell-split.
6. Benchmark fixed seeds inspected during selection — PARTIALLY ADDRESSED. Gates re-verified after every change; an untouched-seed external simulation remains future work.
7. Extreme single-cell rewiring fraction (miscalibrated null) — RESOLVED with (1): the cell-level 798/800 fraction falls to a calibrated 491/800 at the patient level (limited by 3 normal donors; reported descriptively).
8. Stale downstream artifacts — RESOLVED. Contrasts re-run, aggregate_aggressiveness.py rebuilds the recurrence table, survival_full.py supersedes the legacy survival script, provenance in PROVENANCE.md.
9. No family-wide multiplicity plan — RESOLVED. ANALYSIS_PLAN.md declares families; module_axes.py applies BH across the Aim-1 family.
10. No PH diagnostics / external validation — RESOLVED. Schoenfeld PH tests on every Cox (caught a stage PH violation; stage-stratified sensitivity confirms PT-identity HR 0.66); external CPTAC replication in cptac_external.py.

Artifact Version Warning

The current contrast files were regenerated after adding three-component PCA adjustment and TOST equivalence. However:

• silent_shifters_aggressiveness.csv still contains the older 60-gene set;
• survival_report.txt and cptac_crossomic_validation.csv validate that older set;
• mutation_specificity.csv and mutation caches reflect an older method state; and
• narrative reports mix historical and current results.

These artifacts must be regenerated from one method commit before they are used as a single inferential chain.

Reproducibility

There is no pinned ccRCC environment or lockfile yet. The main dependencies are:

numpy
pandas
scipy
matplotlib
lifelines
h5py
requests
cptac

Example environment:

python3.12 -m venv .venv
source .venv/bin/activate
python -m pip install numpy pandas scipy matplotlib lifelines h5py requests cptac

The system python3 on the audited machine points to Python 3.14 without SciPy. The benchmark was therefore verified in an isolated Python 3.12 environment.

Data downloads are not version-pinned or checksum-verified. Record data-source versions, hashes, method commit, parameter values, and artifact timestamps for any release run.

Running the Analyses

Synthetic benchmark

python methodology_benchmark/benchmark.py
python methodology_benchmark/run_benchmark.py

TCGA and CPTAC

python ccrcc_data/fetch_ccrcc_data.py
python ccrcc_analysis/prep_ccrcc.py

python ccrcc_analysis/analyze_contrast.py grade
python ccrcc_analysis/analyze_contrast.py stage
python ccrcc_analysis/analyze_contrast.py concordant
python ccrcc_analysis/analyze_contrast.py tumor_normal
python ccrcc_analysis/analyze_mutation_contrast.py

There is currently no script that rebuilds the aggressiveness recurrence table from the three current contrast outputs. That aggregation must be added before rerunning survival or legacy CPTAC validation.

GSE159115

python ccrcc_data/fetch_scrna_ccrcc.py GSE159115

Extract the 10x HDF5 files into ccrcc_data/scrna/GSE159115/h5/, then run:

python ccrcc_analysis/prep_scrna.py
python ccrcc_analysis/analyze_scrna.py
python ccrcc_analysis/scrna_neighborhood.py
python ccrcc_analysis/bulk_module_test.py
python ccrcc_analysis/cptac_pt_validate.py

analyze_scrna.py computes at most two missing cached runs per invocation and must be rerun until all pooled and per-patient caches exist.

Repository Layout

methodology_benchmark/      simulator, rewiring method, benchmark, scoring
ccrcc_analysis/             bulk, mutation, single-cell, survival, CPTAC scripts
ccrcc_analysis/artifacts/   tracked result tables, reports, and figures
ccrcc_data/                 data download scripts; large local cohorts are ignored
writeup/                    historical scope report and study-design memo
cancer_texts/               local literature PDFs; ignored by Git

Required Work Before Scientific Release

1. Define one versioned analysis specification and regenerate every downstream artifact from it.
2. Collapse repeated TCGA samples or model patient clustering.
3. Replace cell-level single-cell inference with patient-level pseudobulk, mixed-effects, or patient-stratified methods.
4. Calibrate equivalence margins separately by assay and define multiplicity families in advance.
5. Validate the pooled null and FDR behavior under new simulation families, untouched seeds, non-Gaussian distributions, unequal variances, batch effects, and realistic correlation structures.
6. Add an external ccRCC cohort for network and clinical replication.
7. Add a pinned environment, automated tests, CI, input hashes, and artifact manifests.

The present implementation is best treated as a transparent methodology prototype with a useful synthetic development harness, credible descriptive module findings, and exploratory real-data network hypotheses.