# Introduction to MFA

Kieran R Campbell

#### 2017-3-3

# Contents

# 1 Introduction

`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.

# 2 Installation

`mfa` can be installed in one of two ways:

## 2.1 From Bioconductor

    if (!requireNamespace("BiocManager", quietly=TRUE))
        install.packages("BiocManager")
    BiocManager::install("mfa")
    library(mfa)

## 2.2 From Github

This requires the `devtools` package to be installed first

    install.packages("devtools") # If not already installed
    devtools::install_github("kieranrcampbell/mfa")
    library(mfa)

# 3 An example on synthetic data

## 3.1 Generating synthetic data

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] 9.31 7.22 6.45 17.59 7.43 ...
    ##   ..- attr(*, "dimnames")=List of 2
    ##   .. ..$ : chr [1:100] "cell1" "cell2" "cell3" "cell4" ...
    ##   .. ..$ : chr [1:40] "feature1" "feature2" "feature3" "feature4" ...
    ##  $ branch     : int [1:100] 0 0 0 1 0 0 0 0 0 1 ...
    ##  $ pst        : num [1:100] 0.394 0.164 0.195 0.859 0.265 ...
    ##  $ k          : num [1:40, 1:2] 5.22 -5.26 9.3 -5.97 -8.28 ...
    ##  $ phi        : num [1:40, 1:2] 8.61 8.1 7.26 9.99 5.17 ...
    ##  $ delta      : num [1:40, 1:2] 0.374 0.157 0.291 0.284 0.38 ...
    ##  $ 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()

## 3.2 Calling `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.

## 3.3 MCMC diagnostics

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)

## 3.4 Plotting results

We can extract posterior mean estimates along with credible intervals using the `summary` function:

    ms <- summary(m)
    print(head(ms))

    ## # A tibble: 6 x 5
    ##   pseudotime branch branch_certainty pseudotime_lower pseudotime_upper
    ##        <dbl> <fct>             <dbl>            <dbl>            <dbl>
    ## 1     -0.347 1                 1               -0.513           -0.115
    ## 2     -0.750 2                 1               -1.02            -0.596
    ## 3     -0.694 2                 1               -0.925           -0.487
    ## 4      1.04  1                 1                0.760            1.53
    ## 5     -0.443 2                 0.997           -0.711           -0.243
    ## 6      0.586 2                 1                0.131            0.685

This has six entries:

- `pseudotime` The MAP pseudotime estimate
- `branch` The MAP branch estimate
- `branch_certainty` The proportion of MCMC traces (after burn-in) for which the cell was assigned to the MAP branch
- `pseudotime_lower` and `pseudotime_upper`: the lower and upper 95% highest-probability-density posterior credible intervals

We 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')

## 3.5 Finding genes that bifurcate

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 x 2
    ##   feature  chi_map
    ##   <chr>      <dbl>
    ## 1 feature1   0.803
    ## 2 feature2   0.815
    ## 3 feature3   0.566
    ## 4 feature4   0.927
    ## 5 feature5   1.11
    ## 6 feature6   0.902

# 4 Advanced usage

## 4.1 The `mfa` class

A 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 iterations
- `thin` - the thinning of the MCMC chain
- `burn` - the number of MCMC iterations thrown away as burn-in
- `b` - the number of branches modelled
- `collapse` - whether collapsed Gibbs sampling was implemented
- `N` - the number of cells
- `G` - the number of features (e.g. genes)
- `feature_names` - the names of the features (e.g. genes)
- `cell_names` - the names of the cells

## 4.2 Accessing MCMC traces

MCMC 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] 13.7 11.5 12.3 15.3 15.5 ...
    ##  - 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] 1.31 1.24 1.33 1.32 1.32 ...

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] 1.086 -0.696 0.727 -0.908 -0.965 ...

This returns a gene-by-branch matrix of posterior estimates.

# 5 Technical

    sessionInfo()

    ## R version 4.1.3 (2022-03-10)
    ## Platform: x86_64-w64-mingw32/x64 (64-bit)
    ## Running under: Windows Server 2012 R2 x64 (build 9600)
    ##
    ## Matrix products: default
    ##
    ## locale:
    ## [1] LC_COLLATE=C
    ## [2] LC_CTYPE=English_United States.1252
    ## [3] LC_MONETARY=English_United States.1252
    ## [4] LC_NUMERIC=C
    ## [5] LC_TIME=English_United States.1252
    ##
    ## attached base packages:
    ## [1] stats     graphics  grDevices utils     datasets  methods   base
    ##
    ## other attached packages:
    ## [1] dplyr_1.0.8      ggplot2_3.3.5    mfa_1.16.3       BiocStyle_2.22.0
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] Rcpp_1.0.8.3        ape_5.6-2           lattice_0.20-45
    ##  [4] tidyr_1.2.0         corpcor_1.6.10      assertthat_0.2.1
    ##  [7] digest_0.6.29       utf8_1.2.2          plyr_1.8.6
    ## [10] R6_2.5.1            MatrixModels_0.5-0  evaluate_0.15
    ## [13] coda_0.19-4         highr_0.9           pillar_1.7.0
    ## [16] rlang_1.0.2         SparseM_1.81        cubature_2.0.4.2
    ## [19] jquerylib_0.1.4     magick_2.7.3        Matrix_1.4-0
    ## [22] rmarkdown_2.13      labeling_0.4.2      MCMCglmm_2.33
    ## [25] stringr_1.4.0       munsell_0.5.0       compiler_4.1.3
    ## [28] xfun_0.30           pkgconfig_2.0.3     BiocGenerics_0.40.0
    ## [31] ggmcmc_1.5.1.1      mcmc_0.9-7          htmltools_0.5.2
    ## [34] tidyselect_1.1.2    tibble_3.1.6        tensorA_0.36.2
    ## [37] bookdown_0.25       codetools_0.2-18    reshape_0.8.8
    ## [40] fansi_1.0.2         crayon_1.5.0        withr_2.5.0
    ## [43] MASS_7.3-55         grid_4.1.3          GGally_2.1.2
    ## [46] nlme_3.1-155        jsonlite_1.8.0      gtable_0.3.0
    ## [49] lifecycle_1.0.1     DBI_1.1.2           magrittr_2.0.2
    ## [52] scales_1.1.1        cli_3.2.0           stringi_1.7.6
    ## [55] farver_2.1.0        bslib_0.3.1         ellipsis_0.3.2
    ## [58] generics_0.1.2      vctrs_0.3.8         RColorBrewer_1.1-2
    ## [61] tools_4.1.3         Biobase_2.54.0      glue_1.6.2
    ## [64] purrr_0.3.4         parallel_4.1.3      fastmap_1.1.0
    ## [67] yaml_2.3.5          colorspace_2.0-3    BiocManager_1.30.16
    ## [70] knitr_1.37          sass_0.4.0          quantreg_5.88
    ## [73] MCMCpack_1.6-1

# References

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.