BACKGROUND
in a flow cytometer a number of characteristics for each individual cell in a sample of ∼  <dig> to ∼  <dig> cells are quantified as they pass through the cytometer in a fluid stream.
the data that are obtained are most often summarized by grouping cells into cell populations; properties of these cell populations are used in many clinical applications—for example monitoring hiv infection and diagnosing blood cancers—and in many branches of medical research  <cit> .
defining the cell populations based on the measured characteristics is in state-of-the-art analyses still done manually by trained operators looking at two-dimensional projections of the data.
the importance of automated methods has risen along with an increase of the dimension of typical flow cytometry data sets due to developments in flow cytometry technology  <cit>  and the emergence of studies with large numbers of flow cytometry samples  <cit> .
furthermore, manual so called gating of cell populations is a subjective process where operators have to take more or less arbitrary decisions for example when there are overlapping populations  <cit> .
automatic cell population identification is hard since flow cytometry measurements are not absolute, while at the same time different samples cannot be directly compared due to technical variation—especially apparent when samples are analyzed at different laboratories  <cit> —and intrinsic biological variation within and between subjects.
despite this, research into automated population identification methods has focused on individual or pooled flow cytometry samples, sometimes attempting to align data at first through normalization procedures  <cit> .
automated methods with the aim to replace manual gating must be able to treat multiple samples jointly and take variation between samples into account, while at the same time make it possible for the user to monitor that variation so that it is not too high for the application at hand.
for example it needs to be decided if a shift in location of a population in a sample can be seen as technical variation and accepted or if the changed marker expression means that it is a different cell phenotype.
these kinds of methods also need to be able to take prior information into account—in manual gating the experience of the operator can be necessary to define a population.
we have developed bayesflow, a method which models variation in cell population location as well as shape, can include prior information for example about cell population location, and gives a result that can be assessed in compact and comprehensive visualizations.
partitioning the cell measurements in a sample into cell populations is essentially a clustering problem.
in the context of flow cytometry data analysis clustering is called automated gating, as opposed to the manual gating performed by operators.
model-based clustering using mixture models has been the most used approach for automated gating .
mixture models are very well suited to describe flow cytometry data because they have a natural biological interpretation based on the cell populations.
examples of other approaches that have been used for automated gating are grid based density clustering  <cit> , spectral clustering  <cit> , hierarchical clustering  <cit>  and k-means clustering  <cit> .
an evaluation of a wide range of automated gating methods was performed in the flowcap i challenge  <cit> .
the discrepancy with manual gating was often quite large even for the best methods, with average f-measures around  <dig>  for both completely automated and manually tuned methods.
large discrepancies between manual and automatically gated samples can be acceptable since the arbitrary decisions taken in manual gating means that the gates could just as well have been set another way.
however, it is important that the gating is consistent between samples so that they can be compared against each other.
joint identification of cell populations in a collection of samples can be accomplished by pooling the samples  <cit>  or matching populations identified separately in the samples  <cit> .
however, in the first approach no variation between samples is taken into account and in the second approach no information is shared between samples.
recently a third approach has been explored, where a bayesian hierarchical model is used to share information between samples while at the same time allowing for variation.
this was first utilized for flow cytometry gating by cron et al.
<cit> , with a hierarchical dirichlet process model with fixed locations and shapes of cell populations.
an extension of this model, also modeling variation in cell population locations has been used to create aspire, a method for anomalous sample detection  <cit> .
bayesflow follows this third approach, but use a differently structured model than what has been used previously, favoring explicit modeling instead of implicit, parametric instead of non-parametric .
this follows the philosophy that mathematical models can never perfectly fit reality, thus it is important to be able to convey the constructed model and its parameters and in what ways it simplifies the data.
for example, in addition to variation in location bayesflow explicitly models variation in cell population shape, whereas aspire models shape variations implicitly by combining gaussian components with the same shape.
this means that an aberrant shape variation of a cell population in a sample can be detected in bayesflow by examining the parameters of the model, which is not possible in aspire.
perhaps more importantly, bayesflow gives a parsimonious model which much fewer parameters—each individual parameter for the components in bayesflow can be assessed through compact visualizations and thus undesired behaviors can be detected and corrected for by change of setup.
moreover, a restriction in aspire which is avoided by bayesflow is that the variation of component location within and between samples is connected to the shape of the components.
in bayesflow, the cells in a sample are clustered using a multivariate gaussian mixture model , where k components describe true and artificial cell populations and one component describes outliers.
artificial cell populations are measurements that cluster together and behave otherwise like real cell populations, but arise for example from dead cells, non-specific binding of markers or doublets; doublets are pairs or groups of cells that pass through the flow cytometer at the same time.
measurements which are not clearly grouped but spread out over the measurement space, for example due to measurement noise, are modeled as outliers.
for each component not representing outliers its mean and covariance matrix is linked to a latent cluster which collects corresponding components across all samples.
in practice this is done by assuming a normal prior for the means and an inverse wishart prior for the covariance matrices of the components linked to a given latent cluster.
the parameters of sample and latent components are jointly estimated by markov chain monte carlo  sampling.
the variation in location and shape between corresponding mixture components across samples is controlled by the priors on parameters of the latent clusters.
the location of component means and shape of components can also be restricted if there is prior information supporting this.
to allow for that flow cytometry data frequently have missing cell populations, we include the possibility that not all components are present in every sample.
a challenge that has to be addressed when analyzing flow cytometry data is that cell populations can be skewed and/or have heavy tails and are then not well described by a single gaussian component  <cit> .
to handle this we use multiple components to model such populations, an approach that have often been employed for flow cytometry data  <cit>  and has the further advantage that the number of cell populations can be automatically detected.
we merge gaussian components into super components with a procedure based on a systematic study of methods for merging mixture components  <cit> .
results from the mcmc sampling and subsequent merging are evaluated in a number of quality tests.
this is a crucial step since what is deemed as a good clustering is application dependent.
in some settings a given amount of variation in location or shape is expected from biological or technical reasons, whereas in others the same variation would indicate a different population.
this also means that it is necessary for the user to choose prior parameters for their application.
to simplify this process we have derived parametrizations so that the same value of the parameters gives a similar effect of the prior on data sets of different sizes.
we verified the ability of the sampling scheme to recover model parameters by fitting the model to a small three-dimensional synthetic data set with  <dig>  million cells in total and a large synthetic data set with in total  <dig> million cells in  <dig> dimensions.
then we applied bayesflow to one of the datasets in the flowcap i challenge, the gvhd dataset, which contains samples from patients who have had organ transplants and might have early signs of graft-versus-host disease.
we show that bayesflow does not only give a result which has the same degree of accordance with manual gating as the best performing methods in flowcap i—which is much higher than what is obtained for other methods based on joint gating with bayesian hierarchical models—it does also give a more similar treatment of different samples than manual gating and the best methods from flowcap i. finally we applied bayesflow, aspire  <cit>  and hdpgmm  <cit>  to a data set with replicated samples from four healthy individuals.
the ratio between intra-donor technical variation and inter-donor biological variation was similar between bayesflow and hdpgmm, which was lower than for aspire.
moreover, bayesflow was the only of the three methods which gave cell populations with clear expression patterns.
methods
model
let yij denote vector valued measurement number i in sample j. here i∈{ <dig> …,nj}, where nj is the number of cells in sample j, and j∈{ <dig> …,j}, where j is the number of samples.
we let the dimension of the observations be denoted d. with k mixture components describing cell populations the probability density for cell measurement i of a flow cytometry sample j is modeled as
f=∑k=1kπjknyij;μjk,Σjk+πj0nyij;μj <dig> Σj <dig>
where n denotes the probability density function of the normal distribution with mean μ and covariance matrix Σ evaluated at y. to facilitate interpretation, the number k should be chosen as small as possible, given that the model pass quality requirements .
the last component represents outliers and its parameters μj0=μ <dig> and Σj0=Σ <dig> are identical across samples.
outliers are often modeled by a uniform density over the measurement space  <cit> ; however due to the curse of dimensionality  <cit> , this is not well behaved when we have more than a few dimensions, in which case a gaussian should perform better.
noise coming from for example dead cells can also be captured in artificial cell populations, and can be excluded in downstream analyses based on the expression patterns.
the vector πj={πj <dig> …,πjk} contains the mixing proportions, i.e.
the proportion of cells described by the component.
to connect cell populations between samples we use a latent layer, assuming that for a given k each μjk and Σjk is drawn from a normal and an inverse wishart distribution respectively.
specifically, in our model, for k= <dig> …,k,
μjk|θk,Σθk∼n,Σjk|Ψk,νk∼iw
where θk, Σθk, Ψk and νk are hyper-parameters describing latent cluster k. these parameters describes the variability between flow cytometry samples, in contrast to μjk,Σjk which describe the distribution of cell measurements within a sample.
the normal and inverse wishart distributions are conjugate priors to the mean and the covariance respectively of the normal distribution, enabling efficient sampling, however they are not jointly conjugate.
we call θk and Ψk/ the latent cluster mean and latent cluster covariance matrix respectively, since they are the a priori expected values of μjk and Σjk.
for the hyper-parameters describing the latent clusters and the mixing proportions we use the following prior distributions:
θk|tk,sk∼n,πj∼d,Σθk|qk,nθk∼iw,νk|λk∼exp,Ψk|hk,nΨk∼w,
where w denotes the wishart distribution and d denotes the dirichlet distribution, which is conjugate prior to the multinomial distribution.
for each νk we assign a exponential prior on the positive natural numbers.
the complete structure of the model is displayed through a directed acyclic graph  in fig.
<dig>
fig.
<dig> directed acyclic graph describing the bayesian hierarchical model.
square boxes indicate that the values are known
the parameters tk and sk define the prior belief of the locations of the latent means θk, whereas the parameters qk and nθk control the spread of mixture component means within a latent cluster and are hence important to control the variation across samples.
a large nθk along with a small qk forces the μjk together; it makes large deviations between Σθk and qk unlikely.
the parameters hk and nΨk control the expected values and the variation of latent covariance matrices as well as the variation among mixture component covariance matrices in a latent cluster.
if nΨk is large each Σjk will be close to Ψk/ for any k, since a high nΨk makes high νk more probable.
finally, to simplify sampling from the posterior distribution of the parameters, we add an component assignment variable xij∈{ <dig> ,…,k} describing which component yij is drawn from.
to comply with , the a priori uncertainty of component membership is modeled by xij∼ mult, where mult denotes the multinomial distribution.
the resulting posterior distribution of all the parameters, denoted jointly by Θ, and x given the data y is given in the additional file 1: section a. in section b we describe the markov chain monte carlo  sampling scheme used to generate posteriors for our model parameters.
the computational bottleneck of the sampling scheme is the sampling of x, with a computational complexity bounded by o. to handle high dimensions diagonal covariance matrices can be used instead, in which case the complexity is bounded by o. however, for datasets with more than  <dig> dimensions the mathematical feasibility of using gaussian mixture models without any prior dimension reduction needs to be seriously considered first, due to the curse of dimensionality  <cit> .
absent components
in some flow cytometry data sets not all cell populations are present in all samples.
in our model this corresponds to that πjk= <dig> for some .
however, mixture component parameters for empty clusters will still affect the mixing of the mcmc for the parameters of the latent cluster.
it can also happen that if a cluster is empty that the mixture component moves and split a neighboring cluster in two.
to avoid this in such data sets we extend the model by introducing a variable zj∈{ <dig> }k that says which components are active in sample j. this has the further advantage that when sampling from the posterior distribution of the model we get the probability for each cluster that it is present in a sample.
we impose a prior on zj which is proportional to exp−cs∑k=1kzji∑k=1kzj> <dig> where i denotes the indicator function and cs> <dig>  the prior makes the model prefer fewer activated clusters so that if there is a very small cluster the likelihood will be larger if it is inactivated, which prevents spurious clusters.
the strength of this prior can be adjusted to the expected size of the smallest clusters.
the changes to  required by this extension are straightforward but inference of the model becomes a bit more involved since removing components reduces the dimension of the model.
to accommodate for this we have included a reversible jump step in our sampling algorithm.
details are given in the additional file 1: section b.
merging latent clusters
to determine the “correct” number of clusters in a data set directly from the data is an ill-defined problem, since what should be considered to be a separate cluster depends on the interpretation of the data.
nevertheless, there are many different criteria which can be used to guide the decision about the number of populations  <cit> .
we use overlap between components—measured by bhattacharyya distance—and unimodality of the resulting super clusters—measured by hartigan’s dip test  <cit> —to determine which latent clusters to merge and to indicate our confidence in the mergers.
in an evaluation of criteria for merging gaussian components to represent more complex distributions, the bhattacharyya distance performed well  <cit> .
bhattacharyya distance merges clusters according to a pattern-based cluster concept as opposed to a modality-based concept  <cit> .
with a pattern-based cluster concept a small dense cluster inside a sparse cluster—for example a well specified cell population inside a region with sparse outliers—will be considered to be different clusters.
this would not be the case for the modality-based cluster concept as long as the generating probability density is unimodal.
the bhattacharyya distance between n and n is
dbhat=1/8·⊤Σ¯−1+1/2·log|Σ¯|/|Σ2||Σ2|,
where Σ¯=/ <dig>  <cit> .
in order to measure bhattacharyya distance between mixtures of gaussian distributions, which is necessary for deciding if super clusters should be merged with other clusters, we approximate each mixture with a gaussian distribution.
the means and the covariance matrices are estimated using a soft clustering of the data points inferred from the sampling of xij, detailed in the additional file 1: section c.
however, it is not obvious how to set a threshold for dbhat, since the appropriate threshold depends on the distribution of the data  <cit> , which is unknown.
because of this we use a low soft threshold d <dig> and a high hard threshold d <dig>  two clusters closer to each other than d <dig> are always merged, two clusters whose distance is between d <dig> and d <dig> are only merged if they fulfill an additional criterion based on hartigan’s dip test for unimodality.
unimodality is an appealing heuristic for defining cell populations, and it has frequently been used for automated gating  <cit> .
it has two main limitations.
the first one, that populations intuitively should be separate if they have very different densities—even when they overlap so that their combined distribution is unimodal—can be bypassed by combining unimodality with a pattern-based merging criterion such as bhattacharyya distance.
the second one, that it is difficult to determine if a multi-dimensional empirical distribution is multimodal, is usually handled by considering one-dimensional projections  <cit> .
this is the approach we take here, using hartigan’s dip test of unimodality for each of the projections onto the coordinate axes where bhattacharyya overlap is low, and for the projection onto fisher’s discriminant coordinate.
if for a proposed merger, any of these projections is found to be multimodal, the clusters are not merged.
further details of the merging procedure are given in the additional file 1: section c.
quality control
to verify that the output of bayesflow fulfills the user’s requirements, a number of checks are performed:
convergence of the mcmc sampler is established by viewing trace plots of sampled parameters.
to ensure that variation of the two different populations are not confused with each other, we require that the bhattacharyya distance as well as the euclidean distance from each sample component to its corresponding latent component should be smaller than these distances to any other latent component which does not belong to the same super cluster.
to ensure that the obtained clusters should not be divided further, hartigan’s dip test is computed for the projections onto the coordinate axes of all super clusters.
projections which have a dip test p-value below  <dig> —the threshold for merging components —are visualized using histograms of quantiles of the weighted data belonging to the cluster.
to ensure that the model fits the data reasonably well, samples from the posterior predictive is compared to the true data in one- and two-dimensional histograms.
to ensure that there are no outliers among the cluster centers, the centers for each cluster are plotted together along one dimension.
additionally, to detect components with aberrant shapes, the eigenvectors corresponding to the largest eigenvalues, multiplied with the corresponding eigenvalues, can be viewed.
if any of the quality criteria is not met, the simulation should be rerun, either using the same or different parameters.
even if the same parameters are used a different result can be obtained due to randomness in the initialization.
experiments
simulated data
in order to verify that the proposed sampling scheme can find the correct model parameters, the mcmc algorithm was applied to two simulated datasets.
the first dataset was three-dimensional, which enables direct visual evaluation.
it had four latent clusters across eighty artificial flow cytometry samples; each sample had  <dig>  cells giving a total of  <dig>  million cells.
one of the latent clusters was present only in eight samples and another one was present in  <dig> samples, so that the ability to find rare cell populations was tested.
moreover, the cluster which was present in only eight samples contained only  <dig> % of the total number of cells, thus also the ability to find small cell populations was tested.
the parameters and the algorithm used for generating the data are given in the additional file 1: section d. <dig>
the second data set was designed to test the ability to handle large data.
it was eight-dimensional, with eleven latent clusters and  <dig> artificial flow cytometry samples.
each sample had measurements of  <dig>  cells, giving a total of  <dig> million cells.
four of the eleven clusters were missing in half of the samples.
additional file  <dig> contains a python script and data for regenerating this data set.
prior parameters and initial values for the mcmc sampler are given in the additional file 1: section d. <dig>  all priors were chosen to be non-informative.
the outlier component was not used for inference in the small dataset, but it was used for the large dataset.
the mcmc sampler ran first for a number of burn-in iterations, then the posterior distribution was explored in a number of production iterations.
during the production iterations, apart from sampling parameters of the model, a value of y was also drawn, i.e.
a sample from the posterior predictive.
for the first synthetic data set  <dig>  burn-in and  <dig>  production iterations were used.
for the second, larger, data set we used  <dig>  burn-in iterations and  <dig>  production iterations.
for the second data set the mcmc sampler was run on amazon cloud, using  <dig> cores.
each iteration took on average one second, so that about  <dig>  h was needed in total.
flow cytometry data
we analyze two flow cytometry data sets with bayesflow: the data set gvhd from the flowcap i challenge—with four markers,  <dig> samples and approximately  <dig>  cells per sample—and a data set obtained from the r package healthyflowdata  <cit>  with technical replicates of pbmc samples from healthy donors—in total  <dig> samples with approximately  <dig>  cells, also measured with four markers.
in the gvhd dataset we can compare the gating obtained from bayesflow with manual gating provided from flowcap as well as automated gating from a wide range of other methods.
in healthyflowdata we can instead compare gating between technical replicates to see if samples are treated in a consistent manner.
for the healtyflowdata dataset we used an exploratory approach with non-informative priors.
we ran multiple simulations and gradually increased the number of components until we passed the quality criteria described under quality control; we finally arrived at using k= <dig> components.
for the gvhd data set we started with an exploratory approach and gradually increased the number of components, but in the quality checks we noted one population in one of the samples which was very hard to capture.
then we decided to use an informative approach for this population.
using a scatter plot, fig.
<dig>  we set boundaries for this population in the dimensions given by the cd <dig> and the cd8b marker and computed its mean and empirical s covariance matrix.
we used the mean to set an informative prior for θk and the mean and the empirical covariance to initialize the component.
prior parameters in both the informative and non-informative case are described in the additional file 1: section e. <dig>
bayesflow applies three data preprocessing steps: 1) data points with extreme values in at least one dimension  are removed.
such data points can lead to components with singular covariance matrices, and a well designed flow cytometry experiment should not have significant populations with such values.
2) the data is scaled using the  <dig> and  <dig> % percentiles q <dig>  and q <dig>  of the pooled data, with the same scaling for all samples, so that q <dig> = <dig> and q <dig> = <dig> for each marker for the pooled data.
this is done in order to be able to set informative priors in an intuitive way.
3) before testing which components should be merged, a very small amount of noise is added to the data .
this is since the discreteness of the original flow cytometry measurements can lead to a striped pattern in the flow cytometry data  <cit>  and also when it is not visible to the human eye it disturbs the dip test.
after preprocessing, parameters for the mcmc sampler were initialized by running the em algorithm on the pooled data, followed by the initialization scheme used for the large synthetic dataset, detailed in the additional file 1: section d. <dig>  we ran  <dig>  burn-in iterations and  <dig>  production iterations of the mcmc sampler for both experiments.
the burn-in period consisted of five phases: in the first phase, the priors on variation in location and shape were modified to force clusters together.
before the second phase, priors parameters were set to normal again.
after the second phase, components which were considered to be outliers were turned off.
they were forced to stay off during a short third phase, but from the forth phase and onwards components were allowed to turn on and off.
label switching was allowed during the initial four phases in order to escape non-desired local minima, but then disallowed.
the values of parameters controlling the simulation during the burn-in and production period are given in additional file 1: table s <dig>
we also applied the two other joint gating methods based on bayesian hierarchical models: aspire  <cit>  and hdpgmm  <cit> .
for aspire parameters were chosen according to the strategy recommended by dundar et al.
<cit> ; details are given in the additional file 1: section e. <dig>  for each run we used in total  <dig>  iterations, of which  <dig>  were set as burn in iterations.
for hdpgmm default parameters were used, with a burn-in period of  <dig>  iterations and a production period of  <dig> iterations.
we ran bayesflow and aspire on a  <dig>  ghz quad core cpu.
a bayesflow run took  <dig>  h for the gvhd dataset and  <dig>  h for healthyflowdata.
aspire took in total  <dig>  h for the gvhd dataset and  <dig>  h for healthyflowdata per run.
four runs of aspire was needed to determine the κi parameters.
hdpgmm was run on a dual core gpu.
it needed  <dig>  h for the gvhd dataset and approximately  <dig> h for the healthyflowdata dataset.
RESULTS
simulated data
we begin by analyzing the smaller data set.
in fig.
<dig> we show univariate and bivariate histograms of all synthetic cell measurements pooled together, as well as the corresponding histograms of the data from a single flow cytometry sample where all four clusters are present.
note that the data when pooled together has a complicated density, as it is in fact a mixture of  <dig> multivariate normal densities.
fig.
<dig>
a one and two dimensional histograms for one synthetic flow cytometry sample containing  <dig>  data points; b histograms of  <dig>  data points drawn uniformly from the pooled data from the synthetic data experiment
in fig.
<dig> we show the same univariate and bivariate histograms, but this time with samples from the posterior predictive distribution of y. from the synthetic cell measurements generated from the inferred models of the datasets it is clear that the inferred models are accurate and capture the variation across samples, which a model only of pooled data cannot do.
fig.
<dig>
a one and two dimensional histograms of  <dig>  posterior draws of y for the flow cytometry sample displayed in fig.
2
a; b histograms of  <dig>  posterior draws of y drawn uniformly from all the flow cytometry samples, thus matching fig.
2
b
fig.
<dig> bayesflow component parameter representations of inferred latent clusters  and mixture components  together with histograms of real data  and synthetic data generated from the model  for healthyflowdata.
the center of each ellipse is the mean and each semi-axis is an eigenvector with length given by the corresponding eigenvalue of the projected covariance matrix.
for the latent clusters the parameters Ψk) are shown, for the mixture components the parameters  are shown.
each component or cluster is depicted with the same color as in fig.
5; different shades of same color corresponds to latent clusters that have been merged
fig.
<dig> summary statistics of inferred cell populations in bayesflow, aspire and hdpgmm, ordered by population size.
for hdpgmm, the six largest components after merging are shown, the remaining components have together at most  <dig>  of the cells in a sample.
the noise component in bayesflow has at most  <dig>  of the cells in a sample.
a locations μ
jk of mixture components that represent each population, in each sample, cf.
fig.
<dig>  b box plots of the soft clusters in the pooled data, cf.
fig.
<dig>  c population proportions across flow cytometry samples
fig.
<dig> cell population which is hard to detect in the gvhd dataset
fig.
<dig> the posterior mean of the mixture component centers, μ
jk , and the true cluster centers  in the small synthetic data experiment
finally in figs.
<dig> and  <dig>  the marginal posterior distributions of the latent cluster parameters θk and Ψk, subtracted by their true values, are presented.
in fig.
<dig> the dot represents the difference between the median of posterior distribution and the true value of each θk.
the vertical lines represent the  <dig>  and  <dig>  % quantiles.
fig.
<dig> displays results for each latent covariance matrix Ψk/ in the same way.
from figs.
<dig> and  <dig> we see that the true parameters of both the means and the covariances are all between the  <dig>  and  <dig>  % quantiles of the posterior distribution.
fig.
<dig> the difference between the true value of each entry in each θ
k and the approximated marginal posterior distribution generated by the mcmc sampler in the small synthetic data experiment.
the black dot represents the median and the vertical line goes between the  <dig>  and  <dig>  % quantiles.
the light gray horizontal line is the  <dig> line
fig.
<dig> the difference between the true value of each of the entries in Ψ
k/ and the approximated marginal posterior distribution generated by the mcmc sampler in the synthetic data experiment.
the black dot shows the median, and the black vertical line goes between the  <dig>  and  <dig>  % quantiles.
the light gray horizontal line is the  <dig> line
the true and estimated cluster centers of the eight-dimensional data set cannot be displayed efficiently with just three dimensions at hand, but a three-dimensional projection is shown in fig.
<dig>  the average error in euclidean distance in the full eight-dimensional space is  <dig> , which can be compared to the average error had the latent mean across samples been used, namely  <dig> , which is the best that could have been obtained from a model not including variation between samples.
the outlier component was used for inference in the results presented here, but omitting it has very small effect.
fig.
<dig> the posterior mean of the mixture component centers, μ
jk , and the true cluster centers  in the large synthetic data experiment for the first three dimensions
in fig.
<dig>  we show the posterior distribution of the latent cluster means where again the dot represents the difference between the median of posterior distribution and the true value of each θk.
the vertical lines are the  <dig>  and  <dig>  % quantiles.
the posterior samples have been divided by the standard deviation of the true θk so that the scales across the clusters are equal.
some of the credibility intervals do not contain zero, but this is explained when studying the intervals that would have been obtained if the true μk were used , since they are almost identical.
fig.
<dig> the difference between the true value of each entry in each θ
k and the approximated marginal posterior distribution generated by the mcmc sampler in the large synthetic data experiment.
the black dot represents the median and the vertical line goes between the  <dig>  and  <dig>  % quantiles.
to get the axis on the same scale for all the clusters, they are scaled by the standard deviation of μ
k. the light gray horizontal line is the  <dig> line.
the red dot and lines is the same however where one uses the true μ
k to estimate θ
k, rather then the μ
k obtained by taking the posterior means of the mixtures
we thus see that cluster centers and credibility intervals for latent clusters are captured well in both synthetic data sets.
flow cytometry data
gvhd
for the analysis of the gvhd dataset we did twelve runs of bayesflow in the informed setup described above.
seven were excluded due to confusion between populations, i.e.
at least one sample component was closest to the wrong latent component; of the remaining five, one more run was excluded since it has not converged, and another two because of multimodal clusters.
this leaved two runs that passed the quality control.
additional file 1: figs.
s <dig> and s <dig> show trace plots and projections of clusters with high dip test values respectively.
table  <dig> accordance with manual gating for gvhd data set.
for hdpgmm we also report the result when components are merged according to our merging procedure.
when this procedure is applied to the results obtained by aspire, no components are merged, i.e.
the original result is identical to what is obtained after mergeing
flowmeans
<dig>
<dig>
<dig>
samspectral
<dig>
<dig>
<dig>
ensemble flowcap
<dig>
one of the two bayesflow runs has the highest accordance with manual gating, the other one is on par with flowmeans and samspectral, which is considerably higher than aspire and hdpgmm.
however, as can be seen in fig.
<dig>  the gating of different samples is arguably most consistent for bayesflow as compared to manual gating, flowmeans and samspectral.
fig.
<dig> gated events according to four methods  of the twelve samples in the gvhd dataset, projected onto the two first dimensions.
for bayesflow, the run with least accordance with manual gating, run  <dig>  is shown.
similar plots for aspire and hdpgmm as well as bayesflow run  <dig> are shown in the additional file 1: figure s6
to get a further understanding of the variability between samples in bayesflow, summary statistics for the obtained components and cell populations are shown in fig.
<dig>
fig.
<dig> summary statistics of the six cell populations obtained by bayesflow  in the dataset gvhd.
the outlier component has at most  <dig>  of the cells in a sample.
a each panel displays the locations μ
jk of all mixture components that represent the population, across all samples.
different shades of a color represent different latent components k. b box plots of the soft clusters in the pooled data.
the boxes go between the quantiles q
km, <dig>  and q
km, <dig> , the whiskers extend to q
km, <dig>  and q
km, <dig> .
the α-quantile for  component k in dimension m, q
km,α, is here defined as qkm,α=mini′j′{yi′j′m:α<∑ij:yijm<yi′j′mwijk}.
c population proportions in each of the twelve flow cytometry samples
healthyflowdata
we did  <dig> runs of bayesflow with k= <dig>  ten of these were excluded due to confusion between populations, moreover two runs were excluded since they had clusters with clearly multimodal distributions.
for the six runs that passed the quality control, 3– <dig> components were turned off across all samples; they are excluded from visualizations.
additional file 1: figs.
s <dig>  s <dig> and s <dig> show trace plots, projections of clusters with high dip test values and eigenvectors of covariance matrices respectively.
in fig.
<dig> we visualize model fit and inter-sample variation for the first of the six runs that passed the quality control by plotting latent and sample components as well as histograms of real data and synthetic data generated from the model, for two different samples and for the pooled data.
we can thus see how shape variations are captured by the model.
the output of bayesflow, aspire and hdpgmm can be compared in fig.
<dig>  the merging procedure we used for bayesflow has been applied for both aspire and hdpgmm, however for aspire no components were merged by this.
in bayesflow each of the populations correspond to clear expression patterns, which is not the case for the other methods.
for example the first population is clearly cd4+cd8- t-cells whereas for both aspire and hdpgmm this population contains both components which are cd8- and components which are cd8+.
we also compare intra-donor variation of cell population size to inter-donor variation for the six bayesflow runs, as well as for aspire and hdpgmm in fig.
<dig>  for aspire there are inter-donor distances which are clearly smaller than some intra-donor distances, which is not the case for bayesflow and hdpgmm.
fig.
<dig> distances within  and between  donors as measured by ℓ
<dig> distance between vectors of population sizes.
for the six bayesflow runs and hdpgmm there is very little or no overlap between within-donor and between-donor distances, whereas for aspire there is clear overlap
discussion
from different runs of bayesflow we can get different representations of data, as in the case of the gvhd dataset.
this is because with highly overlapping populations there might be multiple models representing the data equally well.
but since all samples are gated jointly in every run, the gated populations can still be compared across samples.
the user might have a preference for one representation or the other though, and informative priors can be used to guide bayesflow to a preferred representation.
bayesflow is not aimed at discovery of rare cell populations, but it can be used together with an algorithm specifically designed for detecting rare cell populations in a sample, such as swift  <cit> , and then use informative priors to find how this population occurs across an entire set of samples, in a similar way as was done in the gvhd dataset.
how much clusters should be merged is a decision that needs to be taken by the interpreter of the data.
in some settings one might want to be restrictive with merging and then use higher thresholds.
in others one might want additional mergers after viewing joint one-dimensional projections of the clusters.
the bayesflow pipeline does not in itself include any compensation or any of the non-linear transformations which are often used for flow cytometry data, such as logicle.
compensation is a linear transformation and gaussian mixture models are invariant under linear transformations, so they perform equally well on uncompensated and compensated data.
non-linear transformations such as logicle can make gaussian populations non-gaussian, which makes inference harder.
the flow cytometry data we used for the experiments had already been compensated, the healthyflowdata data set had also been transformed with an asinh transform; details are given in the additional file 1: section e. <dig>
bayesflow finds a joint representation of an entire set of samples.
in order for this representation to be reasonable there has to be sufficient correspondences between samples.
even if for a data set with very little correspondences a joint model could be obtained by using a very large number of components, it would hard to gain any insights from such a model.
in such a case an entirely computational pipeline without the cell population identification step would be preferred.
bayesflow can be computationally intensive if many runs are needed to pass the quality control.
for these cases it would be desirable to complement bayesflow e.g.
with initialization methods that would allow passing the quality control more often, so that few runs in bayesflow would be needed.
fast initialization methods and early quality checks aiming at this would therefore be of interest for the community and is something that we propose for further study.
CONCLUSIONS
in this paper we have presented a new bayesian hierarchical model designed for joint cell population identification in many flow cytometry samples.
the model captures the variability in shapes and locations of the populations between the samples and we have demonstrated its use in an exploratory as well as in a partly informed setting with some prior information.
we showed that for synthetic datasets generated from the model, the parameters were recovered with high accuracy through a mcmc sampling scheme.
the model was then applied to a real flow cytometry data set where a manual gating was available, and it was shown to have very high accordance with manual gating as compared to other automated gating methods, while at the same time the gating was more consistent across samples than either the manual gating or other automated gating methods.
when applied to another flow cytometry data set with technical replicates of blood from healthy donors, bayesflow gave a parsimonious representation of the data, which enables visualization and monitoring of its parameters.
the obtained cell populations had clear expression patterns as opposed to the clusters obtained by aspire and hdpgmm, where for example cd4+cd8- t-cells where in the same cluster as cd4+cd8+ t-cells.
the population sizes obtained by bayesflow and hdpgmm respectively had lower intra-donor variation compared to inter-donor variation than what was obtained from aspire.
many approaches of automated gating of multiple flow cytometry samples in parallel have been aimed at finding features of the data so that either samples can be classified into groups, e.g.
cancer or normal, or they can be used to predict an outcome such as expected time to progression of disease.
features are often designed based on characteristics of cell populations, but usually not so much attention has been given to ensure that they represent actual cell populations.
bayesflow takes the opposite approach and gives a representation of the data according to cell populations, with the same cell populations across the entire set of samples .
the advantages to this approach are among others that the result is directly biologically interpretable and that a rich output is given which can be explored in many different ways which are familiar to someone who is used to manual gating.
in this way we can join the objectivity and ability to work in high dimensions and with many samples of automated gating with the flexibility in interpretation of manual gating.
additional files
additional file  <dig>
supplementary material.
the supplementary material contains the posterior in bayesflow, the mcmc sampling scheme, additional details on the merging of components, information about the data generation, priors and initialization for the synthetic data example; parameters used for aspire, additional details on healthyflowdata, the priors and the initialization procedure used when studying this data set and further results pertaining to the real flow cytometry data set, including fitting gaussian mixture models to individual samples of healthyflowdata with the em algorithm and scatter plots of gvhd for aspire, hdpgmm and bayesflow run  <dig>
additional file  <dig>
data generation files.
a python script for generating the large synthetic dataset, along with means, covariances and weights needed for this.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
kj, jw and mf conceived and planned the study.
jw and kj designed the statistical model and the inference procedure.
kj, jw and mf designed the experiments.
jw and kj implemented bayesflow and ran the experiments.
kj and jw wrote the article with the help of mf.
all authors read and approved the final version of the manuscript.

 in a flow cytometer a number of characteristics for each individual cell in a sample of ∼ to ∼ cells are quantified as they pass through the cytometer in a fluid stream. the data that are obtained are most often summarized by grouping cells into cell populations; properties of these cell populations are used in many clinical applications—for example monitoring hiv infection and diagnosing blood cancers—and in many branches of medical research. defining the cell populations based on the measured characteristics is in state-of-the-art analyses still done manually by trained operators looking at two-dimensional projections of the data. the importance of automated methods has risen along with an increase of the dimension of typical flow cytometry data sets due to developments in flow cytometry technology and the emergence of studies with large numbers of flow cytometry samples. furthermore, manual so called gating of cell populations is a subjective process where operators have to take more or less arbitrary decisions for example when there are overlapping populations. automatic cell population identification is hard since flow cytometry measurements are not absolute, while at the same time different samples cannot be directly compared due to technical variation—especially apparent when samples are analyzed at different laboratories —and intrinsic biological variation within and between subjects. despite this, research into automated population identification methods has focused on individual or pooled flow cytometry samples, sometimes attempting to align data at first through normalization procedures. automated methods with the aim to replace manual gating must be able to treat multiple samples jointly and take variation between samples into account, while at the same time make it possible for the user to monitor that variation so that it is not too high for the application at hand. for example it needs to be decided if a shift in location of a population in a sample can be seen as technical variation and accepted or if the changed marker expression means that it is a different cell phenotype. these kinds of methods also need to be able to take prior information into account—in manual gating the experience of the operator can be necessary to define a population. we have developed bayesflow, a method which models variation in cell population location as well as shape, can include prior information for example about cell population location, and gives a result that can be assessed in compact and comprehensive visualizations. partitioning the cell measurements in a sample into cell populations is essentially a clustering problem. in the context of flow cytometry data analysis clustering is called automated gating, as opposed to the manual gating performed by operators. model-based clustering using mixture models has been the most used approach for automated gating. mixture models are very well suited to describe flow cytometry data because they have a natural biological interpretation based on the cell populations. examples of other approaches that have been used for automated gating are grid based density clustering, spectral clustering, hierarchical clustering and k-means clustering. an evaluation of a wide range of automated gating methods was performed in the flowcap i challenge. the discrepancy with manual gating was often quite large even for the best methods, with average f-measures around for both completely automated and manually tuned methods. large discrepancies between manual and automatically gated samples can be acceptable since the arbitrary decisions taken in manual gating means that the gates could just as well have been set another way. however, it is important that the gating is consistent between samples so that they can be compared against each other. joint identification of cell populations in a collection of samples can be accomplished by pooling the samples or matching populations identified separately in the samples. however, in the first approach no variation between samples is taken into account and in the second approach no information is shared between samples. recently a third approach has been explored, where a bayesian hierarchical model is used to share information between samples while at the same time allowing for variation. this was first utilized for flow cytometry gating by cron et al., with a hierarchical dirichlet process model with fixed locations and shapes of cell populations. an extension of this model, also modeling variation in cell population locations has been used to create aspire, a method for anomalous sample detection. bayesflow follows this third approach, but use a differently structured model than what has been used previously, favoring explicit modeling instead of implicit, parametric instead of non-parametric. this follows the philosophy that mathematical models can never perfectly fit reality, thus it is important to be able to convey the constructed model and its parameters and in what ways it simplifies the data. for example, in addition to variation in location bayesflow explicitly models variation in cell population shape, whereas aspire models shape variations implicitly by combining gaussian components with the same shape. this means that an aberrant shape variation of a cell population in a sample can be detected in bayesflow by examining the parameters of the model, which is not possible in aspire. perhaps more importantly, bayesflow gives a parsimonious model which much fewer parameters—each individual parameter for the components in bayesflow can be assessed through compact visualizations and thus undesired behaviors can be detected and corrected for by change of setup. moreover, a restriction in aspire which is avoided by bayesflow is that the variation of component location within and between samples is connected to the shape of the components. in bayesflow, the cells in a sample are clustered using a multivariate gaussian mixture model, where k components describe true and artificial cell populations and one component describes outliers. artificial cell populations are measurements that cluster together and behave otherwise like real cell populations, but arise for example from dead cells, non-specific binding of markers or doublets; doublets are pairs or groups of cells that pass through the flow cytometer at the same time. measurements which are not clearly grouped but spread out over the measurement space, for example due to measurement noise, are modeled as outliers. for each component not representing outliers its mean and covariance matrix is linked to a latent cluster which collects corresponding components across all samples. in practice this is done by assuming a normal prior for the means and an inverse wishart prior for the covariance matrices of the components linked to a given latent cluster. the parameters of sample and latent components are jointly estimated by markov chain monte carlo sampling. the variation in location and shape between corresponding mixture components across samples is controlled by the priors on parameters of the latent clusters. the location of component means and shape of components can also be restricted if there is prior information supporting this. to allow for that flow cytometry data frequently have missing cell populations, we include the possibility that not all components are present in every sample. a challenge that has to be addressed when analyzing flow cytometry data is that cell populations can be skewed and/or have heavy tails and are then not well described by a single gaussian component. to handle this we use multiple components to model such populations, an approach that have often been employed for flow cytometry data and has the further advantage that the number of cell populations can be automatically detected. we merge gaussian components into super components with a procedure based on a systematic study of methods for merging mixture components. results from the mcmc sampling and subsequent merging are evaluated in a number of quality tests. this is a crucial step since what is deemed as a good clustering is application dependent. in some settings a given amount of variation in location or shape is expected from biological or technical reasons, whereas in others the same variation would indicate a different population. this also means that it is necessary for the user to choose prior parameters for their application. to simplify this process we have derived parametrizations so that the same value of the parameters gives a similar effect of the prior on data sets of different sizes. we verified the ability of the sampling scheme to recover model parameters by fitting the model to a small three-dimensional synthetic data set with million cells in total and a large synthetic data set with in total million cells in dimensions. then we applied bayesflow to one of the datasets in the flowcap i challenge, the gvhd dataset, which contains samples from patients who have had organ transplants and might have early signs of graft-versus-host disease. we show that bayesflow does not only give a result which has the same degree of accordance with manual gating as the best performing methods in flowcap i—which is much higher than what is obtained for other methods based on joint gating with bayesian hierarchical models—it does also give a more similar treatment of different samples than manual gating and the best methods from flowcap i. finally we applied bayesflow, aspire and hdpgmm to a data set with replicated samples from four healthy individuals. the ratio between intra-donor technical variation and inter-donor biological variation was similar between bayesflow and hdpgmm, which was lower than for aspire. moreover, bayesflow was the only of the three methods which gave cell populations with clear expression patterns. methods model let yij denote vector valued measurement number i in sample j. here i∈{ …,nj}, where nj is the number of cells in sample j, and j∈{ …,j}, where j is the number of samples. we let the dimension of the observations be denoted d. with k mixture components describing cell populations the probability density for cell measurement i of a flow cytometry sample j is modeled as f=∑k=1kπjknyij;μjk,Σjk+πj0nyij;μj Σj where n denotes the probability density function of the normal distribution with mean μ and covariance matrix Σ evaluated at y. to facilitate interpretation, the number k should be chosen as small as possible, given that the model pass quality requirements. the last component represents outliers and its parameters μj0=μ and Σj0=Σ are identical across samples. outliers are often modeled by a uniform density over the measurement space; however due to the curse of dimensionality, this is not well behaved when we have more than a few dimensions, in which case a gaussian should perform better. noise coming from for example dead cells can also be captured in artificial cell populations, and can be excluded in downstream analyses based on the expression patterns. the vector πj={πj …,πjk} contains the mixing proportions, i.e. the proportion of cells described by the component. to connect cell populations between samples we use a latent layer, assuming that for a given k each μjk and Σjk is drawn from a normal and an inverse wishart distribution respectively. specifically, in our model, for k= …,k, μjk|θk,Σθk∼n,Σjk|Ψk,νk∼iw where θk, Σθk, Ψk and νk are hyper-parameters describing latent cluster k. these parameters describes the variability between flow cytometry samples, in contrast to μjk,Σjk which describe the distribution of cell measurements within a sample. the normal and inverse wishart distributions are conjugate priors to the mean and the covariance respectively of the normal distribution, enabling efficient sampling, however they are not jointly conjugate. we call θk and Ψk/ the latent cluster mean and latent cluster covariance matrix respectively, since they are the a priori expected values of μjk and Σjk. for the hyper-parameters describing the latent clusters and the mixing proportions we use the following prior distributions: θk|tk,sk∼n,πj∼d,Σθk|qk,nθk∼iw,νk|λk∼exp,Ψk|hk,nΨk∼w, where w denotes the wishart distribution and d denotes the dirichlet distribution, which is conjugate prior to the multinomial distribution. for each νk we assign a exponential prior on the positive natural numbers. the complete structure of the model is displayed through a directed acyclic graph in fig. fig. directed acyclic graph describing the bayesian hierarchical model. square boxes indicate that the values are known the parameters tk and sk define the prior belief of the locations of the latent means θk, whereas the parameters qk and nθk control the spread of mixture component means within a latent cluster and are hence important to control the variation across samples. a large nθk along with a small qk forces the μjk together; it makes large deviations between Σθk and qk unlikely. the parameters hk and nΨk control the expected values and the variation of latent covariance matrices as well as the variation among mixture component covariance matrices in a latent cluster. if nΨk is large each Σjk will be close to Ψk/ for any k, since a high nΨk makes high νk more probable. finally, to simplify sampling from the posterior distribution of the parameters, we add an component assignment variable xij∈{,…,k} describing which component yij is drawn from. to comply with, the a priori uncertainty of component membership is modeled by xij∼ mult, where mult denotes the multinomial distribution. the resulting posterior distribution of all the parameters, denoted jointly by Θ, and x given the data y is given in the additional file 1: section a. in section b we describe the markov chain monte carlo sampling scheme used to generate posteriors for our model parameters. the computational bottleneck of the sampling scheme is the sampling of x, with a computational complexity bounded by o. to handle high dimensions diagonal covariance matrices can be used instead, in which case the complexity is bounded by o. however, for datasets with more than dimensions the mathematical feasibility of using gaussian mixture models without any prior dimension reduction needs to be seriously considered first, due to the curse of dimensionality. absent components in some flow cytometry data sets not all cell populations are present in all samples. in our model this corresponds to that πjk= for some. however, mixture component parameters for empty clusters will still affect the mixing of the mcmc for the parameters of the latent cluster. it can also happen that if a cluster is empty that the mixture component moves and split a neighboring cluster in two. to avoid this in such data sets we extend the model by introducing a variable zj∈{ }k that says which components are active in sample j. this has the further advantage that when sampling from the posterior distribution of the model we get the probability for each cluster that it is present in a sample. we impose a prior on zj which is proportional to exp−cs∑k=1kzji∑k=1kzj> where i denotes the indicator function and cs> the prior makes the model prefer fewer activated clusters so that if there is a very small cluster the likelihood will be larger if it is inactivated, which prevents spurious clusters. the strength of this prior can be adjusted to the expected size of the smallest clusters. the changes to required by this extension are straightforward but inference of the model becomes a bit more involved since removing components reduces the dimension of the model. to accommodate for this we have included a reversible jump step in our sampling algorithm. details are given in the additional file 1: section b. merging latent clusters to determine the “correct” number of clusters in a data set directly from the data is an ill-defined problem, since what should be considered to be a separate cluster depends on the interpretation of the data. nevertheless, there are many different criteria which can be used to guide the decision about the number of populations. we use overlap between components—measured by bhattacharyya distance—and unimodality of the resulting super clusters—measured by hartigan’s dip test —to determine which latent clusters to merge and to indicate our confidence in the mergers. in an evaluation of criteria for merging gaussian components to represent more complex distributions, the bhattacharyya distance performed well. bhattacharyya distance merges clusters according to a pattern-based cluster concept as opposed to a modality-based concept. with a pattern-based cluster concept a small dense cluster inside a sparse cluster—for example a well specified cell population inside a region with sparse outliers—will be considered to be different clusters. this would not be the case for the modality-based cluster concept as long as the generating probability density is unimodal. the bhattacharyya distance between n and n is dbhat=1/8·⊤Σ¯−1+1/2·log|Σ¯|/|Σ2||Σ2|, where Σ¯=/. in order to measure bhattacharyya distance between mixtures of gaussian distributions, which is necessary for deciding if super clusters should be merged with other clusters, we approximate each mixture with a gaussian distribution. the means and the covariance matrices are estimated using a soft clustering of the data points inferred from the sampling of xij, detailed in the additional file 1: section c. however, it is not obvious how to set a threshold for dbhat, since the appropriate threshold depends on the distribution of the data, which is unknown. because of this we use a low soft threshold d and a high hard threshold d two clusters closer to each other than d are always merged, two clusters whose distance is between d and d are only merged if they fulfill an additional criterion based on hartigan’s dip test for unimodality. unimodality is an appealing heuristic for defining cell populations, and it has frequently been used for automated gating. it has two main limitations. the first one, that populations intuitively should be separate if they have very different densities—even when they overlap so that their combined distribution is unimodal—can be bypassed by combining unimodality with a pattern-based merging criterion such as bhattacharyya distance. the second one, that it is difficult to determine if a multi-dimensional empirical distribution is multimodal, is usually handled by considering one-dimensional projections. this is the approach we take here, using hartigan’s dip test of unimodality for each of the projections onto the coordinate axes where bhattacharyya overlap is low, and for the projection onto fisher’s discriminant coordinate. if for a proposed merger, any of these projections is found to be multimodal, the clusters are not merged. further details of the merging procedure are given in the additional file 1: section c. quality control to verify that the output of bayesflow fulfills the user’s requirements, a number of checks are performed: convergence of the mcmc sampler is established by viewing trace plots of sampled parameters. to ensure that variation of the two different populations are not confused with each other, we require that the bhattacharyya distance as well as the euclidean distance from each sample component to its corresponding latent component should be smaller than these distances to any other latent component which does not belong to the same super cluster. to ensure that the obtained clusters should not be divided further, hartigan’s dip test is computed for the projections onto the coordinate axes of all super clusters. projections which have a dip test p-value below —the threshold for merging components —are visualized using histograms of quantiles of the weighted data belonging to the cluster. to ensure that the model fits the data reasonably well, samples from the posterior predictive is compared to the true data in one- and two-dimensional histograms. to ensure that there are no outliers among the cluster centers, the centers for each cluster are plotted together along one dimension. additionally, to detect components with aberrant shapes, the eigenvectors corresponding to the largest eigenvalues, multiplied with the corresponding eigenvalues, can be viewed. if any of the quality criteria is not met, the simulation should be rerun, either using the same or different parameters. even if the same parameters are used a different result can be obtained due to randomness in the initialization. experiments simulated data in order to verify that the proposed sampling scheme can find the correct model parameters, the mcmc algorithm was applied to two simulated datasets. the first dataset was three-dimensional, which enables direct visual evaluation. it had four latent clusters across eighty artificial flow cytometry samples; each sample had cells giving a total of million cells. one of the latent clusters was present only in eight samples and another one was present in samples, so that the ability to find rare cell populations was tested. moreover, the cluster which was present in only eight samples contained only % of the total number of cells, thus also the ability to find small cell populations was tested. the parameters and the algorithm used for generating the data are given in the additional file 1: section d. the second data set was designed to test the ability to handle large data. it was eight-dimensional, with eleven latent clusters and artificial flow cytometry samples. each sample had measurements of cells, giving a total of million cells. four of the eleven clusters were missing in half of the samples. additional file contains a python script and data for regenerating this data set. prior parameters and initial values for the mcmc sampler are given in the additional file 1: section d. all priors were chosen to be non-informative. the outlier component was not used for inference in the small dataset, but it was used for the large dataset. the mcmc sampler ran first for a number of burn-in iterations, then the posterior distribution was explored in a number of production iterations. during the production iterations, apart from sampling parameters of the model, a value of y was also drawn, i.e. a sample from the posterior predictive. for the first synthetic data set burn-in and production iterations were used. for the second, larger, data set we used burn-in iterations and production iterations. for the second data set the mcmc sampler was run on amazon cloud, using cores. each iteration took on average one second, so that about h was needed in total. flow cytometry data we analyze two flow cytometry data sets with bayesflow: the data set gvhd from the flowcap i challenge—with four markers, samples and approximately cells per sample—and a data set obtained from the r package healthyflowdata with technical replicates of pbmc samples from healthy donors—in total samples with approximately cells, also measured with four markers. in the gvhd dataset we can compare the gating obtained from bayesflow with manual gating provided from flowcap as well as automated gating from a wide range of other methods. in healthyflowdata we can instead compare gating between technical replicates to see if samples are treated in a consistent manner. for the healtyflowdata dataset we used an exploratory approach with non-informative priors. we ran multiple simulations and gradually increased the number of components until we passed the quality criteria described under quality control; we finally arrived at using k= components. for the gvhd data set we started with an exploratory approach and gradually increased the number of components, but in the quality checks we noted one population in one of the samples which was very hard to capture. then we decided to use an informative approach for this population. using a scatter plot, fig. we set boundaries for this population in the dimensions given by the cd and the cd8b marker and computed its mean and empirical s covariance matrix. we used the mean to set an informative prior for θk and the mean and the empirical covariance to initialize the component. prior parameters in both the informative and non-informative case are described in the additional file 1: section e. bayesflow applies three data preprocessing steps: 1) data points with extreme values in at least one dimension are removed. such data points can lead to components with singular covariance matrices, and a well designed flow cytometry experiment should not have significant populations with such values. 2) the data is scaled using the and % percentiles q and q of the pooled data, with the same scaling for all samples, so that q = and q = for each marker for the pooled data. this is done in order to be able to set informative priors in an intuitive way. 3) before testing which components should be merged, a very small amount of noise is added to the data. this is since the discreteness of the original flow cytometry measurements can lead to a striped pattern in the flow cytometry data and also when it is not visible to the human eye it disturbs the dip test. after preprocessing, parameters for the mcmc sampler were initialized by running the em algorithm on the pooled data, followed by the initialization scheme used for the large synthetic dataset, detailed in the additional file 1: section d. we ran burn-in iterations and production iterations of the mcmc sampler for both experiments. the burn-in period consisted of five phases: in the first phase, the priors on variation in location and shape were modified to force clusters together. before the second phase, priors parameters were set to normal again. after the second phase, components which were considered to be outliers were turned off. they were forced to stay off during a short third phase, but from the forth phase and onwards components were allowed to turn on and off. label switching was allowed during the initial four phases in order to escape non-desired local minima, but then disallowed. the values of parameters controlling the simulation during the burn-in and production period are given in additional file 1: table s we also applied the two other joint gating methods based on bayesian hierarchical models: aspire and hdpgmm. for aspire parameters were chosen according to the strategy recommended by dundar et al.; details are given in the additional file 1: section e. for each run we used in total iterations, of which were set as burn in iterations. for hdpgmm default parameters were used, with a burn-in period of iterations and a production period of iterations. we ran bayesflow and aspire on a ghz quad core cpu. a bayesflow run took h for the gvhd dataset and h for healthyflowdata. aspire took in total h for the gvhd dataset and h for healthyflowdata per run. four runs of aspire was needed to determine the κi parameters. hdpgmm was run on a dual core gpu. it needed h for the gvhd dataset and approximately h for the healthyflowdata dataset.  simulated data we begin by analyzing the smaller data set. in fig. we show univariate and bivariate histograms of all synthetic cell measurements pooled together, as well as the corresponding histograms of the data from a single flow cytometry sample where all four clusters are present. note that the data when pooled together has a complicated density, as it is in fact a mixture of multivariate normal densities. fig. a one and two dimensional histograms for one synthetic flow cytometry sample containing data points; b histograms of data points drawn uniformly from the pooled data from the synthetic data experiment in fig. we show the same univariate and bivariate histograms, but this time with samples from the posterior predictive distribution of y. from the synthetic cell measurements generated from the inferred models of the datasets it is clear that the inferred models are accurate and capture the variation across samples, which a model only of pooled data cannot do. fig. a one and two dimensional histograms of posterior draws of y for the flow cytometry sample displayed in fig. 2 a; b histograms of posterior draws of y drawn uniformly from all the flow cytometry samples, thus matching fig. 2 b fig. bayesflow component parameter representations of inferred latent clusters and mixture components together with histograms of real data and synthetic data generated from the model for healthyflowdata. the center of each ellipse is the mean and each semi-axis is an eigenvector with length given by the corresponding eigenvalue of the projected covariance matrix. for the latent clusters the parameters Ψk) are shown, for the mixture components the parameters are shown. each component or cluster is depicted with the same color as in fig. 5; different shades of same color corresponds to latent clusters that have been merged fig. summary statistics of inferred cell populations in bayesflow, aspire and hdpgmm, ordered by population size. for hdpgmm, the six largest components after merging are shown, the remaining components have together at most of the cells in a sample. the noise component in bayesflow has at most of the cells in a sample. a locations μ jk of mixture components that represent each population, in each sample, cf. fig. b box plots of the soft clusters in the pooled data, cf. fig. c population proportions across flow cytometry samples fig. cell population which is hard to detect in the gvhd dataset fig. the posterior mean of the mixture component centers, μ jk, and the true cluster centers in the small synthetic data experiment finally in figs. and the marginal posterior distributions of the latent cluster parameters θk and Ψk, subtracted by their true values, are presented. in fig. the dot represents the difference between the median of posterior distribution and the true value of each θk. the vertical lines represent the and % quantiles. fig. displays results for each latent covariance matrix Ψk/ in the same way. from figs. and we see that the true parameters of both the means and the covariances are all between the and % quantiles of the posterior distribution. fig. the difference between the true value of each entry in each θ k and the approximated marginal posterior distribution generated by the mcmc sampler in the small synthetic data experiment. the black dot represents the median and the vertical line goes between the and % quantiles. the light gray horizontal line is the line fig. the difference between the true value of each of the entries in Ψ k/ and the approximated marginal posterior distribution generated by the mcmc sampler in the synthetic data experiment. the black dot shows the median, and the black vertical line goes between the and % quantiles. the light gray horizontal line is the line the true and estimated cluster centers of the eight-dimensional data set cannot be displayed efficiently with just three dimensions at hand, but a three-dimensional projection is shown in fig. the average error in euclidean distance in the full eight-dimensional space is, which can be compared to the average error had the latent mean across samples been used, namely, which is the best that could have been obtained from a model not including variation between samples. the outlier component was used for inference in the results presented here, but omitting it has very small effect. fig. the posterior mean of the mixture component centers, μ jk, and the true cluster centers in the large synthetic data experiment for the first three dimensions in fig. we show the posterior distribution of the latent cluster means where again the dot represents the difference between the median of posterior distribution and the true value of each θk. the vertical lines are the and % quantiles. the posterior samples have been divided by the standard deviation of the true θk so that the scales across the clusters are equal. some of the credibility intervals do not contain zero, but this is explained when studying the intervals that would have been obtained if the true μk were used, since they are almost identical. fig. the difference between the true value of each entry in each θ k and the approximated marginal posterior distribution generated by the mcmc sampler in the large synthetic data experiment. the black dot represents the median and the vertical line goes between the and % quantiles. to get the axis on the same scale for all the clusters, they are scaled by the standard deviation of μ k. the light gray horizontal line is the line. the red dot and lines is the same however where one uses the true μ k to estimate θ k, rather then the μ k obtained by taking the posterior means of the mixtures we thus see that cluster centers and credibility intervals for latent clusters are captured well in both synthetic data sets. flow cytometry data gvhd for the analysis of the gvhd dataset we did twelve runs of bayesflow in the informed setup described above. seven were excluded due to confusion between populations, i.e. at least one sample component was closest to the wrong latent component; of the remaining five, one more run was excluded since it has not converged, and another two because of multimodal clusters. this leaved two runs that passed the quality control. additional file 1: figs. s and s show trace plots and projections of clusters with high dip test values respectively. table accordance with manual gating for gvhd data set. for hdpgmm we also report the result when components are merged according to our merging procedure. when this procedure is applied to the results obtained by aspire, no components are merged, i.e. the original result is identical to what is obtained after mergeing flowmeans samspectral ensemble flowcap one of the two bayesflow runs has the highest accordance with manual gating, the other one is on par with flowmeans and samspectral, which is considerably higher than aspire and hdpgmm. however, as can be seen in fig. the gating of different samples is arguably most consistent for bayesflow as compared to manual gating, flowmeans and samspectral. fig. gated events according to four methods of the twelve samples in the gvhd dataset, projected onto the two first dimensions. for bayesflow, the run with least accordance with manual gating, run is shown. similar plots for aspire and hdpgmm as well as bayesflow run are shown in the additional file 1: figure s6 to get a further understanding of the variability between samples in bayesflow, summary statistics for the obtained components and cell populations are shown in fig. fig. summary statistics of the six cell populations obtained by bayesflow in the dataset gvhd. the outlier component has at most of the cells in a sample. a each panel displays the locations μ jk of all mixture components that represent the population, across all samples. different shades of a color represent different latent components k. b box plots of the soft clusters in the pooled data. the boxes go between the quantiles q km, and q km,, the whiskers extend to q km, and q km,. the α-quantile for component k in dimension m, q km,α, is here defined as qkm,α=mini′j′{yi′j′m:α<∑ij:yijm<yi′j′mwijk}. c population proportions in each of the twelve flow cytometry samples healthyflowdata we did runs of bayesflow with k= ten of these were excluded due to confusion between populations, moreover two runs were excluded since they had clusters with clearly multimodal distributions. for the six runs that passed the quality control, 3– components were turned off across all samples; they are excluded from visualizations. additional file 1: figs. s s and s show trace plots, projections of clusters with high dip test values and eigenvectors of covariance matrices respectively. in fig. we visualize model fit and inter-sample variation for the first of the six runs that passed the quality control by plotting latent and sample components as well as histograms of real data and synthetic data generated from the model, for two different samples and for the pooled data. we can thus see how shape variations are captured by the model. the output of bayesflow, aspire and hdpgmm can be compared in fig. the merging procedure we used for bayesflow has been applied for both aspire and hdpgmm, however for aspire no components were merged by this. in bayesflow each of the populations correspond to clear expression patterns, which is not the case for the other methods. for example the first population is clearly cd4+cd8- t-cells whereas for both aspire and hdpgmm this population contains both components which are cd8- and components which are cd8+. we also compare intra-donor variation of cell population size to inter-donor variation for the six bayesflow runs, as well as for aspire and hdpgmm in fig. for aspire there are inter-donor distances which are clearly smaller than some intra-donor distances, which is not the case for bayesflow and hdpgmm. fig. distances within and between donors as measured by ℓ distance between vectors of population sizes. for the six bayesflow runs and hdpgmm there is very little or no overlap between within-donor and between-donor distances, whereas for aspire there is clear overlap discussion from different runs of bayesflow we can get different representations of data, as in the case of the gvhd dataset. this is because with highly overlapping populations there might be multiple models representing the data equally well. but since all samples are gated jointly in every run, the gated populations can still be compared across samples. the user might have a preference for one representation or the other though, and informative priors can be used to guide bayesflow to a preferred representation. bayesflow is not aimed at discovery of rare cell populations, but it can be used together with an algorithm specifically designed for detecting rare cell populations in a sample, such as swift, and then use informative priors to find how this population occurs across an entire set of samples, in a similar way as was done in the gvhd dataset. how much clusters should be merged is a decision that needs to be taken by the interpreter of the data. in some settings one might want to be restrictive with merging and then use higher thresholds. in others one might want additional mergers after viewing joint one-dimensional projections of the clusters. the bayesflow pipeline does not in itself include any compensation or any of the non-linear transformations which are often used for flow cytometry data, such as logicle. compensation is a linear transformation and gaussian mixture models are invariant under linear transformations, so they perform equally well on uncompensated and compensated data. non-linear transformations such as logicle can make gaussian populations non-gaussian, which makes inference harder. the flow cytometry data we used for the experiments had already been compensated, the healthyflowdata data set had also been transformed with an asinh transform; details are given in the additional file 1: section e. bayesflow finds a joint representation of an entire set of samples. in order for this representation to be reasonable there has to be sufficient correspondences between samples. even if for a data set with very little correspondences a joint model could be obtained by using a very large number of components, it would hard to gain any insights from such a model. in such a case an entirely computational pipeline without the cell population identification step would be preferred. bayesflow can be computationally intensive if many runs are needed to pass the quality control. for these cases it would be desirable to complement bayesflow e.g. with initialization methods that would allow passing the quality control more often, so that few runs in bayesflow would be needed. fast initialization methods and early quality checks aiming at this would therefore be of interest for the community and is something that we propose for further study.  in this paper we have presented a new bayesian hierarchical model designed for joint cell population identification in many flow cytometry samples. the model captures the variability in shapes and locations of the populations between the samples and we have demonstrated its use in an exploratory as well as in a partly informed setting with some prior information. we showed that for synthetic datasets generated from the model, the parameters were recovered with high accuracy through a mcmc sampling scheme. the model was then applied to a real flow cytometry data set where a manual gating was available, and it was shown to have very high accordance with manual gating as compared to other automated gating methods, while at the same time the gating was more consistent across samples than either the manual gating or other automated gating methods. when applied to another flow cytometry data set with technical replicates of blood from healthy donors, bayesflow gave a parsimonious representation of the data, which enables visualization and monitoring of its parameters. the obtained cell populations had clear expression patterns as opposed to the clusters obtained by aspire and hdpgmm, where for example cd4+cd8- t-cells where in the same cluster as cd4+cd8+ t-cells. the population sizes obtained by bayesflow and hdpgmm respectively had lower intra-donor variation compared to inter-donor variation than what was obtained from aspire. many approaches of automated gating of multiple flow cytometry samples in parallel have been aimed at finding features of the data so that either samples can be classified into groups, e.g. cancer or normal, or they can be used to predict an outcome such as expected time to progression of disease. features are often designed based on characteristics of cell populations, but usually not so much attention has been given to ensure that they represent actual cell populations. bayesflow takes the opposite approach and gives a representation of the data according to cell populations, with the same cell populations across the entire set of samples. the advantages to this approach are among others that the result is directly biologically interpretable and that a rich output is given which can be explored in many different ways which are familiar to someone who is used to manual gating. in this way we can join the objectivity and ability to work in high dimensions and with many samples of automated gating with the flexibility in interpretation of manual gating. additional files additional file supplementary material. the supplementary material contains the posterior in bayesflow, the mcmc sampling scheme, additional details on the merging of components, information about the data generation, priors and initialization for the synthetic data example; parameters used for aspire, additional details on healthyflowdata, the priors and the initialization procedure used when studying this data set and further results pertaining to the real flow cytometry data set, including fitting gaussian mixture models to individual samples of healthyflowdata with the em algorithm and scatter plots of gvhd for aspire, hdpgmm and bayesflow run additional file data generation files. a python script for generating the large synthetic dataset, along with means, covariances and weights needed for this. competing interests the authors declare that they have no competing interests. authors’ contributions kj, jw and mf conceived and planned the study. jw and kj designed the statistical model and the inference procedure. kj, jw and mf designed the experiments. jw and kj implemented bayesflow and ran the experiments. kj and jw wrote the article with the help of mf. all authors read and approved the final version of the manuscript. 
the second real data set contains replicated flow cytometry measurements of samples from healthy individuals. bayesflow is applied to two synthetic and two real flow cytometry data sets. flow cytometry is a widespread single-cell measurement technology with a multitude of clinical and research applications. interpretation of flow cytometry data is hard; the instrumentation is delicate and can not render absolute measurements, hence samples can only be interpreted in relation to each other while at the same time comparisons are confounded by inter-sample variation. despite this, most automated flow cytometry data analysis methods either treat samples individually or ignore the variation by for example pooling the data. we introduce bayesflow, a pipeline for latent modeling of flow cytometry cell populations built upon a bayesian hierarchical model. bayesian hierarchical modelsflow cytometrymodel-based clusteringams subject classification for the first real data set, taken from the flowcap i challenge, bayesflow does not only give a gating which would place it among the top performers in flowcap i for this dataset, it also gives a more consistent treatment of different samples than either manual gating or other automated gating methods. expert knowledge can be incorporated through informative priors and the results can be supervised through compact and comprehensive visualizations. bayesflow gives here cell populations with clear expression patterns and small technical intra-donor variation as compared to biological inter-donor variation. as opposed to other joint gating methods, effort is put at ensuring that the obtained partition of the data corresponds to actual cell populations, and the result is therefore directly biologically interpretable. a key requirement for models that include multiple samples is the ability to visualize and assess inferred variation, since what could be technical variation in one setting would be different phenotypes in another. modeling latent relations between samples through bayesflow enables a systematic analysis of inter-sample variation.
 flow cytometry is a widespread single-cell measurement technology with a multitude of clinical and research applications. interpretation of flow cytometry data is hard; the instrumentation is delicate and can not render absolute measurements, hence samples can only be interpreted in relation to each other while at the same time comparisons are confounded by inter-sample variation. despite this, most automated flow cytometry data analysis methods either treat samples individually or ignore the variation by for example pooling the data. a key requirement for models that include multiple samples is the ability to visualize and assess inferred variation, since what could be technical variation in one setting would be different phenotypes in another.  we introduce bayesflow, a pipeline for latent modeling of flow cytometry cell populations built upon a bayesian hierarchical model. the model systematizes variation in location as well as shape. expert knowledge can be incorporated through informative priors and the results can be supervised through compact and comprehensive visualizations. bayesflow is applied to two synthetic and two real flow cytometry data sets. for the first real data set, taken from the flowcap i challenge, bayesflow does not only give a gating which would place it among the top performers in flowcap i for this dataset, it also gives a more consistent treatment of different samples than either manual gating or other automated gating methods. the second real data set contains replicated flow cytometry measurements of samples from healthy individuals. bayesflow gives here cell populations with clear expression patterns and small technical intra-donor variation as compared to biological inter-donor variation.  modeling latent relations between samples through bayesflow enables a systematic analysis of inter-sample variation. as opposed to other joint gating methods, effort is put at ensuring that the obtained partition of the data corresponds to actual cell populations, and the result is therefore directly biologically interpretable. bayesflow is freely available at github. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords bayesian hierarchical modelsflow cytometrymodel-based clusteringams subject classification primary 62p10secondary 62f1568u99http://dx.doi.org/ /501100004063knut och alice wallenbergs stiftelsekungliga vetenskapsakademien, gs magnussons fondissue-copyright-statement© the author 2016 