feat(jump_rr): add compound_no_source7 dataset for production paper alignment

[shntnu]

Context

The JUMP production paper (jump-cellpainting/2024_Chandrasekaran_Production) uses compound profiles that exclude Source 7 (which made batch correction difficult). These profiles live in cpg0042 on S3, while jump_rr currently uses cpg0016 (all sources).

Per discussion in broadinstitute/jump_hub#91 (comment), we need to add compound_no_source7 as an additional dataset so jump_rr stays consistent with the production paper without breaking existing datasets.

What to add

A single new dataset compound_no_source7 pointing to:

https://cellpainting-gallery.s3.amazonaws.com/cpg0042-chandrasekaran-jump/source_all/workspace/profiles_assembled/compound_no_source7/v1.0/profiles_var_mad_int_featselect_harmony.parquet

SHA256: 8e1e5d9e50c8c7c95ed406981b02adb57c7ff9d253192ed52e661115379be6f1

This file is used for both matching (cosine similarity) and feature analysis — no separate _interpretable variant needed, because Harmony correction is per-feature (each corrected feature maps 1:1 to its original CellProfiler feature).

Feature name mapping

The harmony file has generic column names (X_1 through X_758). The CellProfiler names can be recovered positionally from the companion file:

https://cellpainting-gallery.s3.amazonaws.com/cpg0042-chandrasekaran-jump/source_all/workspace/profiles_assembled/compound_no_source7/v1.0/profiles_var_mad_int_featselect.parquet

Both files have exactly 758 features in the same order. Mapping: X_1 → Cells_AreaShape_Compactness, X_2 → Cells_AreaShape_Extent, etc. Only the parquet schema (column names) needs to be read — not the data (~2.7 GB).

Suggested implementation plan per Claude Opus (⚠️ unverified⚠️ )

Files to modify (all under libs/jump_rr/):

File	Change
src/jump_rr/datasets.py	Add COMPOUND_DATASETS constant, _LOCAL_URLS override dict, hash, get_babel_type() helper, early-return in get_profiles_url()
src/jump_rr/feature_mapping.py (new)	get_harmony_feature_mapping() — schema-only read of featselect file, positional rename dict. rename_harmony_features(df) — applies the rename.
src/jump_rr/consensus.py	get_range(): change dataset == "compound" → dataset in COMPOUND_DATASETS
src/jump_rr/metadata.py	write_metadata(): handle compound_no_source7 in prefix selection (line 94) and skip _interpretable suffix for source URL (lines 118-120)
src/jump_rr/calculate_matches.py	Add ("compound_no_source7", 10) to loop. Change 3× dset != "compound" → dset not in COMPOUND_DATASETS. Use get_babel_type(dset) for get_external_mappers().
src/jump_rr/calculate_features.py	Add ("compound_no_source7", 10, 50) to loop. Rename features after loading. Use get_babel_type() for add_pert_type() and get_external_mappers().
src/jump_rr/galleries.py	Add to loop. Change dset != "compound" → dset not in COMPOUND_DATASETS. Use get_babel_type() for mappers.
tests/test_data_sources.py	Add "compound_no_source7" to parametrize lists

Key design decisions:

• COMPOUND_DATASETS = ("compound", "compound_no_source7") — single constant for all compound-type branching (GPU skip, mapper selection, image site range)
• get_babel_type(dset) — maps dataset name to broad-babel plate_type (returns "compound" for both compound variants)
• _LOCAL_URLS — bypass profile_index.json lookup for datasets not in the GitHub manifest
• Feature rename uses pl.scan_parquet(url).collect_schema().names() (reads only parquet footer, not data)

No changes needed to: significance.py, replicability.py, parse_features.py, generate_databases.sh

Verification

1. uv run pytest tests/ — existing tests pass + new URL test
2. uv run ruff check src/ && uv run ruff format src/
3. Smoke test: download dataset and verify feature rename produces CellProfiler names

[shntnu]

Verification: Harmony per-feature correction and feature mapping

The implementation plan claims that Harmony-corrected features map 1:1 to their original CellProfiler features, so positional renaming (X_1 → Cells_AreaShape_Compactness, etc.) is valid. Here's a reproducible verification.

Background: Harmony's per-feature correction means the correction for feature $j$ depends only on the original values of feature $j$ (the cross-feature influence is indirect, through cluster assignment only).

Schema verification (fast, no data download)

# uv run python3
import polars as pl

harmony_url = "https://cellpainting-gallery.s3.amazonaws.com/cpg0042-chandrasekaran-jump/source_all/workspace/profiles_assembled/compound_no_source7/v1.0/profiles_var_mad_int_featselect_harmony.parquet"
featselect_url = "https://cellpainting-gallery.s3.amazonaws.com/cpg0042-chandrasekaran-jump/source_all/workspace/profiles_assembled/compound_no_source7/v1.0/profiles_var_mad_int_featselect.parquet"

harm_cols = pl.scan_parquet(harmony_url).collect_schema().names()
feat_cols = pl.scan_parquet(featselect_url).collect_schema().names()

harm_feats = [c for c in harm_cols if not c.startswith("Metadata_")]
feat_feats = [c for c in feat_cols if not c.startswith("Metadata_")]

assert len(harm_feats) == len(feat_feats) == 758
assert all(c.startswith("X_") for c in harm_feats)
assert harm_feats == [f"X_{i}" for i in range(1, 759)]
assert all(c.startswith(("Cells_", "Cytoplasm_", "Nuclei_")) for c in feat_feats)
print("Schema checks passed: 758 features, X_1..X_758 ↔ CellProfiler names")

Statistical verification (downloads ~column subset of both files)

Three tests confirm Harmony's per-feature correction on the actual data:

# uv run --with scikit-learn python3
import polars as pl
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score
from scipy.stats import ttest_rel

harmony_url = "https://cellpainting-gallery.s3.amazonaws.com/cpg0042-chandrasekaran-jump/source_all/workspace/profiles_assembled/compound_no_source7/v1.0/profiles_var_mad_int_featselect_harmony.parquet"
featselect_url = "https://cellpainting-gallery.s3.amazonaws.com/cpg0042-chandrasekaran-jump/source_all/workspace/profiles_assembled/compound_no_source7/v1.0/profiles_var_mad_int_featselect.parquet"

harm_schema = pl.scan_parquet(harmony_url).collect_schema()
feat_schema = pl.scan_parquet(featselect_url).collect_schema()
harm_feats = [c for c in harm_schema.names() if not c.startswith("Metadata_")]
feat_feats = [c for c in feat_schema.names() if not c.startswith("Metadata_")]
meta_cols = [c for c in harm_schema.names() if c.startswith("Metadata_")]

# Sample 20 evenly-spaced features
n_test = 20
indices = np.linspace(0, len(harm_feats) - 1, n_test, dtype=int)
test_harm_cols = [harm_feats[i] for i in indices]
test_feat_cols = [feat_feats[i] for i in indices]

harm_df = pl.scan_parquet(harmony_url).select(meta_cols + test_harm_cols).collect()
feat_df = pl.scan_parquet(featselect_url).select(meta_cols + test_feat_cols).collect()
joined = harm_df.join(feat_df, on=meta_cols, how="inner")

H = joined.select(test_harm_cols).to_numpy()
O = joined.select(test_feat_cols).to_numpy()
mask = ~(np.isnan(H).any(axis=1) | np.isnan(O).any(axis=1))
H, O = H[mask], O[mask]

# TEST 1: Same-feature correlation
same_corrs = [np.corrcoef(O[:, j], H[:, j])[0, 1] for j in range(n_test)]
print(f"Test 1 — Same-feature correlation: mean={np.mean(same_corrs):.4f}, min={np.min(same_corrs):.4f}")

# TEST 2: Cross-feature correlation preservation
rng = np.random.default_rng(42)
pairs = rng.choice(n_test, size=(50, 2), replace=True)
pairs = pairs[pairs[:, 0] != pairs[:, 1]]
oo = [np.corrcoef(O[:, i], O[:, j])[0, 1] for i, j in pairs]
oh = [np.corrcoef(O[:, i], H[:, j])[0, 1] for i, j in pairs]
diffs = np.array(oh) - np.array(oo)
t_stat, p_val = ttest_rel(oo, oh)
print(f"Test 2 — Cross-feature corr preservation: mean|Δ|={np.mean(np.abs(diffs)):.4f}, paired t p={p_val:.4f}")

# TEST 3: Feature leakage regression
idx = rng.choice(H.shape[0], size=min(10000, H.shape[0]), replace=False)
Hs, Os = H[idx], O[idx]
for j in range(n_test):
    r2_single = r2_score(Hs[:, j], LinearRegression().fit(Os[:, j:j+1], Hs[:, j]).predict(Os[:, j:j+1]))
    r2_full = r2_score(Hs[:, j], LinearRegression().fit(Os, Hs[:, j]).predict(Os))
    print(f"  {test_feat_cols[j][:45]:45s} R²_single={r2_single:.4f}  R²_full={r2_full:.4f}  Δ={r2_full-r2_single:.4f}")

Output

Test 1 — Same-feature correlation: mean=0.9465, min=0.8730
Test 2 — Cross-feature corr preservation: mean|Δ|=0.0178, paired t p=0.4876
  Cells_AreaShape_Compactness                R²_single=0.9031  R²_full=0.9100  Δ=0.0069
  Cells_Correlation_K_Mito_DNA               R²_single=0.8861  R²_full=0.8902  Δ=0.0042
  Cells_Intensity_MADIntensity_DNA           R²_single=0.8829  R²_full=0.8870  Δ=0.0041
  Cells_RadialDistribution_MeanFrac_Mito_3   R²_single=0.8946  R²_full=0.9031  Δ=0.0085
  Cells_Texture_Correlation_AGP_10_01_256    R²_single=0.9170  R²_full=0.9217  Δ=0.0046
  Cells_Texture_InfoMeas1_Mito_10_03_256     R²_single=0.9033  R²_full=0.9095  Δ=0.0063
  Cytoplasm_AreaShape_Zernike_8_0            R²_single=0.9019  R²_full=0.9064  Δ=0.0045
  Cytoplasm_Correlation_RWC_RNA_DNA          R²_single=0.7622  R²_full=0.7858  Δ=0.0237
  Cytoplasm_Granularity_6_RNA                R²_single=0.9302  R²_full=0.9340  Δ=0.0038
  Cytoplasm_RadialDistribution_FracAtD_mit   R²_single=0.9349  R²_full=0.9388  Δ=0.0038
  Cytoplasm_RadialDistribution_RadialCV_mi   R²_single=0.9228  R²_full=0.9286  Δ=0.0058
  Cytoplasm_Texture_Correlation_RNA_10_02_   R²_single=0.9081  R²_full=0.9143  Δ=0.0062
  Cytoplasm_Texture_InfoMeas1_Mito_5_01_25   R²_single=0.8864  R²_full=0.9000  Δ=0.0136
  Nuclei_AreaShape_Zernike_7_3               R²_single=0.9354  R²_full=0.9393  Δ=0.0039
  Nuclei_Granularity_12_Mito                 R²_single=0.8978  R²_full=0.9075  Δ=0.0096
  Nuclei_Granularity_9_RNA                   R²_single=0.8979  R²_full=0.9042  Δ=0.0062
  Nuclei_RadialDistribution_MeanFrac_DNA_4   R²_single=0.8603  R²_full=0.8654  Δ=0.0052
  Nuclei_Texture_Correlation_AGP_3_02_256    R²_single=0.8448  R²_full=0.8588  Δ=0.0140
  Nuclei_Texture_DifferenceEntropy_ER_3_00   R²_single=0.9041  R²_full=0.9125  Δ=0.0083
  Nuclei_Texture_SumVariance_ER_10_01_256    R²_single=0.9349  R²_full=0.9373  Δ=0.0024

Interpretation

Test	Result	What it means
Same-feature corr	mean r = 0.95	Harmony-corrected features are highly correlated with their originals (not 1.0 because the correction varies by cluster)
Cross-feature corr preservation	mean |Δ| = 0.018, p = 0.49	No significant change in cross-feature correlations → features are not mixed
Feature leakage regression	Δ R² = 0.002–0.024	Adding all other original features barely improves prediction of a harmony feature → correction is per-feature

The small residual Δ in Test 3 is consistent with Harmony's indirect cross-feature influence through cluster assignment (details). The positional renaming X_j → CellProfiler name is valid.

[shntnu]

@AnneCarpenter FYI the harmony corrected (yet interpretable) profiles will be addressed via this repo.

But I will defer to @afermg to decide how/whether JUMP Hub should point to these profiles; these profiles are very paper-specific (see notes above), so adding this to the Hub might get messy (and Alán has been carefully curating the Hub to keep it usable / uncluttered)

Alán - happy to chat about this in person during checkin on Tue

[afermg]

Yeah, we should probably have a chat, since this seems like major changes and a lot of new information that may overload users.

[afermg]

We discussed this in person and came up with a short and long-term plan:

Long term:

The act of generating tables for exploration by biologists will remain a common product of processing images into profiles. It would be useful to have a framework to produce parquet and metadata files that leverage Datasette just as JUMP_rr does. If we formalize it into a library with an API it means that we can repurpose this in future studies and generate a corpus of interactive tables for a variety of datasets, with a common interface and exploration workflow. Thus we want to eventually make it a library with well-defined interfaces, but some assumptions that this code makes must be reworked and some columns may be skipped.

This issue though pertains the specific need of this "no_source7" table for JUMP production, and refactoring the code into a formal library may extend beyond the required timeline.

Short term

We want to reuse most of the code, but the main changes are:
(1) Different source profile
(2) Statistics may change if we remove source 7, so reusing the current p-values pulled from Niranj's table may not be appropriate (since we should also recalculate p-values for these new profiles).

The easiest way to solve short term goals while keeping the path open for (1) is developing a new branch that will NOT be merged until JUMP-rr becomes a library that can produce all these tables.

I think that's all:

• Build this as a new branch (ideally minimising the number of changes necessary, I think just changing the input profiles and adding mapping should suffice).
• We will refactor JUMP-rr into a library in the future, and cover this use-case amongst a few others.

Let me know if anything is unclear

Tagging @jfredinh, as possible owner of the branch

[afermg]

@shntnu @jfredinh Following up on my previous comment. Just to confirm a couple of things that are not written anywhere, I have a couple of questions:

• Are we moving to use the un-PCA'd versions of the profiles?

• Are we adding those profiles to the main manifest (either in a branch or straight to main)? These are paper-specific while JUMP_rr is general across JUMP. That's why we had the previous discussion.
  If I update the JUMP_rr table using those right now, they won't match our "canonical" public versions.

• If we are switching, we still need the p-values (since Niranj's won't suffice anymore since the input data is different). Where can I find those? They should be an artifact of running the pipeline for the whole dataset

(2) Statistics may change if we remove source 7, so reusing the current p-values pulled from Niranj's table may not be appropriate (since we should also recalculate p-values for these new profiles).

As a side-note: @shntnu 's modified version of JUMP_rr introduced some cognitive noise, we need to agree on which the canonical version is (unless we are doing this exclusively for JUMP production and it won't be used anywhere else) to be coherent with the ourselves of two years ago.

[jfredinh]

I would assume using the same features as used in the final production paper would be ideal, @shntnu to direct us to where we can find those.
I know we can fetch them from the production repo, but want to make sure we use the correct ones.

[afermg]

(Code ready to use on https://github.com/broadinstitute/monorepo/tree/update-jump-rr-manifest, and profiles on https://zenodo.org/records/20496083), but will wait until we agree on the decisions for this to update broad links and whatnot.

[afermg]

pinging @shntnu as a reminder

[AnneCarpenter]

Pinging @shntnu on this.

[shntnu]

@afermg

Ignore what this below file is (for now), but just focus on this block:

https://github.com/broadinstitute/jump_production/blob/8106ad684382a7a0a6a58de5629b4249818cc5da/data/processed/dataset_reference/dataset_reference.json#L65-L71

That is the profile used in the paper.

Would you be able to work backwards from that and figure out else is needed here?

If you want to dig into that file, this is the notebook that produces - and displays - that JSON index

https://github.com/broadinstitute/jump_production/blob/main/notebooks/dataset_reference.py

uv run --with marimo marimo run notebooks/dataset_reference.py --sandbox

It pulls data from S3 so it should just work

[afermg]

@shntnu Thanks for the link. That doesn't quite answer any of my questions though:

• Is this set of profiles un-PCA'd?
• Does it match the current official manifest for JUMP? (the one in the datasets repo)
• Since it is only one parquet, is it all the modalities (ORF,CRISPR,compound) together?
• Where are the associated p-values that I need to generate the JUMP_rr tables?

I tried to read the script to see if it answered any questions, but the verbosity makes this tricky (and it seems to just generate the json you also linked). Thus I threw a bot at it and I got:

• features are X_0000, thus probably not unPCA'd
• JCP ids start 1 or 0, thus this only contains compounds.

I got a few of the questions sorted out by the bot, with me validating some things, could you confirm this summary and weight in on whether or not it should be in the official manifest?

┌────────────────────────────────┬─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┐
│ Question                       │ Answer                                                                                                                                                      │
├────────────────────────────────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ Which profile did the paper    │ compound_no_source7.parquet, sourced from profiles_var_mad_int_featselect_harmony.parquet.                                                                  │
│ use?                           │                                                                                                                                                             │
├────────────────────────────────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ Is it PCA-based?               │ No. It is the recipe-Harmony profile with 758 per-feature dimensions, not the separate _rsc PCA+Harmony variant.                                            │
├────────────────────────────────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ Can it support interpretable   │ Yes. Harmony columns are X_1…X_758; the pre-Harmony parquet has 758 CellProfiler names in the same order. The production code itself uses this positional   │
│ features?                      │ mapping. Thus the same Harmony profile can drive matches and feature tables after renaming.                                                                 │
├────────────────────────────────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ Were activity p-values         │ Yes. The production pipeline ran copairs on compound_no_source7 / activity_no_target2 / all_sources / default. It produced per-compound mAP, raw p-values,  │
│ recalculated?                  │ corrected p-values, and activity flags.                                                                                                                     │
├────────────────────────────────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ Should it enter the public     │ Still unanswered. The notebook’s manifest only controls the production report. The production repo has its own profile index, but                           │
│ JUMP manifest?                 │ jump-cellpainting/datasets does not currently list compound_no_source7.                                                                                     │
└────────────────────────────────┴─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘