Skip to content

fix: select highest-variance features in pgx.FindClusters - #513

Merged
phisanti merged 2 commits into
edgyfrom
fix/resort-clustering-features
Sep 2, 2026
Merged

fix: select highest-variance features in pgx.FindClusters#513
phisanti merged 2 commits into
edgyfrom
fix/resort-clustering-features

Conversation

@phisanti

Copy link
Copy Markdown
Contributor

While clustering the sepsis multicohort dataset (Tsalik EL, Langley RJ, Dinwiddie DL, et al. An integrated transcriptome and expressed variant analysis of sepsis survival and death. Genome Med. 2014;6(11):111. doi:10.1186/s13073-014-0111-5) I got three "endotypes" that split by cohort instead of by biology, and reran the same call with a different seed and got a different answer. That sent me to the backend.

pgx.FindClusters() had three problems stacked on top of each other:

  1. R/pgx-cluster.R:275 sorted by SD ascending, so top.sd picked the lowest-variance features. On CPTAC_GBM (a dataset from axel) it selected genes with SD 0.004-0.118 against a dataset max of 3.83.
  2. :285 used out$v from irlba directly. Those are unit-norm, so a component holding 46% of the variance got the same weight as pure noise. This one cancelled out any benefit from fixing Packages that are required but cant be installed #1.
  3. :296 called kmeans() with no nstart, so a single random init decided the result.

Fix is 3 lines. But I'd like a second opinion before merging:

  • This changes output. pgx.FindClusters() has no callers in playbase; it's used in omicsplayground for the timeseries gene modules (board.timeseries/R/board_server.R:78).
  • nstart = 20 makes kmeans ~20x slower. Fine at these sizes, but say so if it's called on something much bigger.
  • I left npca <- min(npca, dim(X) - 1) at :281 alone. It lets npca reach
    n_samples - 1, which makes irlba warn about computing too many singular
    values. Harmless once fix issue #232 on OPG #2 is in, but flagging it.

Two related things I hit and did not touch here, happy to open separate issues: find.clusters is hardcoded FALSE at :216 so pgx$samples$cluster is never populated by the current pipeline.

- order(rowSds(X)) sorted ascending, so top.sd selected the LOWEST
  variance features; on CPTAC_GBM this picked genes with SD 0.004-0.118
  against a dataset max of 3.83, making clusters indistinguishable from noise
- use PC scores (v * d) instead of the unit-norm right singular vectors,
  which weighted a 46% variance component the same as pure noise and
  cancelled out any benefit from correct feature selection
- add nstart = 20 to kmeans so runs are reproducible across seeds; adjusted
  Rand index against known sample labels goes from ~0 to 0.92-1.00 and
  reseed stability from 0.01-0.20 to 1.00 on five local datasets
@phisanti
phisanti requested a review from zitoa August 25, 2026 11:52
@ivokwee

ivokwee commented Sep 2, 2026

Copy link
Copy Markdown
Member

@phisanti Approved — CI is green and the branch merges cleanly into `edgy`, so this is ready to merge. 🚀

(Heads-up: the three `.Rd` files in this PR are just roxygen doc-regeneration and already match what is on `edgy` — no behavioral change, so nothing to worry about there.)

@phisanti
phisanti merged commit 5f7f542 into edgy Sep 2, 2026
1 check passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants