Spatial genomics maps the structure, nature and evolution of cancer clones

Genome sequencing of cancers often reveals mosaics of different subclones present in the same tumour1–3. Although these are believed to arise according to the principles of somatic evolution, the exact spatial growth patterns and underlying mechanisms remain elusive4,5. Here, to address this need, we developed a workflow that generates detailed quantitative maps of genetic subclone composition across whole-tumour sections. These provide the basis for studying clonal growth patterns, and the histological characteristics, microanatomy and microenvironmental composition of each clone. The approach rests on whole-genome sequencing, followed by highly multiplexed base-specific in situ sequencing, single-cell resolved transcriptomics and dedicated algorithms to link these layers. Applying the base-specific in situ sequencing workflow to eight tissue sections from two multifocal primary breast cancers revealed intricate subclonal growth patterns that were validated by microdissection. In a case of ductal carcinoma in situ, polyclonal neoplastic expansions occurred at the macroscopic scale but segregated within microanatomical structures. Across the stages of ductal carcinoma in situ, invasive cancer and lymph node metastasis, subclone territories are shown to exhibit distinct transcriptional and histological features and cellular microenvironments. These results provide examples of the benefits afforded by spatial genomics for deciphering the mechanisms underlying cancer evolution and microenvironmental ecology.


Statistics and Reproducibility
BaSISS fields   Tissue blocks were first sliced to obtain sufficient nucleic acids for bulk whole genome sequencing (WGS) and RNAseq analysis. Subsequent serial 10um sections were obtained from the same tissue block and orientation for base specific in-situ sequencing (BaSISS), in-situ sequencing (ISS), histological assessment and immunohistochemistry.

Histopathological assessment
For each tissue block a frozen section was prepared for pathological review using standard haematoxylin and eosin (H&E) staining. H&E review and annotation was performed by experienced breast histopathologists. Annotation of more than 300 microregions was performed without any prior knowledge of the genetic clonal structure (Supplementary Table 2). For DCIS annotation, regions were pre-selected for annotation to enrich for those that are predicted to be contained within a myoepithelial bordered space (ducts/acini or lobules) based on the H&E and/ or ISS histological images ("is single") (Extended Data Fig. 4a; Supplementary Table 2). Additional regions were selected to ensure comprehensive review of the features of the tumour. Tissue sections were fresh frozen so some features such as mitoses and chromatin structure were difficult to assess and not presented in the main analyses.
Consequently, nuclear grade in DCIS is largely based on nuclear size. Nuclear vacuoles were defined as present or absent based on a 20% cut-off. Differentiating between a TDLU and a duct is not always obvious on a single section -larger DCIS close to terminal structures were marked as a TDLU.
Structures that are not branched and are more peripheral are defined as ducts. Lobule contours as indicated in Fig. 4a were predicted based upon spatial proximity, shared surrounding structures and/or interconnections between glandular structures observed across the available z-stack of serial tissue sections. In P1-ER1 and P1-ER2 microregions were either discrete regions such as DCIS or hyperplastic ducts while similarly sized regions of invasive tissue were demarcated manually to permit sampling across the entire tissue surface area (Extended Data Fig. 5b).
In both P2-LN1 and P2-TN1 the histopathologist independently reported the existence of more than one growth pattern. In P2-LN1 the pathologist defined the selected areas according to these growth patterns, aiming to generate similar sized regions (Fig. 5a) Bulk sequence data metrics are reported in Supplementary Table 3.

Somatic Mutation Identification
Genome-wide somatic substitutions were called using CaVEMan (Cancer Variants Through

Genomic Data
Subclone composition was determined from multi-region WGS somatic substitution and copy number data as reported in the original publications 4,5 . The approach detects a finite number of subclones each with defined parameters including the number of mutations within that subclone and the fraction of cells from each sample that contain those subclone specific mutations -the results for each sample are summarised in Supplementary Table 3.

Cancer Cell Fractions
For bulk genomic sequence data the variant allele fraction (VAF) of each mutation, i, is calculated from the number of reads reporting the mutant rmut and the number of reference, or wild-type reads rwt is calculated as However, the VAF cannot be used directly to infer subclonal architecture because in addition to reflecting the proportion of cells that carry that mutation it is also influenced by the tumour purity and the local copy number -with the latter encompassing both the total number of alleles and the number of alleles that carry the mutation 6 . Consequently, mutations from the same subclone may have different VAFs reflecting their different copy number states, also termed their "multiplicity". To allow a comparison of mutations across related samples with different copy number profiles and tumour purity levels we derive the mutation copy number. For mutation i, the mutation copy number ni is obtained as Where ρ is tumour purity, and nlocus,t and nlocus,n are the locus specific total copy numbers in the tumour and normal sample respectively. Subclonal copy-number changes and tumour purity are defined using the Battenberg algorithm as previously described (https://github.com/cancerit/cgpBattenberg).
Confidence intervals for mutation copy numbers are generated using bootstrap resampling of mutant and wildtype reads from each mutant locus (10,000 times) and applying the above formula to resampled reads.
For mutations that are present on multiple alleles arising through duplication (multiplicity ≠ 1) or in regions of subclonal copy number change, the mutation copy number may not accurately reflect the fraction of cells in which the mutation resides. To group mutations by the number of cells containing those mutations, the number of chromosomes bearing a mutation, nchr, must be determined. For each mutation within a region of amplified major copy number C, the observed number of mutant and wildtype reads are compared to the expected VAFi that would result from the mutation being present in 1,2…C chromosome copies also allowing for non-integer values in regions with subclonal copy number values, assuming a binomial distribution. The value of nchr is determined to be that with the maximum likelihood. The fraction of cells bearing mutation i, termed the cancer cell fraction or CCF can then be calculated as Applying this approach to cancers that are believed to derive from a common ancestor, one would therefore expect to identify some clonal mutations that are present in all cancer cells and at a certain multiplicity result in a CCF around 1. the ith mutation is drawn from a binomial distribution as follows

Multidimensional Mutation Clustering
Where Ni is the number of mutant reads and ςi is the expected fraction of reads that would report a mutation if it was present in 100% of cells (at that specific locus given the copy number state and tumour purity). The fraction of tumour cells carrying mutation i, πi, is modelled as arising from a Dirichlet process (DP) with a given probability distribution P0 and a dispersion parameter α. The approach allows co-estimation of the number of cellular populations and their properties using the stick-breaking representation of the Dirichlet process as described in Dentro et al. 6,8 To accommodate multi-region samples, the model is extended into multiple dimensions as previously described 5,8 . Essentially, the number of mutant reads obtained from different tumour region samples are modelled as independent binomial distributions with independent π drawn with a Dirichlet process.
Gibbs sampling is used to estimate posterior distributions of the features of interest. The Markov chain is run for 10,000 iterations (excluding the first 2,000 iterations). The median densities for CCFs is estimated using a Gaussian kernel implemented in R (stats and kernSmooth libraries). The clustering procedure is implemented in the DPClust software package (available here: https://github.com/Wedge-Oxford/dpclust). Version 2.2.8 was used to obtain the analysis presented.
In P2 (samples P2-TN1, P2-TN2, P2-LN1) the final DPClust solution identified a single subclone within the lymph node sample (P2-LN1) however there was evidence supporting the possible existence of a subclone albeit with lower confidence, understood in the context of low tumour purity (~11%). This is most convincing from copy number analysis that identified multiple sub-chromosomal copy number changes (chromosome 1,9,11,12,20) in addition to a 50% subclone at the telomeric portion of 17q in association with the HER2 amplifying breakage fusion bridge event ( Figure 5B). This was further supported by a build-up of posterior density in P2-LN1, indicating that across some iterations DPClust has explored the presence of a subclone in this sample was explored, although ultimately it decided the subclone wasn't required to explain the observed point mutation data.

Principles of Phylogenetic Tree Construction
To infer the evolutionary relationship between subclones, phylogenetic trees were constructed using the "pigeonhole principle" which states that if there are n objects that are placed in m containers if n>m then at least one container must contain more than one object. Extending this to subclonal reconstruction we can appreciate that the sum of subclone CCFs cannot exceed the CCF of their ancestor. For example, if there are 3 subclones with CCFs of 1, 0.9 and 0.5 then there must be a linear evolution -as 1 + 0.9 > 1, CCF subclone 0.9 must be a descendant of the CCF subclone 1; furthermore, as 0.9 + 0.5 > 1, CCF subclone 0.5 must be a descendant of the CCF subclone 0.9. In contrast one may see the possibility of branching evolution where the sum of subclone populations is less than 1. For example, 3 subclones with CCFs of 1, 0.5 and 0.4 could arise either through linear evolution: 1 > 0.5 > 0.4 or along separate ancestral lines because 0.5 + 0.4 < 1. When applying these principles across multi-region samples the compatible underlying tree structures is usually greatly restricted because of the requirement that the same tree structure must be compatible across all samples. Applying these principles one can identify the phylogenetic tree structure(s) most compatible with the underlying data.

Padlock probe design
The padlock probe consists of a single stranded DNA oligo with two target recognition arms in the 3' and 5' end of the oligo, enabling circularization of the padlock probe upon target hybridization. Each target recognition arm typically has a length of 15-25 nucleotides. The two arms are interspaced with a linker sequence comprising a 20-nucleotide anchor primer sequence and a target-specific 4 nucleotide barcode followed by a 5 nucleotide stabilising sequence for sequencing-by-ligation. The used barcodes were selected in such a way that they differ in at least two positions.
In total, three padlock probe gene panels were used, marker/oncology gene-, mutation-(one specifically designed for each case) and immune-expression panels (Supplementary Table 1). Padlock probes were ordered as ultramer DNA oligos from Integrated DNA Technologies (IDT, Leuven, Belgium) with 5'phosphorylation modification and were reconstituted in Tris/EDTA buffer.

Mutation panel designs
Two separate case-specific mutation panels were designed as detailed in Supplementary Table 1. For wildtype-and mutation-specific padlock probes, the sequence upstream of the mutated site resembles the 3' arm of the padlock probe, with the wildtype or mutant nucleotide at the 3' end, whereas the sequence downstream of the mutated site resembles the 5' arm. The target recognition arm lengths were adjusted to have a similar melting temperature of sim ~55°C. For each mutation site, one wild-type and one mutation-specific padlock probe was designed. For P1, in addition to the anchor and barcode sequences a 20-nucleotide sequence used for a hybridization cycle was added in the padlock probe linker sequence. Specific primers were used for the in situ reverse transcription and were designed as the reverse complement sequence of the 5' target recognition arm. Primers were ordered as DNA oligos from Integrated DNA Technologies (IDT) and were reconstituted in Tris/EDTA buffer. In P1, we included 3 mutations that were not assigned to a cluster/ branch by DPclust of WGS data. In accordance with this, BaSISS signals reporting these mutations were found to be non-specific and uninformative, probably relating to them existing at loci with variable copy number states, and they were dropped from further downstream analyses.

Immune panel design
Large scale probe design was facilitated using an in-house Python software package as described previously (Ref PMID 31740815) which utilises ClustalW and BLAST+ to ensure probe specificity.
Each padlock probe of the immune panel was designed to contain two 20 nucleotide long target recognition arms. Only target fragments with melting temperature between 65°C and 75°C were considered. Probes were selected aiming to obtain a distribution along the whole length of the transcript.
Overall, five padlock probes were selected per target gene with the exception of HLA-DRB1, where only two specific probes could be designed. For the T cell specific genes CD8A, FOXP3, EOMES, CD4 and IFNG 20 probes were selected in order to increase the detection efficiency for these target genes.
All padlock probe sequences of the immune panel are shown in Supplementary Table 1.
A combination of random decamer primers (IDT, Leuven Belgium) and specific primers were used for in situ reverse transcription. Specific primers were designed to hybridise to the mRNA 15-20 nucleotides downstream of the target sequence. Primers were ordered as DNA oligos from Integrated DNA Technologies (IDT) and were reconstituted in Tris/EDTA buffer.
Target genes for the immune panel were selected to cover a broad range of immune cell subtype markers with special emphasis on T cell subsets and their regulation (Supplementary Table 1).

Oncology gene panel design
The oncology gene panel has previously been published and includes genes involved in proliferation, EMT, invasiveness, stemness, angiogenesis as well as genes for breast cancer subtyping and oncotypeDX recurrence scoring 9 . The target recognition arms were designed to capture most splice variants of the gene transcripts and blasted to confirm their specificity. Each target recognition sequence had a GC content of 50-55% and a melting temperature of ~55°C. In this older design, the panel included one padlock probe per gene target and the barcodes used were only differing in one position. For the marker gene panel, random decamer primers were used for in-situ reverse transcription (IDT, Leuven Belgium).

In-situ sequencing (ISS)
ISS was performed as described by Ke et al. 10 and modified according to Svedlund et al. 9 , and was used to spatially resolve oncology panel-, mutation panel-and immune panel-gene expression profiles on consecutive sections from the breast tumour tissue blocks.
The ISS library preparation and sequencing is described in detail at protocols.io.bb2giqbw, the steps in the protocol that were modified for breast cancer tissues are indicated below. In brief, library preparation included fixation of tissue sections with 4% PFA for 30 min (step 2) followed by permeabilized with

Image registration
The registration across imaging cycles was performed in two steps: affine registration on DAPI channel and subsequently local warping on anchor channel. For both steps we used algorithms provided in libraries OpenCV-contrib (version 4.3.0) 12 and scikit-image (version 0.17) 13 . In all imaging cycles, before the registration, both DAPI and anchor channels were maximum intensity projected across the image z-stack.
During the affine registration step, we coarsely align images of all cycles to the first one based on the DAPI channel. Firstly, we detect key points in the images of each cycle using the FAST feature detector.
Secondly, for each key point, its surrounding area is described with histograms of oriented gradients using the DAISY feature descriptor. After that, using the key points and their descriptors, the FLANN-based matcher finds correspondences between pairs of key points from reference and moving images and filters out unreliable points. Lastly, the remaining key points are processed using the RANSACbased algorithm that aligns them and estimates affine transformation parameters with 4 degrees of freedom.
The second registration step aligns imaging cycles sequentially using the anchor channel (fluorophore Cy7). We applied the Farneback optical flow algorithm to achieve more accurate registration by warping the images locally, so that RNA spots of different channels can be better aligned despite the presence of nuclei swelling, imperfect stitching and sample distortion. In both steps, we optimised the algorithm by performing computation on the tiled images to reduce memory consumption and accelerate the transformation parameters estimation.

Serial tissue image signal alignment
Analysis of sample P2-LN required IHC signal projection performed on a consecutive slide back to BaSISS slide. To achieve this, we performed a spline-based elastic registration implemented in ImageJ package UnwarpJ 14 .

Downstream Analysis
Allele co-occurrence analysis The relative co-occurrence of BaSISS alleles was calculated using the M-function 16 at a 20µm radius. Inference and confidence interval estimation was performed with dbmss R package 17 . For LCM allele co-occurrence a simple Pearson correlation of variant allele frequencies between regions was used.

Nuclei segmentation
A two-stage pipeline from Caicedo, Goodman et al. 18

Single cell typing
The decoded ISS signals do not directly reveal which cells they derive from require further processing to recover these associations. To do so, ISS data was overlaid with the nuclei segmentation masks and signals were assigned to the nearest nuclei. A conservative distance cut-off of 5µm was used to decrease the chance of missannoation which led to the ~30% of signal loss.
Since ISS data is sparse and restricted to a limited number of targeted genes, a previously published scRNA breast cancer atlas 19 was used for the label transfer. The differences in gene detection rates between scRNA reference and ISS data combined with the low signal per nucleus counts, limited the scope of applicable computational approaches. Therefore, the label transfer was based on the most informative marker genes and was performed in a hierarchical manner as follows.
Single cell annotation data often has hierarchical structure reflecting the hierarchical nature of the cell differentiation. We could leverage this structure as a way of regularisation for the cell type prediction model where the fine cell subtyping (e.g. CD4 T-cell) is dependent on the top annotation level (e.g. Tcell). A hierarchical multinomial logistic regression model trained on the reference scRNA dataset restricted to the genes available in the ISS panel would encode the gene importance for the cell type characterisation on a respective level of hierarchy. We used this approach to assist marker gene selection and arrived with the panel of genes for broad (Immune, Epithelial and Stromal) and specific (B-cells, Fig. 2b,c). The code is available at github.com/dissatisfaction-ai/scHierarchy.

Myeloid, T-cells, Fibroblasts + PVL and Endothelial) cell type levels (Extended Data
The final cell type assignment was performed in 2 steps. At first iteration, nuclei were classified into 3 broad categories ('epithelial', 'immune', 'stromal') if they had any corresponding marker genes assigned. At the second iteration, nuclei with 'immune' and 'stromal' assignments were further subdivided into specific groups based on the presence of the corresponding markers. The identity of nuclei that did not have any marker genes in proximity or had a contradictory assignment was considered 'unknown'.

Digital IHC based classification
To quantitatively describe DAB-stained areas, Haematoxylin-Eosin-DAB colour deconvolution was used 20 . The implementation from the scikit-learn Python package was used.

Spatial mapping of cancer clones
Raw BaSISS data informs us of allelic variant's spatial distributions. However, this data cannot be directly interpreted in terms of cancer clones. To make sense of BaSISS data we designed a Bayesian model that aims to decompose local BaSISS signal counts into clone densities and corresponding genotypes.

Core model
The BaSISS protocol records a series of fluorescent spots, which are decoded into a class of different barcodes corresponding to each targeted allele. This information produces a table of tuples where ! is the allele of spot and ! and ! are its respective two-dimensional coordinates.
The counts of all probes can be represented by a three-dimensional array ∈ |7|×|9|×|:| , where refers to the allele and and are coarse grained coordinates on a regular grid of dimensions | | × | |.
The grid size was chosen to be 108.8µm considering a trade-off between data sparsity, precision and computational cost.
The essential idea is that the expected number of BaSISS signals is decomposed into maps of | | distinct clones ∈ |*|×|9|×|:| each with a distinct genotype ∈ |7|×|*| , The genotype matrix contains the number of each allelic copy in each clone . is a matrix representation of the underlying phylogenetic tree and provides the instructions of the allelic configuration in each branch of the tree.
The maps provide the relative prevalence of each clone in a given area of the grid. As should replace a spatially continuous map, it is modelled by softmax transformed two dimensional latent Gaussian processes.

Latent Gaussian Process modelling of spatial clone maps
Depending on the choice of targeted alleles, signal counts in matrix may appear rather sparse. For a robust clone prevalence map inference, we use Gaussian processes as it allows us to naturally utilise spatial distance information. Each clone is modelled as independent latent two-dimensional Gaussian process with = 0 and covariance matrix with RBF kernel over grid coordinates. Bandwidth was empirically chosen for sample P1 and adjusted for sample P2 according to the median signal/nuclei ratio.
Each Gaussian process was mapped to simplex space using a softmax (multidimensional logistic) transformation with temperature parameter = 2 which represents our weak prior belief that clones could cover the whole range of [0,1] prevalence.

Sources of noise
To fit the above model to data a number of sources of noise need to be considered. The exact nature of and statistical models for the different noise terms will be discussed in the next subsections.

Cellular density
The prevalence of clones is most insightful when expressed as the density of cells per area. A DAPI stained image of nuclei is obtained as part of the BaSISS experimental workflow and utilised to identify nuclei. We use a neural network ensemble model to segment these images and store nuclei counts in a two-dimensional array ∈ |9|×|:| . However, these counts might not always correctly quantify the true cell densities due to possible segmentation errors and the fact that only a fraction of cells might be captured on the slide. Therefore, we treat the nuclei counts 9,: as a Poisson distributed values with mean 9,: representing cell densities which in turn has a weak Gamma prior.

Differential probe specificity
BaSISS' ability to detect allelic variants relies on a multistage process of reverse transcription, padlock probe annealing, amplification and fluorophore detection. The resulting intensity will be the product of these steps and differ between each probe/allele. As only the product of all steps can be identified, we model the variant detection rate ∈ # |7| for each allele 7 ∼ ( = 0.5, = 1).

Allelic confusion
The difference between wild type and mutated allelic variants is often just a single nucleotide. This increases chances for the padlock probe to anneal to the wrong allelic variant resulting in allelic confusion. To model this confusion, we design a sparse transition matrix ∈ (0,1) |7|,|7| populated with transition probabilities { 7,7′ } for each allele. We put a strong regularising Beta prior on 7,7 to shift probability density close to 1 and prevent it going below 0.75 to avoid non-identifiability.

Clone allelic variations
Although the mean clonal expression level is registered in detection rate matrix , it is reasonable to assume that different clones may have some level of diversity in expression of some genes. We model this deviation as a matrix ∈ # |*|×|7| with a log-normally distributed prior such as [ ] = 1.

Homogeneous and inhomogeneous background adjustments
Raw BaSISS data comes from biological images which have additional sources of variability. The first one is the homogeneous additive effect ∈ # |7| which increases the expected value of a particular probe detection globally on the slide.
The second one is a spatially inhomogeneous background which adjusts for local base signal detection variability. In practice it's hard to reliably model this type of background without a harsh regularisation due to its flexibility. We encode this inhomogeneous background as additional pseudo-clones which have corresponding relative prevalence and pseudo-genotype ∈ # |@|×|7| . Prevalence is modelled together with the real clones as a two-dimensional latent Gaussian process. The pseudo-genotype is sampled from a Beta prior which support is stretched to match median copy number among the loci.
We use | | = 1, but it can be increased for a more flexible background effect.
The over-dispersion 7 is sampled from Gamma distribution where distribution with mean 100 and a standard deviation of 10 7 ∼ ( = 100, = 10).

Optional inputs
Depending on which additional data accompanies investigated samples we may want to include them into the proposed Bayesian framework. Currently, two options are implemented: 1. VAFs from adjacent slices

IHC based cell type counts
VAFs from adjacent slices WGS data contains information of the mutated allele frequencies on the whole slide. We can use this data to increase the accuracy of our inference and stabilise the overall prevalence of clones. To pass it into the model we calculate the expected aggregate frequencies of each mutated variant on the slide, A7BCBB ∈ [0,1] |3| , by summation of × product over spatial dimensions.
. We then construct a Beta pseudo-likelihood of A7BCBB with and parameters proportional to the number of mutated and wild-type reads in the WGS experiment. We use the parameter = 10 to inflate uncertainty in WGS data as it comes from proximal but not exactly the same slide, we are observing in To adjust for a disproportionately higher impact of spatial information on the total likelihood we multiply VAF log-likelihood by the number of tiles on the slide.

IHC based cell type counts
IHC data reveals part of cell type heterogeneity. In our study we mainly worked with CD45 IHC staining which indicates the location of lymphocytes. As we know that none of the lymphocytes should contribute to cancer clone weights, we could further regularise clone weights. This is done by

Multi-sample version
In many cases some of the clones belonging to the same phylogenetic tree might not exist or exist only in a small proportion on most of the slides. Having a multi-sample generalisation of the proposed model helps us to correctly infer parameters associated to such clones. Essentially, the only parameters which are shared between all slides are: • genotype matrix , • probe confusion transition matrix • mean probe detection rate and • overdispersed sampling fluctuations .
All the other parameters become slide specific, i.e., cellular density , differential probe specificity , clone allelic variations and homogeneous and inhomogeneous background adjustments .
To allow some level of slide specific probe efficiency variation we multiply mean probe detection rate by a slide specific probe deviation matrix ∈ # |I|×|7| with a log-normally distributed prior such as

Inference
Variational Bayesian Inference is used to approximate the posterior. We use the mean-field version of Automatic Differentiation Variational Inference (ADVI) implementation from pymc3 package 21 .
Briefly, within the ADVI framework, the posterior distribution over unknown parameters is approximated by an appropriately transformed multivariate normal distribution with a diagonal covariance matrix. Inference is achieved by minimising the KL divergence between the variational approximation and the true posterior distribution which is equivalent to the maximisation of the evidence lower bound (ELBO).
Training is stopped when the ELBO stops increasing passing at least 15,000 iterations of training using the ADAM optimiser with a learning rate of 0.01. The marginalised variational posterior of the parameters was used in the downstream analysis. Our experiments showed that several initialisations may be needed to arrive at an optimal and solution. This is mostly due to the stochastic choice of the initial values but may also be influenced by numerical artefacts arising during influence with pymc3.

Model limitations and assumptions
Quantitative estimation of the spatial cellular composition of genetically distinct clones based on allelic information from in-situ RNA sequencing is a challenging problem. The aforementioned model offers a solution that deals with the multiple sources of biological and technical noise to produce clonal maps which have been orthogonally validated (see validation section). However, it comes at the cost of introduced assumptions and limitations that make the solution possible.

Spatial correlation
A Gaussian process prior is used to deal with the sparsity of the ISS data. The assumption is that the clonal composition changes smoothly across space. The biological reasoning behind it is based on the fact that breast cancer tissue exhibits some structural integrity, so when cells divide they remain close to each other. However, it's important to note that the model doesn't take into account the existence of detached compartments such as ducts. As this assumption breaks on the interface between compartmentalised regions, the inferred clonal composition on the borders may appear slightly more intermixed than it is in reality.

Copy number variation and expression
The BaSISS protocol detects RNA, however cancer clones are determined by their genotypes. To be able to conduct inference in relation to genetic clones, a relation between genotype and transcription needs to be stated. We assume a positive linear relationship between copy number of an allelic variant and its expression (and detection), with a possibility of a clone specific deviation. Previous studies suggest that the correlation between CNVs and expression seems to hold for many genes 22,23 .

Probe affinities
A further technical reason that weakens the link between genotype and raw BaSISS abundances stems from the fact that ISS uses molecular probes which bind to cDNA with variable affinity. This can lead to some ambiguity as probe affinities need to be estimated from the observed data, which requires additional information. This is especially the case if there are few clone defining variants available (otherwise poorly and strongly binding probes may average out). To assist estimation of the allele and clone specific detection rate, it is instrumental to use prior information on the clonal composition, at least at the level of the whole tissue section. Given the first step in the BaSISS workflow is to perform subclone discovery by WGS, prior knowledge of clone composition exists in most cases, even considering the possibility of sampling bias.

Homogeneity of clonal expression
Another assumption is that alleles belonging to a particular clone are expressed homogeneously across the slide. Violations of this expectation could be partly compensated by: 1. Multiple alleles defining clonal branches. If only a few alleles are highly variable, the rest could compensate for that. Overdispersion parameters would automatically decrease the influence of alleles with higher heterogeneity than expected.
2. Inhomogeneous background. Part of the heterogeneity could be captured as the background 'residual' clone.

Phylogenetic tree
The clonal genotype matrix / phylogeny is a discrete model parameter, such that the model doesn't have a way to adjust it directly. While it is theoretically possible to select an optimal tree solution based on BaSISS data, we observed that in practice different tree solutions provide very similar ELBO values, likely due to the aforementioned unknown parameters which lead to a practical unidentifiability.
Therefore, attention should be paid during WGS mutation clustering and copy number estimation, as errors in these quantities may lead to unrealistic clonal abundances and growth patterns. However, assessment of the model residuals could indicate problematic alleles and indicate if the genotype matrix is feasible.

Confidence intervals of the inference
Mean-field variational approximation allows relatively fast inference but has several limitations on its own. Variational approximation of the posterior is estimated using the reverse-KL divergence ( || ). If the true posterior is multimodal, the variational approximation is known to fit only one of the modes of , yielding only one of the potential solutions (the likely combination of the parameter values) and missing other potential solutions. In addition, mean-field approximation has a tendency to underestimate the parameter variance 24 .

Clone neighbourhood analysis
Knowing the spatial distribution of clonal densities, one may characterise the clones phenotypically based on additional spatially matched data. These include: • histopathological phenotype, • ISS expression signals assigned to nuclei, • nuclei cell types and • IHC staining.
These modalities of spatially resolved information are located on several consecutive slides of tissue which could be distorted from one layer to another. To integrate information among all the modalities, we annotated structurally similar areas of breast tissue on each slide which should be in physical proximity on the z-stack (see Histopathological assessment).
For differential phenotype analysis, highly clonal areas were selected based on the applicable threshold, CCFclone > 0.7 for P1, and CCFclone > 0.15 and 0.05 for P2-TN1 and P2-LN1 respectively (P2). The lower threshold for the P2-TN1 and P2-LN1 is due to a presence of infiltrations with the cells of a nonepithelial origin.

Cell type specific expression
To measure cell type specific differential expression associated with clonal regions, we record ISS signals which belong to the nuclei classified as the cell type of interest. In this case, the distribution of observed ISS expression signals "J of genes across regions is modelled as a result of a clone specific expression rate (J , However, the sum-to-one constraint of the cell type frequency, requires a softmax (logit) link function and Multinomial likelihood for the cell types: ", = (, + ", , ).

IHC data
For each region, the number of IHC-positive and IHC-negative nuclei was counted using QuPath 1 .
Then, IHC-positive and negative nuclei were treated as 2 distinct cell types, so the Cell type composition model would be applicable.

Inference and statistical testing
The models were implemented in numpyro 25 and HMC was used for the inference of the parameters.
Differentially expressed genes were determined through pairwise comparisons. The difference between clone specific expression rates (J was calculated for each clone pair. As HMC produces empirical posterior distributions with finite samples the distinction of small differences is challenging and probabilities/quantiles are limited to 1/n where n is the HMC sample size and many. genome library was also generated from whole blood derived DNA and sequenced in the same way. As described above for bulk WGS analysis, point mutations were called using CaVEMan and copy number using ASCAT and Battenberg algorithms and mutation clusters identified using DPclust v2.2.8 (Supplementary Data Table 5).

Mutation timing estimates
We report that there is evidence of early divergence of the subclones detected in P1 (prior to 50% of in m time). Evolutionary timing estimates can be sought for P1 but are not attempted for P2 due to relatively low purity. Based on the assumption of a relatively constant mutation rate within a given

Supplementary Discussion Notes
A key observation from these data is the variable relationship between genotype and phenotype. We examined three primary breast tumours with intermixed DCIS and invasive disease, and in each case, a clone was identified that existed simultaneously in the two distinct histological states. In contrast, albeit in a single examined case, we observe remarkable consistency between genetic lineage and DCIS grade spanning centimetres of tissue. While occasional examples of 'upgrading' occur, which could be accounted for by unsampled genetic evolution, we find no evidence of transition to a lower grade. The lack of plasticity within the grade phenotype space is aligned with studies that report that DCIS and invasive cancers are usually of the same grade 28,29 . Our findings could reflect stable differentiation states acquired by divergent cells, prior to DCIS onset 30 .
Another key message, that echoes those of a recent spatial study in melanoma 31 , is that measuring entire tissue sections rather than small fields of view is advantageous for understanding biology.
Firstly, this approach reveals that patterns of cancer growth and organisation can vary across different scales. For example, in P1-D1, P1-D2 and P1-D3, at macroscopic scale we see DCIS clones co-exist in multiple lobules, across an entire breast lobe, while at the microscopic level they are mainly segregated, distending individual acini or forming abutting populations. There are various possible explanations for the observed DCIS growth patterns such as clonal co-operation or genetic drift 32,33 .
A number of recent mathematical modelling studies have emphasised the importance of tissue architecture in shaping evolution and in the reported case it is certainly plausible that altered rates of fixation and absorption arising from the physical bottlenecks in the different sized ducts and ductules could account for the observed genetic patterns [34][35][36] . Secondly, within a lymph node the panoramic view allows one to appreciate that subclones repeatedly foster the same ecosystems -one that is immune depleted within the sinuses while the other clustered around B-cell aggregates consistent with a clone specific interactions with the adaptive immune response 37 . This way of examining tissues naturally leads to a multitude of questions that each warrant follow up in its own extended cohort. The tools we have developed and approaches demonstrated in this study will provide others with a framework for achieving this.
Particular advantages of the technology are that it is capable of interrogating very large tissue sections and it is comparably cheap, unlike solely relying on sequencing based methods 38 . However, of course the de novo mutation detection step must also be considered in terms of the complexity and cost of the entire workflow. The granularity of clones mapped by BaSISS depends on this critical phase and can be adapted according to the precise research question. Here we used multi-region WGS (~40 fold coverage per sample) to identify and map relatively broad subclone populations but the spatial genomics approach could be applied to single genomes, high depth targeted capture approaches or more detailed phylogenies obtained through LCM or even single cell approaches. In theory, the approach also holds the potential to create three dimensional genomic tomographs by aligning consecutive tissue sections and this will be particularly important for examining branched anatomical structures such as breast, prostate and lung.