Feature Flow

Feature Flow is a Snakemake workflow for reducing high-dimensional feature tables into interpretable modules and testing their associations with phenotypes. It retains parallel results for the original features, connecting dimension reduction with feature-level biological follow-up.

The workflow is designed for heterogeneous quantitative features, including gene measurements, immune-repertoire summaries, cell proportions, and cell-communication scores. It provides a tracked analysis from preprocessing and module discovery through statistical testing, optional module annotation, result summarisation, and visualisation.

Overview

Feature Flow supports three input modes: a single feature table, a batch of dataset directories, or multiple tables combined before analysis.

The full workflow performs the following stages:

1. Input validation and combination — validates inputs and optionally row-binds or column-binds compatible feature tables.
2. Feature assessment and preprocessing — records normality, filters features, optionally imputes missing values, benchmarks imputation, and can assess distribution shifts.
3. Module reduction — constructs hierarchical-clustering modules and module eigengenes using all samples or a configured discovery subset.
4. Module annotation — optionally describes gene-like modules using pathway enrichment or names other modules using a generic LLM-backed annotator.
5. Module-set enrichment — optionally tests module members against configured GMT gene sets.
6. Association testing — fits configured parametric and rank-based models for module eigengenes and original member features.
7. Result summarisation — combines model outputs, adjusts p-values, records assumption checks, and selects final statistical fields.
8. Visualisation — produces module-level and original-feature plots with phenotype groups and available annotations.

Requirements

• Snakemake 8 or later
• Python 3
• R and the packages used by the workflow scripts
• The Snakemake SLURM executor plugin for cluster execution
• OpenAI credentials when generic module description is enabled
• Network access to g:Profiler when its annotation backend is enabled

R is invoked through the rscript_bin configured for each run. Run the workflow from the project root so its R environment and relative input, result, and log paths resolve consistently.

Shared R utilities

Reusable R functions are pinned as the Git submodule workflow/scripts/utils. Include it when cloning Feature Flow:

git clone --recurse-submodules https://github.com/liezeltamon/feature-flow.git

For an existing clone, including Feature Flow nested inside another project, initialise all recorded submodules recursively:

git submodule update --init --recursive

Utility development remains in the canonical utils repository. After a utility change is committed and pushed there, update Feature Flow's pinned revision deliberately:

git -C workflow/scripts/utils fetch origin main
git -C workflow/scripts/utils checkout origin/main
git add workflow/scripts/utils
git commit -m "build: update shared utilities"

Development can also start inside workflow/scripts/utils: switch that submodule checkout to a branch, commit and push to the utils remote, then record the new submodule revision in Feature Flow as above.

Each host repository records its own tested utility revision, so the code is shared without silently changing existing workflow runs.

Input data

Data contract

• Feature tables must have the configured sample_key as the first column.
• Sample IDs stay as explicit columns, not row names.
• Numeric columns after sample_key are treated as features.
• Non-numeric feature-table columns after sample_key are ignored.
• Metadata tables provide grouping variables, individual IDs, batch variables, subset keys, and other model covariates.
• Feature-table sample IDs must be present in metadata; metadata is aligned to feature-table order.
• Batch mode expects immediate dataset directories containing:
  • bulk_x_features.csv
  • bulk_metadata.csv
• Combine mode supports:
  • row_bind: same features, different samples
  • column_bind: same samples, different features
• Contract-style feature IDs are kept intact as the official feature names.

Configuration

The default example is config/config.yaml. Additional examples for batch, gene-like, row-bind, and column-bind inputs are available under config/examples/.

Configuration controls preprocessing thresholds, module discovery subsets, annotation backends, model methods, fixed and random effects, contrasts, member-testing scopes, result selection, and plotting.

univariate.member_test_scopes supports independent original-feature analyses over named sets of phenotype groups. Each scope changes the samples used for fitting and can define its own minimum complete-sample requirement, custom contrasts, and global tests. Configurations without explicit scopes use the compatibility scope all_groups.

Core components

• workflow/Snakefile

  • Defines the full workflow graph and final targets.
  • Supports single, batch, and combine input modes.
  • Discovers dataset IDs from input directories in batch mode.
  • Sends module-level testing to univariate_test/ and summarise_tests/.
  • Sends original/member-feature testing to scope-specific univariate_test_members/ and summarise_tests_members/ outputs.
  • Fits each configured member-test scope independently; a scope changes the samples used for fitting, not merely which statistics are displayed.
  • Chooses the canonical module table for downstream testing:
    • gene-like renamed module table when gene module description is enabled
    • generic-renamed module table when generic description is enabled
    • raw modules_hc table otherwise

• workflow/scripts/io_helpers.R

  • Provides shared feature-table and metadata validation.
  • Preserves sample_key as character data so IDs are not silently converted.
  • Enforces sample-key presence and uniqueness.
  • Extracts numeric feature columns and records ignored non-numeric columns.
  • Aligns metadata rows to feature-table sample order.

• workflow/scripts/combine_tables.R

  • Combines multiple input datasets before downstream processing.
  • Uses row_bind when tables share feature columns and contain different samples.
  • Uses column_bind when tables share samples and contain different feature columns.
  • Writes combined feature table, combined metadata table, and an input manifest.

• workflow/scripts/test_normality_wrapper.R

  • Runs feature-level normality tests on numeric feature columns.
  • Writes normality_test_results.csv and normality_summary.csv.
  • Used as an annotation source for summaries.
  • Does not decide whether final univariate p-values come from parametric or nonparametric tests.

• workflow/scripts/preprocess_impute.R

  • Filters features by missingness and optional minimum unique values.
  • Optionally imputes missing values.
  • Optionally benchmarks imputation error across configured methods and missingness frequencies.
  • Optionally filters features by imputation benchmark performance.
  • Optionally applies distribution-shift filtering.
  • Writes processed bulk_x_features.csv, benchmark.rds, and preprocess_summary.csv.

• workflow/scripts/evaluate_preprocessing.R

  • Aggregates preprocessing benchmark outputs across datasets.
  • Produces imputation-error threshold plots.
  • Produces feature-retention summaries.
  • Uses the configured error metric, thresholds, and benchmark directory.

• workflow/scripts/modules_hc.R

  • Reduces preprocessed features into hierarchical-clustering modules.
  • Uses configured subset samples for module discovery.
  • Writes module eigengene tables, module member tables, and eigengene RDS output.
  • Uses repeated subsampling only when imputation was actually applied.
  • Can filter features by subset-specific unique values before module reduction.

• workflow/scripts/describe_module_genes.R

  • Annotates gene-like modules using g:Profiler GO:BP enrichment.
  • Gene-like datasets are identified by dataset IDs containing var_genes.
  • Derives enrichment IDs from module member feature IDs using the configured separator.
  • Prefers significant pathway labels and falls back to suggestive labels when needed.
  • Writes renamed module tables, module_annotations.csv, and per-module audit files.

• workflow/scripts/describe_module_generic.py

  • Names modules from top-loading module members using a generic OpenAI-backed descriptor.
  • Intended for non-gene or mixed feature modules where pathway enrichment is not appropriate.
  • Writes renamed module table and module_annotations.csv.
  • Writes per-module prompt, raw response, parsed response, and meaning text for auditability.

• workflow/scripts/module_set_enrichment.R

  • Tests whether module members are enriched for configured GMT gene sets.
  • Supports module-member feature IDs with contract-style separators.
  • Can use only important-loading members or all members.
  • Writes mHG and ORA enrichment tables per dataset and subset.

• workflow/scripts/univariate_test.R

  • Fits per-feature models for each configured method, usually lmer and rank_transform.
  • Writes result_matrices.rds and failed_features.csv.
  • Stores effect-size, p-value, global p-value, assumption-test, and sample-count vectors.
  • Module mode uses --missing_response_policy fail.
  • Member mode uses --missing_response_policy drop_feature_sample.
  • Member mode drops NA samples separately per feature and applies each configured scope's minimum complete-sample requirement.
  • univariate.member_test_scopes supports independent fits over named group sets; scope-specific custom contrasts and global tests are validated before fitting.
  • Feature values are scaled before testing, preserving existing pipeline behavior.
  • Uses an internal response column for modeling; output feature names remain unchanged.

• workflow/scripts/summarise_tests.R

  • Converts method-specific result_matrices.rds files into wide summary CSVs.
  • Writes one per-source summary_df.csv and one aggregate summary_df.csv.
  • Adjusts contrast p-values within each feature column.
  • Chooses final parametric vs nonparametric fields using model_assumptions_met.
  • model_assumptions_met == TRUE means final fields use the parametric method.
  • model_assumptions_met == FALSE means final fields use the nonparametric method.
  • model_assumptions_met == NA leaves final fields as NA.
  • Joins is_normal from normality outputs when available.
  • Reads module-annotation specifications from the Snakemake-generated module_annotation_manifest.csv, avoiding command-line length limits for large batch runs. The legacy --module_annotation_paths argument remains supported for smaller manual runs.

• workflow/scripts/plot_features.R

  • Plots module-level summaries and original feature distributions.
  • Uses metadata for group labels and sample covariates.
  • Uses module annotations when available.
  • Writes plot outputs under plot_features/<dataset_id>/<subset_name>/.

• docs/slurm.md

  • Documents general SLURM execution, targeted runs, and recovery entry points.
  • Provides neutral resource variables and copyable command examples.

Output

• Resolved configuration and original-feature assessment:

<results_dir>/resolved_dataset_config.csv

<results_dir>/test_normality/<dataset_id>/
  normality_test_results.csv
  normality_summary.csv

• Preprocessing:

<results_dir>/preprocess_impute/<dataset_id>/
  bulk_x_features.csv
  benchmark.rds
  preprocess_summary.csv

• Aggregate preprocessing evaluation:

<results_dir>/evaluate_preprocessing/<error_metric>/
  error_distribution.pdf
  error_threshold*.pdf
  feature_retention_summary*.csv/pdf

• Module reduction:

<results_dir>/modules_hc/<dataset_id>/<subset_name>/
  bulk_x_features.csv
  module_eigengenes_hc.rds
  module_members/

<results_dir>/test_normality_modules/<dataset_id>/<subset_name>/
  normality_test_results.csv
  normality_summary.csv

<results_dir>/feature_distribution_summary_modules/<dataset_id>/<subset_name>/
  feature_distribution_summary.csv
  feature_distribution_overview.csv

• Module description:

<results_dir>/describe_module_genes/<dataset_id>/<subset_name>/
  bulk_x_features.csv
  module_annotations.csv
  modules/

<results_dir>/describe_module_genes_ora/<dataset_id>/<subset_name>/
  bulk_x_features.csv
  module_annotations.csv
  modules/

<results_dir>/describe_module_generic/<dataset_id>/<subset_name>/
  bulk_x_features.csv
  module_annotations.csv
  modules/

• Module-set enrichment:

<results_dir>/module_set_enrichment/<dataset_id>/<subset_name>/
  mhg_enrichment.csv
  ora_enrichment.csv

• Module-level univariate testing:

<results_dir>/univariate_test/<dataset_id>/<method>/
  result_matrices.rds
  failed_features.csv

<results_dir>/summarise_tests/<dataset_id>/
  summary_df.csv

<results_dir>/summarise_tests/
  summary_df.csv

• Member/original-feature univariate testing:

<results_dir>/univariate_test_members/<member_scope>/<dataset_id>/<method>/
  result_matrices.rds
  failed_features.csv

<results_dir>/feature_distribution_summary/<member_scope>/<dataset_id>/
  feature_distribution_summary.csv
  feature_distribution_overview.csv

<results_dir>/summarise_tests_members/<member_scope>/<dataset_id>/
  summary_df.csv

<results_dir>/summarise_tests_members/<member_scope>/
  summary_df.csv

<results_dir>/summarise_tests_available/
  summary_df.csv
  summarise_status.csv

<results_dir>/summarise_tests_members_available/<member_scope>/
  summary_df.csv
  summarise_status.csv

• Visualisation:

<results_dir>/plot_features/<dataset_id>/<subset_name>/

<results_dir>/plot_features_members/<dataset_id>/

• Logs:

<logs_dir>/<rule_name>/
<logs_dir>/<rule_name>.<details>.log

Usage

Run Feature Flow from the project root. Start by inspecting the DAG with the default example configuration:

snakemake \
  --snakefile workflow/Snakefile \
  --configfile config/config.yaml \
  --cores 1 \
  --dry-run

Remove --dry-run to execute the workflow:

snakemake \
  --snakefile workflow/Snakefile \
  --configfile config/config.yaml \
  --cores 1

To use another configuration:

snakemake \
  --snakefile workflow/Snakefile \
  --configfile path/to/config.yaml \
  --cores 8

Input, output, and log paths in the configuration are interpreted relative to the working directory. Workflow-owned rules and scripts are resolved relative to the Feature Flow source tree, allowing the workflow to be used from a standalone checkout or mounted as a submodule.

SLURM execution

docs/slurm.md documents full, targeted, recovery, and dry-run execution. Review the config path, targets, resources, and credentials before using an example.

Recovery workflows

The main workflow/Snakefile plans the complete workflow. Additional entry points support partially completed production runs:

• workflow/Snakefile.univariate_available retries module and original-feature tests whose inputs are available;
• workflow/Snakefile.describe_enrichment_available completes available module descriptions and enrichment analyses;
• workflow/Snakefile.summarise_available creates summaries from results that completed;
• workflow/Snakefile.test provides a reduced workflow entry point for testing.

These entry points are useful with --keep-going, where an unusable small-data branch can fail without stopping independent valid datasets.

Repository layout

feature-flow/
  config/
    config.yaml
    examples/
  docs/
    slurm.md
  workflow/
    Snakefile
    rules/
    scripts/
      utils/  # pinned shared R utility submodule
  .test/
    data/
    tests/

Authors

• Liezel Tamon
  • University of Oxford
  • ORCID

References

Köster, J., Mölder, F., Jablonski, K. P., Letcher, B., Hall, M. B., Tomkins-Tinch, C. H., Sochat, V., Forster, J., Lee, S., Twardziok, S. O., Kanitz, A., Wilm, A., Holtgrewe, M., Rahmann, S., & Nahnsen, S. Sustainable data analysis with Snakemake. F1000Research, 10:33, 2021. https://doi.org/10.12688/f1000research.29032.2

Bashford-Rogers Lab. vdjremix. https://github.com/Bashford-Rogers-lab/vdjremix