Add function for feature prioritization - #17
Conversation
This function computes feature improvement by reading and processing feature data from parquet files, applying various transformations and calculations to derive insights on feature importance and contributions.
Added computeFeatureScore function to calculate feature scores based on various metrics and cluster information.
…ion and introducing wide table functions for global score calculations
…update parameters for clarity
… feature and cluster analysis
AbhirupaGhosh
changed the title
|
This script introduces a workflow for prioritizing molecular features associated with resistance across drugs and drug classes. This workflow identifies features that are consistently important across repeated model fits (seeds), making the results less noisy and more interpretable. The workflow:
The run would be |
|
Minor suggestions for now or later:
@AbhirupaGhosh, can you clarify the following:
|
Standardize feature column names across the codebase (Variable -> variable, Importance -> importance, Sign -> sign) and update all uses (mutate, group_by, joins, transmute, distinct). Rename summariseFeatureAcrossSeeds -> summariseFeaturesAcrossSeeds, buildSignalNetwork -> buildFeatureNetwork, and plotSignalNetworkD3 -> plotFeatureNetworkD3. Adjust join keys (variable == feature), grouping and summarise calls, examples, and roxygen docs to match the new names. Includes minor formatting and whitespace cleanups. These changes unify naming and clarify the feature/cluster network API.
jananiravi
left a comment
jananiravi
left a comment
There was a problem hiding this comment.
Added a few minor fixes inline + some comments for clarification.
All other scripts are named do_something. so 'prioritize_features' or 'rescore_features' or something like that?
Note from Slack
In summariseClusters() and buildFeatureNetwork(), when a top-ranked domain/COG feature gets joined to cluster_feature, it fans out into every cluster carrying that domain — including clusters with no real connection to the resistance signal, just shared domain architecture.
Concretely:
- A single promiscuous COG (fan-out up to 189) inflates frequency/n_variables/cluster_mean_rank_score for up to 189 different clusters in summariseClusters(), none of which may be the true driver.
- In buildFeatureNetwork(), feature_cluster_edges will make that one feature a hub node connected to up to ~189 cluster nodes, which will dominate/clutter plotFeatureNetworkD3()'s output and overstate specificity.
So, is this fan-out intentional (i.e., "show all plausible clusters a domain could implicate"), or should there be a cap/weighting (e.g., down-weight edges by 1/n_clusters, or exclude domains above some fan-out threshold) so ubiquitous COGs/Pfams don't swamp the network? This is a modeling decision, not a coding bug.
| protein_scale_realization <- ranked_features |> | ||
| dplyr::filter(!is.na(cluster), shuffled == FALSE) |> | ||
| top_clusters <- top_features |> | ||
| dplyr::left_join(cluster_feature, by = dplyr::join_by(Variable == feature)) |> |
There was a problem hiding this comment.
In dplyr::left_join(top_features, cluster_feature, by = dplyr::join_by(variable == :
Detected an unexpected many-to-many relationship between `x` and `y`.
ℹ Row 1 of `x` matches multiple rows in `y`.
ℹ Row 54683 of `y` matches multiple rows in `x`.
ℹ If a many-to-many relationship is expected, set `relationship = "many-to-many"` to silence this warning.
|
|
||
| cluster_table <- top_clusters |> | ||
| dplyr::mutate(model_id = make_model_id(drug_label, drug_or_class)) |> | ||
| dplyr::left_join( |
There was a problem hiding this comment.
In dplyr::left_join(feature_table, cluster_feature, by = dplyr::join_by(variable == :
Detected an unexpected many-to-many relationship between `x` and `y`.
ℹ Row 1 of `x` matches multiple rows in `y`.
ℹ Row 151601 of `y` matches multiple rows in `x`.
ℹ If a many-to-many relationship is expected, set `relationship = "many-to-many"` to silence this warning.
| dplyr::ungroup() |> | ||
| dplyr::group_by(drug_label, drug_or_class) |> | ||
| dplyr::filter( | ||
| # n_seeds == max(n_seeds), |
There was a problem hiding this comment.
intentionally left out? what's the consistency across seeds you're looking for?
There was a problem hiding this comment.
yes, I am using the threshold_sd_rank as a proxy for consistency.
max(n_seeds) is too stringent
| ) |> | ||
| dplyr::mutate( | ||
| contribution = Importance / sum(Importance, na.rm = TRUE), | ||
| contribution = importance / sum(importance, na.rm = TRUE), |
There was a problem hiding this comment.
cumulative impact first then rank?
There was a problem hiding this comment.
ranking based on the contribution or importance will give the same. I kept contribution if we want to care about contribution %tile.
| mean_rank_score = mean(rank_score, na.rm = TRUE), | ||
| best_rank = min(rank, na.rm = TRUE), | ||
| rank_consistent = dplyr::n_distinct(rank) == 1, | ||
| rank_sd = sd(rank, na.rm = TRUE), |
There was a problem hiding this comment.
rank_score, mean_rank_score are normalized but not this one? are these comparable?
There was a problem hiding this comment.
I used rank_sd because it is easy to dry run with whole numbers.
| contribution = Importance / sum(Importance, na.rm = TRUE), | ||
| contribution = importance / sum(importance, na.rm = TRUE), | ||
| rank = dplyr::dense_rank(dplyr::desc(contribution)), | ||
| n_features = dplyr::n(), |
There was a problem hiding this comment.
Is this to get totals by feature scale? If so, group_by then n() within summarize?
|
@AbhirupaGhosh @epbrenner A few more follow-up questions while reworking the scoring/reconciliation logic in Model performance for a drop-poor-models filter (what @AbhirupaGhosh mentioned yesterday). We're adding a default arg to exclude features from poorly-performing models (MCC-based, cutoff = 0, i.e., above-random) before pooling anything into prioritization — a model with weak MCC shouldn't have its "top features" treated the same as a model that's actually predictive. (one more Q from Claude) Context for why this is important: we noticed thoughts/comments? |
(autogen) Add igraph and networkD3 to DESCRIPTION and export new utilities in NAMESPACE. Improve scoreFeaturesWithinSeed: add core_contribution_threshold and exclude_feature_types, handle ties, compute cumulative contribution and in_core. Use median_rank_score across summariseFeaturesAcrossSeeds/topFeatures and weight clusters by cluster membership in summariseClusters. Update buildFeatureNetwork to use median-based scores, adjust edge weights, and add network building/plotting helpers. Comment out removed demo data docs in R/data.R and add many generated .Rd man pages for new/updated functions to unblock documentation generation.
jananiravi
left a comment
jananiravi
left a comment
There was a problem hiding this comment.
@AbhirupaGhosh, a few notes based on what's done and what's pending
Many-to-many join warnings — fixed. Added relationship = "many-to-many" to both joins (summariseClusters() and buildFeatureNetwork()'s feature_cluster_edges) now that the fan-out is an acknowledged, intentional part of the design rather than an accidental warning.
Domain/COG fan-out (Slack thread) — resolved via down-weighting rather than a cap (needs checking). When a feature maps to N clusters, its rank-score contribution to each cluster (and its edge weight in the network plot) is now divided by N, so a promiscuous COG/domain/etc. no longer inflates every cluster it touches equally, while still being visible everywhere it genuinely maps. cluster_max_rank_score is deliberately left unweighted — it represents "the single best feature backing this cluster," and specificity shouldn't matter for that statistic the way it does for the mean/median. Worth noting for anyone reviewing: this affects more than domains/COGs — we checked, and every feature type except proteins (which are 1:1 by construction) has real fan-out, with struct showing the most individual variables affected before it was pulled out of scoring entirely (see below).
Sign filtering (sign == "NEG") — left as-is this round; still provisional, now documented inline at the filter itself (not just in this thread) so it doesn't read as an unexplained hardcode to a future reader. plot_ml.R:1163 already treats NEG as -Importance elsewhere in the package, which is useful context but doesn't resolve the open question of whether NEG-only is the right selection criterion — that's being tracked separately, not folded into this round.
Naming convention — the original comment was about the .R file name (feature_rescoring.R vs. a do_something.R pattern), not function casing. Checked: the repo's file naming is actually mixed — run_ML.R, run_Fisher_tests.R, merge_ml_results.R, plot_ml.R, prep_ml.R are verb-first, but core_ml.R is topic-first like this one. Leaving the file name as-is; if it gets renamed, that should be a separate commit from you to preserve blame/provenance rather than bundled into this PR's changes.
One more, not from this thread but adjacent: plotTopClusters() in plot_ml.R:1133-1180 independently reimplements the same Variable cleanup regex (strip .NCBIFAM, leading X, domain suffix, fig.→fig|) that now lives in scoreFeaturesWithinSeed(). Flagging as a follow-up — consolidating into a shared internal helper touches a second file outside this PR's scope, so not fixing it here.
Also added a comment for you in data.R.
A few more things, since two open pieces would benefit from your read given they touch your original contribution/rank_score design directly. Flagging them as things to give a shot at, not decisions already made on your behalf.
Struct is now fully out of scoring (already implemented) — which we have discussed at length. scoreFeaturesWithinSeed() drops it via exclude_feature_types = "struct" before any rank/contribution math runs. Reasoning: its IDs are composite (three joined gene identifiers) rather than single entities like the other five scales, and its candidate-variable count dwarfs the others by orders of magnitude — pooling it in was comparing a compound signal against five primary ones on an incomparable scale. It's reserved for post-hoc biological annotation once top clusters are identified (via the existing cluster_feature mapping, which it already joins into cleanly), not for scoring/ranking/thresholding. Flag if this doesn't match what you had in mind for it.
Bin/count reconciliation — want to take a crack at this one first before I dig in? We (Claude and I) checked empirically: n_features (row count in the top-features table) between binary and counts encoding of the same drug/feature_type can differ by up to 33× (e.g. one args case: binary=232, counts=7) — not a small effect. That traces back to extractTopFeats()'s cumulative-importance cutoff already having run per-model before this data reaches feature_rescoring.R, so n_features here reflects how concentrated each model's own importance distribution happened to be, not the size of the candidate feature space. Given that, averaging two already-normalized subtype-level scores isn't safe. The fix we'd suggest: pool the raw seed-level rows from both subtypes for a given (species, drug_label, drug_or_class, feature_type, variable) and recompute contribution/rank/rank_score fresh against the combined set, rather than averaging two subtype-level medians — slots in between scoreFeaturesWithinSeed() and summariseFeaturesAcrossSeeds(). If you know of cases where binary/counts runs used different n_top_feats/prop_vi_top_feats settings, that'd confirm whether the 33× spread is a real model-shape signal or a pipeline-config artifact — worth checking either way, since it changes how much to trust n_features as an input.
Contribution normalization + priority_score — also yours to take a shot at if you want it. contribution (share of a seed's total importance) has the same cross-group scale problem as n_features above: its typical magnitude is roughly 1/n_features, so even within the five remaining scales, args contributions run ~30× larger than proteins contributions purely from group size, not signal strength. Before contribution can safely combine with rank_score into a single priority_score, it needs a group-size-normalized variant — one option is contribution_enrichment = contribution * n_features (an enrichment ratio, 1.0 = average share, higher = concentrated signal), computed once bin/count reconciliation above lands. So, a simple product may miss this num_features magnitude that may play a role. cum_contrib/in_core (cumulative-contribution core-membership, with tied blocks handled so a threshold never splits a tie -- per Monday discussion) are already in and checked against real edge cases — worth a look since the actual priority_score formula (how rank_score and contribution combine) is your call to make. @epbrenner and I can look at your implementation as soon as you add them.
Still blocking -- need from you, separately: the MCC-based model-drop filter needs to know where the merged performance parquet lives (or whether to run buildPerfPq() fresh) — see earlier note (Slack discussion) on that.
Go through my recent commits to catch anything else I may have missed here.
2 participants
This function computes feature improvement by reading and processing feature data from parquet files, applying various transformations and calculations to derive insights on feature importance and contributions.
Description
What kind of change(s) are included?
Checklist
Please ensure that all boxes are checked before indicating that this pull request is ready for review.