From Bioconductor

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

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)

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] 13.5 11.9 11.1 2.8 15.5 ...
##   ..- 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 1 1 0 0 0 0 0 0 ...
##  $ pst        : num [1:100] 0.654 0.589 0.644 0.283 0.907 ...
##  $ k          : num [1:40, 1:2] 6.93 -8.43 -7.42 9.14 -7.49 ...
##  $ phi        : num [1:40, 1:2] 7.17 8.09 6.91 9.6 5.5 ...
##  $ delta      : num [1:40, 1:2] 0.383 0.0128 0.3851 0.4239 0.3837 ...
##  $ 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()

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.

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)

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 × 5
##   pseudotime branch branch_certainty pseudotime_lower pseudotime_upper
##        <dbl> <fct>             <dbl>            <dbl>            <dbl>
## 1     -1.17  1                     1           -1.35            -0.973
## 2     -0.854 1                     1           -0.990           -0.611
## 3     -0.366 2                     1           -0.464           -0.104
## 4      0.819 1                     1            0.462            0.979
## 5     -1.68  1                     1           -1.95            -1.55 
## 6     -0.926 1                     1           -1.14            -0.752

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

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 × 2
##   feature  chi_map
##   <chr>      <dbl>
## 1 feature1   1.54 
## 2 feature2   0.619
## 3 feature3   1.01 
## 4 feature4   0.530
## 5 feature5   1.38 
## 6 feature6   0.677

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

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 τ is just an interation-by-gene matrix:

str(m$traces$tau_trace)

##  num [1:1000, 1:40] 13.1 11.9 11.3 13.5 14.1 ...
##  - 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.984 -1.06 -1.072 -1.01 -1.081 ...

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.96 0.726 0.953 -0.97 0.967 ...

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