mfa
is an R package for fitting a Bayesian mixture of factor analysers to infer developmental trajectories with bifurcations from single-cell gene expression data. It is able to jointly infer pseudotimes, branching, and genes differentially regulated across branches using a generative, Bayesian hierarchical model. Inference is performed using fast Gibbs sampling.
mfa
can be installed in one of two ways:
if (!requireNamespace("BiocManager", quietly=TRUE))
install.packages("BiocManager")
BiocManager::install("mfa")
library(mfa)
This requires the devtools
package to be installed first
install.packages("devtools") # If not already installed
devtools::install_github("kieranrcampbell/mfa")
library(mfa)
We first create some synthetic data for 100 cells and 40 genes calling the mfa
function create_synthetic
. This returns a list with gene expression, pseudotime, branch allocation, and various parameter estimates:
synth <- create_synthetic(C = 100, G = 40)
print(str(synth))
## List of 7
## $ X : num [1:100, 1:40] 7.671 0 0.315 0 2.85 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:100] "cell1" "cell2" "cell3" "cell4" ...
## .. ..$ : chr [1:40] "feature1" "feature2" "feature3" "feature4" ...
## $ branch : int [1:100] 1 0 1 0 0 1 1 0 0 1 ...
## $ pst : num [1:100] 0.311 0.671 0.999 0.807 0.459 ...
## $ k : num [1:40, 1:2] -9.6 7.79 -6.86 6 -9.44 ...
## $ phi : num [1:40, 1:2] 7.94 5.08 5.05 8.76 8.4 ...
## $ delta : num [1:40, 1:2] 0.311 0.247 0.176 0.26 0.028 ...
## $ p_transient: num 0
## NULL
We can then PCA and put into a tidy format:
df_synth <- as_data_frame(prcomp(synth$X)$x[,1:2]) %>%
mutate(pseudotime = synth$pst,
branch = factor(synth$branch))
and have a look at a PCA representation, coloured by both pseudotime and branch allocation:
ggplot(df_synth, aes(x = PC1, y = PC2, color = pseudotime)) + geom_point()
ggplot(df_synth, aes(x = PC1, y = PC2, color = branch)) + geom_point()
mfa
The input to mfa
is either an ExpressionSet
(e.g. from using the package Scater) or a cell-by-gene expression matrix. If an ExpressionSet
is provided then the values in the exprs
slot are used for gene expression.
We invoke mfa
with a call to the mfa(...)
function. Depending on the size of the dataset and number of MCMC iterations used, this may take some time:
m <- mfa(synth$X)
print(m)
## MFA fit with
## 100 cells and 40 genes
## ( 2000 iterations )
Particular care must be paid to the initialisation of the pseudotimes: by default they are initialised to the first principal component, though if the researcher suspects (based on plotting marker genes) that the trajectory corresponds to a different PC, this can be set using the pc_initialise
argument.
As in any MCMC analysis, basic care is needed to make sure the samples have converged to something resembling the stationary distribution (see e.g. Cowles and Carlin (1996) for a full discussion).
For a quick summary of these, mfa
provides two functions: plot_mfa_trace
and plot_mfa_autocorr
for quick plotting of the trace and autocorrelation of the posterior log-likelihood:
plot_mfa_trace(m)
plot_mfa_autocorr(m)
We can extract posterior mean estimates along with credible intervals using the summary
function:
ms <- summary(m)
print(head(ms))
## # A tibble: 6 × 5
## pseudotime branch branch_certainty pseudotime_lower pseudotime_upper
## <dbl> <fct> <dbl> <dbl> <dbl>
## 1 -0.472 1 1 -0.740 -0.280
## 2 0.851 2 1 0.645 1.13
## 3 0.981 1 1 0.732 1.50
## 4 1.01 2 1 0.821 1.38
## 5 0.130 2 1 0.00336 0.381
## 6 1.02 1 1 0.672 1.41
This has six entries:
pseudotime
The MAP pseudotime estimatebranch
The MAP branch estimatebranch_certainty
The proportion of MCMC traces (after burn-in) for which the cell was assigned to the MAP branchpseudotime_lower
and pseudotime_upper
: the lower and upper 95% highest-probability-density posterior credible intervalsWe can compare the inferred pseudotimes to the true values:
qplot(synth$pst, ms$pseudotime, color = factor(synth$branch)) +
xlab('True pseudotime') + ylab('Inferred pseudotime') +
scale_color_discrete(name = 'True\nbranch')
And we can equivalently plot the PCA representation coloured by MAP branch:
mutate(df_synth, inferred_branch = ms[['branch']]) %>%
ggplot(aes(x = PC1, y = PC2, color = inferred_branch)) +
geom_point() +
scale_color_discrete(name = 'Inferred\nbranch')
A unique part of this model is that through an ARD-like prior structure on the loading matrices we can automatically infer which genes are involved in the bifurcation process. For a quick-and-dirty look we can use the plot_chi
function, where larger values of inverse-chi imply the gene is associated with the bifurcation:
plot_chi(m)
To calculate the MAP values for chi we can call the calculate_chi
function, which returns a data_frame
with the feature names and values:
posterior_chi_df <- calculate_chi(m)
head(posterior_chi_df)
## # A tibble: 6 × 2
## feature chi_map
## <chr> <dbl>
## 1 feature1 0.880
## 2 feature2 0.604
## 3 feature3 0.487
## 4 feature4 1.31
## 5 feature5 0.309
## 6 feature6 0.961
mfa
classA call to mfa(...)
returns an mfa
object that contains all the information about the dataset and the MCMC inference performed. Note that it does not contain a copy of the original data. We can see the structure by calling str
on an mfa
object:
str(m, max.level = 1)
## List of 10
## $ traces :List of 10
## $ iter : num 2000
## $ thin : num 1
## $ burn : num 1000
## $ b : num 2
## $ collapse : logi FALSE
## $ N : int 100
## $ G : int 40
## $ feature_names: chr [1:40] "feature1" "feature2" "feature3" "feature4" ...
## $ cell_names : chr [1:100] "cell1" "cell2" "cell3" "cell4" ...
## - attr(*, "class")= chr "mfa"
This contains the following slots:
traces
- the raw MCMC traces (discussed in following section)iter
- the number of MCMC iterationsthin
- the thinning of the MCMC chainburn
- the number of MCMC iterations thrown away as burn-inb
- the number of branches modelledcollapse
- whether collapsed Gibbs sampling was implementedN
- the number of cellsG
- the number of features (e.g. genes)feature_names
- the names of the features (e.g. genes)cell_names
- the names of the cellsMCMC traces can be accessed through the traces
slot of an mfa
object. This gives a list with an element for each variable, along with the log-likelihood:
print(names(m$traces))
## [1] "tau_trace" "gamma_trace" "pst_trace"
## [4] "theta_trace" "lambda_theta_trace" "chi_trace"
## [7] "eta_trace" "k_trace" "c_trace"
## [10] "lp_trace"
For non-branch-specific variables this is simply a matrix. For example, for the variable \(\tau\) is just an interation-by-gene matrix:
str(m$traces$tau_trace)
## num [1:1000, 1:40] 9.55 7.91 7.59 9.17 11.71 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:40] "tau[1]" "tau[2]" "tau[3]" "tau[4]" ...
We can easily get the posterior mean by calling colMeans
. More fancy posterior density estimation can be perfomed using the MCMCglmm
package, such as posterior.mode(...)
for MAP estimation (though in practice this is often similar to posterior mean). We can estimate posterior intervals using the HPDInterval(...)
function from the coda
package (note that traces must be converted to coda
objects before calling either of these).
Some variables are branch dependent, meaning the traces returned are arrays (or tensors in fashionable speak) that have dimension iteration x gene x branch
. An example is the \(k\) variable:
str(m$traces$k_trace)
## num [1:1000, 1:40, 1:2] -0.829 -0.882 -0.88 -0.74 -0.678 ...
To get posterior means (or modes, or intervals) we then need to use the apply
function to iterate over the branches. To find the posterior means of k
, we then call
pmean_k <- apply(m$traces$k_trace, 3, colMeans)
str(pmean_k)
## num [1:40, 1:2] -0.899 0.81 -0.788 0.986 -0.242 ...
This returns a gene-by-branch matrix of posterior estimates.
sessionInfo()
## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.20-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] dplyr_1.1.4 ggplot2_3.5.1 mfa_1.28.0 BiocStyle_2.34.0
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 tensorA_0.36.2.1 xfun_0.48
## [4] bslib_0.8.0 Biobase_2.66.0 GGally_2.2.1
## [7] lattice_0.22-6 vctrs_0.6.5 tools_4.4.1
## [10] generics_0.1.3 parallel_4.4.1 tibble_3.2.1
## [13] fansi_1.0.6 highr_0.11 pkgconfig_2.0.3
## [16] Matrix_1.7-1 RColorBrewer_1.1-3 ggmcmc_1.5.1.1
## [19] lifecycle_1.0.4 cubature_2.1.1 compiler_4.4.1
## [22] farver_2.1.2 MatrixModels_0.5-3 mcmc_0.9-8
## [25] munsell_0.5.1 tinytex_0.53 codetools_0.2-20
## [28] SparseM_1.84-2 quantreg_5.99 htmltools_0.5.8.1
## [31] sass_0.4.9 yaml_2.3.10 pillar_1.9.0
## [34] jquerylib_0.1.4 tidyr_1.3.1 MASS_7.3-61
## [37] cachem_1.1.0 magick_2.8.5 nlme_3.1-166
## [40] ggstats_0.7.0 tidyselect_1.2.1 digest_0.6.37
## [43] purrr_1.0.2 bookdown_0.41 labeling_0.4.3
## [46] splines_4.4.1 fastmap_1.2.0 grid_4.4.1
## [49] colorspace_2.1-1 cli_3.6.3 magrittr_2.0.3
## [52] survival_3.7-0 utf8_1.2.4 ape_5.8
## [55] corpcor_1.6.10 withr_3.0.2 scales_1.3.0
## [58] MCMCglmm_2.36 rmarkdown_2.28 coda_0.19-4.1
## [61] evaluate_1.0.1 knitr_1.48 rlang_1.1.4
## [64] MCMCpack_1.7-1 Rcpp_1.0.13 glue_1.8.0
## [67] BiocManager_1.30.25 BiocGenerics_0.52.0 jsonlite_1.8.9
## [70] R6_2.5.1 plyr_1.8.9
Cowles, Mary Kathryn, and Bradley P Carlin. 1996. “Markov Chain Monte Carlo Convergence Diagnostics: A Comparative Review.” Journal of the American Statistical Association 91 (434): 883–904.