{'fraction_non_alphanumeric': 0.03068577658414657, 'fraction_numerical': 0.00922700561065759, 'mean_word_length': 4.613972234661889, 'pattern_counts': {'":': 1, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Coarsened mixtures of hierarchical skew normal kernels for flow cytometry analyses

Flow cytometry (FCM) is the standard multi-parameter assay used to measure single cell phenotype and functionality. It is commonly used to quantify the relative frequencies of cell subsets in blood and disaggregated tissues. A typical analysis of FCM data involves cell classification - the identification of cell subgroups in the sample - and comparisons of the cell subgroups across samples. While modern experiments often necessitate the collection and processing of samples in multiple batches, analysis of FCM data across batches is challenging because the locations in the marker space of cell subsets may vary across samples. Differences across samples may occur because of true biological variation or technical reasons such as antibody lot effects or instrument optics. An important step in comparative analyses of multi-sample FCM data is cross-sample calibration, whose goal is to align cell subsets across multiple samples in the presence of variations in locations, so that variation due to technical reasons is minimized and true biological variation can be meaningfully compared. We introduce a Bayesian nonparametric hierarchical modeling approach for accomplishing calibration and cell classification simultaneously in a unified probabilistic manner. Three important features of our method make it particularly effective for analyzing multi-sample FCM data: a nonparametric mixture avoids prespecifying the number of cell clusters; the hierarchical skew normal kernels allow flexibility in the shapes of the cell subsets and cross-sample variation in their locations; and finally the"coarsening"strategy makes inference robust to small departures from the model, a feature that becomes crucial with massive numbers of observations such as those encountered in FCM data. We demonstrate the merits of our approach in simulated examples and carry out a case study in the analysis of two FCM data sets.


Introduction
Flow cytometry (FCM) is a standard biological assay used to measure single cell features, referred to as "markers", and is commonly used to quantify the relative frequencies of cell subsets in blood and disaggregated tissues. These features may be general physical and chemical properties of a cell or the amounts of specific biological macromolecules on the cell surface or within a cell. For example, general markers include forward scatter (FSC) which reflects cell size, side scatter (SSC) which reflects the granularity of a cell, and viability probes which are enhanced in dying or dead cells. Most markers correspond to labeled probes to cell surface or intracellular proteins. For example, a fluorescent-labeled antibody can be used to indicate the density of CD3, a protein associated with the T cell receptor whose presence is indicative of a T cell. An alternative technology based on time-of-flight spectroscopy (mass cytometry) uses heavy metal isotopes instead of fluorescent dyes to label probes, but the data generated is very similar and we treat data from flow and mass cytometry interchangeably in this manuscript.
A typical modern cytometer can provide researchers with measurements on dozens of such markers at the single-cell level. Most commonly, these high-dimensional data are visualized and thereby analyzed via scatter plots of pairs of margins. Such data often involve multiple samples collected for the purpose of comparative analysis: the comparison of cell subset features across different samples. However, variations in cell subset locations (in the marker space) frequently occur due to uncontrolled technical reasons unrelated to the underlying biological differences. In particular, technical differences often make it extremely challenging to compare samples processed in different laboratories. A prerequisite for proper analysis of such data is cross-sample calibration-aligning cell subsets across multiple samples-so that cell subset properties can be meaningfully compared.
A comparative analysis of FCM data generally involves cell classification-i.e., the iden-tification of cell subsets (or clusters) in the samples. The standard approach is through manual gating, during which a human expert uses sequences of two-dimensional windows (the so-called "gates") to filter for cells of interest on two-dimensional scatter plots of the data.
Because manual gating on paired scatters becomes increasingly laborious as the dimensionality of the data grows, automated methods for cell subset identification from FCM data are becoming increasingly necessary. However, none of the existing automated classification methods based on clustering performs calibration in-house and therefore cannot be directly applied to multi-sample studies.
It is worth noting that consideration of cross-sample variability, which may be construed as a step toward calibration, appears in a number of previous works. For example, Cron et al. (2013) developed Bayesian hierarchical modeling extensions to the Dirichlet Process Mixture (DPM) of Gaussian kernels, in which each sample has its own set of weights, thus accounting for sample variability in subset sizes. Dundar et al. (2014) further expanded the approach to also model sample variability in kernel parameters by adding a Dirichlet process mixture (DPM) prior on the means of the Gaussian kernels. Lee et al. (2015) proposed a hierarchical structure where each sample is a finite mixture model of multivariate skew-t distributions with random effects on the location parameters. Soriano and Ma (2019) suggest using a DPM of hierarchical Gaussian kernels to allow for variability in cell subset locations across samples.
Building on these works, we continue to adopt the nonparametric mixture framework and allow all of the aforementioned features of cross-sample variation in our model. In particular, we explicitly model sample variability in relative frequencies by allowing each sample its own set of weights and incorporate a hierarchical multivariate skew normal kernels to characterize both the flexible shapes of the cell subsets as well as their variation across samples. This hierarchical modeling approach gives rise to a natural scheme for performing calibration, as explained in Section 2.3. In addition, we incorporate into our inference recipe a posterior "coarsening" strategy (Miller and Dunson, 2018), which is a general and principled approach for allowing small deviations from a statistical model through likelihood tempering in a Bayesian model. With this strategy, our method becomes robust to shape-misspecification in the kernels and mitigates the undesirable issue of DPMs in producing a diverging number of clusters as sample size grows, thereby increasing the efficacy of our model in classification and calibration.
Before describing our method in detail, we provide a quick review of some other previous statistical approaches to analyzing FCM data as well as some related literature in Bayesian modeling which will place our method in context. Murphy (1985) applied K-means cluster analysis to FCM data. Generally, K-means classification tends to generate clusters of similar spatial dimension, and lack statistical interpretation. Bakker Schut et al. (1993) perform cluster analysis using K-means initialized with a large number of seed points, followed by a modified nearest neighbor technique to reduce the large number of subclusters.
This method caters to symmetric clusters. Boedigheimer and Ferbas (2008) applied finite Gaussian mixtures to FCM data, and adopted a frequentist estimation strategy based on the expectation-maximization (EM) algorithm. Chan et al. (2008) fit the data with a finite mixture of multivariate Gaussians using standard conjugate Bayesian analysis and Gibbs sampling for inference. Methods that model FCM data with Gaussian kernels suffer from the obvious weakness that clusters are clearly asymmetric. Malsiner-Walli et al. (2017) offered a Bayesian model that allows a finite mixture of mixtures of Gaussian kernels and demonstrate this approach on FCM data. This method automates the selection of the number of clusters, and allows asymmetric clusters. In the frequentist literature, Pyne et al. (2010) formulated the FLAME framework that models FCM data with a finite mixture model of skew-t distributions. O'Hagan et al. (2016) suggested using the multivariate normal inverse Gaussian distribution kernels in the context of finite mixture models. Lo et al. (2008) proposed a finite mixture of t distributions with a Box-Cox transformation in order to reduce asymmetry. This method is data-and variable-dependent, which makes full automation of the process difficult. Arellano-Valle et al. (2009) discussed Bayesian mixtures of multivariate skew normal kernels, but did not provide a unified approach that handles all three parameters of the skew normal distribution. Frühwirth-Schnatter and Pyne (2010) developed a Bayesian, fully conjugate multivariate finite mixture model with multivariate skew normal and skew-t distributions. As is further discussed in Section 2.4 this formulation is constrained by the strong conjugate prior structure and does not allow for intuitive treatment of calibration for FCM data. Hejblum et al. (2017) used the Bayesian formulation of Frühwirth-Schnatter and  that also accommodates dependencies within the data using a sequential sampler.
Our approach addresses several important difficulties faced by these existing approaches.
Methods based on finite mixture models require determining the number of clusters in advance (and more so in a hierarchical setting that requires determining the number of clusters for each sample) thus treating it as a tuning parameter. These approaches run into obstacles when new data is introduced after the selection of the number of clusters: if a new sample reflects biological variation in the sense that a new cluster is formed in addition to the clusters in the control sample, this cluster cannot be identified as one when the number of clusters is fixed to that of the control sample. Thus, sophisticated analyses involving comparison of data with biological variation across batches are severely limited. In contrast to finite mixtures, an infinite mixture such as a DPM avoids the presetting of the number of clusters.
However, a common issue that arises with DPMs is that the number of estimated clusters always grows to infinity as the number of observations grow, even as the underlying structure of the data remains the same. A common way to handle this difficulty is to use the symmetric approximation to DPMs (Ishwaran and James, 2001) and cap the number of clusters at a relatively low number. With data sets as large as a typical FCM sample, this essentially reduces DPMs to finite mixture models and the maximal number of clusters will often just be the estimated number of clusters. In contrast, we take a model robustification approach to address the issue by incorporating a strategy called model "coarsening" (Miller and Dunson, 2018). In addition, existing mixture model based approaches using Gaussian and t kernels almost invariably require post hoc merging of clusters due to the massive sample size of FCM data. The need for such is mitigated in our approach by the flexibility of the skew kernels as well as the coarsening. Finally, when estimating mixtures with multivariate skew kernels, frequentist estimation with EM or Bayesian approaches with conjugate priors and Gibbs sampling face difficulties in overcoming the multimodality of the likelihood function of those kernels. We design a hybrid sampler that embeds a Population Monte Carlo (PMC) step into a Gibbs sampler, which reduces the risk of the sampler to be caught in a local mode and allows us to choose and sample from more reasonable though non-conjugate priors effectively.
The rest of the paper is organized as follows. Section 2 provides the details of our Bayesian hierarchical model, the coarsening strategy, as well as the inference recipe for achieving classification and calibration jointly. In Section 3 we provide numerical examples on simulated data that demonstrate the efficacy of our model. In Section 4 we carry out a case study on a 6-dimensional FCM data set as well as a 19-dimensional mass cytometry data set.

A Bayesian nonparametric hierarchical model
The complexity of the observed structures in FCM data requires flexible statistical models to characterize their key features. Our model captures those by a DPM with several modeling choices that allow sample variability in weights and kernel locations as well as flexibility in the kernel shapes. We choose to work with DPMs as they do not limit the number of clusters a priori. Since FCM data rarely involves symmetric clusters we choose to work with skew normal (SN) kernels. These are characterized by a location parameter ξ, a scale parameter Σ, and a skew parameter α. (See Supplementary Materials Section S1 for further details about the multivariate SN distribution.) To allow the clusters to vary in sizes across samples, we endow each sample with its own set of mixing weights. In order to allow differences in the location of the clusters, we posit a hierarchical structure for the SN kernels. Figure 1 presents a full graphical view of our model.
We next describe the model components in detail. i the observation number within a sample j; k index is for cluster labels. y i,j is the ith observation in sample j; T i,j the cluster assignment random variable specifying which cluster y i,j belongs to. Each z i,j is associated with y i,j and is used in augmenting the likelihood. ξ j,k is the sample and cluster specific location parameter. ω j,k s are cluster and sample specific DPM weights. α k and Σ k are, respectively, the skewness vector and covariance matrix for each cluster, they are shared among all samples. ξ 0,k and E k are, respectively, the mean and the covariance of the mutlivariate normal distribution from which the sample specific means are drawn.
In the following, we assume that the data are drawn from J samples: y i,j ind " F j , i " 1, . . . , n j and j " 1, . . . , J such that n j is the number of observations in sample j P t1, . . . , Ju and n " ř J j"1 n j is the total number of observations across all samples.
Cluster assignment and weights. Let K be a countable set of cluster labels shared over all samples. For each i " 1, . . . , n j and j " 1, . . . , J, let T i,j be the latent cluster assignment variable for observation i in sample j, such that T i,j " k if and only if y i,j belongs to cluster k for some k P K. As discussed above, we allow each sample to have its own set of weights, ω j,k :" PpT i,j " kq. The subscript j indicates that cluster sizes may vary across samples.
Note that this also allows a cluster to appear in some samples but not in others; the latter samples will simply have very small weights for that cluster. We assign an independent GEMpηq prior (Ewens, 1990) on the sample specific weights. For each j P t1, . . . , Ju, let tω j,k u kPK iid " GEMpηq.
In order to learn η, we also assign it a prior, η " Gammapa η , b η q.
Hierarchical multivariate skew normal kernels. Assume now that each F j has the density f j p¨q " ÿ kPN ω j,k¨g p¨| λ j,k q where g is the density function of the multivariate SN distribution, and each λ j,k is comprised of the three parameters of the SN distribution: λ j,k " tξ j,k , Σ k , α k u. In particular, for each cluster k we assume that the scale parameter Σ k and the skew parameter α k are the same across samples. We allow sample variability in the sample-specific location parameter ξ j,k by assuming that they are distributed around a grand "centroid" cluster mean ξ 0,k following a multivariate Gaussian with covariance E k . Thus we have the following hierarchical kernel for generating an observation in each cluster.
We will further discuss the justification for only letting the location parameter vary between samples in Section 2.2.
We assign multivariate Gaussian and inverse Wishart prior for the means and covariances of the kernels: and follow Liseo and Parisi (2013) and Parisi and Liseo (2018) in the assignment of a prior for the skewness parameter: where δ is a transformed version of α. The priors for δ and Σ induce priors on an alternative equivalent parametrization in terms of ψ and G, which are derived from δ and Σ by multiplying the following Jacobian term, defined separately for each cluster k P K: |J k rpξ 0,k , Σ k , δ k q Ñ pξ 0,k , G k , ψ k qs| " p ź j"1 pG k pj, jq`ψ k pjq 2 q´1 2 .
In practice, we will infer on ψ k and G k and then transform them back to the original parameters. (See Supplementary Materials Section S1 for further details on the alternative parameterizations of the multivariate SN distribution and its computational advantages.) We further assign an inverse-Wishart prior to the covariance of the normal distribution of the cluster means around the grand mean: This completes the specification of our hierarchical model. Figure 1 provides a graphical model representation of the full hierarchical model while Figure 2 illustrates the structure of the data generated from this model.

Coarsening
We use the "coarsening" strategy introduced in Miller and Dunson (2018) to make inference robust to model-misspecifications. This is particularly relevant in the context of FCM datawhere clusters with parametric distributions are not likely to perfectly fit and the massive sample size makes the resulting inference sensitive to such model misspecifications.
Generally, for a model family indexed by some parameter θ, one defines an "idealized distribution" as a member in the model family that we use to represent the data generative mechanism. (Here it is the full hierarchical model presented above.) The coarsening approach assumes the existence of unobserved "ideal data", say Y 1 , . . . , Y n from this distribution. The observed data, say y 1 , . . . , y n , are assumed to be drawn from a true distribution which is in a R-neighborhood (under some discrepancy measure d) of the idealized distribution. When the observed and ideal distributions are the same, Bayesian inference is conducted on the standard posterior distribution, ppθ | Y 1 " y 1 , . . . , Y n " y n q. When they differ, Bayesian inference is performed on the "coarsened posterior": ppθ | dptY 1 , . . . , Y n u, ty 1 , . . . , y n uq ă Rq.
That is, the "posterior" is computed conditional on the event that the empirical distributions of the observed and ideal data are within a R-ball defined by a discrepancy measure d. R can be taken as a random variable, and assigned a prior. When R " exppγq and d is the Kullback-Leibler divergence (or d n a consistent estimator of it), the coarsened posterior is approximated by: πpθ | d n ptY 1 , . . . , Y n u, ty 1 , . . . , y n uq ă Rq 9 " πpθq n ź i"1 p θ py i q ζ where 9 " means "approximately proportional to", ζ " γ γ`n , and πpθq is the prior on θ. Given some ζ P r0, 1s, ś n i"1 p θ py i q ζ is referred to as the "power likelihood" and πpθq ś n i"1 p θ py i q ζ is referred to as the "power posterior". This form is easily implemented within the context of Gibbs sampling for mixture models, and we do so for our model. In our context, we consider the entire hierarchical DPM model as the idealized distribution and treat ζ as a tuning parameter. Generally, the type of deviations that the coarsening procedure tolerates depends on the discrepancy measure. The Kullback-Leibler divergence is sensitive to changes in the shape of a distribution (as opposed to, say, moving probability mass to outlying regions).
This choice allows the hierarchical kernel in our model to focus on the location alone, which simplifies the model and the corresponding sampling algorithm.

Cross-sample calibration
To calibrate multiple samples, we shift each observation in a cluster in a sample by the difference between the grand "centroid" mean and the sample-specific mean of that cluster.
That is, when T i,j " k we compute a corrected value for each observation by adjusting for the shift in the mean: y i,j " y i,j´p pξ j,k`ωk δ k a 2 {πq´pξ 0,k`ωk δ k a 2 {πqq " y i,j´p ξ j,k´ξ0,k q whereỹ i,j is the shifted observation corresponding to the original observation y i,j . In the above, ω and δ are alternative parameterizations to Σ and α. (See Supplementary Materials Section S1 for further details on the multivariate SN distribution.) To incorporate the posterior uncertainty in the cluster assignment, we follow the technique suggested in Soriano and Ma (2017) for calibration by integrating out the cluster assignment variables T i,j s: Erỹ i,j | ty i,j us " y i,j´E rξ j,k´ξ0,k | ty i,j us « y i,j´1 N N ÿ n"1 pξ j,T pnq i,j´ξ 0,T pnq i,j q.
A desirable byproduct of this integrating-out strategy is that through it we also bypass potential label switching issue, which will be discussed further in the next subsection.

Computational strategies
The multimodality of the SN likelihood poses difficulties to EM estimation in frequentist settings and to Gibbs samplers in Bayesian settings. For example, Frühwirth-Schnatter and Pyne (2010) offer a conjugate structure and a Gibbs sampler in order to perform Bayesian estimation of the parameters of the SN distribution. Liseo and Parisi (2013) demonstrate how this approach may fail when multimodality arises in the likelihood even in a simple setting. Our experimentation of the conjugate prior approach also suggested that a large amount of data is required in order to overwhelm the prior and allow for correct inference on highly skewed kernels. Furthermore, the conjugate prior structure entails the elicitation of a joint prior for the location parameter ξ and the skewness parameter ψ, which in the current context of application to FCM data will make tasks such as cross-sample calibration difficult.
Liseo and Parisi suggest utilizing the Population Monte Carlo (PMC) approach to tackle the problem of multimodality while allowing prior modeling on the location parameter indecently from those of the other SN parameters-that is, ξ K K pΣ, δq a priori. This strategy particularly suits our need as it allows us to select flexible and intuitive priors that offer a simple hierarchical structure for the location parameter as in our model. Because the resulting priors are not conjugate, we cannot use a vanilla Gibbs sampler. (In contrast, Hejblum et al. (2017), Dundar et al. (2014), and Soriano and Ma (2019) all constrained their models so that the prior structure will be conjugate and allow blocked Gibbs sampling.) We thus construct a hybrid "Gibbs-PMC" sampler that utilizes PMC moves for the SN parameters while using Gibbs moves on the rest of the model.
A "Gibbs-PMC' hybrid Sampler. Our sampler uses PMC moves for the SN parameters ξ j,k , G k , ψ k , ξ 0,k , E k and z i,j , and then given summary statistics (mean) of all particles for these parameters; it uses Gibbs moves for the other parameters parameters π j,k and T , and Metropolis-Hastings moves for just the single hyperparameter a. The full sampling algorithm is as follows. • Step 0: Initialize a population of M particles ξ 1:M j,k , G 1:M k , ψ 1:M k , ξ 1:M 0,k , E 1:M k and z 1:M i,j . • Step t ą 0: -Let n k be the number of observations in cluster k and n j,k the number of observations for sample j in cluster k.
-Update the Dirichlet pseudo-count parameter a using a Metropolis-Hastings step with the proposal: a˚| a " Gammapa 2¨a 0 , a¨a 0 q where a 0 is calibrated during the burn-in iterations.
-Sample from the full conditional of the mixture weights.  here q pmq is the joint proposal for each particle.
-Scale the t pmq u to sum to 1. -Merge clusters for which the Kullback-Leibler divergence is smaller than a preset threshold.
The full conditional and proposal distributions used in the sampler are described in Supplementary Materials Section S2.
In implementing the sampler, we utilize the finite-dimensional symmetric Dirichlet distribution approximation (Ishwaran and James, 2001) to the Dirichlet process mixture. With this approximation, the J infinite sequences of mixture weights tω j,k u kPK are replaced for each j " 1, . . . , J by: tω j,k u kPK " Dirichletpη{K, . . . , η{Kq where we need to choose the maximal number of clusters K. A known issue of this approximation for DPM when applied to large data sets is that the posterior will concentrates on a diverging number of clusters and K is often the de facto number of estimated clusters. In our approach, however, this issue is addressed with the coarsening strategy-we can set K to be very large while the estimated number of clusters will still be much smaller than K even for large and noisy data sets such as those from FCM.
A common problem in inference algorithms for mixture models that utilize a cluster assignment variable is known as label switching. Since the values the cluster assignment variable assigns (the labels) to the different clusters are exchangeable, the prior and posterior distributions for the parameters of the mixture components are symmetric with respect to permutations of the labels. This problem does not arise in our framework when calibration is the sole purpose as our calibration strategy integrates out the different labeling scenarios. For cell classification, on the other hand, label switching needs to be addressed and many possible strategies are available. Our software implementation handles this issue post hoc using the method of Cron and West (2011). A coherent classification is maintained by choosing a reference classification, which we take from the last iteration of the MCMC chain. Then, for each saved classification a cost matrix is computed (based only on the values of the cluster assignment variable) and minimized by selecting a permutation on its columns using the Hungarian algorithm of Munkres (1957). The resulting minimizing permutation is then applied to the cluster labels.

Simulation study
We demonstrate our method using two-dimensional simulated data for ease of visualization.
To examine the robustness of our method to model misspecification and the effects of coarsening, we apply our method to two different data sets. The first is simulated exactly from our hierarchical model, i.e., the "ideal" distribution. The second data set is generated by "distorting" the first, "ideal", data set to induce model misspecification. In the distorted samples each cluster is narrowed asymmetrically (see supplement for source code for generating these data). Figure 3 presents both the "ideal" data (top row) and the "distorted" data (bottom row).

Sample #1
"Ideal" data Sample #2 "Ideal" data Sample #3 "Ideal" data Distorted data Distorted data Distorted data Figure 3: Simulated data. Top row: an "ideal" from our generative model in 2 dimensions with 3 clusters in 3 samples. Bottom row: a distorted version of the "ideal" data. The colors correspond to true cluster labels. The circles designate true sample-specific cluster means. The crossed circles designate true grand means across the samples.
We apply COMIX with a maximal number of 150 clusters to each data set. In order to examine the effects of coarsening, we show model estimates and calibration results with coarsening (ζ " 0.2) and without it (ζ " 1). We have found in many numerical studies that ζ " 0.2 provides robust results for classification and calibration in many contexts where the observed data deviate from the theoretical skew kernels. As such we adopt 0.2 as a default value in our software. In applications we recommend the user to still carry out a context-specific sensitivity analysis for ζ to ensure that an appropriate value is selected. We provide such an example in our later case studies. The classification and calibration results for the data drawn from the generative model are shown in Figure 4 and the results for the distorted data are shown in Figure 5. In Figure 4, where the data is perfectly drawn from the generative model, each true cluster is identified as such both with and without coarsening. However, in Figure 5, where the distorted data are not exactly from SN kernels as our generative model prescribes, the true clusters are broken to several clusters without coarsening (ζ " 1). Coarsening helps identify and align the distorted clusters correctly.  Figure 4: Calibrated samples for the "ideal" data from our generative model with coarsening (top row, ζ " 0.2) and without (bottom row, ζ " 1). Dark grey: the observations in the "ideal" data prior to calibration. The colors indicate estimated cluster label assignment. The crossed circles indicate estimated grand cluster means.  Figure 5: Calibrated samples for the distorted data with coarsening (top row, ζ " 0.2) and without (bottom row, ζ " 1). Dark grey: the observations in the "distorted" data prior to calibration. The colors indicate estimated cluster label assignment. The crossed circles indicate estimated grand cluster means.
4 Case study: the EQAPOL data 4.1 A 6-dimensional data set We apply our method to FCM samples processed in 10 different laboratories participating in the External Quality Assurance Program Oversight Laboratory (EQAPOL) program. In this study, peripheral blood samples from the same subject was sent to different laboratories together with lyophilized antibody panels for processing. Because the samples are from the same individual and subjected to the same sample preparation protocol, the analytical results should be almost identical in principle. We demonstrate our calibration method to the 10 samples with six features (FSC, SSC, CD3, CD4, CD8, multiplexed IFN-γ + IL-2). Figure 6 shows the t-SNE plots for the raw and calibrated data for the 10 labs. As in the previous section, we show the results for ζ " 0.2 and ζ " 1. The aligned data are much more similar across samples, and both large and small clusters are properly aligned. Specifically, when ζ " 0.2, in all samples the number of estimated clusters is 11, in accordance with the underlying assumption for the batch control study that all samples are the same and thus have the same number of clusters. When ζ " 1 the number of estimated clusters in the different samples varies between 13 to 15. This suggests that coarsening here is indeed useful in addressing model misspecification in the context of FCM data. We carry out a sensitivity analysis on the value of ζ and provide the details in Supplementary Materials Section S3.
Finally, Figure 7 presents the marginal densities of the six markers before and after calibration.

Raw Data
Figure 6: t-SNE plots for the 6-dimensional data. Top row: raw data. Middle row: calibrated data for four samples with ζ " 0.2. Bottom row: calibrated data for four samples with ζ " 1.
Color coding corresponds to estimated cluster assignment variable. The three labs highlighted with red boundary demonstrate where the classification and calibration for ζ " 1 is inferior to that of ζ " 0.2.

A 19-dimensional data set
We further apply our method to a publicly available higher dimensional data set collected using a mass cytometer. We pre-processed the data with an Arcsinh transform using per channel parameters in the binary FCS file. Following Figure 2B in Kleinsteuber et al. (2016) we remove from the samples a "spike-in" batch control based on expression levels of the CD45 barcodes. We then extract the 19 markers described in Kleinsteuber et al. (2016) and demonstrate the calibration and classification results for Patient #1 batch control across the 3 experimental settings. Figure 8 shows t-SNE plots for the raw and calibrated data the of 3 samples for ζ " 1 and ζ " 0.2. Here too, the calibrated data is very similar across the samples, attaining our goal. In this case, we get that for ζ " 1 the number of estimated clusters is 11 across all samples, whereas for ζ " 0.2 it is 10. Figure 9 shows some marginal scatter plots of the raw and calibrated data for the two settings. Although the results are largely similar, here too the smaller number of estimated clusters due to coarsening are beneficial, as the lowest cluster in the CD45RO-Perforin plot of the left sample for ζ " 1 appear to be artificial and uninformative.
The results for these two data sets show that COMIX works well for calibration of flow and mass cytometry data from 6 to 19 dimensions, which span the range of dimensions used in the vast majority of such experimental data sets. Calibration is broadly useful not just for multi-center studies, but also for studies across batches of data as cytometer performance characteristics, antibody lots, and sample preparation often vary over time, necessitating time-consuming and error-prone manual adjustment of gates across batches in manual analysis, and decreasing the robustness of automated methods that ignore the need for calibration. Hence, COMIX has broad application for all but the smallest flow and mass cytometry experimental data sets.

Discussion
We have presented a principled probabilistic approach for calibrating and classifying multisample FCM data. Our approach utilizes a flexible Bayesian nonparametric mixture model with multivariate SN kernels to incorporate the key features of FCM data are potentially high-dimensional and have massive sample sizes, and incorporate the "coarsening" strategy to make inference robust to model-misspecification. Moreover, we constructed a Gibbs-PMC ' '  '  Calibrated Data, ζ = 1 Figure 9: Marginal scatter plots for some markers and all samples in high dimensional FCM data. Blue: raw data. Red: calibrated data for ζ " 0.2 and ζ " 1.
hybrid sampler, which embeds PMC move for the SN parameters into a Gibbs sampler, thereby addressing the multimodality of the posterior on the kernel parameters.
As described in our first case study, we applied the method to data from the EQAPOL program. This program is funded by National Institute of Allergy and Infectious Disease (NI-AID) to improve the immune monitoring proficiency of national and international laboratories participating in NIAID-funded HIV clinical research. The current bottleneck for the EQAPOL flow cytometry proficiency testing program is the need for the central laboratory to identify deviations of the site gating strategy from a reference or consensus strategy for the purpose of feedback and remediation of site analysts. Currently, this is achieved by centralized manual re-analysis of the data sets from every participating site, an extremely labor-intensive and time-consuming activity. We are currently developing a semi-supervised processing pipeline that combines the calibration of data sets using COMIX with programmatic extraction of the gating strategy from FlowJo (the standard software used for manual analysis) workspace XML files to automatically report deviations of the gating strategy from the consensus or reference. We expect that this semi-supervised processing will eliminate the need for manual re-analysis but still provide the EQAPOL program with the information needed to provide site remediation effectively and allow program expansion to include more sites by removing the manual re-analysis bottleneck.
While our method is motivated by and developed for the purpose of analyzing multi-sample FCM data, the modeling and inference techniques used are generally applicable to other multi-sample data that can be effectively modeled by mixtures. In particular, the idea of adopting a flexible kernel, allowing hierarchical structure on the kernel, and incorporating coarsening to further robustify the method can all be readily applied to other types of data.

S2 Full Conditionals and MCMC Proposal
The MCMC sampling is based on the following full conditional distributions, where ζ is the tuning parameter of the power likelihood.
3. Latent random effects z i,j : f pz i,j | T i,j " k, . . .q " ' % ϕpz i,j´mj,k piq, v k q z i,j ě 0 ϕpz i,j`mj,k piq, v k q z i,j ą 0 Where v k " p1`ζ¨ψ T k G´1 k ψ k q´1 and m j,k piq is the ith element of the vector m j,k " v k pζ¨py k ,j´1 n b ξ j,k q T G´1 k ψ k q where y k ,j is the vector of all n j observations in sample j and cluster k.

EQAPOL data: sensitivity analysis
We perform a sensitivity analysis to the values of ζ for the 6-dimensional EQAPOL data set.
We fit the COMIX model for each ζ in t0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1u. In Figures S1 and S2 we show marginal density plots for the raw data (top row) and the calibrated data (each row corresponds to a different ζ). Overall, the calibrated samples are much more aligned compared to the raw data. However, for several of the ζ values we highlighted with red circles some anomalies in the calibration. In addition to the above figures, Table S1 presents the number of estimated clusters (number of unique entries in the estimated cluster assignment variable) for each sample and value of ζ. As our data here is a batch control, we may assume that the true number of clusters is the same across all samples. Indeed, the samples for which the estimated number of clusters is the same across all clusters tend to be more aligned and stable compared to the other samples. Based on this sensitivity analysis, values from among t0.1, 0.2, 0.5, 0.9u are our best candidates for the value of ζ. The lower values will prefer larger clusters, while the higher one will cater to smaller ones. The choice of ζ to apply from among those candidates may depend on expert's input as well as the desired emphases for analysis of further data.