Multiplex chromatin interaction analysis by signal processing and statistical algorithms

The single-molecule multiplex chromatin interaction data generated by emerging non-ligation-based 3D genome mapping technologies provide novel insights into high dimensional chromatin organization, yet introduce new computational challenges. We developed MIA-Sig (https://github.com/TheJacksonLaboratory/mia-sig.git), an algorithmic framework to de-noise the data, assess the statistical significance of chromatin complexes, and identify topological domains and inter-domain contacts. On chromatin immunoprecipitation (ChIP)-enriched data, MIA-Sig can clearly distinguish the protein-associated interactions from the non-specific topological domains.

transcriptomics data are not apt for multiplex chromatin interactions data since: 1) the signal for chromatin interactions is point data (fragment is captured or not captured) rather than continuously valued data (gene expression level), and 2) multiplex chromatin interaction data are inherently sparser than the single-cell transcriptomics data, due to the lack of cell barcodes.
Therefore, we devised a distance test with an entropy filter based on the biological knowledge that most meaningful chromatin interactions occur in a certain distance range, while those outside the range are likely noise (Lajoie et al., 2015). By converting the distances between fragments into a probability vector, we compute the normalized Shannon entropy (Shannon, 1948), ranging from 0 to 1. If a droplet contains a single complex, the fragments are presumably close and equally spaced, resulting in high entropy close to 1. In the case of a doublet, two independent complexes would be separated by a single large distance, resulting in low entropy close to 0, which can then be separated into two singlets (Figure 1b). To identify significant chromatin complexes, a resampling-based distance test is applied before and after the entropy filter (Figure 1c; Methods). As a result, we retained 55,995 statistically significant complexes in the Drosophila S2 ChIA-Drop data out of 3,075,926 putative complexes (Supplementary Figure 1). Filtering to retain significant complexes preserves the TADs along the diagonal of the 2D heat maps, while reducing the off-diagonal noise (Figure 1d; visualization through Juicebox (Durand et al., 2016)). A shift in distance distributions supports that meaningful interactions are captured within 10 kb and 1 Mb (Figure 1e; Supplementary  Figure 2). Of the significant chromatin complexes, 15,055 (27%) were from the entropy filtering step that resolved doublets and triplets (Supplementary Figure 3).
From the significant complexes, it is desirable to automatically call TADs for downstream analyses. Many TAD calling algorithms exist for Hi-C data (Zufferey et al., 2018), yet all are based on pairwise contacts. To retain multiplexity information, we developed an algorithm to call TADs directly from ChIA-Drop data (Methods). The idea is to convert complexes into 1D signal track, then apply wavelet transformation (Mallat, 1989) to smooth the signal while retaining clear change points (Supplementary Figure 4). MIA-Sig called 335 TADs with wider range of sizes than 513 TADs called by pairwise 'Insulation Score' ('InS') approach (Supplementary Figure 5). Compared to InS TADs, the MIA-Sig TADs are less likely to overlap active regions characterized by high H3K27ac and low H3K27me3 (Figure 1f), which are known to be the gaps between TADs in Drosophila (Rowley et al., 2017). This pattern is observed genome-wide: MIA-Sig TADs have higher inactive mark (H3K27me3) in than InS TADs, and MIA-Sig gaps have higher active mark (H3K27ac) than InS gaps (Supplementary Figure 6). Most interactions occur within a single TAD, but 23% of significant complexes also cross 2 or more TADs, consistent with previous findings (Paulsen et al., 2019). Thus, we identify frequent interactions involving multiple TADs by counting the occurrences and performing a binomial test (Methods). A set of TADs with frequent contacts are ultimately assigned low p-values (Supplementary Figure 7), which can guide the researchers to perform validation experiments.
Similar to ChIA-PET, ChIA-Drop can also enrich chromatin complexes involving a specific protein, such as RNAPII or CTCF. We implemented an enrichment test to estimate the significance of binding intensity of observed chromatin complexes (Figure 2a; Methods) compared to an empirical null distribution (Supplementary Figure 8). As a result, we retain significant complexes with their fragments in highly enriched domains characterized by high RNA-seq expression and H3K27ac signal with abundant RNAPII ChIA-PET loops (Figure 2b). Genome-wide patterns confirm that significant complexes are biased towards active regions, whereas insignificant complexes are not (Supplementary Figure 9). Moreover, significant complexes have higher median H3K27ac signals and lower median H3K27me3 signals than insignificant complexes (Figure 2c-d). A detailed view around a few genes shows that significant complexes are more likely to retain promoter-centric interactions than insignificant complexes (Figure 2d; visualization through ChIA-View ). This pattern is prevalent genome-wide, with 69% of significant complexes containing at least 1 promoter compared to only 30% of insignificant complexes (Figure 2f). Notably, significant complexes are most likely to capture one active promoter and one or more non-promoters-possibly enhancers-while insignificant complexes are prone to detect interactions among nonpromoters (Supplementary Figure 10). Among the promoter-involving fragments, those in significant complexes have higher median gene expression than those in insignificant ones.
As with many experimental protocols, the chromatin immunoprecipitation step is not 100% efficient and typically yields a 20-40% efficiency rate (Tang et al., 2015). Thus, we take advantage of the fact that enriched ChIA-Drop data sets also contain some background signal for chromatin complexes that did not specifically involve the protein of interest, similar to nonenriched ChIA-Drop data. Through the MIA-Sig enrichment test on RNAPII ChIA-Drop data, we can extract the non-enriched complexes from the insignificant complexes, which approximately emulate the ChIA-Drop data (Figure 2g).
Given that the frontier of the nuclear architecture field is now single-cell and singlemolecule 3D genome mapping, it is imperative to develop algorithms to analyze data from these novel experimental protocols. We have presented an approach to the imminent problem of extracting statistically significant complexes from noisy signals, calling TADs, and identifying frequent inter-TAD contacts (Figure 3). In addition, we offer a practical strategy to extract nonenriched ChIA-Drop from RNAPII ChIA-Drop. As a publicly available software package, MIA-Sig provides a valuable algorithmic framework for multiplex chromatin interaction data to be utilized by the broader scientific community. MIA-Sig nonetheless has its own drawbacks. One key assumption in the distance test is that a fragment far from the other fragments is likely a droplet contamination resulting in a doublet, a behavior yet to be confirmed experimentally and statistically. In performing the enrichment test for RNAPII ChIA-Drop data, we do not use a background distribution model and instead draw an empirical null distribution via random sampling. A disadvantage of this approach is the computational cost, which can be demanding for large human datasets.
We envision that MIA-Sig will be broadly applicable to any type of multiplex chromatin interaction data ranging from ChIA-Drop, SPRITE, to GAM, under the aforementioned assumptions and with modifications (Supplementary Note 1). Here, we focused on the Drosophila ChIA-Drop and RNAPII ChIA-Drop data as a proof of concept and demonstrated that MIA-Sig filters and retains only the highly informative complexes. We anticipate that researchers will be able to utilize this refined data to dissect heterogeneous population within TADs and to identify interactions involving multiple regulatory elements. Figure 1: Performance of MIA-Sig on ChIA-Drop data. (a) ChIA-Drop experiments are designed to encapsulate each chromatin complex in a droplet, but the encapsulation is a random process and sometimes results in more than one complex in a droplet (multiplets). (b) MIA-Sig aims to detect multiplets by computing the normalized Shannon entropy "#$% (Methods). It separates a complex at the largest distance if "#$% is smaller than a threshold, which is 0.7 in this example. (c) Summary statistics of the distance test indicate that entropy filter resolves around 500,000 doublets and 85,000 triplets, from which 15,055 complexes pass the second distance test. (d) 2D heatmap comparison of original (bottom triangle) and significant (upper triangle) complexes demonstrates that MIA-Sig removes off-diagonal noise. (e) Empirical cumulative distribution function for the neighboring distances of original and significant complexes (two-sided Kolmogorov-Smirnov test statistic= 0.47, p-value < 2.2 × 10 /01 ). (f) Comparison of TADs called from significant putative complexes (MIA-Sig) and from enumerating all pairs of fragments (Insulation score). MIA-Sig more specifically separates active regions (high H3K27ac and low H3K27me3) rather than assigning them to TADs. (a) MIA-Sig performs an enrichment test on RNAPII-enriched ChIA-Drop data by retaining complexes with fragments in strong binding regions, which also correspond to RNAPII ChIA-PET peaks. (b) The significant complexes are pronounced in regions with high level of transcription, abundant loops, and active histone mark; insignificant complexes tend to be in inactive regions. (c) Log of H3K27ac signal for fragments in significant and insignificant complexes (one-sided Mann-Whitney U test, p-value < 2.2 × 10 /01 ). (d) Log of H3K27me3 signal for fragments in significant and insignificant complexes (one-sided Mann-Whitney U test, p-value < 2.2 × 10 /01 ). (e) Fragment coverage profile of significant complexes is similar to that of RNAPII ChIA-PET, with 45 promoter-centric multiplex interactions (green: non-promoters, light green: promoters). By contrast, insignificant complexes do not show any strong binding peaks in coverage, and 91 multiplex interactions are non-specific (turquoise: non-promoters, light turquoise: promoters). (f) Genome-wide, significant complexes have higher proportion of active promoter fragments than insignificant complexes do (two-sided K-S test statistic= 0.39, pvalue < 2.2 × 10 /01 ). (g) The insignificant RNAPII ChIA-Drop complexes from the enrichment test are comparable to the significant ChIA-Drop complexes from the distance test. TADs (black lines) are called by MIA-Sig on the latter complexes. ChIA-Drop putative complexes from ChIA-DropBox pipeline are inputs to the distance test, which assigns p-values to each complex to quantify its significance. Refined complexes enable TAD calling directly on the multiplex data. A binomial test identifies frequent contacts among multiple TADs. RNAPII-enriched ChIA-Drop putative complexes are assigned significance according to the level of enrichment.

Data availability
The

Notation
An input dataset contains a set of chromatin complexes, each with 2 or more fragments. Let % be the set of fragments contained in the -th 'observed complex' (OC), for ∈ {1,2, … , }, and = | % | is the size of the set denoting the number of fragments in a complex. Each fragment is subscripted by the complex index and superscripted by the fragment index, and encodes the genomic location of its origin expressed as a triplet of chromosome, start, and end position. The distance between fragments % B and % C is start( % C )-end( % B ), and neighboring (fragment-to-fragment; F2F) distances are encoded in a vector ]. Finally, we can partition complexes 0 , F , … , [ into ] , where is the number of fragments in a complex ( W belongs to J since it has 3 fragments).

Distance test for non-enriched multiplex chromatin interactions data Empirical null distribution and first distance test.
Assuming that complexes are independent of chromosome, we perform the distance test separately for each chromosome. Motivated by the fact that each fragment class ] has distinct distributions in F2F distances, we construct the expected null background distribution by randomly rewiring fragments. Specifically, all neighboring distances EFE ( % ) for ∈ {1,2, … , } are placed in a bucket . For each observed ] , we randomly draw − 1 elements (with replacement) from to create 100,000 'expected complexes' (EC) b ] for ∈ {1,2, … ,100000} and store them in ] ′. Note that since we only care about distance between fragments, we can assume that every fragment starts at (chr, 1, 500) and each fragment is of equal length. In practice, we store minimum information to save compute memory (implementation details below). For each % in ] , we compare its total F2F distance to total F2F distance in ] e and record the proportion of expected complexes that have shorter distances than the observed complexes as the estimated 'raw p-value'. Formally, for a % ∈ ] , where 1 { * } is an indicator function. Assuming that complexes in each fragment class are independent, we subsequently separate the raw p-values by ] and adjust them for multiple hypothesis testing using Benjamini-Hochberg method (Benjamini, Hochberg, 1995) with false discovery rate (FDR) of 0.1. The complexes with adjusted p-value ≤ 0.1 are considered to be statistically significant and are classified as 'pass1' ( ],tBuu0 ). Of those insignificant complexes with adjusted p-value > 0.1, we 'fail1' ( ],vBwx0 ) those with 2 fragments ( % ∈ F with adjusted p-value > 0.1) and treat others in a separate category called 'defer' ( ],Tyv ). These 'deferred' complexes are passed onto the entropy filter to correct for droplet contamination. Entropy filter. Some complexes in the 'deferred' category may be due to the experimental noise that can be computationally detected. Specifically, this step aims to computationally correct for the undesired phenomenon of a droplet containing more than one chromatin complex (referred to as 'doublet' for two, 'triplet' for three, and 'multiplet' for 2 or more). In single-cell RNA-seq (scRNA-seq; single-cell transcriptome) experiments, the outcome of a doublet would be a vector of real numbers indicating average expression of the two cells. By contrast, ChIA-Drop data only provide binary values indicating if a fragment was captured or not, with a variable number of fragments. Therefore, the effect of two complexes accidentally being encapsulated in a single droplet would be a large distance in the data. This assumption is based on the observation from Hi-C and ChIA-PET data analysis that true interactions occur within certain range of genomic span. Our goal is to identify complexes with one dominating distance between fragments. Using the probability vector of the neighboring distance, we quantify the likelihood of a dominating event. Formally, for an observed complex % with fragments and EFE ( % ) = [ 0 , F , … , "/0 ], we compute the normalized Shannon entropy (Shannon, 1948) The normalization factor log F ( − 1) ensures that "#$% ( ) ∈ [0,1] for any probability vector . Generally, "#$% is small when only one or two of w of is large, in which case we presume that a complex is a multiplet and need to separate into singlets. For each observed complexes in the 'deferred' category, we compare its normalized Shannon entropy to the average normalized Shannon entropy of the expected complexes in the corresponding class; if the former is smaller, then we separate the observed complex at the longest distance interaction. In other words, for % ∈ ],Tyv , if where ( % ƒ , % ƒ"0 ) = max EFE ( % ). Furthermore, if the second largest distance is at least 0 of the largest distance, we also separate at the second longest distance. is a variable parameter and we set it to 2 in our analyses; the larger the , the more likelihood of a 'second cut' (implying a triplet . We adjust raw p-values using Benjamini-Hochberg method with false discovery rate (FDR) of 0.1. The complexes with adjusted p-value ≤ 0.1 are classified as 'pass2' ( ],tBuuF ) and others are 'fail2' ( ],vBwxF ). A diagram of the distance test is illustrated in Supplementary Figure 1a. Implementation, results, and analysis. MIA-Sig takes putative chromatin complexes as the input, which are results of the ChIA-DropBox  data processing and visualization pipeline. The 'distance test' python (v3.6) script encompasses all parts using the following packages: numpy, random, statsmodels, itertools, os, sys. We used the parameters --gen dm3 --fdr 0.1 --cef 2 --sz 100000 to run the script on GSM3347523 dataset, which used 1.8 Gigabyte of memory and 13 minutes of cpu time. To save memory, we store minimal information for the null: total distance for expected complexes, and their mean entropy for each fragment class. Two runs with the same parameters should yield identical results because we set seeds in the construction step for the expected complexes. By saving the first 1000 expected complexes for each class in a chromosome, we can compare our expected null model to the biological null model, which is the 'pure DNA' described in . Plotting the neighboring distances, we observed that both the computational null and pure DNA are unimodal with peaks between 1 Mbps and 10 Mbps for all classes (Supplementary Figure 1b). After confirming that our expected complexes do emulate long-range noise, we obtained detailed statistics of each step resulting in 55,995 significant complexes (Supplementary Figure 1c). Complexes in each of the 'original', 'significant' ('pass1'+'pass2'), and 'insignificant' ('fail1'+ 'fail2') categories are converted into a .short format by enumerating over all pairs of fragments in a complex. Three .short files are then converted into .hic files via juicer (v1.7.5) to be visualized in juicebox. A 5 Mbps window on chr3L shows that the original data exhibit both the signal and noise, which are separated by MIA-Sig into significant and insignificant, respectively (Supplementary Figure 2a). The original observed complexes have a bimodal distribution for high fragment classes, which is a distinct behavior from the null distribution (Supplementary Figures 1b, 2b). The density plot further support that significant complexes retained short distances or a mix of short and long distances. By contrast, insignificant complexes are only comprised of unimodal long distances (Supplementary Figure 2b). Consistent with an observation that high-fragment complexes contribute to the structure more than the low-fragment complexes , MIA-Sig assessed the majority of high-fragment complexes as significant (Supplementary  Figure 2c). We next investigated the effects of the entropy filter, which was designed to remove doublets and triplets. Of the 1,452,878 complexes in the deferred category ranging from = 3 to = 8, MIA-Sig identified 60% (869,065) to be singlets, 34% (498,291) to be doublets, and 6% (85,522) to be triplets, yielding 548,342 singletons ( 0 ) and 1,573,871 complexes ( ‹F ) (Supplementary Figure 3). For each class, singletons had the highest normalized Shannon entropy, followed by doublets and triplets. The entropy filter step allowed MIA-Sig to identify additional 15,055 complexes as significant, which amounts to 27% of the total significant complexes.

TAD calling for non-enriched multiplex chromatin interactions data
Generating 1D signal track. Existing TAD calling algorithms for pairwise Hi-C data generally fall into two categories: 1) signal segmentation after conversion from 2D contact maps into 1D tracks measuring interaction intensities along the genome, 2) community detection directly on the 2D heatmap by treating each bin as a node on an undirected graph. We take the first approach and convert our complexes into 1D signal track. A conventional pairwise approach would enumerate over all pairs of fragments in a complex and record their spans. However, multi-fragment complexes may over-contribute since the number of pairs grows quadratically: z " F { = "("/0) F , where is the number of fragments in a complex. Instead, we allow each complex to only contribute linearly in by recording its span weighted by . More precisely, coordinates are (chrom( % 0 ), start( % 0 ), end( % " ), ) for an % with fragments. We finally obtain a 'weighted complex span coverage' by accumulating the coordinates over all given complexes. Smoothing and segmentation. Our next task is to segment the 1D track into regions with high signal and annotate them as TADs. In an ideal case, we can achieve this goal by computing the slope of the signal s and by recording critical points where the slope is 0. However, our signal has a basepair resolution and thus is not smooth, resulting in too many false critical points. A common way to smooth the signal is by a moving average window, but using a large window size would lose the resolution and yield TADs with fuzzy boundaries. Moreover, due to the inherent nature of TAD sizes, a window size parameter optimal in one region may not be optimal in another . We avoid this parameter tuning step by instead applying a discrete wavelet transformation, which decomposes signal into high-frequency component and lowfrequency component (Mallat, 1989). Of note, the low-frequency component generally retains the smoothed version of the signal without affecting the shape, which is helpful for us to find accurate TAD boundaries (Supplementary Figure 4). Using this 'smoothed' signal, we compute the slope and finetune TAD coordinates. Implementation, results, and analysis. The 'tad calling' python (v3.6) script encompasses all parts using the following packages: numpy, os, scipy, pywt, itertools, sys. We used the parameters --cat PASS --bs 1000 -sp drosophila --r dm3 to run the script on significant complexes from the distance test of GSM3347523 dataset, which used 84 Megabyte of memory and 1 minute of cpu time. Before generating the 1D signal track, we separate two fragments if they are more than 100 Kb apart, based on the upper range of general TAD sizes by organisms (Dekker and Heard, 2015). Coverage was generated by BEDtools (Quinlan et al., 2010) using 'genomecov' function and the coverage is binned into 1 kb windows via 'makewindows' and 'map' commands. Signal smoothing was done by pywt package using the parameters 'bior1.1' for wavelet function and '3' for the level. MIA-Sig called 335 TADs over the 6 chromosomes, with a median size of 200 Kb (Supplementary Figure 5a). For a comparison, we also tested insulation score as follows: .hic file (of all pairs of fragments) are converted into contact matrices via Juicer's 'dump' function with a dense matrix option (-d) in the Juicer tool (v1.7.5); insulation score script (https://github.com/dekkerlab/cworld-dekker/tree/master/scripts/perl) is executed with 100Kb insulation square size, 100Kb delta window size for 10Kb resolution contact maps with balanced normalization. Insulation score (InS) called 513 TADs with a median of 150Kb, and did not call any TADs larger than 500 Kb (Supplementary Figure 5b,c). When we examined the gaps (defined as regions between two TADs, if any), MIA-Sig also had wider size range than InS (Supplementary Figure  5d,e). For each TAD called by MIA-Sig and InS, we compute the total H3K27me3 signal and plot the genome-wide behavior (Supplementary Figure 6a). Overall, MIA-Sig has higher inactive signal in TADs than InS. The gap regions in Drosophila are known to be transcriptionally active and should positively correlate with H3K27ac signal. We confirm that MIA-Sig has slightly higher median active signal than InS (Supplementary Figure 6b). Note that we did not perform any normalization by region size because both algorithms segment the genome into either a TAD or a gap, so the region size should also be a feature. Histone marks provide biological evidence that MIA-Sig TADs are inactive and gaps are active, but ChIA-Drop fragment counts provide a direct measure of TAD and gap intensities. Using the BEDtools command 'intersect -c', we count the number of fragments in each region. MIA-Sig generally captured more fragments in TADs than InS did (Supplementary Figure 6c), and less fragments in gaps than InS (Supplementary Figure 6d). Finally, we annotate each fragment in significant and insignificant complexes as 'TAD' or 'gap' as called by MIA-Sig. For each complex, we count the number of TADs with at least 2 fragments within each TAD. Only 5% of the insignificant complexes had fragments in 1 or 2 TADs, and the rest were not contributing to the TAD structure ( Supplementary  Figure 7a), validating the observation from 2D heatmaps. By contrast, only 26% of the significant complexes were not in TADs, a majority (51.3%) in intra-TAD interactions, and many (23%) connected two or more TADs. By observing that 12,884 complexes involve 2 to 21 TADs, we next sought to characterize if multiple complexes connect the same set of TADs

Inter-TAD binomial test for non-enriched multiplex chromatin interaction data Motivation and intuition.
Our goal is to evaluate the statistical significance of these TAD combinations based on the frequency of occurrence measured by the number of complexes therein. The problem is simple for a pair of TADs: we may treat a TAD as a ChIA-PET loop anchor and apply tools based on hypergeometric test. However, our data are now multi-dimensional. For instance, suppose that there are 5 TADs, and 5 combinations 'A-C', 'B-C', 'B-C-D', 'A-B-E', 'A-D-E' (Supplementary Figure 7b). The pair 'B-C' appears 4 times on its own, but also appears 3 times as a part of the triple 'B-C-D'. Moreover, some parts of a combination may appear elsewhere with the same number of TADs: given 'B-C-D' and 'A-C-D', 'C-D' appears twice. Therefore, we propose a counting scheme based on occurrence of 'expanded pairs'. Methods. The notations used defined in this section are independent from those in other sections. We let the th combination be w = { w 0 , w F , … , w • }, where each w " ∈ { 0 , F , … , [ }, = | w | is the number of TADs involved, and we partition each w into the same class ] if | w | = . All pairs of TADs in w are in ( w ) = {{ , }: ≠ , for , ∈ w } and | ( w )| = "("/0) F . For each w , we record the number of pairs in the same class as ( w ) = ™ ™ 1 jqšB(›R oe ) j∈šB(•) •∈ž Ÿ and the number of exact appearance in higher class as ( w ) = ∑ 1 ›R oe ⊂j j∉ž Ÿ . Using these two numbers, we compute the appearance of 'pairs' in same class and higher class: ; the alternative hypothesis is that the observed probability is greater than the expected probability . A detailed example is provided using same notations (Supplementary Figure 7b). Implementation, results, and analysis. A python script 'inter-TAD binomial test' implements the method using packages numpy, itertools, scipy, statsmodels, os, and sys. Of 6,861 unique combinations involving 2 to 21 TADs, 915 (13%) were identified as statistically significant. An example illustrates that a pair of TADs with strong signal in the heatmap and many complexes in the linear view has lower pvalue than that with weak signal (Supplementary Figure 7c). Here, we assumed that frequency of interactions between TADs is independent from their distance and sizes, and we also did not distinguish contacts with 2 fragments from those with 10 fragments. These parameters may be incorporated in the future version.

Enrichment test for RNAPII-enriched multiplex chromatin interaction data
Motivation. The above sections are designed to analyze non-specific multiplex interaction data analogous to the Hi-C data. With an additional step of chromatin immunoprecipitation, protein-enriched multiplex data reveal protein-specific interactions similar to the population average ChIA-PET loops. In a typical ChIA-PET analysis, loops anchored in strong binding peaks are considered to be more reliable than those with weak or no peaks. Extending this notion to the multiplex data, we developed an enrichment test for RNAPII ChIA-Drop data. Our end goal is to retain complexes with fragments in strong binding peaks. One approach is to call peaks and only keep complexes that overlap the peak regions. However, peak calling algorithms have their own model assumptions that may not hold for ChIA-Drop data. Even with accurate peak regions, assigning statistical significance to each complex is not a trivial problem since the null distribution is unclear. Thus, we take an alternative-inevitably the computationally expensiveapproach by sampling the background null distribution for each complex. Statistical test. The idea is to take the observed complex and place it on a random location of the same chromosome and compare the mean coverage between the observed and the expected. Through many rounds of re-sampling, we obtain the p-value by counting the number of occurrences in which the expected coverage exceeds the observed coverage (Supplementary Figure 8a). More precisely, for an observed complex % = { % 0 , % F , … , % " }, we randomly draw an integer ∈ {1, … , length( ℎ ) − start( % 0 )} and the shift = start( % 0 ) − . The first expected complex is then 0 Repeating this process 10,000 times, we obtain b % for ∈ {1, … ,10000}. We can then compute the raw p-value of the th observed complex as: between coordinates and . Raw p-values are separated by chromosomes and are adjusted via the Benjamini-Hochberg method with false discovery rate (FDR) of 0.1. The complexes with adjusted pvalue ≤ 0.1 are considered to be statistically significant and are classified as 'pass'; others are considered insignificant or 'fail'. Implementation, results, and analysis. A python script 'enrichment test' utilizes packages numpy, random, statsmodels, os, sys. GSM3347525 RNAPII ChIA-Drop data are pre-processed to exclude fragments mapped to the repetitive regions in the genome (dm3.rmsk.bed), and 769,803 complexes remain as 'GSM3347525NR'. The most time-consuming part of the algorithm is to obtain the fragment coverage at a given location, since we need to search for a start and end index in a bedgraph or a bigwig file. With at least 769803 * 2 * 10000 = 1.54 × 10 0Y operations, we realized that python implementations of exact search would be intractable. As means to reduce the runtime, we store the bedgraph file into bins of size 10 bps and store only the 4 th column 'value'. The solution then turns into a simple lookup operation, yielding an approximation that is close to the exact solution. Our code is 'parallelized' by chromosome, each using around 5 hours cpu time and 230 Megabytes of memory (Supplementary Figure 8b). MIA-Sig identified 190,226 complexes (24.7%) as statistically significant. We ensure that our empirical null distribution does behave randomly by comparing the enrichment scores of the observed complexes in chr2L with those of 1,000 expected complexes generated for each observed complex (Supplementary Figure 8c). Zooming in further, we note that the histogram of the observed is shifted to the right of the histogram of the expected null (Supplementary Figure 8d). Using the active and inactive regions defined in , we count the number of fragments therein for significant and insignificant complexes (Supplementary Figure 9a). For each active and inactive region, we compute the number of significant complexes fragments and their log10 values are plotted (Supplementary Figure 9b); K-S test supports that significant complexes are indeed more likely to be in active regions than inactive regions. By contrast, insignificant complexes have no bias towards or against active regions (Supplementary Figure 9c). We define a gene promoter as +/-1kb from the transcriptional start site (TSS) annotated by UCSC genome browser. Note that typically +/-250 bp is used for drosophila, but we extend it to accommodate ChIA-Drop protocol-specific features. A gene is active (6466 genes) if total RNA-seq level is greater than 5, and inactive (8874 genes) otherwise. A fragment is 'active promoter' if it overlaps the promoter of an active gene. In general, significant complexes have higher proportion of promoter fragments than insignificant complexes (Supplementary Figure 9d), and the skew is more pronounced for active promoters (Supplementary Figure 9e). Inactive promoters serve as a control, in which both significant and insignificant complexes display similar patterns in the number of inactive promoter fragments (Supplementary Figure 9f,g). Supplementary Figure 3: Effects of the entropy filter. (a) A box plot of normalized Shannon entropy is provided for sub-complexes after filter. Of complexes with three fragments (in J ), 499,613 are identified as 'singlets' due to high entropy, and 284,540 are considered to be 'doublets' due to low entropy. x-axis is labled with the size of sub-complexes, with counts in a parenthesis. For instance, a 5-fragment complex can be split into a singleton (1) and a 4-fragment complex (4), or into a complex (2) and another complex (3); alternatively, a triplet has three sub-complexes (1,1,3 and 1,2,2). A general trend is that entropy is highest for those without any splits, lowest for a doublet with a singleton, and increases as the size of sub-complexes balance to be roughly equal.

Supplementary Note
Supplementary Note 1. Extending MIA-Sig to GAM and SPRITE data This note discusses the technical principles of the current three experimental protocols (GAM, SPRITE, and ChIA-Drop) for capturing multiplex interaction data, and the applicability of MIA-Sig to all multiplex data generated by different methods including GAM and SPRITE.