Introduction to bionmf

Overview

bionmf fits non-negative matrix factorization (NMF) to a gene-by-sample expression matrix. The model is

V≈WH,W≥0,H≥0. V \approx W H, \qquad W \ge 0,\ H \ge 0.

  • Rows of VV are genes; columns are samples.
  • Columns of WW are latent gene programs.
  • Rows of HH are sample loadings on those programs.

The fit uses Lee–Seung multiplicative updates under the Frobenius loss, implemented in base R.

library(bionmf)

Load the bundled example

The package ships a simulated matrix (50 genes ×\times 100 samples, with 3 distinct expression patterns) as inst/extdata/expression_matrix.csv. Use expression_matrix_path() so the file is found after install, without hard-coded working-directory paths.

path <- expression_matrix_path()
V <- load_expression_matrix(path)
dim(V)
#> [1]  50 100

load_expression_matrix() checks that values are non-missing and non-negative, which NMF requires.

Fit NMF

Choose a factorization rank (number of programs), a maximum number of iterations, and a relative-error tolerance. A seed makes the random initialization reproducible.

fit <- run_nmf(
  V,
  rank = 3L,
  max_iter = 200L,
  tol = 1e-4,
  seed = 42L
)

fit$iterations
#> [1] 128

The result is a list with factor matrices and the error trajectory:

  • fit$W — genes ×\times factors
  • fit$H — factors ×\times samples
  • fit$error_history — Frobenius error after each iteration
  • fit$iterations — number of updates performed (may stop early)
dim(fit$W)
#> [1] 50  3
dim(fit$H)
#> [1]   3 100

Diagnostics

Reconstruction quality is summarized by the Frobenius residual ∥V−WH∥F\|V - WH\|_F and by the fraction of variance explained.

reconstruction_error(V, fit$W, fit$H)
#> [1] 20.79332
explain_variance(V, fit$W, fit$H)
#> [1] 0.9962696

Plots

The convergence curve should decrease and then flatten when the relative change in error falls below tol.

plot_error_history(fit$error_history)

Heatmaps of WW and HH show which genes load on each program and how samples use those programs.

plot_nmf_heatmaps(fit$W, fit$H)

Using your own data

Pass any gene-by-sample CSV with gene IDs in the first column:

V <- load_expression_matrix("path/to/your_expression.csv")
fit <- run_nmf(V, rank = 4L, seed = 1L)

See ?run_nmf, ?load_expression_matrix, and ?explain_variance for argument details.