# Transcriptional coregulation in cis around a contact insulation site revealed by single-molecule microscopy

Maciej A Kerlin, Ilham Aboulfath-Ladid, Julia Roensch, Chloe Jaubert, Aude Battistella, Kyra J E Borgman, Antoine Coulon

► To cite this version:

Maciej A Kerlin, Ilham Aboulfath-Ladid, Julia Roensch, Chloe Jaubert, Aude Battistella, et al.. Transcriptional coregulation in cis around a contact insulation site revealed by single-molecule microscopy. 2025. hal-05403707

HAL Id: hal-05403707

https://hal.science/hal-05403707v1

Preprint submitted on 8 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Transcriptional coregulation *in cis* around a contact insulation site revealed by single-molecule microscopy

Maciej A. Kerlin<sup>1,2,#</sup>, Ilham Aboulfath-Ladid<sup>1</sup>, Julia Roensch<sup>1,2</sup>,
Chloé Jaubert<sup>1,2,§</sup>, Aude Battistella<sup>2¶</sup>, Kyra J. E. Borgman<sup>1,2,‡,*</sup>, Antoine Coulon<sup>1,2,*</sup>

1. Institut Curie, PSL Research University, Sorbonne Université, CNRS UMR3664, Laboratoire Dynamique du Noyau, 75005 Paris, France.

2. Institut Curie, PSL Research University, Sorbonne Université, CNRS UMR168, Laboratoire Physique des Cellules et Cancer, 75005 Paris, France.

# Present address: Center for Integrative Genomics (CIG), University of Lausanne (UNIL), Lausanne, Switzerland

§ Present address: L'Atelier des Jours à Venir, 64210 Bidart, France.

¶ Present address: Ecole Normale Supérieure, IBENS, CNRS UMR8197, Inserm U1024, 75005 Paris, France.

‡ Present address: Institut Curie, Single-Cell Initiative, Cytometry Platform, Curie CoreTech, 75005 Paris, France.

* Correspondance: kyra.borgman@curie.fr, antoine.coulon@curie.fr

**Article type:** Short report (~1500 words, 3 main figures)
Abstract

Chromosome conformation in mammals is closely related to gene regulation. Within topologically associating domains, where genomic contacts are enriched, genes tend to show correlated expression across tissues and conditions, suggesting domain-wide mechanisms coregulating multiple genes, such as enhancer sharing or local histone mark spreading. At the single-cell level, where transcription occurs in sporadic bursts, transcriptional coordination has been observed between proximal genes, but how the local folding of mammalian chromosomes influences gene coregulation *in cis* at individual alleles remains unclear. Using single-molecule microscopy, we imaged nascent transcription from three adjacent genes located around a strong contact insulation site at the *FOS* locus, during the estrogen response in human breast cancer cells. To interpret this data, we developed two new analysis approaches to dissect the sources of (co)variation in gene activities: one to separate allele-extrinsic, allele-intrinsic, and gene-autonomous components; and another one to quantify the contributions of burst co-occurrence and burst size correlations. We find that transcriptional variability is largely gene-autonomous, yet correlations between genes display distinct patterns and occur almost exclusively *in cis*. Correlations are stronger between proximal and less insulated genes. However, unexpectedly, substantial correlations also occur across the strong insulation site and, under certain conditions, two proximal genes on the same side can exhibit uncorrelated burst occurrences. By disentangling burst co-occurrence from burst size correlations, we reveal transcriptional patterns suggesting two distinct coregulatory mechanisms influenced by local chromosome folding.

Introduction

The mammalian genome is highly organized in both space and sequence. At the sub-megabase scale, chromosomes fold into domains of enriched contacts, called topologically associating domains¹,². These domains have been shown to constrain gene regulation by enhancers³,⁴ and have been proposed to function as insulated regulatory units⁵,⁶, though some evidence nuances this view⁷,⁸. Still, the tendency of genes within these domains to show correlated expression across tissues and conditions¹,⁹⁻¹¹—especially when functionally related¹²—suggests the existence of domain-wide coregulatory mechanisms, such as enhancer sharing or histone mark spreading.

Transcriptional correlations provide a powerful tool to study gene coregulation at different scales, from population-level analyses to single-cell studies¹³,¹⁴. At the single-cell level, where transcription occurs in sporadic bursts, proximity between genes has long been linked to transcriptional correlations¹⁵. More recently, correlation patterns between genes have been used to investigate local regulatory mechanisms in mammals¹⁶,¹⁷, flies¹⁸,¹⁹, and yeast²⁰,²¹. However, how the local folding of mammalian chromosomes shapes gene coregulation *in cis* at individual alleles remains unclear.

Here, we combine single-molecule imaging with new data analysis approaches to investigate transcriptional correlation patterns among three adjacent human genes flanking a strong contact insulation site. We show that transcriptional coordination within a contact domain occurs at the level of individual alleles. However, we also uncover unexpected features that challenge a strict role for genome folding in transcriptional coordination and reveal correlation patterns suggesting two distinct coregulatory mechanisms, whose origins we speculate on.
# Results and Discussion

## Population-level data identify local, contact-dependent correlations in the transcriptional response of genes to estrogen

We first used public GRO-seq data²² to quantify, genome-wide, the local correlations between the (population-averaged) transcriptional response of genes to estradiol (E2) in MCF7 cells. We found that genes within a few hundred kilobases tend to have similar temporal profiles of transcriptional response (Figure 1A), unless strongly insulated from one another as measured by Hi-C²³ (Figure 1B). This suggests local coregulatory mechanisms linked to 3D organization (itself shaped by genomic distance and chromosome folding mechanisms), consistent with our current understanding of enhancer-gene control³. Yet, it remains unclear whether such coregulation between estrogen-responsive genes occurs locally *in cis* and how it manifests through transcriptional bursting.

## Single-molecule imaging of the joint activities of three genes at individual alleles with combined DNA-RNA FISH

To address this, we chose three adjacent genes: the proto-oncogene *FOS*, separated from *JDP2* and *BATF* by a strong contact insulation site (Figures 2A, S1A). All three genes are upregulated following E2 exposure, while showing basal transcriptional activity in the absence of E2 (Figures 1A, S1B). The genomic proximity of these genes (< 250 kb), their strong functional relationships (AP-1 family transcription factors²⁴), and the abundance of estrogen receptor α (ERα) binding sites in the vicinity (Figure S1B) hints towards the existence of local mechanisms that coregulate them.

We combined single-molecule RNA FISH with oligonucleotide-based DNA FISH to co-image the transcriptional activities of *FOS*, *JDP2*, and *BATF*, at each allele of the locus in MCF7 cells (Figure 2B). With an oligonucleotide library²⁵ spanning the whole locus, we used the DNA FISH channel to identify individual alleles (which copy number substantially varied between nuclei; Figure S1D-E). We observed bright RNA FISH spots at the identified DNA loci, corresponding to transcription sites, from which we estimated the number of nascent RNAs using the brightness of single RNAs as a reference (Figures S1F-G). With this approach, we calculated the distributions of nascent RNA counts and the fractions of active alleles for each gene following treatment with E2 or vehicle control (DMSO) for 40 min (Figure 2C), while maintaining cells in G1 using palbociclib (Figure S1C). We confirmed that all three genes go from basal to upregulated transcription upon E2 exposure (Figure 2C) and undergo transcriptional bursting (Figure S1H). The resulting joint distributions of nascent RNAs at individual alleles for all three gene pairs (Figure 2D) are the basis for analyzing how transcription is locally coordinated across the locus.

## Two mathematical approaches to separate sources of (co)variations in gene activity

Transcriptional correlations between genes were broken down into multiple additive components using two strategies (Figures 3A-B and S2).

First, we developed an approach to separate out global effects due to cell-to-cell variability (e.g. cell cycle, cell-wise differences in E2 responsiveness...). Our method extends the classical *intrinsic/extrinsic noise* framework²⁶, generalizing it to variable allele copy numbers and multiple genes per allele, yielding *allele-intrinsic* and *allele-extrinsic* pairwise correlations between genes (Figure 3A). It relies on computing covariance matrices from both allele-level measurements and their cell-wise sums, then using arithmetic to extract covariances across all possible cell-wise *trans* permutations and isolate *cis*-level effects—i.e. what would be observed
in the absence of cell-to-cell variability (see Methods; Figure S2A-B). While applied here between *single genes at individual alleles* within *single cells* of a *population* (Figure 3C), our generic framework can apply recursively —e.g. from single-cell and/or spatial-omics data (Figure S2C-D)— opening avenues for multiscale analyses of heterogeneities and correlation patterns¹⁴.

Second, to further dissect transcriptional correlations, we reasoned that, just as bursting is characterized by burst frequency (how often bursts occur) and burst size (number of RNAs per burst), correlations between genes may arise from burst co-occurrence (e.g. partially synchronized burst events, or simply covarying burst frequencies) and/or burst size correlations (e.g. correlated burst magnitudes, even if occurrence may be uncorrelated). To separate these effects, we combine the *law of total covariance* with a perturbation approach to break down the covariances into several components: one stemming from the contingency table of the burst/no-burst binarized data, while the rest characterizes correlations involving burst sizes (Figure S2E).

## Transcriptional fluctuations are predominantly gene-autonomous

Before interpreting the transcriptional correlation patterns between *FOS*, *JDP2*, and *BATF*, an important observation is that these correlations occur almost exclusively *in cis* and are rather moderate (Figure 2A). Indeed, the transcriptional variability of each gene considered individually is dominated by allele-intrinsic –or *cis*– effects (with only minor contributions from allele-extrinsic –or *trans*– effects; Figure S3) and even *in cis*, the covariances between genes remain moderate compared to the variances (Figure S3). Together, these findings indicate that, in our experiments (i.e. on populations of cells exposed to the same E2 treatment and arrested in G1), transcriptional fluctuations are largely gene-autonomous.

## The most proximal and least insulated gene pair exhibits the strongest *cis* correlations

Transcriptional correlations *in cis*, although moderate, offer a basis to study local mechanisms that coregulate transcription. Across all three gene pairs and in both uninduced (DMSO) and induced (E2) conditions, we find that the strongest *cis* correlation occurs in the most proximal and least insulated gene pair (*JDP2-BATF*; Figure 3A). This observation, made here for the first time at the single-allele level, aligns with our population-level analysis (Figure 1) and with previous reports at the population¹,⁹ and single-cell¹² levels, and support the possibility of domain-wide *cis*-coregulatory mechanisms relying on 3D chromosome organization (e.g. enhancer sharing, contact-dependent histone mark spreading).

In addition, while we cannot disentangle the contributions of 1D genomic distance and 3D contacts (since the least insulated genes are also the closest), correlation strength does not appear to be strictly determined by gene order nor by genomic distance, in line with earlier observations at larger scales¹⁷.

## Transcriptional correlation patterns nuance the role of local genomic contacts in gene coregulation

Although our results and many observations in the field¹,³,¹²,¹⁷ suggest 3D chromosome folding as a dominant factor in gene (co)regulation, two unexpected observations nuance this view.

First, we find that substantial transcriptional correlations can occur across a strong contact insulation site. Specifically, correlations between genes on either side of the contact insulation site (*FOS-JDP2* and *FOS-BATF*) are only ~3-fold weaker than for genes on the same side (*JDP2-BATF*) (Figure 3B), despite the ~7:1 intra-to-inter-domain contact ratio between the two domains delimited by the insulation site (Figure SA1). This indicates that contact insulation
does not make genes operate in isolation, but rather only moderately uncouple their activities, consistent with the idea that regulatory mechanisms may bypass contact insulation sites⁷.

Conversely, the two genes on the same side of the contact insulation site (JDP2-BATF) can, under certain conditions, show a surprising lack of transcriptional coordination. Indeed, separating the contributions of burst co-occurrence and burst size correlations reveals that the former is negligible in uninduced (DMSO) conditions for all three gene pairs (Figure 3B). Among other implications discussed below, this finding shows that JDP2 and BATF, two adjacent and proximal genes (< 100 kb apart) without a strong contact insulation site between them, can exhibit mostly uncorrelated burst occurrences.

## Burst occurrence and size statistics outline two distinct mechanisms of local gene coregulation

Transcriptional correlations *in cis* are present in uninduced conditions and increase upon E2 induction (Figure 3A). Our analysis disentangling burst co-occurrence from burst size correlations reveals a remarkably clear explanation for this behavior: burst co-occurrence arises almost exclusively upon E2 induction, while burst size correlations account for all basal pre-induction correlations and remain largely unchanged by induction (Figure 3B). This suggests two distinct coregulatory mechanisms, both of which appear shaped by local chromosome folding (strongest correlations in the most proximal and least insulated gene pair; Figure 3B).

The molecular origins of these two modes of coregulation remain speculative. Enhancers having been shown to modulate burst frequency³,¹⁴,¹⁶ and to coordinate burst timing in flies¹⁸,¹⁹, we propose that the burst co-occurrence we observe, which is almost entirely induction-dependent, may stem from the activity of the many ERα-bound enhancers at the locus (Figure S1A). This contrasts with an earlier report of disfavored co-bursting between genes sharing a single enhancer¹⁶, possibly reflecting here a collective behavior of many enhancers. In turn, burst size correlations reveal preexisting allele-to-allele variations that remain mostly unchanged during the 40 min of E2 induction. Burst sizes having been associated to histone marks, we propose that the burst size correlation we observed may originate from contact-dependent local spreading of histone marks²⁷,²⁸, hence constrained by genome folding and contact insulation patterns.

Further work will be needed to elucidate the mechanisms underlying such local transcriptional coordination between genes.
<figure><img src="image_2.png" /><figcaption>**Figure 1** – Population-level data suggest local transcriptional coordination in response to estrogen. (A) From GRO-seq measurements²² of nascent transcription over time in response to estradiol (E2) in MCF7 cells, the similarity of transcriptional time profiles between genes is assessed using the Pearson correlation coefficient (i.e. insensitive to differences in offset or scale of the time profiles) and plotted against their genomic distance as a two-dimensional histogram. (B) Same analysis using only gene pairs separated by a strong contact insulation site as measured by Hi-C²³.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2** – Measuring nascent transcription of three adjacent genes at individual alleles. (A) Hi-C contact map and insulation score (low values indicate strong contact insulation) at our genomic locus of interest. (B) MCF7 cell imaged in five colors after combined oligo-based DNA FISH and single-molecule RNA FISH (gray: DAPI). Five allele loci are visible in the DNA channel. In the RNA channels, bright spots colocalize with allele DNA loci, revealing ongoing bursts of nascent RNAs, while dim spots correspond to single RNAs (see Methods and Figure S1D-I). Shown are maximum-intensity projections of band-pass filtered images. Scale bar: 5 µm. (C) Quantification of nascent transcription of *FOS*, *JDP2*, and *BATF*, measured by FISH at individual alleles, following a 40-min treatment with 100 nM estradiol (E2) or vehicle control (DMSO), represented as number of nascent RNAs (top) and fraction of active alleles (bottom). (D) Same data, represented as joint distributions of number of nascent RNAs at individual alleles, for all three gene pairs, in uninduced (DMSO) and induced (E2) conditions.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 3** – Multilevel decomposition of transcriptional correlations. (A) For each gene pairs, the Pearson correlation coefficient between nascent RNAs at individual alleles is separated into two additive components, isolating the contribution of *trans* and *cis* effects. (B) The latter is further separated into two additive components, quantifying correlations due to burst co-occurrence and burst size correlations, respectively. Averages of 5 replicates (dots) are shown as bars, with standard errors as error bars. See Methods for statistical tests and *p*-values. (C) Our data and computational approaches yield a multiscale understanding of transcriptional correlation patterns, from population-level heterogeneity to allele-specific coupling and gene bursting dynamics.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure S1 – Local genomic context and transcription at the FOS-JDP2-BATF locus.** Genomic context of our locus of interest described by publicly available data from MCF7 cells: (A) Hi-C contact map<sup>23</sup> and an insulation score (low score indicates strong insulation), (B) ChIP-seq for CTCF and ERα on E2 treated cells<sup>29,30</sup>, and GRO-seq over time after E2 treatment<sup>22</sup>. (C) Histograms of EdU signal measured in MCF7 cells, after a 48 hrs treatment with different drug combinations, followed by 8 hrs of EdU incorporation. Percentages of arrested cells are indicated. (D) Distribution of fluorescence intensity of spots detected within nuclei of MCF7 cell, in the DNA FISH channel, using 12,005 fluorescent oligonucleotide probes labeling our locus of interest. Vertical line: intensity threshold above which a spot is considered to be an allele locus. (E) Histogram of the number of allele loci per nucleus. (F) Distributions of pairwise distances (left) between allele loci or (right) between allele loci and RNA FISH spots in each of the three RNA FISH channels (FOS, JDP2, and BATF). Vertical line: 650 nm. (G) Distributions of fluorescence intensity of FOS RNA spots, separating spots within 650 nm of an allele locus (yellow) from the rest (gray). The latter being considered single RNAs, the peak of their distribution is used, here and thereafter, to normalize fluorescence intensities. Vertical line: intensity threshold for an RNA spot to be unambiguously brighter than a single RNA, used to classify alleles a transcriptionally active. (H) Distributions of the number of nascent RNAs at individual alleles for all three genes, in uninduced and induced conditions. Black curve: expected distributions if the genes were not undergoing transcriptional bursting, i.e. producing individual RNAs as random and uncorrelated events.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure S2 – Mathematical approaches to decompose transcriptional correlations.** (A) Considering measurements of nascent RNAs on two genes at individual alleles in single cells, the covariance $\sigma_{alleles}^2$ of such measurements and the covariance $\sigma_{cells}^2$ of their sums or averages within cells can be combined to deduce the covariance $\sigma_{trans}^2$ across all the possible *trans* permutations. Intuitively, $\sigma_{cells}^2$ contains the products of the cell-wise sums (i.e. $X_1Y_1 + X_2Y_2 + \cdots$) which, when expanded in terms of single-allele measurement (e.g. $(x_1 + x_2 + x_3)(y_1 + y_2 + y_3) + (x_4 + x_5)(y_4 + y_5) + \cdots$), contains all intra-cellular pairs (i.e. $x_1y_1 + x_1y_2 + x_1y_3 + x_2y_1 + x_2y_2 + x_2y_3 + x_3y_1 + x_3y_2 + x_3y_3 + x_4y_4 + x_4y_5 + x_5y_4 + x_5y_5 + \cdots$), out of which the *cis* terms can be subtracted out using $\sigma_{alleles}^2$ (i.e. containing only the *cis* terms, $x_1y_1 + x_2y_2 + x_3y_3 + x_4y_4 + x_5y_5 + \cdots$). See Methods. (B) The *trans* covariance $\sigma_{trans}^2$ captures extrinsic factors (cell-to-cell differences in cell cycle, availability of machinery...) and equals the allele-level covariance $\sigma_{alleles}^2$ under the null hypothesis that genes X and Y at the same alleles do not correlate more or less than any other (non-*cis*) pair within cells. The residual $\sigma_{cis}^2 = \sigma_{alleles}^2 - \sigma_{trans}^2$ hence captures intrinsic variability, i.e. effects that are specifically occurring in *cis*. Illustrated are the two extreme scenarios where the variability is dominated by extrinsic factors (left, i.e. $\sigma_{alleles}^2 \approx \sigma_{trans}^2$) or by intrinsic factors (right, i.e. $\sigma_{alleles}^2 \approx \sigma_{cis}^2$). (C) Our general approach can be used recursively on data sets with several hierarchical levels: the cell-wise aggregated data (sums or averages of allele-level data) can itself be aggregate further at the cell-type level, then at the tissue levels... etc, each time decomposing further the origin of the variability. (D) Potential application in a broader context includes single-cell and/or spatial multi-omics data, which may include time (development, response to stimuli...) and experimental conditions, as well as other ‘omics’ measurements than gene transcription (e.g. scATAC, scCUT&Tag, ...). (E) The joint distribution of nascent RNA counts from two genes can be separated into four quadrants, indicating for each allele whether either gene is bursting or not. The total covariance of RNA counts can be written as the sum of the covariance within quadrants and the covariance of the average RNA counts within quadrants, which can themselves be expressed as a rescaled version of the contingency table and deviation terms. Separating out the contribution of the contingency table from the covariance of RNA counts characterizes burst co-occurrence, while all the other terms involve burst sizes.</figcaption></figure>
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Figure S3 – Covariance matrices.** (A) The covariances calculated across all the alleles measured from a cell population can be separated into two components, respectively capturing *trans* effects (i.e. cell-to-cell heterogeneity) and *cis* effects (i.e. (co)variations that are not explained by cell-wise effects). See Methods and Figure S2A-B. Each component corresponds to a covariance matrix, the values of which are shown in panel (B) with the same color code. (B) *Trans* and *cis* components of the variances and covariances of nascent RNA counts, in uninduced (DMSO) and induced (E2) conditions. Each bar is the average of n = 5 replicates, shown as points. Numbers below the bars indicate the contribution of the *cis* and *trans* to the total (co)variance (e.g. $\sigma_{\text{trans}}^{(\text{FOS,FOS})^2}/\sigma_{\text{alleles}}^{(\text{FOS,FOS})^2} = \sigma_{\text{trans}}^{(\text{FOS,FOS})^2}/(\sigma_{\text{trans}}^{(\text{FOS,FOS})^2} + \sigma_{\text{cis}}^{(\text{FOS,FOS})^2}) = 4.0\%$ in DMSO conditions).</figcaption></figure>
# Material and Methods

## Cell lines and cell culture

Human breast cancer cells MCF7 (ATCC) were grown at 37°C with 5% CO₂ in DMEM medium (Gibco, ref. 41966052) supplemented with 10% FBS (Gibco, ref. 10500064) and 5% Penicillin/Streptomycin (thereafter referred to as *regular culture medium*) and split twice a week.

For FISH sample preparation, 200k cells were plated in 12-well plate on 18 mm glass coverslips (Karl Hecht, ref. 11817093) coated with 20 µg of human fibronectin (Merck, ref. FC010) in 1 mL of *regular culture medium* and allowed to attach for 48 hours. To eliminate estrogen activity, we prepared a culture medium aiming to minimize amounts of molecules with estrogen-like activity (i.e. phenol red and estrogens from the FBS). This medium, thereafter referred to as *estrogen-free medium*, is made of phenol red-free DMEM (Gibco, ref. 21063029) supplemented with 10% charcoal-stripped FBS (PanBiotech, ref. P30-2302). After the cells have attached to the coverslips in *regular culture medium*, cells were (step 1) washed twice in 1X PBS (Gibco, ref. 14190250) and once in phenol red-free DMEM; then (step 2) cultured in *estrogen-free medium* for two hours, with the medium replaced once during this period; then (step 3) placed in 1mL of *estrogen-free medium* containing 1 µM palbociclib (Sigma, ref. PZ0383) and cultured for 48 hrs. Finally (step 4) 1µL of DMSO containing 17-β-estradiol (E2, Sigma-Aldrich E2758) or not (vehicle control, labeled DMSO) was added 40 min before fixation. The final concentration of E2 in the medium was 100 nM.

Estrogen-deprived MCF7 cells are known to accumulate at the G1/S transition and to enter S-phase upon E2 induction³¹,³². To avoid confounding effect due to cell cycle progression in our study, we used palbociclib, a CDK4/6 inhibitor that blocks cells at the G1/S transition³³, the efficacy of which we tested in our MCF7 cells by quantifying DNA synthesis by EdU incorporation³⁴ over a period of 8 hrs, following 48 hrs in *estrogen-free medium* as in step 3 above but with or without 1 µM palbociclib and supplemented or not with 100 nM E2 (Figure S1C). This data confirms the effect of estrogen deprivation on cell cycle progression (marked increase of EdU-negative cells in estrogen deprived conditions) and the efficacy of palbociclib whether cells are estrogen deprived or not (0% of cells with EdU signal), consistent with previous reports³⁵. Palbociclib was used as described above (steps 1–4) for all the data generated in this study, in which we observe induction of the studied estrogen-responsive genes (Figure 2C).

The experiments were performed on 5 different single-cell-derived populations of MCF7 cells and treated as replicates. To obtain the single-cell-derived populations, we separated individual cells from the population using flow cytometry and plated them individually in 96-well plates. Cells were then cultured in *regular culture medium* under standard conditions and upon reaching confluency replated to progressively larger culture dishes, until 15 cm plate. Cells were then collected and frozen at -80°C and transferred to liquid nitrogen for the long-term storage.

## RNA FISH probes

Custom RNA FISH probe sets of 20-nucleotide-long oligonucleotides were designed against *FOS*, *JDP2* and *BATF*, using the online Stellaris® Probe designer software (LGC, Biosearch Technologies) and ordered as DNA oligonucleotides conjugated respectively with CAL Fluor® Red 610, Quasar®670, and Quasar®570 fluorescent dyes, from Biosearch Technologies.
**FOS (mRNA):**
cagccactgcttttataaca agatgCGgttggagtacgag cagtcttggcttcctagatg aacatcatCGtgtCGgttag tcgtagtctgCgttgaagc
gagtgtgtagtaagagaggct catgctggagaagaggctg ggaatgaagttggcactgga ctggTCgagatggCagtga cgaagaGGgtgagggctct tgctgatCGgttcacaacg
ttctcatcttcctaggtgtc aatctCGgttCTgcaagaCag gcatcatCaggatcttCgga cagacatctcttCTggaag ccagtcagatCagaGGgagcag
ctgatgctcttcgacaggttc caggaagtcatcaagaGGgct agataggtCCatgtCTgcca agtCTgCTgcatagaaggaC agtgCagCGagTCaggaagg tacaGgaaACCctCTaggaa
gtttCacGCacagataaGgt gatgCTTtCAgTCCTgag tgaaggttCaggtCCACAC gtttgtCataCTtCCCTgaaCTcagctCTCTgaagtgCact ggtCAcaACatgCTactaCT
ccccatagattagtttaatgc gcaccaggtaattCCAataa caacttggCCaggaCacagc gtcgtatcttttCTtagtat tCAcaacaatgcatgatCagt aatgTCagaACaCttCagC
**tccacatgtCAaaagacCTC gtcGCcatTCaacttaaATgc**

**JDP2 (introns 1+2):**
ctgatgCTCCTCCAattaag agggCTtaaaggacGTTgg gataCGCGaATgCagttCT gataCTCggttaAAACTCtCC gtcTgCaaAagCTtagCagTG caagttgagCTggagCTgaa ttttCAACAACACCTCCACA
tcCTCAaAGGCataggagac gactgCCAagaCaaGCAca aCAagCAGCAGgtCTCACAC ccagCAGCCAGataaTgtaa ccACACagCCAGataaTgtaa tCAgAAATggagCTgCTCAAC cCAAGagCTgagTgaGGGC
atccagaAGGCGagaAGCAa gagaatCACACCAGgatGAC tgGagaATCACACCACaAGC gTCCACagagaAGgtCAGaaA agaATCGgagTgagAgTG cGGTCTaAggttCAAggttt
cctgTGAAGCCAACTTGAag ccaggaataAgttggaAGCC ataGGCTggaATCCAAGaAGC ttGGGaaAGAgGGCTattgC cCTggttGtCCatACggaAA TACccACTtCAAggaAACTg
accaAGGAACCAgtagCAacg aCAgCTgatATCaccagC agaACACCTCAGgtCTCt tttCACCTCAAgTTatCA gCTtagATgCTCgttATCag gACCTattCCTgtagCAaaA
tagtggCagaTAgagagCTC aACTCTgagCTggagCTgga atGTGGgttGtCAgTgtCT gaaAGttCCCTgaaATaACTCC taACTgTgCTCCTACCTaTG gAgTTTTattTTggTgtCCT
aAggtGaaATTTGGTtCCTC aATCtCTtCCatCTaCattt atttGtaACCTgtttggc cTAgATCattCCAGagAggk gTTAATgTACTCagATGGCC tCCatCGGACaaATCTCTCt
aAAATCCACatGTGGTgCT gACCTgaggAggaATTCCTT gATgagAgCTgCAGgtgCAA tCTtCCACgAAAgAAATCCT gATagttgCTCtaATaTCCC
caccAgTGCCCGagAAAGAA gACagGCTTCtCagAAACg tCAgCTtTGgCCTgCTgCT gATgCAACACTTCgCGCC tACCACCTCAACACAGCTTgC cCTTgCAtgTGAagACagAA
gCagAGGGCCAAAGAcaag gCTtCAAgTTACagACAGTgg cagCTactggCTgCTaATgAg gAGggAggACgATCAaaAT gacaACACCTgCAGTgCAGA
cagCTCAgTgCGGTACag gCCTCAAAACgCTTgTta gagggAggACgATCAaaAT gacaACACCTgCAGTgCAGA tCACACACCTgCAGTgCAGA aAgCATTGGAggTCGGTTC
cgcCTTTaggtTTtagTAC tCagCTgCAAgTTCTTCC cagCTTtCAGACCATagatt agATCTtCAGACCATagattC aggCAGtttGCCatATgCAA tGCAcaATATGAAAATCCgCT aAgTgAgACCCggCTgTACa
gaAgGGGTatAggagCgAT gTCAACAGCTCTGGTTAG tCAgTGAACACTTCgCGCC tAGgAgACgAgTgAggCTgCA aAgAACACCCTgCAGAg gCCTCTgCAAgGCAaaACAC
tgaATTgCTCCatgCGTTg tAAATCtCCCTaTTtCTTT tCAcACAAACACCTgTCgtgA gACatCAGCACAaAGCACTT aACACtgATgtgTTgACCATg agCAGttGGCACCAAAACag

**BATF (intron 2):**
tCAaAGCTatggTCAAgGG catACACATCCCACACaAG cCTtCTtaAAgCACTCtCAc tagaaAgGCCatttCTCTgct tCGtCACTCGCgagACaATGG
aAgGGATgCgACTatGTgT aAgagGAatCTggTACagCG attgCATctttCAtagACCC tGGTgATCATCCCAgataAg tGGgACTagAggCAACTCA gACCACTCAagACagAggTAa
ggaACTCgACttCagACgCT CCTAgAACACCATCTAggG tATgAGCAATGGGGAgCAag tGTAAGAgAggTAAgTGGT aACCAAGACCACTAgGGGGG
gaATTtCTTgagCAgAgGG ttGCTAggCTattAgAggTAa CTTattgATTTTTGGTgCCA atCATgCTaATCTCCAgAT gATgCATTatttCCCCCTAAA cCTCTgtttCCTtgTAAgCT
aatCagAgCTgttGGTCag gTGTgCAATgagGTgAagC gACTgttGAgCACagTTAg tAGGGCACCAAAACTtGAC gCTTAAggtGGTAAACCAgA gTTGTCTGCCCTAAAaAAT
agCTGCCAACTATACCTAC tgTgCACCCTCATTaATAC attCagTGCTgAAAgCCTTG gATgAGCTgACTCTTTGGC ccTATATgCACatGACTTG gACTATatGACTTGGGGCT
tgtgtCacGATgtGGGGAg agaAGCAGTGGAggAgAA taADAgACgATgCTgTGGTC gACCTACTCATTGGCGAAA ataAgttGCattGGCAGCag cCTtCActTGaAGCCAgAAa
ttgAACACgTTgCTgCCAa agCTaATTTgATgTAggGCA atAggtgattgTCatCTCGG tttCACCTgAATCTCAGaca gCAagCATTGGaaAGAggTCA

## Design and production of iFISH/Oligopaints library

A library of 12,005 oligonucleotides homogenously covering a 1.17 Mb genomic region encompassing our locus of interest (chr14:75,302,947-76,473,730, hg19 coordinates) was obtained from the 'Full 40mers' database of the iFISH platform²⁵. To avoid the binding of DNA FISH probes to RNAs from our 3 genes of interest, all probes are (+)-strand sequences, hence targeting the (-)-strand, corresponding to the template (non-coding) strand of *FOS*, *JDP2* and *BATF*. The 40nt-long oligonucleotide sequences were extended according to Bintu *et al.*³⁶ by adding (i) on either sides, 20nt-long sequences for secondary targeting (not used in this publication), then (ii) in 5', a sequence (CATCAACGCCACGATCAGCT) to be targeted by a primer for the retro-transcription (RT) reaction, and finally (iii) on either sides, two sequences (CGGCTCGCAGCGTGTAAACG in 5' and CGGTGCATTTGCGGGAAGAC in 3') to be targeted by two PCR primers for amplification of the whole library. The resulting library of 140nt-long sequences was ordered from GenScript as a Classic Oligo Pool-12K oligonucleotide pool.

From the ordered oligonucleotide pool, DNA FISH probes were then produced as described in Beliveau *et al.*³⁷ with small modifications in order to make the whole library of primary DNA FISH probes directly fluorescent. Briefly, the ordered oligonucleotide pool was first PCR amplified with Phusion Hot-start Master Mix (New England BioLabs), using forward primer CGGCTCGCAGCGTGTAAACG and reverse primer *TAATACTGACTCACTATAGGGGTCTTCCCGCAAATGCACCG*, where the italicized part is a T7 promoter sequence. The PCR product was purified using DNA Clean & Concentrator-5 (DCC-25) kit (Zymo research). The resulting double-stranded DNA template was then *in vitro* transcribed using HiScribe T7 Quick High Yield RNA Synthesis Kit (New England BioLabs), producing a reverse-complement single-stranded RNA copy of the original library. This was then retro-transcribed with Maxima H Minus RT Transcriptase (ThermoFisher Scientific) using as RT primer a DNA oligonucleotide ATCGACCCGGCATCAACGCCACGATCAGCT ordered from IDT DNA coupled in 5' to a 6-FAM (fluorescein) fluorophore. RNAs were finally degraded from the resulting RNA-DNA duplex using a NaOH-EDTA solution for 10' at 95°C and the resulting fluorescently-labeled single-stranded DNA oligonucleotides were purified and concentrated using DNA Clean & Concentrator-25 (DCC-25) kit (Zymo research). All steps
involving RNA were performed under RNase-free working conditions and with the addition of RNasin Plus RNase Inhibitor (Promega).

## Combined RNA-DNA FISH

The cells were crosslinked for 10 minutes at room temperature (RT) using Fixation Solution (4% paraformaldehyde (EMS, ref. 15710)/1X PBS (Gibco, ref. 14190250)) and subsequently washed twice for 10 min with 1X PBS. The cells were then permeabilized for 20 min at RT with Permeabilization Buffer (0.5% Triton X-100 (EUROMEDEX, ref. 2000-C)/1X PBS), washed three times for 1 min with 1X PBS and then treated with 0.1N HCl for 5 min. Next, the coverslips were incubated in the dark for 2 days in prehybridization buffer (2X SSC (EUROMEDEX, ref. EU0300-A)/50% formamide (Ambion, ref. AM9342)) at RT. After prehybridization, each coverslip was gently placed on microscope slides with a 30 µL drop of Primary Hybridization Buffer containing amplified primary iFISH/Oligopaints probes (2X SSC/50% formamide/10% dextran sulfate (Sigma, ref. D8906-50G)/0.2 µM of primary probes). The coverslips were then sealed with fixogum (MP Biomedicals, ref. 11430456) and left for 30 min in the dark at RT for the fixogum to solidify. Afterwards, the whole slides were placed on a heating block at 80°C for 3 minutes to denature DNA. Following denaturation, the slides were placed in a humid chamber and incubated overnight at 37°C in the dark. After hybridization of the primary probes, fixogum was gently removed and the coverslips were washed three times for 10 min with 2X SSCT (2xSSC/0.2% Tween) at 37°C, then twice for 7 min in 0.2X SSCT (0.2X SSC/0.2% Tween) at 65°C and finally briefly rinsed once with 4X SSCT (4X SSC/0.2% Tween) at RT and then twice with 2X SSC at RT. Subsequently, the coverslips were placed on microscope slides with a 25 µL drops of the Secondary Hybridization Buffer containing secondary smRNA FISH probes, as well as fluorescent readout probes aiming to visualize a small 20kb DNA region at the contact insulation site, although we did not use it in this publication (2X SSC/30% formamide/10% dextran sulfate/0.1 µM of each of the smRNA FISH probes/1 µM of secondary readout probes) in the same manner as for the primary hybridization. The coverslips were then sealed with fixogum, placed in a humid chamber and incubated overnight at 30°C in the dark. After secondary hybridization, the fixogum was gently removed and the coverslips were washed three times for 10 min in Wash Buffer (2X SSC/30% formamide) at 30°C, rinsed twice with 2xSSC at RT and then stained with 5 ng/mL DAPI in 2X SSC for 10 min at 30°C. The coverslips were then rinsed twice with 2X SSC, mounted on a microscope slide with Prolong Gold (Fisher Scientific, ref. P36934) and left overnight in the dark at RT to dry. After preparation, the slides were stored in the dark at -20°C.

## Microscopy and image quantifications

The FISH samples were imaged on an epifluorescence microscope made of a Nikon Ti-2 stand with a 60x 1.4NA CFI Plan Apochromat Lambda D oil immersion objective, a SpectraX (Lumencor) light source, a H117 XY stage (Prior), a Nano-ZL100 stage piezo (Mad City Labs), Semrock filters and dichroics, and an ORCA-Flash4.0 V3 sCMOS camera (Hamamatsu), and controlled by the MicroManager software³⁸. Images were analyzed using custom-written analysis scripts (https://github.com/CoulonLab/FISHingRod). Nuclei were first segmented based on DAPI staining, using a threshold-based approach. Improperly segmented objects were excluded first based on their area, then through manual curation using only the DAPI channel. Fluorescent spots were detected in the DNA and RNA FISH channels. The 3D spot coordinates in each channel were corrected for chromatic shift using an affine transform (translation and rescaling) determined using images of broad-spectrum fluorescence beads (TetraSpeck, Invitrogen, ref. T14792) acquired on the same microscope.
Using, for the DNA FISH, a large library of 12,005 oligonucleotides that are directly fluorescent (no secondary oligonucleotide for labeling) permitted non-ambiguous identification of the alleles (fluorescence spot intensities of most alleles are well separated from background; Figure S1D). MCF7 cells having an instable karyotype, the copy number of alleles, even from a single-cell derived population, is variable (Figure S1E). We discarded nuclei with less than 2 alleles from our analysis, as well as nuclei with a high fluorescence background around RNA FISH spots. Plotting the 3D distances between the allele DNA loci and the detected RNA spots in each channel revealed a clear enrichment up to 650 nm (Figure S1F), suggesting that these correspond to nascent RNAs. Indeed, using this as a cutoff distance, RNA spots away from any allele loci showed a narrow distribution of fluorescence intensities, as expected for individual RNAs, while the fluorescence intensity distribution of RNA spots close to an allele locus showed a long tail towards larger values, as expected for genes with bursts of nascent RNAs (Figure S1G).

For each RNA channel, the single-RNA spot intensity was used to normalize fluorescence intensity values, yielding absolute numbers of RNAs. These normalization values showed little variation between experiments. When the spot intensity distribution lacked a clear single-RNA peak (e.g. due to high background), we used the value obtained in other experiments for the same channel. For each allele DNA locus, the amount of nascent RNAs was determined, for each channel, from the brightest RNA spot within 650 nm or assigned to 0 when no RNA spot was found within this distance (Figure 2C, top). Note that, as with most methods in the field, the resulting nascent RNA counts are underestimated, since RNAs become fully fluorescent only after the targeted region is fully transcribed and, for intronic probes, as long as it is unspliced. Furthermore, single nascent RNAs are indistinguishable from already-released single RNAs that happen to be within 650 nm of an allele DNA locus. Hence, we can only accurately determine the presence and number of nascent RNAs when the fluorescence signal exceeds a certain threshold (e.g. 2 single-RNA-equivalent fluorescence units for *FOS*; Figure S1G, dashed line). When binarizing our measurements to determine the fraction of active/bursting alleles (Figure 2C, bottom), we selected thresholds as low as possible while maintaining low false-positive calls: namely $c_{FOS} = 2$, $c_{JDP2} = 2.5$, and $c_{BATF} = 2$ in single-RNA-equivalent fluorescence units.
# Analysis of transcriptional correlations

## Population-level correlation patterns from public GRO-seq and Hi-C data

We used Global Run-On sequencing (GRO-seq) data set from Hah et al.²², mapping transcriptionally engaged RNA polymerases in MCF7 cells upon induction with 100 nM estradiol (E2). From the provided table (Gene Expression Omnibus GSE27463) quantifying GRO-seq reads in all RefSeq-annotated genes, we note $X_k$ and $Y_k$ the measurements obtained on two genes X and Y in experimental condition k (i.e. at 0', 10', 40' or 160' after induction with E2) and $X = \sum_k X_k/N$ and $Y = \sum_k Y_k/N$ their means, where $N = 4$ is the number of experimental conditions.

The covariance and Pearson correlation, across experimental conditions, between the population-level measurement of the transcriptional activities of genes X and Y are

$$ \sigma_{\text{pop.}}^{(X,Y)} = \frac{1}{N} \sum_k (X_k - \bar{X})(Y_k - \bar{Y}) \quad \text{and} \quad r_{\text{pop.}}^{(X,Y)} = \frac{\sigma_{\text{pop.}}^{(X,Y)} \sigma_{\text{pop.}}^{(Y,Y)}}{\sigma_{\text{pop.}}^{(X,X)} \sigma_{\text{pop.}}^{(Y,Y)}} $$

This correlation measures how similar the *shapes* of the time profiles of the transcriptional response of two genes are, up to arbitrary offset and rescaling of either time profile —that is, independently of basal transcription levels and of the magnitude of the E2-induced response. In other words, the correlation between ($X_{0'}$, $X_{10'}$, $X_{40'}$, $X_{160'}$) and ($Y_{0'}$, $Y_{10'}$, $Y_{40'}$, $Y_{160'}$) is the same as between ($aX_{0'} + b$, $aX_{10'} + b$, $aX_{40'} + b$, $aX_{160'} + b$) and ($Y_{0'}$, $Y_{10'}$, $Y_{40'}$, $Y_{160'}$) for any value of $a > 0$ and $b$.

Pearson correlation values close to 1 indicate similar transcriptional time profiles (e.g. the activity of two putative genes A and B both increasing after 10', plateauing until 40', then reverting at 160', even if A starts with a high basal level and shows small E2-induced variations while B does not). Values close to 0 indicate no similarity (e.g. gene A vs a gene C which activity keeps increasing). Values close to -1 indicate opposite transcriptional time profiles (e.g. gene A vs a gene D which would be downregulated after 10' and until 40', then reverting at 160').

Figure 1A presents the 2D histogram, for all gene pairs genome-wide within 1 Mb of each other, of the Pearson correlation between their transcriptional response to E2 and the genomic distance between their promoters. Figure 1B presents the same analysis but only for the top 20% most insulated gene pairs, as calculated from the Hi-C data of Rodriguez et al.²³ obtained in MCF7 cells (Gene Expression Omnibus GSE121443). For this, we used the cooler package³⁹ to generate 1 kb- and 10 kb-resolution matrices for each experiment, pooled them together, and calculated an insulation score at each genomic position by summing the contact frequency values observed between two 140kb-regions upstream and downstream of the position of interest, leaving a 10kb-gap between them. A low value of the insulation score indicates a strong contact insulation.
**Approach 1 – Separating cis (allele-intrinsic) and trans (allele-extrinsic) components in transcriptional correlations from RNA FISH data**

– Principle –

Consider we have nascent RNA measurements at a total of *n* alleles from a given experiment. Let (*x*<sub>*i*</sub>, *y*<sub>*i*</sub>) be the *n* measurements from genes X and Y at each individual allele *i*, and their means be *x* = Σ<sub>*i*</sub> *x*<sub>*i*</sub>/n and *y* = Σ<sub>*i*</sub> *y*<sub>*i*</sub>/n (unless otherwise noted, sums on *i* are taken over [1, *n*]).

The covariance across alleles of the allele-level measurements is

$$
\sigma_{\text{alleles}}^2 = \frac{1}{n} \sum_{i} (x_i - \bar{x})(y_i - \bar{y})
$$

Let us consider the *n* alleles to be partitioned among *N* cells. Let *n*<sub>*j*</sub> denote the number of alleles in cell *j*, so that *n* = Σ<sub>*j*</sub> *n*<sub>*j*</sub> (sums on *j* are taken over [1, *N*]). Let us note the cell-wise summed measurements as *X*<sub>*j*</sub> = Σ<sub>*i*∈Ω<sub>*j*</sub></sub> *x*<sub>*i*</sub> and *Y*<sub>*j*</sub> = Σ<sub>*i*∈Ω<sub>*j*</sub></sub> *y*<sub>*i*</sub> where *i* ∈ Ω<sub>*j*</sub> denotes the alleles *i* in cell *j*. Using the same mean subtraction as above, such that Σ<sub>*i*∈Ω<sub>*j*</sub></sub> (*x*<sub>*i*</sub> − *x*) = *X*<sub>*j*</sub> − *n*<sub>*j*</sub>*x* and Σ<sub>*i*∈Ω<sub>*j*</sub></sub> (*y*<sub>*i*</sub> − *y*) = *Y*<sub>*j*</sub> − *n*<sub>*j*</sub>*y*, the covariance of the cell-wise summed measurements is

$$
\sigma_{\text{cells}}^2 = \frac{1}{N} \sum_j (X_j - n_j x)(Y_j - n_j y)
$$

We note that $N\sigma_{\text{cells}}^2 = \sum_j \sum_{i,i' \in \Omega_j} (x_i - x)(y_{i'} - y)$ enumerates all the possible pairs within cells, including both cis and trans pairs. Hence, the covariance over all possible trans pairs, i.e.

all the possible intra-cellular permutations (explicitly $\sigma_{\text{trans}}^2 = \frac{\sum_j \sum_{i,i' \in \Omega_j; i \neq i'} (x_i - x)(y_{i'} - y)}{\sum_j n_j(n_j - 1)}$) can be simply calculated as

$$
\sigma_{\text{trans}}^2 = \frac{N \sigma_{\text{cells}}^2 - n \sigma_{\text{alleles}}^2}{\sum_j n_j (n_j - 1)}
$$

as illustrated in Figure S2A.

The part of $\sigma_{\text{alleles}}^2$ explained by cell-wise (allele-extrinsic) effects corresponds to $\sigma_{\text{trans}}^2$ and their difference, which we note $\sigma_{\text{cis}}^2$, isolates the part of the covariance originating from local (allele-intrinsic) effects in cis:

$$
\sigma_{\text{alleles}}^2 = \sigma_{\text{trans}}^2 + \sigma_{\text{cis}}^2
$$

If the covariance across alleles were fully explained by *trans* effects, then $\sigma_{\text{alleles}}^2 = \sigma_{\text{trans}}^2$ and $\sigma_{\text{cis}}^2 = 0$ (Figure S2B, left). On the opposite, if alleles behaved independently and there was no cell-to-cell heterogeneity, then $\sigma_{\text{trans}}^2 = 0$ and the measured covariance would fully come from *cis* effects $\sigma_{\text{alleles}}^2 = \sigma_{\text{cis}}^2$ (Figure S2B, right). In the general case, $\sigma_{\text{trans}}^2$ quantifies cell-to-cell heterogeneities or other sources of inter-allelic coupling and $\sigma_{\text{cis}}^2$ is what would be observed if there were no such inter-allelic effects.
## - Remarks -

1. In the case of two identical copies of a gene in each cell, our approach simplifies into the classical extrinsic and intrinsic noise definition²⁶, which decomposes the variance as $\sigma_{gene}^2 = \sigma_{extrinsic}^2 + \sigma_{intrinsic}^2$. Denoting $(x_j, x'_j)$ the two measurements in each cell j, we have $\sigma_{gene}^2 = \frac{\sum_j (x_i - x_j)^2}{N} = \frac{\sum_j (x'_i - x_j)^2}{N}$ and $\sigma_{extrinsic}^2 = \frac{\sum_j (x_i - x_j)(x'_i - x_j)}{N}$ with $\bar{x} = \frac{\sum_j x_i}{N} = \frac{\sum_j x'_i}{N}$.

2. We chose to use cell-wise summed measurements $(X_j, Y_j)$, but this approach can be adapted to use cell-wise averages $\bar{x}_j = \frac{x_i}{n_j}$ and $\bar{y}_j = \frac{y_i}{n_j}$. Namely, we write

$$\tilde{\sigma}_{\text{cells}}^2 = \frac{1}{\sum_j n_j^2} \sum_j n_j^2 (\bar{x}_j - x)(\bar{y}_j - y)$$

(notice the unusual weighting) which only differs from $\sigma_{cells}^2$ by a scaling factor corresponding to the average number of allele pairs per cell (i.e. $\tilde{\sigma}_{cells}^2 \frac{\sum_j n_j^2}{N} = \sigma_{cells}^2$). Then, we get $\sigma_{trans}^2 = (\sum_j n_j^2 \tilde{\sigma}_{cells}^2 - n)/\sum_j n_j(n_j - 1)$ and the rest is identical.

3. Although related, this approach differs from the classical way to separate within-group and between-group effects using the law of total covariance, which reads:

$$\sigma_{\text{alleles}}^2 = \frac{1}{n} \sum_i (x_i - x)(y_i - y) \\ = \frac{1}{n} \sum_j \sum_{i \in \Omega_j} (x_i - x)(y_i - y) \\ = \frac{1}{n} \sum_j \sum_{i \in \Omega_j} (x_i - \bar{x}_j)(y_i - \bar{y}_j) + \frac{1}{n} \sum_j \sum_{i \in \Omega_j} (\bar{x}_j - x)(\bar{y}_j - x)$$

$$\text{Within-group covariance: } \sigma_W^2 = \frac{1}{n} \sum_j \left[ n_j \frac{\sum_{i \in \Omega_j} (x_i - \bar{x}_j)(y_i - \bar{y}_j)}{n_j} \right] + \frac{1}{n} \sum_j n_j (\bar{x}_j - x)(\bar{y}_j - x) \quad \text{Between-group covariance: } \sigma_B^2 = \frac{1}{n} \sum_j n_j (\bar{x}_j - x)(\bar{y}_j - x) \quad \text{Covariance of cell-wise averages}$$

Notice how the weighting of the terms in the sum differs between $\sigma_B^2$ and $\tilde{\sigma}_{cells}^2$ above. Substituting $\bar{x}_j = \frac{\sum_{i \in \Omega_j} x_i}{n_j}$ and $\bar{y}_j = \frac{\sum_{i \in \Omega_j} y_i}{n_j}$ leads to $\sigma_B^2 = \frac{1}{n} \sum_j \frac{1}{n_j} \sum_{i, i' \in \Omega_j} (x_i - x)(y_i' - y)$ and $\tilde{\sigma}_{cells}^2 = \frac{1}{\sum_j n_j^2} \sum_j \sum_{i, i' \in \Omega_j} (x_i - x)(y_i' - y)$. Both expressions enumerate all the possible intra-cellular pairs, but the former gives them unequal weighting between cells, making it impossible to subtract the cis pairs using $\sigma_{alleles}^2$.

Two useful limit cases to consider are:

* If all cells have the same number of alleles, then $\sigma_B^2 = \tilde{\sigma}_{cells}^2$. In this case, we get

$$\sigma_{\text{cis}}^2 = \frac{n}{n-N} \sigma_W^2 \quad \text{and} \quad \sigma_{\text{trans}}^2 = \sigma_B^2 - \frac{N}{n-N} \sigma_W^2$$

* In addition, for large groups ($n/N \gg 1$, which does not hold here for alleles in cells, but may be valid in other situations; see Generalization below), we get

$$\sigma_{\text{cis}}^2 \approx \sigma_W^2 \quad \text{and} \quad \sigma_{\text{trans}}^2 \approx \sigma_B^2 \quad (\text{unless } \sigma_{\text{trans}}^2 \ll \sigma_{\text{cis}}^2)$$
## Application

Applying this approach to different gene pairs, e.g. X, Y, Z, yields covariance matrices:

$$ \frac{\begin{bmatrix} \sigma_{\text{alleles}}^{(X,X)} & \sigma_{\text{alleles}}^{(X,Y)} & \sigma_{\text{alleles}}^{(X,Z)} \\ \sigma_{\text{alleles}}^{(Y,X)} & \sigma_{\text{alleles}}^{(Y,Y)} & \sigma_{\text{alleles}}^{(Y,Z)} \\ \sigma_{\text{alleles}}^{(Z,X)} & \sigma_{\text{alleles}}^{(Z,Y)} & \sigma_{\text{alleles}}^{(Z,Z)} \end{bmatrix}}{\sigma_{\text{alleles}}^2} = \frac{\begin{bmatrix} \sigma_{\text{trans}}^{(X,X)} & \sigma_{\text{trans}}^{(X,Y)} & \sigma_{\text{trans}}^{(X,Z)} \\ \sigma_{\text{trans}}^{(Y,X)} & \sigma_{\text{trans}}^{(Y,Y)} & \sigma_{\text{trans}}^{(Y,Z)} \\ \sigma_{\text{trans}}^{(Z,X)} & \sigma_{\text{trans}}^{(Z,Y)} & \sigma_{\text{trans}}^{(Z,Z)} \end{bmatrix}}{\sigma_{\text{trans}}^2} + \frac{\begin{bmatrix} \sigma_{\text{cis}}^{(X,X)} & \sigma_{\text{cis}}^{(X,Y)} & \sigma_{\text{cis}}^{(X,Z)} \\ \sigma_{\text{cis}}^{(Y,X)} & \sigma_{\text{cis}}^{(Y,Y)} & \sigma_{\text{cis}}^{(Y,Z)} \\ \sigma_{\text{cis}}^{(Z,X)} & \sigma_{\text{cis}}^{(Z,Y)} & \sigma_{\text{cis}}^{(Z,Z)} \end{bmatrix}}{\sigma_{\text{cis}}^2} $$

Figure S3 shows the values of covariance matrices $\sigma_{\text{trans}}^2$ and $\sigma_{\text{cis}}^2$ built from the nascent RNA FISH data on *FOS*, *JDP2* and *BATF*, in both uninduced (DMSO) and induced (E2) experimental conditions. For each gene and each condition, the cis-to-trans ratio is displayed (e.g. in E2 condition, $\sigma_{\text{cis}}^{(JDP2,JDP2)^2}/\sigma_{\text{trans}}^{(JDP2,JDP2)^2} = 89.8$), as well as the percentage of its cis variance that is explained by its covariance with either of the other two genes (e.g. in E2 condition, $\sigma_{\text{cis}}^{(FOS,JDP2)^2}/\sigma_{\text{cis}}^{(JDP2,JDP2)^2} = 4.2\%$ and $\sigma_{\text{cis}}^{(JDP2,BATF)^2}/\sigma_{\text{cis}}^{(JDP2,JDP2)^2} = 15.8\%$).

The Pearson correlation across alleles, between any gene pair X and Y, can be decomposed as

$$ r_{\text{alleles}}^{(X,Y)} = \frac{\text{Pearson correlation coefficient}}{\frac{\sigma_{\text{alleles}}^{(X,Y)^2}}{\sigma_{\text{alleles}}^{(X,X)} \sigma_{\text{alleles}}^{(Y,Y)}}} = \frac{\text{Trans component:}}{\frac{r_{\text{trans}}^{(X,Y)}}{\sigma_{\text{alleles}}^{(X,Y)^2}}} + \frac{\text{Cis component:}}{\frac{r_{\text{cis}}^{(X,Y)}}{\sigma_{\text{alleles}}^{(X,Y)^2}}} + \frac{\sigma_{\text{trans}}^{(X,Y)^2}}{\sigma_{\text{alleles}}^{(X,X)} \sigma_{\text{alleles}}^{(Y,Y)}} + \frac{\sigma_{\text{cis}}^{(X,Y)^2}}{\sigma_{\text{alleles}}^{(X,Y)^2} \sigma_{\text{alleles}}^{(Y,Y)}} $$

which is what Figure 3A represents.

Note that the denominators are the same so that $r_{\text{trans}}^{(X,Y)}$ and $r_{\text{cis}}^{(X,Y)}$ represent two additive components of $r_{\text{alleles}}^{(X,Y)}$. Hence, $r_{\text{trans}}^{(X,Y)}$ and $r_{\text{cis}}^{(X,Y)}$ are not the Pearson correlation coefficients of the cis and trans components (i.e. calculated from the covariance matrices $\sigma_{\text{trans}}^2$ and $\sigma_{\text{cis}}^2$). However, since most of the variability occurs in cis in our case ($\sigma_{\text{trans}}^2 \ll \sigma_{\text{cis}}^2$, Figure S3B), the cis and total covariances are very similar ($\sigma_{\text{alleles}}^2 \approx \sigma_{\text{cis}}^2$), so that what is shown on Figure 3A (i.e. the cis component of the Pearson correlation, $r_{\text{cis}}^{(X,Y)}$ above) is very close to the Pearson correlation of the cis component of the variability $\sigma_{\text{cis}}^{(X,Y)^2}/(\sigma_{\text{cis}}^{(X,X)} \sigma_{\text{cis}}^{(Y,Y)})$.

## Generalization

Our approach can be applied recursively to experimental data where sample elements are organized hierarchically. In such cases, the approach is applied iteratively, starting from the lowest level, and each iteration uses the aggregated data from the previous iteration as its new elementary data.

As illustrated in Figure S2C, after applying our approach to nascent RNA counts from individual alleles (level 0) and their cell-wise aggregated values (level 1; i.e. summed or averaged within each cell), if the data came for instance from a spatial transcriptomics experiment on a tissue section with various cell types and tissue areas, a second iteration could be applied to the cell-level measurements (level 1) and their aggregated values over cells of the same cell type (level 2). A third iteration could then be applied to cell-type-level measurements (level 2) and their aggregated values over various tissue areas (level 3), and so on.
Beyond this example, our approach is not limited to FISH-based transcriptomics data. It may also be well-suited to analyze data from single-cell and spatial multi-omics technologies⁴⁰, including both microscopy- and sequencing-based methods, since the resulting data is often intrinsically multi-level and hierarchical. (Note that allele-level measurements are not always possible with these approaches, in which case the lowest level in the hierarchy can be directly the single-cell level). In this context (Figure S2D), our approach can be readily extended by incorporating:

(i) at the lower level in the data hierarchy: alongside or instead of gene transcription, other single-cell ‘omics’ measurements, such as DNA accessibility by single-cell ATAC-seq, or protein occupancy or histone marks by single-cell CUT&Tag;

(ii) as extra higher levels in the hierarchy: beyond spatial organization (groups and subgroups of cells, as in Figure S2C), temporal aspects (e.g. developmental time in an organism, response to treatment or perturbation) and experimental conditions (e.g. drug treatments).

## Approach 2 – Separating the contributions of burst co-occurrence and burst size correlations

### – Principle –

Consider the contingency table obtained by classifying genes X and Y at each allele as undergoing or not a transcriptional burst. This equates to binarizing the data into $(\dot{x}_i, \dot{y}_i)$ as

$$ \dot{x}_i = \begin{cases} 0 & \text{if } x_i < c_X \\ 1 & \text{if } x_i \ge c_X \end{cases} \quad \text{and} \quad \dot{y}_i = \begin{cases} 0 & \text{if } y_i < c_Y \\ 1 & \text{if } y_i \ge c_Y \end{cases} $$

where $(c_X, c_Y)$ are the thresholds to consider the genes as undergoing a transcriptional burst.

The associated contingency table is

$$ \begin{array}{c} \text{Gene Y} \\ \left| \begin{array}{c|c} \text{burst} & \begin{bmatrix} \dot{n}_{01} & \dot{n}_{11} \\ \dot{n}_{00} & \dot{n}_{10} \end{bmatrix} \\ \text{no burst} & \text{no burst} \end{array} \right| \\ \hline \text{Gene X} \end{array} $$

where $\dot{n}_{00}, \dot{n}_{01}, \dot{n}_{10}$, and $\dot{n}_{11}$ denote the number of alleles where $(\dot{x}_i, \dot{y}_i)$ are respectively (0,0), (0,1), (1,0), and (1,1) (with $\dot{n}_{00} + \dot{n}_{01} + \dot{n}_{10} + \dot{n}_{11} = n$), corresponding to the partition $\dot{\Omega}_{00}, \dot{\Omega}_{01}, \dot{\Omega}_{10}$, and $\dot{\Omega}_{11}$ of the ensemble of alleles $\Omega$.

This contingency table describes the correlations of the burst occurrences between the two genes, irrespective of burst sizes. The associated covariance is

$$ \dot{\sigma}^2 = \frac{1}{n} \sum_i (\dot{x}_i - \dot{x})(\dot{y}_i - \dot{y}) = \frac{\dot{n}_{00}\dot{n}_{11} - \dot{n}_{01}\dot{n}_{10}}{n^2} $$

with $\dot{x} = \frac{\sum_i \dot{x}_i}{n} = \frac{\dot{n}_{10} + \dot{n}_{11}}{n}$ and $\dot{y} = \frac{\sum_i \dot{y}_i}{n} = \frac{\dot{n}_{01} + \dot{n}_{11}}{n}$ the fraction of bursting alleles for either gene.

What follows shows how the covariance $\sigma_{\text{alleles}}^2$ of the nascent RNA counts $(x_i, y_i)$ can be decomposed as a weighted sum of a contribution from the covariance of the binarized data $\dot{\sigma}^2$, describing correlations of burst occurrences independently of their sizes, and other components, altogether noted $\hat{\sigma}^2$, describing correlations involving burst sizes.
We first describe this approach independently of the previous one to separate cis and trans effect and address in the subsequent subsection how both approaches can be combined. For consistency, we first note here $\sigma^2 = \sigma_{alleles}^2$ the covariance of the non-binarized data $(x_i, y_i)$.

## – Definitions and notation –

Let us introduce the following terminology and notation:

* We call quadrants the four regions of the $(x_i, y_i)$ scatter plot, defined by the thresholds $c_X$ and $c_Y$, corresponding to the partition of alleles into $\dot{\Omega}_{00}$, $\dot{\Omega}_{01}$, $\dot{\Omega}_{10}$, and $\dot{\Omega}_{11}$.

* Quadrant-wise averages: We note $(x_q^+, y_q^+)$ the averages of $(x_i, y_i)$ within each quadrant

$$ \begin{cases} x_q^+ = \sum_{i \in \dot{\Omega}_q} x_i / \dot{n}_q \\ y_q^+ = \sum_{i \in \dot{\Omega}_q} y_i / \dot{n}_q \end{cases} \quad \text{with} \quad q \in \{00, 01, 10, 11\} $$

and use the shorthand $(x_i^+, y_i^+)$ to represent, for each allele $i$, the values of $(x_q^+, y_q^+)$ in its quadrant (i.e. $(x_i^+, y_i^+) \equiv (x_q^+, y_q^+)$ with $q$ such that $i \in \dot{\Omega}_q$). The notation ‘+’ symbolizes the $(x_i, y_i)$ scatter plot being divided into quadrants.

* Half-plane marginal averages are defined as follows:

* $x_0^\|$ and $x_1^\|$ denote the averages of the $x_i$ values that are respectively lower or greater than the threshold $c_X$

$$ x_0^\| = \frac{\sum_{i, x_i < c_X} x_i}{\sum_{i, x_i < c_X} 1} = \frac{\sum_{i \in \dot{\Omega}_{00} \cup \dot{\Omega}_{01}} x_i}{\dot{n}_{00} + \dot{n}_{01}} \quad \text{and} \quad x_1^\| = \frac{\sum_{i, x_i \ge c_X} x_i}{\sum_{i, x_i \ge c_X} 1} = \frac{\sum_{i \in \dot{\Omega}_{10} \cup \dot{\Omega}_{11}} x_i}{\dot{n}_{10} + \dot{n}_{11}} $$

* and $y_0^-$ and $y_1^-$ are the averages of $y_i$ that are respectively lower or greater than $c_Y$

$$ y_0^- = \frac{\sum_{i, y_i < c_Y} y_i}{\sum_{i, y_i < c_Y} 1} = \frac{\sum_{i \in \dot{\Omega}_{00} \cup \dot{\Omega}_{10}} y_i}{\dot{n}_{00} + \dot{n}_{10}} \quad \text{and} \quad y_0^- = \frac{\sum_{i, y_i \ge c_Y} y_i}{\sum_{i, y_i \ge c_Y} 1} = \frac{\sum_{i \in \dot{\Omega}_{01} \cup \dot{\Omega}_{11}} y_i}{\dot{n}_{01} + \dot{n}_{11}} $$

These are marginal quantities since $x_0^\|$ and $x_1^\|$ are defined independently of $y_i$, and $y_0^-$ and $y_1^-$ independently of $x_i$. On the $(x_i, y_i)$ scatter plot, $x_0^\|$ and $x_1^\|$ (resp. $y_0^-$ and $y_1^-$) are the averages within the two half-planes defined the vertical (resp. horizontal) threshold $c_X$ (resp. $c_Y$). As a notation shorthand, the value $x_0^\|$ or $x_1^\|$ (resp. $y_0^-$ or $y_1^-$) corresponding to each allele $i$ and each quadrant $q$ are noted $x_i^\|$ (resp. $y_i^-$) and $x_q^\|$ (resp. $y_q^-$).

* The binarized data $(\dot{x}_i, \dot{y}_i)$, as defined earlier, are 0 and 1 values representing whether $x_i$ (resp. $y_i$) is below (0) or above (1) the threshold $c_X$ (resp. $c_Y$).

To give a few concrete examples: $x_{11}^+$ is the average number of nascent RNAs on gene X when both genes are bursting simultaneously. $y_{01}^+$ is the average number of nascent RNAs on gene Y when gene Y is bursting and gene X is not. $x_1^\|$ is average number of nascent RNAs on gene X when bursting. $\dot{y}_i = 1$ indicates that gene Y is bursting on allele $i$, in which case $y_i^- = y_1^-$.

Note that the raw nascent RNA counts $(x_i, y_i)$, their quadrant-wise averages $(x_i^+, y_i^+)$, and their half-plane marginal averages $(x_i^\|, y_i^\|)$ all have the same average $(x, y)$

$$ x = \frac{\sum_i x_i}{n} = \frac{\sum_i x_i^+}{n} = \frac{\sum_i x_i^\|}{n} = \frac{\sum_q \dot{n}_q x_q^+}{n} = \frac{\sum_q \dot{n}_q x_q^\|}{n} \quad \text{and} \quad y = \frac{\sum_i y_i}{n} = \frac{\sum_i y_i^+}{n} = \frac{\sum_i y_i^\|}{n} = \frac{\sum_q \dot{n}_q y_q^+}{n} = \frac{\sum_q \dot{n}_q y_q^\|}{n} $$
and that $x = x_0^\| (1 - \dot{x}) + x_1^\| \dot{x}$ and $y = y_0^- (1 - \dot{y}) + y_1^- \dot{y}$.
## Derivations

Our approach to connect the statistics of the RNA counts and the binarized data consists in expressing the former from the latter, using the three following steps:

$$ 
\begin{align*}
x_i &= (x_i - x_i^+) + \underbrace{x_i^+}_{\to \Delta x_i^+ + x_i^-} + \underbrace{x_i^-}_{\to x_0^+ + (x_1^- - x_0^-) \dot{x}_i} \dot{x}_i \\
y_i &= (y_i - y_i^+) + \Delta y_i^+ + \underbrace{y_0^- + (y_1^- - y_0^-) \dot{y}_i}_{\text{(same as above)}} \quad \text{(same as above)}
\end{align*}
$$

and examining how this breaks down the covariance into meaningful additive components.

Step (i) consists in applying the law of total covariance to the four quadrants:

$$ 
\begin{align*}
\sigma^2 &= \frac{1}{n} \sum_i (x_i - x)(y_i - y) \\
&= \frac{1}{n} \sum_q \sum_{i \in \hat{\Omega}_q} (x_i - x_q^+ + x_q^+ - x)(y_i - y_q^+ + y_q^+ - y) \\
&= \underbrace{\frac{1}{n} \sum_q \left[ \dot{n}_q \frac{\sum_{i \in \hat{\Omega}_q} (x_i - x_q^+)(y_i - y_q^+)}{\dot{n}_q} \right]}_{\text{Average of...}} + \underbrace{\frac{1}{n} \sum_q \dot{n}_q (x_q^+ - x)(y_q^+ - y)}_{\text{Covariance of quadrant-wise averages}}
\end{align*}
$$

The first term quantifies the covariance within each individual quadrant, among which quadrant $\hat{\Omega}_{11}$ describes how burst sizes are covary when both genes are bursting simultaneously, while the other quadrants should have a small contribution since one or both genes have close to zero nascent RNA count.

The second term quantifies how the average nascent RNA counts within quadrants covary across the four quadrants. It differs from the covariance $\dot{\sigma}^2$ of the contingency table because, although the number of elements in each quadrant are the same as in the contingency table

$[\dot{n}_{01} \ \dot{n}_{11}]$, the associated values are not the unit square $\begin{bmatrix} (0,1) & (1,1) \\ (0,0) & (1,0) \end{bmatrix}$, but $\begin{bmatrix} (x_{01}^+, y_{01}^+) & (x_{11}^+, y_{11}^+) \\ (x_{00}^+, y_{00}^+) & (x_{10}^+, y_{10}^+) \end{bmatrix}$, which in the general case is a non-rectangle quadrilateral. If it were a rectangle, then, as we show below, this second term in the equation above would simply be $\dot{\sigma}^2$ multiplied by a prefactor. Hence, the rational for steps (ii) and (iii) is to map this quadrilateral first onto a rectangle, then onto the unit square.

Step (ii) consists in re-expressing the quadrant-wise averages $(x_q^+, y_q^+)$ in terms of deviations $(\Delta x_q^+, \Delta y_q^+)$ from the half-plane marginal averages, noted $\Delta x_q^+ = x_q^+ - x_q^-$ and $\Delta y_q^+ = y_q^+ - y_q^-$.

$$ 
\underbrace{\begin{bmatrix} (x_{01}^+, y_{01}^+) & (x_{11}^+, y_{11}^+) \\ (x_{00}^+, y_{00}^+) & (x_{10}^+, y_{10}^+) \end{bmatrix}}_{\text{Quadrant-wise averages}} = \underbrace{\begin{bmatrix} (x_0^-, y_1^-) & (x_1^-, y_1^-) \\ (x_0^-, y_0^-) & (x_1^-, y_0^-) \end{bmatrix}}_{\text{Half-plane marginal averages}} + \underbrace{\begin{bmatrix} (\Delta x_{01}^+, \Delta y_{01}^+) & (\Delta x_{11}^+, \Delta y_{11}^+) \\ (\Delta x_{00}^+, \Delta y_{00}^+) & (\Delta x_{10}^+, \Delta y_{10}^+) \end{bmatrix}}_{\text{Deviations}}
$$

The deviations terms $(\Delta x_q^+, \Delta y_q^+)$ capture how the average burst size of either gene varies when conditioned on the fact that the other gene is bursting or not. For instance, $\Delta x_{11}^+$ and $\Delta x_{10}^+$ represent by how much the average nascent RNA counts on gene X during a burst (i.e. $x_1^-$) changes knowing whether gene Y is bursting or not (i.e. difference with $x_{11}^+$ or $x_{10}^+$).
The covariance of the quadrant-wise averages (second term in the previous equation) rewrites

$$ \frac{1}{n} \sum_{q} \dot{n}_{q} (x_{q}^{+} - x)(y_{q}^{+} - y) = \frac{1}{n} \sum_{q} \dot{n}_{q} (x_{q}^{\perp} - x)(y_{q}^{\perp} - y) + \frac{1}{n} \sum_{q} \dot{n}_{q} \left[ \Delta x_{q}^{+} \Delta y_{q}^{+} + \Delta x_{q}^{+} (y_{q}^{\perp} - y) + \Delta y_{q}^{+} (x_{q}^{\perp} - x) \right] = \hat{\sigma}^{2} $$

where we factored out into the second term all the contributions of the deviation $(\Delta x_q^+, \Delta y_q^+)$. This second term hence represents how the size of bursts on either gene depends on the presence of a burst on the other gene. Combined with the average covariance within quadrants (first term of the previous equation), this leads to a term we call $\hat{\sigma}^2$, capturing altogether how the size of bursts on either gene depends on the activity of the other gene:

$$ \hat{\sigma}^2 = \frac{1}{n} \sum_q \dot{n}_q \left[ \frac{\sum_{i \in \hat{\Omega}_q} (x_i - x_q^+)(y_i - y_q^+)}{\dot{n}_q} + \Delta x_q^+ \Delta y_q^+ + \Delta x_q^+ (y_q^- - y) + \Delta y_q^+ (x_q^- - x) \right] $$

Finally, the remaining term, labeled $\check{\sigma}^2$ above, is a simple rescaling of $\dot{\sigma}^2$, the covariance of the contingency table. Indeed, applying step (iii) to express by rewriting and substituting the means $(x, y)$ by their expression in terms of $(\dot{x}, \dot{y})$, we obtain

$$ \begin{aligned} \check{\sigma}^2 &= \frac{1}{n} \sum_i (x_i^\perp - x)(y_i^\perp - y) \\ &= (x_1^\perp - x_0^\perp)(y_1^\perp - y_0^\perp) \frac{1}{n} \sum_i (\dot{x}_i - \dot{x})(\dot{y}_i - \dot{y}) \\ &= (x_1^\perp - x_0^\perp)(y_1^\perp - y_0^\perp) \frac{\dot{n}_{00} \dot{n}_{11} - \dot{n}_{01} \dot{n}_{10}}{n^2} \\ &= \hat{\sigma}^2 \text{ (covariance of contingency table)} \end{aligned} $$

Therefore, we obtain

$$ \begin{array}{c} \text{Total covariance} \\ \underline{\hspace{2.5cm}} \\ \sigma^2 \end{array} = \begin{array}{c} \text{Burst co-occurrence} \\ \underline{\hspace{2.5cm}} \\ \check{\sigma}^2 \end{array} + \begin{array}{c} \text{Burst size correlations} \\ \underline{\hspace{2.5cm}} \\ \hat{\sigma}^2 \end{array} $$

The first term on the right-hand side stems from the contingency table and exclusively **quantifies how the occurrence of bursts correlates between genes**, irrespective of their sizes. The second term captures effects related to burst sizes, **quantifying how the size of bursts on either gene correlates with the activity of the other**.
## Application

Practically, $\sigma^2$, $\check{\sigma}^2$, and $\hat{\sigma}^2$ can be easily calculated as follows:

a. Calculate the covariance $\sigma^2$ from the nascent RNA counts $(x_i, y_i)$.

b. Calculate the covariance $\check{\sigma}^2$ from the $(x_i^1, y_i^-)$ values (simply obtained by replacing all the $x_i$ values below $c_X$ by their average, all the $x_i$ values above $c_X$ by their average, all the $y_i$ values below $c_Y$ by their average, all the $y_i$ values above $c_Y$ by their average).

c. Calculate the third covariance as the difference between the first two: $\hat{\sigma}^2 = \sigma^2 - \check{\sigma}^2$.

Correlation matrices are built by considering different gene pairs for X, Y, ... etc.

$$ \underbrace{\begin{bmatrix} \sigma^{(X,X)} & \sigma^{(X,Y)} \\ \sigma^{(Y,X)} & \sigma^{(Y,Y)} \end{bmatrix}}_{\hat{\sigma}^2} = \underbrace{\begin{bmatrix} \check{\sigma}^{(X,X)} & \check{\sigma}^{(X,Y)} \\ \check{\sigma}^{(Y,X)} & \check{\sigma}^{(Y,Y)} \end{bmatrix}}_{\check{\sigma}^2} + \underbrace{\begin{bmatrix} \hat{\sigma}^{(X,Y)} & \hat{\sigma}^{(X,Y)} \\ \hat{\sigma}^{(X,Y)} & \hat{\sigma}^{(X,Y)} \end{bmatrix}}_{\hat{\sigma}^2} $$

Using this decomposition of the covariances, the Pearson correlation can be separated into:

$$ r^{(X,Y)} = \frac{\text{Pearson correlation coefficient}}{\frac{\sigma^{(X,Y)}^2}{\sigma^{(X,X)} \sigma^{(Y,Y)}}} = \frac{\text{Burst co-occurrence:}}{\frac{\check{\sigma}^{(X,Y)}^2}{\sigma^{(X,X)} \sigma^{(Y,Y)}}} + \frac{\text{Burst size correlation:}}{\frac{\hat{\sigma}^{(X,Y)}^2}{\sigma^{(X,X)} \sigma^{(Y,Y)}}} $$

Note that the denominators are the same, so that $\check{\sigma}^{(X,Y)}$ and $\hat{\sigma}^{(X,Y)}$ represent two additive components of $r^{(X,Y)}$. Hence, $\check{\sigma}^{(X,Y)}$ and $\hat{\sigma}^{(X,Y)}$ are not the Pearson correlation coefficients calculated from the covariance matrices $\check{\sigma}^2$ and $\hat{\sigma}^2$.

## Combining both approaches

Our two approaches to separate (i) cis and trans contributions and (ii) burst co-occurrence and burst size correlation effects are orthogonal and can be combined. Indeed, the former

<table><thead><tr><th></th><th>Total =</th><th>Burst co-occurrence</th><th>+</th><th>Burst size correlations</th></tr></thead><tbody><tr><th>Total =</th><td>σ<sup>2</sup><sub>alleles</sub></td><td>σ<sup>2</sup><sub>alleles</sub></td><td></td><td>σ<sup>2</sup><sub>alleles</sub></td></tr><tr><th>trans component + cis component</th><td>σ<sup>2</sup><sub>trans</sub></td><td>σ<sup>2</sup><sub>trans</sub></td><td></td><td>σ<sup>2</sup><sub>trans</sub></td></tr><tr><th></th><td>σ<sup>2</sup><sub>cis</sub></td><td>σ<sup>2</sup><sub>cis</sub></td><td></td><td>σ<sup>2</sup><sub>cis</sub></td></tr></tbody></table>

holds for any allele-level measurements grouped by cells. Hence, we use the exact same equations (page 16) and same categorization of the n alleles into N cells (i.e. ensembles $\Omega_j$) to

a. calculate $\sigma_{alleles}^2$, $\sigma_{trans}^2$, and $\sigma_{cis}^2$ from the $(x_i, y_i)$ values and their cell-wise sums,

b. calculate $\check{\sigma}_{alleles}^2$, $\check{\sigma}_{trans}^2$, and $\check{\sigma}_{cis}^2$ from the $(x_i^1, y_i^-)$ values and their cell-wise sums.

Then, the two approaches are orthogonal in the sense that it is equivalent to:

(c.) apply the equations on page 16 to obtain $\hat{\sigma}_{trans}^2$, and $\hat{\sigma}_{cis}^2$ from the allele-level and cell-level covariances calculated as $\hat{\sigma}_{alleles}^2 = \sigma_{alleles}^2 - \check{\sigma}_{alleles}^2$

c. or, more simply, deduce directly all three covariances $\hat{\sigma}_{alleles}^2$, $\hat{\sigma}_{trans}^2$, and $\hat{\sigma}_{cis}^2$ by subtracting $\check{\sigma}_{alleles}^2$, $\check{\sigma}_{trans}^2$, and $\check{\sigma}_{cis}^2$ from $\sigma_{alleles}^2$, $\sigma_{trans}^2$, and $\sigma_{cis}^2$.

Figure 3B represents the separation of $r_{cis}^{(X,Y)}$ into the two additive sub-components $\check{r}_{cis}^{(X,Y)} = \check{\sigma}_{cis}^{(X,Y)}^2 / (\sigma_{alleles}^{(X,X)} \sigma_{alleles}^{(Y,Y)})$ and $\hat{r}_{cis}^{(X,Y)} = \hat{\sigma}_{cis}^{(X,Y)}^2 / (\sigma_{alleles}^{(X,X)} \sigma_{alleles}^{(Y,Y)})$. Denominators are allele-level standard deviations, since $r_{cis}^{(X,Y)}$ is itself already a component of $r_{alleles}^{(X,Y)} = r_{trans}^{(X,Y)} + r_{cis}^{(X,Y)}$.
## Sample sizes, replicates, and statistical analysis

Both DMSO-control and E2-induced experiments were performed in 5 replicates, using 5 different single-cell-derived MCF7 cell populations (see *Cell lines and cell culture*). The number of analyzed alleles (n) and nuclei (N) for all 5 replicates were n = 6093, 17066, 12637, 13255, 11148 and N = 1160, 2890, 2187, 1835, 1869 in DMSO condition and n = 7005, 7157, 16476, 21658, 8647 and N = 1274, 1523, 2649, 3159, 1439 in E2 condition. Pearson correlation coefficients and their sub-components were calculated independently for each replicate and all statistical tests were then performed across the 5 replicates using either a 2-sample two-sided t-test with Welch correction (to account for unequal variances) or a 1-sample two-sided t-test, as reported in the following table.

<table><thead><tr><th>Test</th><th>Sample 1<br>(5 replicates)</th><th>Sample 2<br>(5 replicates)</th><th>p-value</th><th>Conclusions</th></tr></thead><tbody><tr><td rowspan="6">2-sample<br>two-sided<br>t-test<br>with Welch<br>correction</td><td>r<sub>alleles</sub><sup>(FOS,JDP2)</sup>|<sub>DMSO</sub></td><td>r<sub>cis</sub><sup>(FOS,JDP2)</sup>|<sub>DMSO</sub></td><td>6.84e-01</td><td>n.s.</td></tr><tr><td>r<sub>alleles</sub><sup>(FOS,BATF)</sup>|<sub>DMSO</sub></td><td>r<sub>cis</sub><sup>(FOS,BATF)</sup>|<sub>DMSO</sub></td><td>3.63e-01</td><td>n.s.</td></tr><tr><td>r<sub>alleles</sub><sup>(JDP2,BATF)</sup>|<sub>DMSO</sub></td><td>r<sub>cis</sub><sup>(JDP2,BATF)</sup>|<sub>DMSO</sub></td><td>4.14e-01</td><td>n.s.</td></tr><tr><td>r<sub>alleles</sub><sup>(FOS,JDP2)</sup>|<sub>E2</sub></td><td>r<sub>cis</sub><sup>(FOS,JDP2)</sup>|<sub>E2</sub></td><td>9.95e-01</td><td>n.s.</td></tr><tr><td>r<sub>alleles</sub><sup>(FOS,BATF)</sup>|<sub>E2</sub></td><td>r<sub>cis</sub><sup>(FOS,BATF)</sup>|<sub>E2</sub></td><td>9.18e-02</td><td>n.s.</td></tr><tr><td>r<sub>alleles</sub><sup>(JDP2,BATF)</sup>|<sub>E2</sub></td><td>r<sub>cis</sub><sup>(JDP2,BATF)</sup>|<sub>E2</sub></td><td>5.94e-01</td><td>n.s.</td></tr><tr><td rowspan="6">1-sample<br>two-sided<br>t-test</td><td>r<sub>cis</sub><sup>(FOS,JDP2)</sup>|<sub>DMSO</sub></td><td>—</td><td>2.98e-04</td><td>**</td></tr><tr><td>r<sub>cis</sub><sup>(FOS,JDP2)</sup>|<sub>E2</sub></td><td>—</td><td>4.69e-03</td><td>**</td></tr><tr><td>r<sub>cis</sub><sup>(FOS,BATF)</sup>|<sub>DMSO</sub></td><td>—</td><td>7.12e-03</td><td>**</td></tr><tr><td>r<sub>cis</sub><sup>(FOS,BATF)</sup>|<sub>E2</sub></td><td>—</td><td>1.56e-04</td><td>**</td></tr><tr><td>r<sub>cis</sub><sup>(JDP2,BATF)</sup>|<sub>DMSO</sub></td><td>—</td><td>3.78e-04</td><td>**</td></tr><tr><td>r<sub>cis</sub><sup>(JDP2,BATF)</sup>|<sub>E2</sub></td><td>—</td><td>1.01e-03</td><td>**</td></tr><tr><td rowspan="10">2-sample<br>two-sided<br>t-test<br>with Welch<br>correction</td><td>r<sub>cis</sub><sup>(FOS,JDP2)</sup>|<sub>DMSO</sub></td><td>r<sub>cis</sub><sup>(FOS,JDP2)</sup>|<sub>E2</sub></td><td>4.92e-02</td><td>*</td></tr><tr><td>r<sub>cis</sub><sup>(FOS,BATF)</sup>|<sub>DMSO</sub></td><td>r<sub>cis</sub><sup>(FOS,BATF)</sup>|<sub>E2</sub></td><td>2.41e-04</td><td>**</td></tr><tr><td>r<sub>cis</sub><sup>(JDP2,BATF)</sup>|<sub>DMSO</sub></td><td>r<sub>cis</sub><sup>(JDP2,BATF)</sup>|<sub>E2</sub></td><td>7.28e-03</td><td>**</td></tr><tr><td>r<sub>cis</sub><sup>(FOS,JDP2)</sup>|<sub>DMSO</sub></td><td>r<sub>cis</sub><sup>(JDP2,BATF)</sup>|<sub>DMSO</sub></td><td>1.02e-03</td><td>**</td></tr><tr><td>r<sub>cis</sub><sup>(FOS,BATF)</sup>|<sub>DMSO</sub></td><td>r<sub>cis</sub><sup>(JDP2,BATF)</sup>|<sub>DMSO</sub></td><td>5.85e-04</td><td>**</td></tr><tr><td>r<sub>cis</sub><sup>(FOS,JDP2)</sup>|<sub>E2</sub></td><td>r<sub>cis</sub><sup>(JDP2,BATF)</sup>|<sub>E2</sub></td><td>1.86e-03</td><td>**</td></tr><tr><td>r<sub>cis</sub><sup>(FOS,BATF)</sup>|<sub>E2</sub></td><td>r<sub>cis</sub><sup>(JDP2,BATF)</sup>|<sub>E2</sub></td><td>7.31e-03</td><td>**</td></tr><tr><td>r<sub>cis</sub><sup>(FOS,JDP2)</sup>|<sub>DMSO</sub></td><td>r<sub>cis</sub><sup>(FOS,JDP2)</sup>|<sub>DMSO</sub></td><td>9.51e-01</td><td>n.s.</td></tr><tr><td>r<sub>cis</sub><sup>(FOS,BATF)</sup>|<sub>DMSO</sub></td><td>r<sub>cis</sub><sup>(FOS,BATF)</sup>|<sub>DMSO</sub></td><td>8.43e-01</td><td>n.s.</td></tr><tr><td>r<sub>cis</sub><sup>(JDP2,BATF)</sup>|<sub>DMSO</sub></td><td>r<sub>cis</sub><sup>(JDP2,BATF)</sup>|<sub>DMSO</sub></td><td>4.67e-01</td><td>n.s.</td></tr><tr><td>r<sub>cis</sub><sup>(FOS,JDP2)</sup>|<sub>E2</sub></td><td>r<sub>cis</sub><sup>(FOS,JDP2)</sup>|<sub>E2</sub></td><td>3.29e-01</td><td>n.s.</td></tr><tr><td>r<sub>cis</sub><sup>(FOS,BATF)</sup>|<sub>E2</sub></td><td>r<sub>cis</sub><sup>(FOS,BATF)</sup>|<sub>E2</sub></td><td>1.93e-03</td><td>**</td></tr><tr><td>r<sub>cis</sub><sup>(JDP2,BATF)</sup>|<sub>E2</sub></td><td>r<sub>cis</sub><sup>(JDP2,BATF)</sup>|<sub>E2</sub></td><td>2.21e-02</td><td>*</td></tr></tbody></table>

Red text highlights the differences being tested between sample 1 and sample 2. The background colors match the color code in the figures, denoting the 3 different gene pairs. (n.s. *p* ≥ 0.05; * *p* < 0.05; ** *p* < 0.01)
**Resource availability:** Data and code are available on a Zenodo repository at http://doi.org/10.5281/zenodo.14865160. This includes all the processed data and analysis scripts, and a subset of raw image files (i.e. first images for each replicate) to match the size limitations of Zenodo repositories (50 GB). The whole dataset (1.6 TB) is available upon request. Image analysis scripts for spot detection and nuclear segmentation are available at https://github.com/CoulonLab/FISHingRod.

**Acknowledgments:** The authors wish to thank Magda Bienko for advice on FISH approaches, Vittore F. Scolari for help with Hi-C data processing, and all the members of our research group for discussions and comments. We also thank the Flow Cytometry Core Facility of Institut Curie (CYTPIC), the BMBC platform of UMR168, and the Cell and Tissue Imaging Facility (PICT-IBiSA) of Institut Curie, member of the national infrastructure France-BioImaging (https://ror.org/01y7vt929) supported by the French National Research Agency (ANR-24-INBS-0005 FBI BIOGEN). This work received funding from the Centre National de la Recherche Scientifique (CNRS); the Institut Curie; the ATIP-Avenir program of CNRS and INSERM, the Plan Cancer of the French ministry for research and health; the European Union's Horizon 2020 research and innovation program, through the Marie Sklodowska-Curie grant agreement No 666003 (H2020-MSCA-COFUND-2014 IC-3i PhD program of Institut Curie) and the European Research Council (ERC) (grant agreement 757956, project "4D-GenEx"); the LabEx DEEP (ANR-11-LABX-0044 and ANR-10-IDEX-0001-02); and the program Fondation ARC (grant agreement PJA 20161204869).

**Author contributions:** MAK, KJEB, and AC designed the study. AC obtained funding. MAK and KJEB performed experiments. IAL, JR, CJ, and AB produced essential reagents. MAK, KJEB, and AC developed analysis methods and performed data analysis. MAK, KJEB, and AC wrote the manuscript with inputs from all the co-authors.

**Declaration of interests:** The authors declare no competing interests.

**Declaration of generative AI and AI-assisted technologies:** Generative AI was used occasionally during the preparation of this publication for enhancing clarity. After using this tool, the authors reviewed and edited the resulting text as needed and take full responsibility for the content of the publication.

## References

1. Nora, E. P. et al. Spatial partitioning of the regulatory landscape of the X-inactivation centre. *Nature* **485**, 381–385 (2012).

2. Dixon, J. R. et al. Topological domains in mammalian genomes identified by analysis of chromatin interactions. *Nature* **485**, 376–380 (2012).

3. Zuin, J. et al. Nonlinear control of transcription through enhancer–promoter interactions. *Nature* **604**, 571–577 (2022).

4. Lupiánez, D. G. et al. Disruptions of Topological Chromatin Domains Cause Pathogenic Rewiring of Gene-Enhancer Interactions. *Cell* **161**, 1012–1025 (2015).

5. Hnisz, D., Day, D. S. & Young, R. A. Insulated Neighborhoods: Structural and Functional Units of Mammalian Gene Control. *Cell* **167**, 1188–1200 (2016).

6. Da Costa-Nunes, J. A. & Noordermeer, D. TADs: Dynamic structures to create stable regulatory functions. *Current Opinion in Structural Biology* **81**, 102622 (2023).

7. Chakraborty, S. et al. Enhancer–promoter interactions can bypass CTCF-mediated boundaries and contribute to phenotypic robustness. *Nat Genet* **55**, 280–290 (2023).

8. Despang, A. et al. Functional dissection of the Sox9–Kcnj2 locus identifies nonessential and instructive roles of TAD architecture. *Nat Genet* **51**, 1263–1271 (2019).

9. Symmons, O. et al. Functional and topological characteristics of mammalian regulatory domains. *Genome Research* **24**, 390–400 (2014).
10. Le Dily, F. et al. Distinct structural transitions of chromatin topological domains correlate with coordinated hormone-induced gene regulation. *Genes Dev.* **28**, 2151–2162 (2014).

11. Zhan, Y. et al. Reciprocal insulation analysis of Hi-C data shows that TADs represent a functionally but not structurally privileged scale in the hierarchical folding of chromosomes. *Genome Res* **27**, 479–490 (2017).

12. Tarbier, M. et al. Nuclear gene proximity and protein interactions shape transcript covariations in mammalian single cells. *Nat Commun* **11**, 5445 (2020).

13. De Smet, R. & Marchal, K. Advantages and limitations of current network inference methods. *Nat Rev Microbiol* **8**, 717–729 (2010).

14. Bohrer, C. H., Fursova, N. A. & Larson, D. R. Enhancers: A Focus on Synthetic Biology and Correlated Gene Expression. *ACS Synth. Biol.* **13**, 3093–3108 (2024).

15. Raj, A., Peskin, C. S., Tranchina, D., Vargas, D. Y. & Tyagi, S. Stochastic mRNA Synthesis in Mammalian Cells. *PLoS Biol* **4**, e309 (2006).

16. Bartman, C. R., Hsu, S. C., Hsiung, C. C.-S., Raj, A. & Blobel, G. A. Enhancer Regulation of Transcriptional Bursting Parameters Revealed by Forced Chromatin Looping. *Molecular Cell* **62**, 237–247 (2016).

17. Bohrer, C. H. & Larson, D. R. Synthetic analysis of chromatin tracing and live-cell imaging indicates pervasive spatial coupling between genes. *eLife* **12**, e81861 (2023).

18. Fukaya, T., Lim, B. & Levine, M. Enhancer Control of Transcriptional Bursting. *Cell* **166**, 358–368 (2016).

19. Levo, M. et al. Transcriptional coupling of distant regulatory genes in living embryos. *Nature* **605**, 754–760 (2022).

20. Lenstra, T. L., Coulon, A., Chow, C. C. & Larson, D. R. Single-Molecule Imaging Reveals a Switch between Spurious and Functional ncRNA Transcription. *Molecular Cell* **60**, 597–610 (2015).

21. Patel, H. P. et al. DNA supercoiling restricts the transcriptional bursting of neighboring eukaryotic genes. *Molecular Cell* **83**, 1573–1587.e8 (2023).

22. Hah, N. et al. A Rapid, Extensive, and Transient Transcriptional Response to Estrogen Signaling in Breast Cancer Cells. *Cell* **145**, 622–634 (2011).

23. Rodriguez, J. et al. Intrinsic Dynamics of a Human Gene Reveal the Basis of Expression Heterogeneity. *Cell* **176**, 213-226.e18 (2019).

24. Atsaves, V., Leventaki, V., Rassidakis, G. Z. & Claret, F. X. AP-1 Transcription Factors as Regulators of Immune Responses in Cancer. *Cancers* **11**, 1037 (2019).

25. Gelali, E. et al. iFISH is a publically available resource enabling versatile DNA FISH to study genome architecture. *Nat Commun* **10**, 1636 (2019).

26. Swain, P. S., Elowitz, M. B. & Siggia, E. D. Intrinsic and extrinsic contributions to stochasticity in gene expression. *Proc. Natl. Acad. Sci. U.S.A.* **99**, 12795–12800 (2002).

27. Abdulla, A. Z., Vaillant, C. & Jost, D. Painters in chromatin: a unified quantitative framework to systematically characterize epigenome regulation and memory. *Nucleic Acids Research* **50**, 9083–9104 (2022).

28. Abdulla, A. Z., Salari, H., Tortora, M. M. C., Vaillant, C. & Jost, D. 4D epigenomics: deciphering the coupling between genome folding and epigenomic regulation with biophysical modeling. *Current Opinion in Genetics & Development* **79**, 102033 (2023).

29. Lee, B.-K. et al. Cell-type specific and combinatorial usage of diverse transcription factors revealed by genome-wide binding studies in multiple human cells. *Genome Res.* **22**, 9–24 (2012).

30. Liu, Z. et al. Enhancer Activation Requires trans-Recruitment of a Mega Transcription Factor Complex. *Cell* **159**, 358–373 (2014).

31. Castoria, G. et al. Non-transcriptional action of oestradiol and progestin triggers DNA synthesis. *EMBO J* **18**, 2500–2510 (1999).

32. Hamelers, I. H., Van Schaik, R. F., Sussenbach, J. S. & Steenbergh, P. H. 17beta-Estradiol responsiveness of MCF-7 laboratory strains is dependent on an autocrine signal activating the IGF type I receptor. *Cancer Cell Int* **3**, 10 (2003).

33. Rocca, A., Farolfi, A., Bravaccini, S., Schirone, A. & Amadori, D. Palbociclib (PD 0332991): targeting the cell cycle machinery in breast cancer. *Expert Opinion on Pharmacotherapy* **15**, 407–420 (2014).

34. Salic, A. & Mitchison, T. J. A chemical method for fast and sensitive detection of DNA synthesis in vivo. *Proceedings of the National Academy of Sciences* **105**, 2415–2420 (2008).

35. Lypova, N. et al. Targeting Palbociclib-Resistant Estrogen Receptor-Positive Breast Cancer Cells via Oncolytic Virotherapy. *Cancers (Basel)* **11**, E684 (2019).

36. Bintu, B. et al. Super-resolution chromatin tracing reveals domains and cooperative interactions in single cells. *Science* **362**, eaau1783 (2018).

37. Beliveau, B. J. et al. In situ super-resolution imaging of genomic DNA with OligoSTORM and OligoDNA-PAINT. *Methods Mol Biol* **1663**, 231–252 (2017).

38. Edelstein, A., Amodaj, N., Hoover, K., Vale, R. & Stuurman, N. Computer Control of Microscopes Using µManager. *Current Protocols in Molecular Biology* **92**, (2010).

39. Abdennur, N. & Mirny, L. A. Cooler: scalable storage for Hi-C data and other genomically labeled arrays. *Bioinformatics* **36**, 311–316 (2020).

40. Vandereyken, K., Sifrim, A., Thienpont, B. & Voet, T. Methods and applications for single-cell and spatial multi-omics. *Nat Rev Genet* **24**, 494–515 (2023).