i=1 Pnmols (pu,i ∗ (di − μu)2) i=1

i=1 mu = X

pu,i nmols i=1

The difference of the parameters μ and σ between iterations was used as the convergence criteria, with the convergence threshold of 0.005. After the algorithm converged, the MRF prior was discarded and only the densities of the normal distribution were used to estimate the cell assignment probabilities: pc(di) =

mc∗N (di|μc,σc) mc∗N (di|μc,σc)+mb∗N (di|μb,σb) . That was done because the MRF prior consistently push probabilities to be close to either 0.0 or 1.0, which corresponds to binary ( 4 ) ( 5 ) classification. In contrast, using only normal densities results in more gradual probability values, which can be integrated better with the further probabilistic parts of the algorithm.

When a prior segmentation Lprior was available (e.g. DAPI), it was used as a constraint for the optimization. Given the prior segmentation confidence cprior ∈ [0.0, 1.0], the Expectation step 2 was restricted pc(di) to be greater or equal than cprior for all molecules, assigned to cells in the prior segmentation:

0, this approach converges to the case without a prior segmentation, whereas cprior → 1 will ensure that the molecules assigned to cells in the prior segmentation will be necessarily recognized as intracellular.

Segmentation of cell types. The same MRF formalism can be applied to perform de-novo clustering of molecules by cell type of origin. Moreover, when the scRNA-seq data is available, the approach can match the transcriptional identifies of the clusters identified in the spatial data to those observed in the provided scRNA-seq. To perform the segmentation, we considered gene identity of a transcript as observable data, and modeled the whole dataset as a mixture of Categorical distributions representing different cell types. The number of components ncomps is an input algorithm parameter that must be specified at the start. Then, the basic Expectation step on the iteration t works as the following:

p u ∈ adj(i) = 1 − fMRF zi = u |{zj , wi,j , pc,j }j∈adj(i) = exp

Y fu(gi|W, vu, z\i, pc) = fMRF zi = u |{zj , wi,j , pc,j }j∈adj(i) ∗ Cat(gi|vu) ∗ p u ∈ adj(i) (t) pu,i = p(ut)(gi|W, V, z\i, pc) =

fu(gi|W, vu, z\i, pc)

Psn=co1mps fs(gi|W, vs, z\i, pc)

These formulas differ from the previous case ( 4 ) in several aspects. First, given that gene identities are categorical, and not continuous variables, the two Normal distributions were replaced with ncomps Categorical ones to describe this kind of data. The parameter vu = {vu,q}q∈1:Ngenes here defines expression probabilities per gene for component u. Second, they remove dependency of the density on the segment size mu, as its presence forced the algorithm to prefer components of larger size, and without strong signal from data, the MRF prior tended to eliminate all but the largest component. Such problems did not arise in the previous application, because there two Normal distributions described the observed data much better than one, ensuring that both components (background and signal) were maintained. In contrast, here the whole dataset can be described as a single Categorical distribution with a high quality of fit, pushing the algorithm to converge towards one component. Next, the molecule confidence pc,j , obtained from the background estimation step, is used here for estimation of fMRF . This step is very important, because only intracellular molecules are expected to be clustered over cell types, while background molecules can be arranged in arbitrary patterns. Finally, the updated formulas introduce p u ∈ adj(i) , which ensures continuity of the resulting segmentation: without it, a molecule can be assigned to some cluster even if none of its neighbors belong to that cluster, as fMRF zi = u |{zj , wi,j , pc,j = 0 }j∈adj(i) = exp(0) = 1. Indeed, such ability to re-assign molecules globally can be desirable, as the EM algorithm can get trapped in local minima and is generally sensitive to initialization. And the assignments to the components distant from the local neighborhood can help to reach convergence to a global optimum. However, the probability of assigning a molecule to a cluster that is not connected to it cannot be inferred solely from the density fMRF (zi = u) = 1; it also depends on the scale of the other terms. So, in addition to p u ∈ adj(i) the formulas introduce an additional term, pglobal = 0.05, which defines the probability of assigning the molecule to a component regardless of the connectivity. Then, the Expectation step can be adjusted as:

fu(gi|W, vu, z\i, pc) Cat(gi|vu) p(ut)(gi|W, V, z\i, pc) = (1 − pglobal) ∗ Psn=co1mps fs(gi|W, vs, z\i, pc) + pglobal ∗ Psn=co1mps Cat(gi|vs)) ( 8 ) The downside of the global assignment is the reduced continuity of clusters. To correct for that effect, after the algorithm converged, pglobal was set to 0, which reduces to the basic version of the Expectation step, and the EM iterations were carried out further until convergence to the final result.

The Maximization step fit the parameters for Cat(gi|vu) using the assignment pu,i from the Expectation step: vu,q =

Pi: gi=q (pc,i ∗ pu,i) + 1 Ptn=ge1nes Pj: gj=t (pc,i ∗ pu,j ) + 1 ( 9 ) It is important to note here that the vector vu is a non-normalized probability density, due to the way pseudo-counts were incorporated: instead of adding ngenes to the denominator, 1 was added as it allowed to preserve spiking structure of the sparse distribution.

When prior information about transcriptional composition of cell clusters is available (e.g scRNA-seq cell types), the Maximization step is changed to take it into account. Given prior expression fraction μu,q from the cluster u and the gene q and its standard deviation σu,q, the estimate was adjusted based on the Z-score value ζu,q = vu,q−μu,q as: σu,q 0 vu,q =  vu,q,   

μu,q + sign(ζu,q) ∗  μu,q + sign(ζu,q) ∗ |ζu,q| 4

+ 0.75 p|ζu,q| + 1.5 − ∗ σu,q, √ 3 if |ζu,q| < 1 if 1 ≤ |ζu,q| < 3 ,

( 10 ) ∗ σu,q, if |ζu,q| ≥ 3 This function corresponds to no penalty for all deviation from the mean within one σ, linear penalty for deviations less than 3σ, and a super-linear penalty otherwise (Supplementary Fig. 14a). If the standard deviation was not available, it was set to the mean value by default: σu,q = μu,q. It is important to note here that for a given prior clustering reasonable results can often be obtained by running only the Expectation step without the Maximization, which corresponds to setting σu,q = 0, ∀u, q. Particularly, the results shown on the Supplementary Fig. 5 were obtained with only iterating over the Expectation step.

The algorithm was initialized using a vector of gene frequencies, estimated over the whole dataset, multiplied by uniform noise in [0.95; 1.05]: vuin,qit = |{i: gi=q}|∗Unif(0.95,1.05) . Next, these nmols values were normalized by the sum over all genes: vu,q = probabilities were estimated as pu,i =

Cat(gi|vu)

Psn=co1mps Cat(gi|vs)) maximal change in pu,i between iterations, weighted by pc,i: max 1≤i≤nmols The default threshold for the convergence was set to 0.01. 1≤u≤ncomps vuin,qit

Ptn=g1enes vuin,tit . The initial assignment . Convergence was determined based on the

(t) (t−1) |pu,i − pu,i | ∗ pc,i .

Cell segmentation. Assignment of molecules to a cell of origin is another case of the MRF segmentation. In the most basic form, a cell can be modeled as a Multivariate Normal distribution over positions of molecules within a cell and a Categorical distribution over the cell gene composition. Thus, the density of a cell k in the molecule i is: prior fk(obsi|W, pc, z\i) = fMRF (zi = k|W, pc, z\i) ∗ pk ∗ M vN ormalk(xi) ∗ Catk(gi), (11) where pprior is the prior probability of the component k, which is proportional to the number of k molecules, assigned to this component, fMRF (zi = k|W, pc, z\i) is the MRF term, and M vN ormal() is the multivariate Normal distribution (bivariate for the 2D implementation). However, fitting this model by the EM algorithm does not work well for several reasons. First, the number of components cannot be defined beforehand and has to be inferred by the algorithm. To perform such inference, the algorithm we use a Dirichlet prior over the number of components (the approach called Bayesian Mixture Modeling or BMM). Second, not all of the molecules belong to cells, and filtration of the background molecules by a hard threshold over pc,i does not always work well. So the algorithm was adjusted to deal with raw probabilities pc,i, and a separate background component fbg was added to the mixture. Finally, dense homogeneous regions within a tissue do not have sufficient information to be segmented based on their transcritpional composition, so the model (11) will segment the whole region into a single component. Such situations, however, can be resolved much better by utilizing information about the expected physical size of a cell. This was done by introducing a global scale parameter sglobal that was used as a prior for estimating the covariance matrices of the bivariate Normal distribution. If a prior segmentation is provided, sglobal is inferred based on this data, and otherwise it has to be specified by the user.

Initialization. The algorithm is initialized with a large number of components uniformly distributed over 2D space. Setting the initial number of cells to be much larger than the expected number greatly improved the convergence of the algorithm. The initial cell centers were selected uniformly across 2D space, and each molecule was assigned to the nearest center. The center selection was done using a strategy similar the one used for subsampling of the neighbourhood composition vectors (NCVs): the molecules are ranked by the sum over the x and y coordinates, and then the algorithm selects a subset from this array uniformly across indices. When molecule background confidences were available, only molecules with true signal confidence greater than 0.25 are used. The MRF was initialized using Delaunay triangulation in the same way as described above. However, given that the triangulation uses only the information about spatial positions, the MRF was then further adjusted to reflect the neighbourhood composition similarities, based on NCVs. Specifically, NCVs were estimated for each molecule, and the resulting matrix of NCVs was transformed using PCA. Pairwise Pearson linear correlation of the PC vectors, ρi,j , was estimated for any two molecules i and j that were connected by an edge in the Delaunay graph, and the MRF edge weight was then multiplied by max(ρi,j , 0.01).

Fitting Bayesian Mixture Model. Fitting of the Bayesian Mixture Model was performed by incorporating Stick-breaking process into the Stochastic EM algorithm. The algorithm iterates over the four steps over a pre-defined number of iterations N iter (500 by default, which was enough for convergence on all of the tests): 1. Maximize parameters of the distributions given existing assignment of molecules by cells (Maximization step) 2. Sample empty components from the Dirichlet prior (Distribution Sampling step) 3. Stochastically assign molecules to components given the exiting components (Stochastic

Expectation step) 4. Remove all components that have less than two molecules assigned to them

After finishing the iterations, the algorithm re-estimates molecule assignment by averaging it over the last N est iterations (N iter/10 by default).

Maximization. The Maximization step fits the parameters of the Normal and Categorical distributions based on the current assignment of molecules to components. For the Categorical distribution of the component k, non-normalized probability vk,q of the gene q being expressed was estimated as a fraction of the gene q across the observed molecules (smoothed with pseudo-counts), weighted by the molecule confidence pc,i. To avoid numerical problems, all confidences were restricted with 0.01 from the bottom: p0c,i = max(pc,i, 0.01), vk,q = (Pi(: P(gij=:qj)∈&k(ip∈0ck,j))p+0c,1i)+1 . This procedure is similar to the one employed in estimation of probabilities for the cell type segmentation. However, here hard assignment of cell labels was used for the sake of performance. For the Multivariate Normal distribution, the mean μk was estimated as the weighted average over the positions of the assigned molecules: μk,u = PiP:ii∈: ki∈xki,pu0c∗,pi0c,i , where u ∈ {1, 2}. The weighted Maximal Likelihood Estimator was also used for the covariance matrix: Sk = Pi: iP∈ki(:xii∈∗kxpi0cT,i∗p0c,i) . After estimating the covariance matrix, it was adjusted based on the global scale sglobal. The most popular solution for such adjustment uses the Wishart prior over the covariance matrix, as this prior is conjugate for the Normal distribution. However, the Wishart prior is parametrized with the expected covariance matrix and the number of degrees of freedom, and thus does not allow to explicitly control the magnitude of deviation from the expected covarate matrix. Therefore we instead relied on the non-conjugate Normal-like prior on the eigenvalues of the covariance matrix. This prior was parametrized with the expected size of the eigenvalues sglobal, and their standard deviation σglobal, which can be specified by the user, or set to 0.25 ∗ sglobal by default. Then, the adjustment starts by performing eigen decomposition over Sk to calculate its eigenvalues λk and the eigenvector matrix Qk. Next, the eigenvalues are adjusted using the formula λ0k,u = λk,u + sign(λk,u − sglobal) ∗ q |λk,σug−losbgalolbal| ∗ σglobal, where u ∈ [1, 2] . This transformation corresponds to quadratic penalty over deviation Z-scores Zk,u, reducing the deviation Zk,u ∗ σglobal to pZk,u ∗ σglobal. Next, to account for the components with low number of samples, it is further adjusted as λ0k0,u = m∗σglobal+nk∗λ0k,u , where m is the expected minimal number of molecules per cell m+nk and nk is the number of molecules assigned to the component k. Finally, the adjusted covariance matrix is estimated as Sk0 = Qk ∗ (λ0k0,1)2 0 ! ∗ Qk−1. When the distribution parameters 0 (λ0k0,2)2 are estimated, the component prior probabilities ppkrior are sampled from the Dirichlet distribution, proportionally to the number of molecules per component nk: pprior ∼ Dirichlet (max(n, α)), where α is the Dirichlet Process parameter, set to 0.2 by default. Density of the background component was also estimated during this step, as it is constant and depends only on the parameters of the cell components:

fbg = fbpgosition ∗ fbggene fbpgosition = M vN ormal [3sglobal, 3sglobal]|μ = [0, 0], Σ = fbggene =

1 NXcomps Ncoms i=1 1

X |expressedi| q∈expressedi

Catk(q) ! (12) (13) (14) where expressedi is a set of all genes, expressed in the component i. The position part of this density corresponds to the level of three expected standard deviations (i.e. 3 ∗ sglobal) from the cell center. And the composition part estimates the average expression probability across all cells and all genes expressed in them.

Distribution sampling. After the Maximization step, Nnew comps = βnew ∗ Ncomps new components were sampled from the prior distributions. The parameter βnew was set to 0.3 by default. For the Normal distribution, centers were sampled from all molecule positions with the weights, proportional to the molecule confidences: μ∗ = [xi, yi]; i ∼ Cat(pc). The diagonal covariance matrix was sampled from the global scale prior: S∗ = s012 s022! ; si ∼ N sglobal, σg2lobal . Sampling gene composition parameters requires proper modeling of sparsity of the expression vectors, which varies greatly between protocols and cell types. It is therefore unclear how to capture this with a parametric prior distribution. Instead, the algorithm sampled gene composition parameters uniformly from the existing components: v∗ = vk; k ∼ U nif orm(1 : Ncomps).

Stochastic Expectation. During the Expectation step, the algorithm iterates over all the molecules and stochastically re-estimates their component assignment. The algorithm starts by determining candidate components for each molecule, as well as their MRF prior densities. First, if a component k already has some molecules {i}zi=k assigned to it, it is included as a candidate for all molecules j that are connected to the assigned molecules {i}zi=k: (zi = k) ⇒ (k ∈ candidatesj , ∀j ∈ adj(i)). The MRF prior densities are defined as the weighted sum over all edges coming from the molecules that are assigned to this component: f (zi = k|W, pc) = P(j∈adj(i))&zj=k (pc,j ∗ wi,j ). For every component k that was just sampled from the prior and still has no assigned molecules, the algorithms found the molecule i, closest to the component center μk and included component k as a candidate for all molecules j that are connected to i: i = argmin(|xl − μk|) ⇒ (k ∈ candidatesj , ∀j ∈ ∪(adj(i), i))

l∈1:Nmols In this case, the edge weight was set to 1, which is the maximal possible weight given the scaling procedure. The background component was included as a candidate to all molecules, and its MRF prior probability was set to   fMRF (zi = bg|W, pc, z\i) = max  X j∈adj(i) & zj=bg ((1 − pc,j ) ∗ wi,j ) , 1 .

  It’s important to note that the exponent was not used for the MRF densities here, as it gave visibly worse results on the tests. After the MRF weights are estimated, the full formula for the component density is: fk(obsi|W, pc, z\i) = pc,i ∗ fMRF (zi = k|W, pc, z\i) ∗ ppkrior ∗ M vN ormalk(xi) ∗ Catk(gi) fbg(obsi|W, pc, z\i) = (1 − pc,i)) ∗ fMRF (zi = bg|W, pc, z\i) ∗ fbpgosition ∗ fbggene Finally, assignment for the molecule i is sampled from the estimated densities.

Using molecule cell type segmentation information. To improve gene composition purity of the cell segmentation, results from the cell type segmentation can be utilised here by penalizing assignment of molecules from different clusters to the same cell. For that, the cell type segmentation is ran prior to the cell segmentation, and each molecule has a cluster label assigned to it. By default, number of clusters is set to 4, as almost every dataset has 4 cell types, while having larger number of clusters could lead to over-segmentation. Given a cluster label per molecule, on the maximization step the algorithm estimates the most represented cluster per cell. Then, on the expectation step, densities of the cell components are penalized in case the cluster of the molecule cli does not match to the cluster clk of the cell: fk0 (obsi|W, pc, z\i) = fk(obsi|W, pc, z\i) ∗ γI[clk=cli], where γ is the penalty term (0.25 by default) and I is the indicator function.

Using a prior segmentation. For many datasets, the auxiliary microscopy stains, such as DAPI, can be used to generate a prior segmentation. Such a segmentation can be used to resolve the cases where molecules do not provide sufficient information. The complexities of segmenting and aligning the auxiliary stains, however, mean that the quality of such segmentations can vary. To incorporate this optional information, Baysor can accept prior cell segmentation labels per molecule, as well as the prior segmentation confidence parameter cprior ∈ [0; 1]. At cprior = 0, the prior segmentation would be ignored, whereas cprior = 1 forces the Baysor segmentation not to violate the prior segmentation assignments. Increasing cprior from 0 to 1 gradually increases importance of the prior segmentation. The prior segmentation penalty is evaluated only for the molecules that are assigned to some cell (but not to background) in both Baysor and the prior segmentations. This accounts for the fact that the imaging-based segmentations may miss some cells or portions of cells that can still be deduced from the spatial transcriptomics data. The most obvious example of such situation are the DAPI-based segmentations, which cover only molecules within (15) the cell nuclei, leaving most of the cytoplasm molecules unannotated. Baysor algorithm, therefore, treats background labels (i.e. not within segmented region) within the prior segmentations as "unknown", instead of explicitly assigning these molecules to the background component. The situation in which some non-background molecules from the prior segmentation are recognized as background by Baysor is dealt with during the Background segmentation step. Specifically, if a molecule i is recognized as intracellular in the prior, its confidence can not be less than c2prior: pc,i ≥ c2prior, ∀i : ziprior 6= background. As a result, there are two possible types of contradictions between the two segmentations: (i) multiple Baysor components are present within one prior segment, and (ii) one Baysor cell component touches multiple prior segment.

The penalties are applied during the stage at which the densities of components are estimated for a given molecule, by multiplying fk(obsi|W, pc, z\i) by the penalty term βk,i. If this molecule has a prior assignment to some segment, all but the component with the largest intersection with this prior segment are penalized. This is done using the following procedure. After the Distribution Sampling step, the algorithm estimates number of molecules nske,qg per segment q for each of the cell components k. Then, based on these numbers it estimates the main prior segment per component ςk. This main segment label is used to separate the cases where two Baysor components are present within one prior segment (in this case they would have the same main segment id) from the cases where a Baysor component touches multiple segments (it would have one main segment, but would get penalized for touching all other segments). During the Expectation step for the molecule i, the algorithm finds the component u∗, which has the largest intersection with the segment ziprior = q across all candidate components for any molecule i that has this segment as their main assignment or do not have a main segment at all: u∗ = argmin nske,qg . If any two components have the same k:ςk∈(q,∅) number of molecules for this segment, the component with larger total number of molecules is chosen as the main component. The main component per segment incurs no penalty (βu∗,i = 1), while the others are penalized as the following:

βk,i = pp11 −− ccpprriioorr nnnnqssskuskeee,∗e,qgggq,gq cpcpriroiror∗exp(3∗cprior) iofthςker∈w(isqe, ∅) ( 16 ) where nseq is the total number of molecules per segment q. The first type of penalty prevents two q cells from existing in the same prior segment (over-segmentation problem), while the second type prevents one cell from taking over a big part of a segment assigned to a different cell (overlapping problem). The second penalty type also solves the under-segmentation problem: if a Baysor cell overlaps two prior segments, then as soon as a new component is sampled inside one of these segments, the existing "doublet" cell incurs very large penalty for all other molecules of this segment. These functional form of the penalties, while somewhat arbitrary, has allowed to express a desired curve shape for the penalties (see Supplementary Fig. 14b,c). The penalties are also scaled in [0; 1], which allows to keep the Distribution Sampling step without changes. Freshly sampled cells have no penalty, which corresponds to βk,i = 1.

Complete Baysor workflow. Below is the description of the full Baysor workflow, as implemented by the command-line interface: 1. Read the data frame with information about molecules. Filter out the genes with total number of molecules below threshold if it was specified. 2. If provided, load the prior segmentation mask. Filter out any prior segments that have less than mprior molecules. Estimate the global scale sglobal based on this data. 3. Estimate confidence pc,i per molecule 4. If the specified number of cell clusters (4 by default) is greater than 1, run the cell type segmentation using the confidences, estimated above. Otherwise, assign all molecules to the same cluster. 5. Initialize the Bayesian Mixture Model algorithm 6. Optionally, split data by frames to enable parallel processing 7. Run the BMM algorithm 8. Re-estimate assignment by averaging over the last N est iterations 9. If the data was split by frames, merge them back into one 10. Save the segmentation results Inferring the algorithm parameters. Baysor implementation derives most of the parameter estimates based on the minimal expected number of molecules per cell m, which must be specified by the user, as well as the number of genes measured by the assay (ngenes): • The initial number of cells is set to Ncienlilts

= Nmolecules/m, where Nmolecules is the total number of molecules. • The number of principal components for the PCA transformation of NCVs for adjusting the

MRF is estimated as min (max (ngenes/3, 30) , 100, ngenes) • The nearest neighbor index for estimating distances di during the Background Segmentation is set to max (m/2 + 1, 2) • The number of nearest neighbors for NCV coloring is set to max(ngenes/10), m, 3) • If a segmentation mask was provided, it is used to infer the global scale sglobal and its standard deviation σglobal. To do so, Baysor first approximates cell radii from the number of pixels per prior segment n pix i

pix as ri = ni /π. Then, sglobal = median(ri) and σglobal = 1.4826 ∗ i MAD(ri), which are the median and the adjusted median absolute deviation over the radii.

i • Without a segmentation mask, the global scale sglobal is a required input parameter, and the deviation σglobal is set to 0.25 ∗ sglobal by default.

• The minimal number of molecules per prior segment mprior = max(m/4, 2).

Segmentation parameters. Parameters of the segmentation runs for different datasets are shown in the Supplementary Table 3.

Benchmarks. To evaluate Baysor, we have shown its performance in separating distinct cell types, using Pearson Correlation benchmarks. We have also shown stability of its convergence for different random number generator seeds, as well as runtime profiling.

Correlation benchmarks. As we currently lack independent techniques to establish ground truth on cell segmentations, we developed a benchmark to evaluate the differences between two segmentations. For each source cell c in the segmentation zsrc the benchmark finds a target cell t∗ from the segmentation ztarget that has the largest overlap with c (Fig. 4c). Then it splits the molecules from c over the part that overlaps with t∗ from the molecules in the non-overlapping part. Measuring similarities of gene compositions of these two parts as a Pearson linear correlation (ρc), the benchmark evaluates the distribution of such similarities over all cells from the source segmentation: ρ = {ρc}c∈zsrc . Taking two segmentations z1 and z2 as the input, the benchmark returns two distributions ρ1 and ρ2 for the cases where zsrc = z1, ztarget = z2 and zsrc = z2, ztarget = z1 correspondingly. Based on these two distribution, the segmentation is said to be better if it has higher ρi values. More formally, given two segmentations zsrc and ztarget, the benchmark procedure does the following: 1. All cells i with the total number of molecules below the threshold (nimols < m) are removed 2. Among the remaining cells, a contingency matrix between the two segmentation assignments from both segmentations.

is estimated. 3. Based on the contingency matrix, for each cell c of the segmentation zsrc and cell t from the segmentation ztarget, the overlap fraction fco,vter is estimated. 4. For each cell c ∈ zsrc, the target cell t∗ ∈ ztarget is determined as the cell with the highest molecule overlap fraction fco,vt∗er = |{i:(zsrc,i=c) & (ztarget,i=t∗)}i∈1:Nmols | .

ncmols 5. Only cells from zsrc with the overlap fco,vt∗er ∈ (25%, 75%) are selected for further analysis. 6. For each cell c ∈ zsrc, the molecules from this cell are partitioned into (i) the main overlap lcm,ta∗in = {i : (zsrc,i = c) & (ztarget,i = t∗)}, and (ii) the rest lcr,ets∗t = {i : (zsrc,i = c) & (ztarget,i 6= t∗)}. 7. Expression vectors vcm,ta∗in and vcr,ets∗t are estimated over the molecules from lcm,ta∗in and lcr,ets∗t correspondingly. 8. Pearson linear correlation coefficients ρc between the vectors vcm,ta∗in and vcr,ets∗t are estimated.

If these two sets of molecules were produced by different cell types, the correlation ρc is expected to be low, suggesting that the segmentation zsrc was erroneous for the cell c. However, similarly to the classical hypothesis testing framework, a high correlation ρc does not mean that the segmentation zsrc is correct: it only suggests that the molecules were obtained from the same cell type, and it is still possible that they originated from different cells. Stability benchmarks. To assess the stability of the algorithm convergence, Baysor was run 10 times with different random number generator seeds on the MERFISH dataset. Each run used 400 iterations each, with parameters "scale=6.16" and "min_molecules_per_cell=30". The final assignment was obtained by averaging over the last 100 iterations ("assignment_history_depth=100"). The molecule confidences and clusters were pre-estimated only once, as these procedures are deterministic. Pairwise Adjusted Rand Index and Mutual Information between results of different runs were calculated as stability metrics.

Performance benchmarks. For runtime performance profiling of the molecule clustering (Supplementary Table 1), MERFISH dataset was used with the molecules subset to those with coordinates x < −3300 and y < −3300 (338023 molecules in total). Then, molecule confidence was estimated, using confidence_nn_id=26 parameter setting. The MRF clustering was then ran five times for each value of k ∈ {2, 4, 6, 8, 10}, with parameters max_iters=5000 and n_iters_without_update=100.

To evaluate the performance of the NCV Leiden clustering, Pagoda2 package was used. For each run, the NCV matrix was estimated using 50 nearest neighbors, normalized by the total count of each NCV, and reduced to projections on the top 50 principal components. A k-NN graph was then built using 30 nearest neighbors (using the default cosine similarity distance metric) and leiden community detection algorithm was applied with the default resolution=1.0. For the profiling, the same parameters as specified in the Supplementary Table 3 (Prior=No) were used. Each dataset was segmented 5 times using a single thread. All benchmarks were ran on Lenovo Thinkpad X1 laptop with Intel(R) Core(TM) i7-8850H @ 2.60GHz CPU and 64GB RAM.

Quality metrics. After finishing the segmentation, Baysor extracts the following summary metrics per cell: • Cell area: area of the convex hull around the cell molecules • Density: cell area divided by the number of molecules in cell • Elongation: ratio of the two eigenvalues of the cell covariance matrix • Average confidence: average confidence of the cell molecules scRNA-seq processing. Analyzing NCVs using scRNA-seq tools. To generate clustering and embedding of NCVs (Figure 1b,c) Pagoda2 package (https://github.com/kharchenkolab /pagoda2/) was used. The data was pre-processed using the total count normalization, and then building 50-NN graph over the normalized counts using the default cosine distance metric. Leiden clustering with parameters "resolution=8" and "n.iterations=15" was used for annotation expansion (see below). Visualization of the MERFISH and the osm-FISH datasets was done by building a joint matrix over the paper and the Baysor segmentations, reducing its dimensionality with PCA to 10 top PCs (MultivariateStats.jl package) and then running UMAP embedding (UMAP.jl package) using Euclidean distance and parameters "spread=2.0" and "min_dist=0.1" over the joint matrix. Hierarchical clustering with the Ward linkage method and 70 clusters (Clustering.jl package) was used for the annotation expansion. scRNA-seq Annotation. The annotation for Allen smFISH data was generated using CellAnnotatoR package. The cell type markers were inferred from Mouse VISp scRNA-seq data produced for the SpaceTx consortium by the Allen Brain Institute (personal communication), and then applied it to the corresponding FISH data. For inference, the "merged_cluster_smFISH" annotation shared with the dataset was used. First, the cell type hierarchy was built using "broad_class" as the first level, class prefix from "merged_cluster_smFISH" as the second and the full "merged_cluster_smFISH" labels as the third level of the hierarchy. Second, cell types "CR", "Astro", "Endo", "Macrophage", "Oligo" and "SMC" were removed, as they could not be distinguished by the genes, measured in the Allen smFISH dataset. Next, CellAnnotatoR marker inference procedure was utilized to obtain the markers using this hierarchy, and the resulting marker list was adjusted by hands to improve the quality of the annotation. Finally, CellAnnotatoR was used to apply the identified markers to the NCV data. To annotate MERSISH and osm-FISH data, the marker list was compiled manually using the information published by the protocol authors, and then CellAnnotatoR was applied in the same manner. In all three cases, the resulting annotation was expanded over the clusters found in the previous step (see "Analyzing NCVs using scRNA-seq tools").

Polygon visualization. To visualize the boundaries of the cells, a Kernel Density Estimation (KDE) based algorithm was used. The algorithm builds a grid over 2D space, assigns each grid node to the cell with the highest density of molecules around this node and draws a polygon around this labels on the grid. More formally: 1. A uniform 4-connected grid was created over the 2D space, with the grid step specified as an input parameter. 2. The density of molecules for each cell was estimated over the grid nodes using the KDE implementation in the KernelDensity.jl package. KDE bandwidth - a parameter of the algorithm - was set to the 0.5 * "grid step" by default. To improve runtime performance, only the nodes within three bandwidths of the cell molecules were taken into account. 3. Each node was assigned to the cell with the maximal density in this node. If the maximal density was below threshold (10−5 by default), the node was assigned to the background. 4. For each cell, the graph of boundary nodes were determined as the grid nodes of the cell that were adjacent to the nodes from the other cells or from the background. The edges between these boundary nodes were accepted into the boundary graph. 5. For each cell, a minimal spanning tree was built over its boundary graph using the Kruskal algorithm. 6. For each cell, the longest path in this tree was extracted using the Dijkstra algorithm, and the resulting path was transformed into a polygon by connecting its beginning and its end. DAPI watershed segmentation with ImageJ. To generate prior segmentations using DAPI stains, the Watershed segmentation was performed using ImageJ23 software. Each staining was segmented by the following procedure: 1. The image was converted to 8-bit format ("Image" / "Type" / "8-bit"). 2. Median filter with one pixel radius was applied ("Process" / "Filters" / "Median..."). 3. Auto Threshold with the "Default" method and "Ignore black" option was applied for image binarization ("Image" / "Adjust" / "Auto Threshold") 4. Watershed segmentation was applied ("Process" / "Binary" / "Watershed").

VISP multiplexed smFISH data generation. Multiplexed smFISH data of mouse primary visual cortex (VISp) was generated as part of the SpaceTx consortium. Tissue processing was carried out as previously described27, with some modifications. The description is taken from19.

Silanizationof coverslips (#1.5, Thorlabs CG15KH) was performed by plasma cleaning for 30 min in a Plasma-Prep III (SPI 11050-AB), followed by vapor deposition of 3-aminopropyltriethoxysilane (APES, Sigma A3648) in a vacuum for 10 minutes. Coverslips were then washedin 100% methanol for 2 × 5 minutes, allowed to dry, and stored in a dust-free environment until use.

Fresh-frozen mouse brain tissue was sectioned at 10 μm onto silanized coverslips, let dry for 20 min at -20◦C, then fixed for 15 min at 4◦C in 4% PFA inPBS. Sections were washed 3 × 10 min in PBS, then permeabilized and dehydrated with chilled 100% methanol at -20◦C for 10 min and allowed to dry. Sections were stored at -80 ◦C until use. Frozen sections were rehydrated in 2X SSC (Sigma 20XSSC, 15557036) for 5 min, then treated 10 min with 8% SDS (Sigma 724255) in PBS at room temperature. Sections were washed 5 times in 2X SSC. Sections were then incubated in hybridization buffer (10% Formamide (v/v, Sigma 4650), 10% dextran sulfate (w/v, Sigma D8906), 200μg/mL BSA (ThermoFisher AM2616), 2 mM ribonucleoside vanadyl complex (New England Biolabs S1402S), 1 mg/ml tRNA (Sigma 10109541001) in 2X SSC) for 5 min at 37◦C. Probes were diluted in hybridization buffer at a concentration of 250 nM and hybridized at 37◦C for 2 h. Following hybridization, sections were washed 2 × 10 min at 37◦C in wash buffer (2X SSC, 20% Formamide), and 1 × 10 min in wash buffer with 5 μg/ml DAPI (Sigma 32670), then washed 3 times with 2X SSC. Sections were then imaged in Imaging buffer(20 mM Tris-HCl pH 8, 50 mM NaCl, 0.8% glucose (Sigma G8270), 30 U/ml pyranose oxidase (Sigma P4234), 50 μg/ml catalase (Abcam ab219092). Following imaging, sections were incubated 3 × 10 min in stripping buffer (65% formamide, 2X SSC) at 30◦C to remove hybridization probes from the first round. Sections were then washed in 2X SSC for 3 × 5 min at room temperature before repeating the hybridization procedure.

The multiplexed smFISH image data was collected and processed using methods previously described27, except that the images from different rounds of hybridization were registered in (x,y) based on the DAPI signal. The spot locations and raw data are available on request.

Code availability. Baysor package is available at https://github.com/kharchenkol ab/Baysor. Baysor parameters for different datasets are reported in the Supplementary Table 3. The code to reproduce the results is available at https://github.com/kharchenkolab /BaysorAnalysis/.