GPSeq reveals the radial organization of chromatin in the cell nucleus

With the exception of lamina-associated domains, the radial organization of chromatin in mammalian cells remains largely unexplored. Here, we describe genomic loci positioning by sequencing (GPSeq), a genome-wide method for inferring distances to the nuclear lamina all along the nuclear radius. GPSeq relies on gradual restriction digestion of chromatin from the nuclear lamina towards the nucleus center, followed by sequencing of the generated cut sites. Using GPSeq, we mapped the radial organization of the human genome at 100 kb resolution, which revealed radial patterns of genomic and epigenomic features, gene expression, as well as A/B subcompartments. By combining radial information with chromosome contact frequencies measured by Hi-C, we substantially improved the accuracy of whole-genome structure modeling. Finally, we charted the radial topography of DNA double-strand breaks, germline variants and cancer mutations, and found that they have distinctive radial arrangements in A/B subcompartments. We conclude that GPSeq can reveal fundamental aspects of genome architecture.


Supplementary Figures
Dashed circles: nuclear lamina. Solid circles: nuclear center. Red segments represent repetitive sequences that were masked out before calculating the GPSeq score (see Supplementary Table   10). Sample size information for each plot is available in Supplementary Table 11. All the source data for this figure are from HAP1 cells. Wilcoxon test, two-sided. n, number of genomic windows analyzed. (f) Number of 100 kb genomic windows belonging to different A/B subcompartments across ten concentric nuclear layers. In all the violin plots in the figure, each box spans from the 25th to the 75th percentile and whiskers extend from -1.5×IQR to +1.5×IQR from the closest quartile, where IQR is the inter-quartile range. Dots: outliers (data falling outside whiskers). All the source data for this figure are from HAP1 cells. Figure 7. Radial organization of A/B subcompartments along individual chromosomes. For each chromosome, the distributions of the GPSeq scores calculated in nonoverlapping 100 kb windows are shown separately for each subcompartment. Asterisks indicate the P-value (Wilcoxon test, two-sided) of the corresponding comparison between two different subcompartments (ns: P > 0.05; * P ≤ 0.05; ** P ≤ 0.01; *** P ≤ 0.001; **** P ≤ 0.0001). For exact P and n values see Supplementary Table 11. In all the boxplots, boxes span from the 25th to the 75th percentile and whiskers extend from -1.5×IQR to +1.5×IQR from the closest quartile, where IQR is the inter-quartile range. Dots: outliers (data falling outside whiskers).

Supplementary
All the source data for this figure are from HAP1 cells.  The upper triangle shows the inter-bead 3D distances calculated from 10,000 H(intra)

Supplementary
structures. The bottom triangle shows the KR-normalized Hi-C contact frequency matrix, with each element raised to the power of -0.25. The reported correlation coefficients are for 1 Mb resolution, while for simplicity the plot shows averaged values over 10 Mb genomic windows (points). (b) Correlation between the average inter-bead 3D distance in H(intra) structures and the KR-normalized Hi-C contact frequency, with each element raised to the power of -0.25.
Each dot represents a pair of 10 Mb non-overlapping genomic windows, each obtained by averaging 1 Mb non-overlapping bins. n = 47,531 genomic window pairs (points) were analyzed. (c) Distribution of the average distance from the modeled nuclear surface of 1 Mb beads in H(intra) structures, separately for each chromosome. In all the violin plots, each box spans from the 25th to the 75th percentile, whiskers extend from -1.5×IQR to +1.5×IQR from the closest quartile, where IQR is the inter-quartile range. Dots: outliers (data falling outside whiskers). chr9:22 and chr22:9 are the derivative chromosomes of the t(9;22)(q34;q11. All the source data for this figure are from HAP1 cells.

Supplementary Methods
Information about antibodies, cell lines, data and code availability, and statistics is provided in the Life Sciences Reporting Summary available online. None of the cell lines used is included in the list of misidentified cell lines; therefore we did not authenticate them. In the case of HAP1, we seeded the cells directly onto 22x22 mm coverslips placed in 6-well plates and grew them until they reached ~70% confluency in each well. In the case of GM06990, which grow in suspension, we first centrifuged them for 5 min at 300 g and re-suspended them in 1X PBS. We then dispensed the cell suspension onto 22x22 mm coverslips pre-coated with Poly-L-Lysine (Sigma, cat. no. P8920-100 ml) and placed inside a 6-well plate and incubated for 10 min at room temperature (RT). We fixed the cells in 0.4X PBS (Thermo Fisher Scientific, cat. no. AM9625)/4% paraformaldehyde (EMS, cat. no. 15710) for 10 min at RT, followed by quenching of unreacted paraformaldehyde in 1X PBS/125 mM glycine for 5 min at RT. Subsequently, we washed the cells three times, 5 min each, with 1X PBS/0.05% Triton X-100 at RT and permeabilized them in 1X PBS/0.5% Triton X-100 for 20 min at RT. Following overnight incubation in 1X PBS/20% glycerol at RT, we subjected the cells to four cycles of freeze-and-thaw in liquid nitrogen, and then washed them three times, 5 min each in 1X PBS/0.05% Triton X-100 at RT. Afterwards, we incubated the cells in 0.1 N HCl for 5 min at RT and quickly rinsed them twice in 1X PBS/0.05% Triton X-100 at RT. Lastly, we rinsed the cells in 2X SSC buffer (Thermo Fisher Scientific, cat. no. AM9763) and stored them in 2X SCC/0.05% NaN3 up to one month at 4 °C.

Preparation of cells for YFISH and
Preparation of YFISH and GPSeq adapters. The list of oligonucleotides (oligos) used to prepare YFISH and GPSeq adapters is available in Supplementary Table 1. We purchased individual oligos from Integrated DNA Technologies as desalted oligos at 100 mM concentration in nuclease-free water. UMIs were generated by random incorporation of the four standard dNTPs using the 'Machine mixing' option. We diluted 10 μl of each forward oligo in 90 μl final volume of a mix containing 10 μl of 10X Polynucleotide Kinase (PNK) buffer and 2 μl of T4 PNK (NEB, cat. no. M0201). We incubated the samples for 1 h at 37 °C, after which we added 10 μl of the corresponding antisense oligonucleotides, and incubated the samples for 5 min at 95 °C, followed by gradual cooling down to 25 °C over a period of 45 min (-1.55 °C/min) in a PCR thermocycler. Ready-to-use adapters can be stored at -20 °C up to one year. 3D DNA FISH. We designed and produced all the DNA FISH probes used to validate GPSeq with the iFISH pipeline, which we recently established1. The genomic coordinates of the probes are available in Supplementary Table 3. We performed DNA FISH using a modification of the original 3D DNA FISH protocol2, which we recently described3. Briefly, we first incubated the coverslips with the pre-hybridization buffer (PHB) containing 2X SSC/5Denhardt's solution/50 mM sodium phosphate buffer/1 mM EDTA/100 ng/ml ssDNA/50% formamide, pH 7.5-8.0, for 1 h at 37 °C. During this time, we prepared the first hybridization mix (HM-1) by mixing the single-locus probes (up to 6 single-locus probes together, typically using a final concentration of 3.2 nM per probe) at 1:9 vol./vol. ratio with 1.1X first hybridization buffer (HB-1) containing 2.2X SSC/5.5 Denhardt's solution/55 mM sodium phosphate buffer/1.1 mM EDTA/111 ng/ml ssDNA/55% formamide/11% dextran sulfate, pH 7.5-8.0. We removed the coverslips from PHB and placed each of them on 10 l of HM-1 on a microscope slide. We then sealed the coverslips with fixogum (MP Biomedical, cat. no. 11FIXO0125) and waited until the fixogum solidified. We performed DNA denaturation for 3 min at 75 °C on a heating block. Afterwards, we incubated the samples for 18 h at 37 °C. The next day, we washed the coverslips three times, 10 min each at 37 °C in 2X SSC/0.2% Tween, while shaking, followed by two washes, 7 min each at 60 °C in 0.2X SSC/0.2% Tween pre-warmed at 60 °C, inside a water bath, a brief wash in 4X SSC/0.2% Tween at RT, two brief washes in 2X SSC and one final short wash in 2X SSC/25% formamide. Next, we placed the coverslips on 100 l of the second hybridization mix (HM-2) containing the secondary fluorescently labeled oligos (one color per locus, up to 6 colors together), each at a final 20 nM concentration in 2X SSC/25% formamide/10% Dextran sulfate/1mg/ml E. coli tRNA/0.02% BSA, and incubated the samples for 3 h at 30 °C in a humidity chamber. Afterwards, we washed the coverslips for 1 h at 30 °C in 2X SSC/25% formamide, followed by 30 min at 30 °C in 0.1 ng/l Hoechst 33342 in 2X SSC/25% formamide. We imaged all the samples using wide-field microscopy as described below.

Visualization of DSBs by immunofluorescence (IF). We visualized DSBs in HAP1 cells
using an anti-phospho-Histone H2A.X (Ser139) antibody (Millipore, cat. no. 05-636). We performed IF by adapting a previously published protocol for nuclear proteins4. Briefly, we incubated the samples with the primary antibody diluted 1:200 vol./vol. in a blocking buffer consisting of 1X PBS/0.1% Tween 20/1% BSA, for 1 hour at room temperature (RT). For primary antibody detection, we incubated the samples with the secondary antibody diluted 1:500 vol./vol. in the blocking buffer, for 1 hour at RT. We imaged all the samples using widefield microscopy as described below.
We fixed HAP1 cells following the same protocol used for YFISH and GPSeq. We incubated different samples for 1, 5, 10, 15, 30 min or 1 h in the presence of the primary antibody diluted 1:500 vol./vol. in blocking buffer (1X PBS/0.1% Tween 20/1% BSA). Immediately after incubation, we post-fixed the samples in 1X PBS/4% PFA for 10 min at RT, followed by incubation with the secondary antibody diluted 1:500 vol./vol. in blocking buffer for 1 hour at RT. We imaged all the samples using wide-field microscopy as described below.
Wide-field epifluorescence microscopy. We used wide-field microscopy to image all the YFISH, 3D DNA FISH as well as IF samples. Briefly, we rinsed the samples twice in 2X SSC before mounting them with ProLong Diamond Antifade Mountant (Thermo Fisher Scientific, cat. no. P36965). We imaged all the samples using a 100X 1.45 NA objective mounted on a custom-built Eclipse Ti-E inverted microscope system (Nikon) controlled by the NIS Elements software (Nikon) and equipped with an iXON Ultra 888 ECCD camera (Andor Technology).
For YFISH and 3D DNA FISH, we acquired multiple image stacks per sample, each consisting of 40-110 focal planes spaced 0.2 or 0.3 µm apart. We imaged TetraSpeck Microspheres (0.1 m, fluorescent blue/green/orange and dark red, Thermo Fisher Scientific, cat. no. T7279) before or after each imaging session and used the images to correct for chromatic aberrations and shifts between channels using our in-house suite DOTTER. Microsystems). To block STED laser light in excess, we used a sharp notch filter in front of the ATTO647N detector. We acquired dual-color axial stacks sequentially, frame-by-frame at a scan speed of 400 lines per sec. We tuned the super-resolution STED pixel size based on the depletion power applied laterally and axially (80% laterally and 20% axially). Before the analysis, we deconvolved all the stacks with the Huygens Software (Scientific Volume Imaging).

YFISH image analysis.
We converted raw images from ND2 format (Nikon) to uncompressed TIFF format using the nd2_to_tiff script from the pygpseq Python3 package available on GitHub (https://github.com/ggirelli/pygpseq/). We deconvolved all the channels using the Huygens Professional v17.04 Software (Scientific Volume Imaging) with the following parameters: CMLE algorithm, null background, and signal-to-noise ratio (SNR) of 7, in 50 iterations. To perform 3D deconvolution, we estimated a theoretical point spread function using the same software, considering both our microscope setup and the optical configuration used to acquire the images. We then estimated the background level of both DNA and YFISH channels as the median of the background intensities, and subtracted these values from the original image (we set negative intensity voxels to zero). After deconvolution, we performed 3D segmentation of the nuclei stained with Hoechst 33342, in each field of view, using the tiff_auto3dseg script of the pygpseq package. Briefly, the script combines two binary masks generated with the threshold_otsu and threshold_local methods from the scikit-image.filters package5 with a logical AND operation. Then, the script discards objects touching the XY contour of the image, fills any holes in the masks, and performs a dilate-fill-erode operation.
To identify out-of-focus images, we used the tiff_findoof (v0.3.1) tool from the pygpseq package. The script discards the stacks where the peak of the gradient magnitude of the stack intensity over Z does not fall in a range of 50% of the stack around the mid slice. For each segmented nucleus, we estimated the volume, shape, surface, flattened size (in Z-projection), sum of intensity, average intensity, a shape descriptor, and the center of mass. Specifically, we ran the script using the following parameters: --compressed -rn -a 200 130 130 --an-type 3duy --nuclear-sel -t 10. In particular, the -u parameter allows extracting all the nuclear bounding boxes for further analysis. For each voxel, the pipeline produced all channel images, nuclear masks, and EDT-based matrices with absolute distance from the nuclear lamina and center (defined as the 1% most distant voxels from the lamina). Then, we used the pygpseq-scripts (v0.0.1, DOI: 10.5281/zenodo.3365634) suite of Python3 and R scripts, available at https://github.com/ggirelli/pygpseq-scripts to analyze the extracted nuclear bounding boxes.
First, we ran the extract_nuclear_features.py script to obtain the characteristics of all extracted nuclei, which we fed in turn to the select_nuclei.R script, with parameter -k 2, to select only nuclei from cells in the G1 phase of the cell cycle. To achieve this, we fitted a sum of Gaussians model to both the distribution of nuclear volume (in voxels) and the integral over the full nuclear volume of the DNA staining intensity. If the fitting failed, we instead fit a single Gaussian. We selected nuclei that fell in a range of ±k around the peak of the major Gaussian.
If also the single Gaussian fitting failed, we selected nuclei falling in a range of ±FWHM/2 around the major distribution peak, where FWHM is the full-width-half-maximum. We retained for further analysis only nuclei selected in such a manner. Then, we used the extract_nuclear_vx.py script to obtain the following measurements, for every single voxel of the selected nuclei: intensity from each channel, absolute distance from nuclear lamina, absolute distance from nuclear center, normalized distance from nuclear lamina (defined as the absolute distance from the nuclear lamina divided by the sum of the absolute distance from lamina and center). Next, we fed the voxel data to the extract_condition_profiles.py, extract_nuclear_profiles.py, and extract_nuclear_radii.py scripts. These scripts build radiality profiles from the nuclear periphery inwards for: (1) each condition (e.g., time of digestion) by pooling all the nuclear voxels (extract_condition_profiles.py); (2) each nucleus (extract_nuclear_profiles.py); and (3) single straight-line trajectories, going from the nuclear center of mass (CoM) to the nuclear surface (extract_nuclear_radii.py). Specifically, we ran the extract_nuclear_radii.py script with default parameters, to draw 200 straight trajectories and sample 100 points homogeneously in each of them, for 500 randomly picked nuclei for each condition. The 200 trajectories all depart from the nuclear CoM, are homogeneously spread in space using a spherical Fibonacci lattice, and terminate at their point of intersection with the nuclear surface, determined by a triangular meshing algorithm over the nuclear 3D mask. In all the three cases, we built radiality profiles by first dividing the nuclear radius into 200 bins, and then by assigning the intensity values of the voxels (or sampled points) of interest to the corresponding bins, based on the normalized distance from the nuclear lamina. Then, we calculated the median of each bin and reported it alongside the midpoint of each bin. Finally, we ran the extract_profile_descriptors.R script to retrieve information on the peak, inflection point, and contrast of the profile. Specifically, the script fits a 5th-degree polynomial curve to the profile data and uses the uniroot.all function from the rootSolve R package (v1.7) to identify the peak and inflection point of the profile, and the relative intensities. The profile contrast is then calculated as the ratio of the intensity at the peak over the intensity at the inflection point.
3D DNA FISH image analysis. We deconvolved the DNA staining channel and performed 3D segmentation as described above for YFISH. To identify fluorescence dots corresponding to the detected DNA loci in 3D, and to correct for chromatic aberrations, we used our in-house suite DOTTER. In order to calculate the radial position of the selected FISH dots, we passed the output of DOTTER to the gpseq_fromfish (v7.0.2) script from the pygpseq package. For each FISH dot, we computed its distance from the nuclear edge (contour of segmented DNA stain) as well as from the nuclear center. Then, we normalized the distances dot-wise, by dividing them by the sum of the corresponding distances from the nuclear edge and center. We defined the nuclear center as the set of voxels in the top percentile of the distribution of distances from the nuclear edge. For each voxel, we calculated its 3D distance from the nuclear edge using an anisotropic 3D Euclidean transform on the nuclear mask using the ndimage.morphology.distance_transform_edt script from scipy (v1.1.0). We then normalized the distance of each voxel from the nuclear edge, using the same procedure as for the DNA FISH dots. This analysis was implemented as a snakemake flow6, also available on GitHub (https://github.com/ggirelli/iFISH-singleLocus-analysis).
Pre-processing of GPSeq sequencing data. We demultiplexed raw sequencing data based on the RA5 indexes, either using the BaseSpace Sequence Hub cloud service of Illumina, or manually with bcl2fastq (v2.18). We quality-checked the generated fastq files with fastqc (v0.11.4). We selected the reads containing the full prefix (UMI_barcode_restriction site) using scan_for_matches. Afterwards, we trimmed the reads to remove the prefix (including the restriction site) and aligned the remaining part against the human reference genome (Grch37/hg19 GCA_000001405.1) using bwa-mem. We retained primary alignments with a mapping quality equal to or higher than 30, while we discarded unmapped reads, chimeric reads, and reads mapped to chrY that is not present in HAP1 cells. We retained only the reads whose 5' end is less than 20 bp away from the position of a HindIII or MboI recognition site (RS) in the reference human genome, depending on which enzyme was used. The reason why we do not strictly require all the reads to align exactly to the position of RS is that the T7 polymerase and reverse transcriptase used during the GPSeq library preparation are prone to occasionally skip some bases, leading to the resulting reads to align slightly downstream of the RS. Afterwards, we recovered the UMI sequence of each aligned read and we filtered the reads based on the quality of the UMI sequence using an approach similar to the one used by the fastq_quality_filter tool (http://hannonlab.cshl.edu/fastx_toolkit/index.html). We deduplicated UMI sequences mapped to the same restriction site and generated a BED file containing the genomic coordinate and number of de-duplicated UMIs associated with each restriction site. We performed all the above steps using gpseq-seq-gg (v2.0.2), a bash/Python/R custom-designed pipeline that is available on GitHub (https://github.com/ggirelli/gpseq-seqgg). A newer version of the pipeline (v2.0.3) with improvements to the efficiency of the UMI de-duplication script is available at DOI: 10.5281/zenodo.3264757. This version is optimized to deal with datasets generated with 4-base cutters such as MboI used in this study. Finally, we corrected the BED files generated by gpseq-seq-gg for the presence of the t(9;22) translocation using an ad hoc Python script available at http://github.com/ggirelli/bed-fix-chromrearrangement (v0.0.1, DOI: 10.5281/zenodo.3365906), using the following parameters: -1 chr9:133681295 -2 chr22:23632359.

Repeat masking
We manually annotated a list of masked regions (Supplementary Table 10), using the RepeatMasker track available from the UCSC Table browser for Grch37/hg19. The list includes peri-centromeric and peri-telomeric regions (up to 30 Mb) characterized by abnormally high GPSeq score across all experiments. The largest masked regions correspond to heterochromatic regions of variable length and to the p-arms of acrocentric chromosomes. Moreover, we also included two small regions on the q-arm of chrX and chr5, which yielded abnormally high read counts. When calculating the GPSeq score genome-wide, we masked out any genomic window overlapping with the masked regions using the -M parameter of the gpseqc_estimate script, from the gpseqc Python3 package available on GitHub (https://github.com/ggirelli/gpseqc/).

Effect of resolution and sequencing depth.
To investigate the effect of the resolution (i.e., genomic window size) on the GPSeq score, we calculated the median number of reads, median number of enzyme recognition sites, median number of reads per recognition site, and standard deviation of the number of reads per recognition site, inside genomic windows of 5, 10, 15, 20, 25, 50, 100, 250, 500 kb, and 1, 5, 10 Mb for the longest digestion times (2 hours for HindIII, 30 min for MboI). To assess the effect of the sequencing depth, we iteratively removed 5% of the original de-duplicated UMIs, using an R script and assigning to each de-duplicated UMI the same probability of being removed at any cycle. First, we calculated the same genomic window statistics (window size of 1 Mb, with a step of 100 kb) as described above. Then, we calculated the GPSeq score from the sub-sampled BED files, both in a genome-wide fashion (1 Mb overlapping genomic windows with a step of 100 kb and 100 kb genomic windows) as well as using genomic windows (1 Mb and 100 kb) centered on the midpoint of the DNA FISH probes shown in Supplementary Fig. 1a. Following this, we calculated the Pearson's correlation coefficient between the GPSeq score at different simulated sequencing depths and the median distance from the nuclear lamina measured by DNA FISH. Lastly, we calculated the Spearman's correlation coefficient between the genome-wide GPSeq score of each experiment at different simulated sequencing depths, and all the other experiments at full sequencing depth.

Conversion of GPSeq scores to physical distances.
To convert the GPSeq score into absolute physical distances, we first retrieved the average HAP1 nuclear volume (in voxels) from the nuclear.summaries.csv file automatically generated by the pygpseq pipeline and converted it into cubic nanometers (by multiplying it by our voxel size: 130×130×200 nm). Then, we multiplied the log2 GPSeq score by the absolute average nuclear radius to obtain a transformed GPSeq score (in nm).

GC-content.
We obtained the fraction of guanines and cytosines (GC) in each genomic window (1 Mb and 100 kb) by applying the letterFrequency function from the Biostrings R package (https://www.bioconductor.org/packages/release/bioc/html/Biostrings.html) to the genome sequence provided by the BSgenome.Hsapiens.UCSC.hg19 R package (http://bioconductor.org/packages/release/data/annotation/html/BSgenome.Hsapiens.UCSC.h g19.html). We excluded genomic windows overlapping masked regions listed in Supplementary Table 10 or for which it was not possible to obtain a GPSeq score.

ATAC-seq.
We obtained ATAC-seq data for the HAP1 cell line from the GEO repository GSE111047 (Supplementary Table 5). We used the Parker Lab's ATAC-seq processing pipeline (https://github.com/ParkerLab/ATACseq-Snakemake) to obtain the genome-wide accessibility signal in bedGraph format. To directly compare the GPSeq scores and genome accessibility, we binned the ATAC-seq bedGraph files using the bioTrackBinner R script GPSeq score profiles. We visualized GPSeq score profiles along chromosome ideograms by using the ggkaryo2 R package (v0.0.3) available at http://github.com/ggirelli/ggkaryo2 (DOI: 10.5281/zenodo.3358887). We generated a circular chromosome-wise GPSeq score representation with the same package by reporting the log2 of the GPSeq score on the y-axis, the genomic chromosome coordinate (in bp) on the x-axis arranged circularly using the coord_polar function, and the connecting consecutive data points by segments colored based on the GPSeq score. We marked in red the peri-centromeric regions listed in Supplementary   Table 10.
Pizza-plots. We generated scatter 'pizza-plots' using the ggplot2 R package by setting the yaxis to follow the log2 GPSeq score and reporting the chromosomes on the x-axis arranged circularly using the coord_polar function. We then jittered the points at the same radial position (i.e., same GPSeq score) using the geom_jitter function and coloring the data points based on orthogonal data track scores. We generated tiled pizza-plots by dividing the log2 GPSeq score into 10 equally sized bins from 0 to 1 and assigning data points to each bin in a chromosomewise fashion. Afterwards, we plotted each bin as in scatter pizza-plots but replacing the geom_jitter by the geom_tile function and normalizing the values to have the tiles of each layer (log2 GPSeq score bin) or slice (chromosome) summing up to 100%. To visualize the radial distribution of all chromosomes at once, we first filled a disk with circular sectors, one per chromosome. Then, for each of five layers of equal thickness (not shown) we updated the number of pixels per chromosome to be proportional to the number of genomic sites in that layer in our GPSeq data. The order of chromosomes was arbitrary.

DNA methylation and histone marks. We obtained DNA methylation and histone mark
ChIP-seq and ChIPmentation data from ENCODE (Supplementary Table 5). For DNA methylation, we retrieved the corresponding BED file, while for ChIP-seq datasets we retrieved the corresponding bigWig file of the merged replicates signal ('fold change over control') mapped to the Grch37/hg19 genome assembly. ChIP-seq and ChIPmentation datasets for which signal tracks were not available were processed with a custom pipeline. Briefly, we removed the adapter sequences using TrimGalore (v0.4.4_dev) and discarded reads shorter than 20 nt.
We then aligned the reads to the UCSC hg19/GRCh37 genome assembly using bwa-mem (v0.7.17-r1188) with default options. We filtered out the reads that failed to align and those with MAPQ ≤ 30, and discarded PCR duplicates using the Picard MarkDuplicates (v2.18.11) module. Finally, we generated genome coverage tracks in bigWig format using the bamCoverage module from deeptools (v3.2.1) with --binSize 50. We used these bigWig and BED tracks to calculate the signal at 1 Mb or 100 kb resolution, by applying bioTrackBinner to them. Afterwards, we combined ChIP-seq and DNA methylation binned tracks with the GPSeq score calculated at 1 Mb or 100 kb resolution and averaged the tracks signal across replicates. Finally, we built radiality profiles following the same procedure as described above for DamID.
Gene expression and RNA Pol II occupancy. We obtained HAP1 RNA-seq raw files from the GEO repository GSE95015 and RNA Pol II ChIP-seq data from GSE107599 (Supplementary Table 5). We removed the adapter sequences using TrimGalore (v0. 4.4_dev) and discarded reads shorter than 16 nt. We mapped the reads using Bowtie2 (v2.3.4.1; parameters: --sensitive-local) to a list of human rRNAs and tRNAs retrieved from NCBI (https://www.ncbi.nlm.nih.gov/). We then mapped the reads that failed to align in the previous step to the UCSC hg19/GRCh37 genome assembly by using the Gencode (v27) gene annotation as a reference and the STAR software (v2.6.0c; parameters: --twopassMode Basic --alignSJoverhangMin 8 --alignSJDBoverhangMin 1 --sjdbScore 1 --alignIntronMin 20 --alignIntronMax 1000000 --alignMatesGapMax 100000 --outFilterMultimapNmax 1 --outFilterMismatchNmax 999 --outFilterMismatchNoverReadLmax 0.04). We discarded PCR duplicates using the Picard MarkDuplicates module (v2.18.11). To quantify the RNA-seq signal at the gene level, we used the QoRTs QC module (v1.3.0; parameters: --minMAPQ 255) and the Gencode gene annotation (v27) as a reference. We normalized gene counts to TPM (Transcript Per Million tags). To quantify the RNA-seq signal at 1 Mb or 100 kb resolution, we used the bamCount function from the bamsignals package (v1.12.1) in R. For RNA Pol II ChIPseq, we performed raw data processing and genomic tracks generation in the same way as described for histone marks. Finally, we built radiality profiles following the same procedure as described above for DamID.
Gene set radial enrichment. We obtained the list of homeobox genes from the Homeobox Database8, the list of housekeeping genes from Eisenberg and Levanon9, and the list of hallmark gene sets from MSigDB (http://software.broadinstitute.org/gsea/index.jsp). In one approach, we generated gene set enrichment plots by assigning to each gene the GPSeq score and layer of the 100 kb genomic window that comprised its TSS. Then, for each layer, we calculated the fraction of gene set entries over the fraction of all genes. In the other approach, we compared the GPSeq score distributions and calculated the statistical significance of the shift of each gene set from the background ("all") gene distribution using a two-sided Wilcoxon rank sum test.

TFBS radial distribution. We retrieved the list of transcription factor binding sites (TFBS)
predicted in the UCSC hg19/GRCh37 genome assembly from JASPAR202010. For each TFBS, we counted its occurrence within 1Mb and 100 kb genomic windows. We then calculated the GC-content of each TFBS as the average probability of each base in the motif being either a G or a C. Lastly, we used the ggplot2 script to plot the Pearson's correlation coefficient between the number of TFBS per genomic window and the GPSeq score of the same window, against the GC-content of each TFBS.
Linear models of the GPSeq score. We retrieved chromosome sizes ( ) based on the Grch37/hg19 genome assembly available from the UCSC Table Browser We replaced with different epigenetic mark and replication time data. A complete list of the estimated parameters, β, R2, and model P-values is available in Supplementary 9. We calculated the prediction error (PE) as the square root of the mean square prediction error. Table 5).

Replication timing. We obtained Repli-seq data from ENCODE (Supplementary
We retrieved the bigWig files of the signal in each of the cell cycle phases (G1, S1, S2, S3, S4, G2) and of their wavelet transform. We binned the bigWig files with bioTrackBinner and combined them with the GPSeq score calculated at the same resolution. We used the cell cycle phase-specific signal to color data points in scatter pizza-plots.

Hi-C contacts, A/B compartments and subcompartments.
We downloaded HAP1 Hi-C raw data from the 4D nucleome project web portal (https://www.4dnucleome.org/) (Supplementary Table 5). We processed the data using Juicer11 by aligning the data against the human reference genome (Grch37/hg19 GCA_000001405.1) using bwa-mem. We classified the Hi-C contacts into four groups based on their position within the genome-wide contact frequency matrix. We defined as 'diagonal' and 'near-diagonal' the Hi-C contacts occurring between loci lying inside the same or to linearly adjacent genomic windows, respectively. We classified the remaining contacts as 'intra' or 'inter', depending on whether the genomic loci were located on the same or on different chromosomes. Pearson's and Spearman's correlation coefficients between the full Hi-C distance matrix and the averaged 3D spatial bead distance matrix.

Radial distribution of A/B compartments in chromflock structures.
To assess the polarity of A/B compartments in chromflock structures built using GPSeq and Hi-C data, we calculated the median radial position based on the distance of each bead from the nuclear surface, for all beads belonging to either A or B compartments, and then calculated the difference in the median radial position between A and B compartments for every structure.

Orientation of A/B subcompartments in chromflock structures.
For simplicity, let us first assume that we are interested in studying the orientation of genomic compartments with two different properties, 1 and 2. For each chromosome in each chromflock structure, we first identify the beads with these properties: 1 = { , ∈ 1} and 2 = { , ∈ 2}. We also use the set of beads with any of these two properties, i.e., 1 = 1⋃ 2. We then define the normalized vector, , from the center of mass of 2, ( 2), to ( 1), and name the normalized vector from ( ) to 0 as . We define the orientation, as: i.e., = 1 if beads with property 1 are comprised between ( 2) and 0, and = −1 if ( 2) is between ( 1) and 0. If a chromosome were completely polarized with respect to 1 and 2, we would be able to draw a plane that completely separates 1 and 2. We defined polarity, as the quotient of beads that can be correctly classified as 1 or 2 by a linear classifier.

Radial distribution of gene fusions and mingling.
We retrieved the cancer-related gene fusions from TCGA (https://www.tumorfusions.org/). We overlapped chromflock structure beads and gene fusions using the foverlaps function from the data.table R package, and labeled as "Fusions" the beads overlapping with at least one gene fusion junction in any cancer type and as "Controls" all the remaining beads.

Radial distribution of DSBs.
To assess the radial distribution of endogenous DSBs, we used a published dataset (GSM3444988), which we previously obtained from K562 cells16, by applying a modification of our BLISS method17 (Supplementary Table 5). We computed the number of unique DSB ends per 100 kb genomic window. We defined as 'genic' the genomic windows with their center overlapping on either strand with any of the GENCODE v19 reference gene annotations and considered the remaining windows as 'intergenic'. We normalized the BLISS signal by dividing it by the length of the 'genic' or 'intergenic' portion of the genome in each of ten concentric nuclear layers. To analyze the radial distribution of DSBs in different A/B subcompartments, we assigned the genomic windows to different A/B subcompartments as described above for genic vs. 'intergenic regions, and then averaged the BLISS signal across all the windows in each layer.
Intermingling in chromflock structures. To assess mingling between heterologous chromosomes in the chromflock structures built using both GPSeq and Hi-C data (1 Mb resolution), we first extracted all the neighboring beads, found in a sphere four times the bead radius, for each bead, belonging to a given chromosome, . For every chromosome with beads in , except , we constructed the convex hull, using the chromosome's beads in .
We defined as a mingling event the case in which the centroid of falls inside . We then defined the mingling frequency as the fraction of chromflock structures in a given set, in which a bead mingles.

Supplementary Tables
Supplementary Table 1

Supplementary Note 1. Calculation of the GPSeq score
To convert the sequencing data generated by GPSeq into radiality maps, we have developed a computational framework that allows computing a centrality score, namely the GPSeq score.
Below we describe in detail the centrality scores that we have tested and how the GPSeq score is calculated.

Introduction
A GPSeq experiment consists of digestion times (herein named conditions), which yield 1 , 2 , … , sets of processed sequencing reads (one set for every digestion time). Each condition is characterized by a specific number of de-duplicated reads, ( ), distributed along the genomic recognition sites, , of the restriction enzyme used in the experiment. Let us denote ( , ) as the number of de-duplicated reads mapped to a given recognition site, and ( ) as the number of recognition sites for a given restriction enzyme considered in the condition (for more details on how to consider recognition sites, please refer to the section "Recognition site domain" below). Thus, we have: Let us then represent a genomic window, by the genomic coordinates of the first, and the last, bases included in the window on a specific chromosome, : Each window is characterized by a specific number of de-duplicated reads, ( , ) and of recognition sites, ( , ) for a given restriction enzyme. Each recognition site can be considered as a small genomic window. Thus, the number of reads and recognition sites can be written as: Considering only reads that are mapped to the recognition sites of the restriction enzyme in use, we can define the restriction (digestion) probability as: While the of the same window can be compared across different conditions, it is difficult to use this measure to compare different windows in the same condition, since different genomic windows can contain a different number of recognition sites, which directly affects the restriction probability. In other words, genomic windows with a higher number of recognition sites are expected to show a higher compared to windows with a lower number of recognition sites. To take this into account, we define the restriction probability as the average probability, ( , ) that a single recognition site in a given genomic window is cut: This probability takes ( , ) into consideration and can be compared both across different genomic windows in the same condition as well as across conditions. At the same time, each genomic window is characterized by a different distribution of sequencing reads along its recognition sites. Thus, we can define the mean of the number of reads per site, and the variance, as:

Definition of centrality scores
Initially, we computed the correlation between calculated for genomic windows of 1 Mb and 100 kb centered on the 68 DNA FISH probes shown in Supplementary Fig. 1a, and the median distance from the nuclear lamina measured using the same probes, which showed a reproducible increase in the correlation with increasing digestion times (Supplementary Note 1 Fig. 1a-e). This observation prompted us to ponder on possible ways in which we could combine the different restriction conditions into a single centrality score. Specifically, we formulated two groups of candidate centrality estimates, based either on the read variance, or on the restriction probability, .

Probability-based centrality scores
We first hypothesized that centrally located genomic windows should have higher in samples digested for a longer time, compared to shortly digested samples, while peripheral genomic windows should show the opposite behavior. Therefore, we combined into a probabilitybased centrality measure ( ) by comparing each condition with the immediately preceding one: Additionally, we defined a centrality measure based on the cumulative restriction probability:

Variability-based centrality scores
We then reasoned that peripheral windows should be fully digested at any restriction condition, thus showing a stable read variance, . On the other hand, internal windows would show increasing numbers of reads, and at the same time a higher read variance, . Hence, we defined a variability-based centrality score, ( ) by comparing the last and first condition: Furthermore, we defined similar centrality measures based on the coefficient of variation, ( ) and of dispersion, ( ):

Combination of restriction conditions
In the above formulations, we combined different restriction conditions by comparing each digestion time with the immediately preceding one: we refer to this as adjacent approach. At the same time, we also tested two additional ways of combining different conditions: (i) the fixed approach, which compares each condition with the initial one, ( and 1 ); and (ii) the two-points approach, which only combines the last and initial condition, ( and 1 ). The different formulations, based on these three approaches, are summarized in Supplementary Note 1 Table 1. It is important to note how the adjacent and two-points approaches lead to the same formulation for variability-based centrality estimates, as they are based on differences rather than on ratios.

Fixed approach Two-points approach Adjacent approach
Probability-  Table 1. Different centrality formulations, either probability-or variability-based, using different approaches to combine different restriction conditions of the same GPSeq experiment.
One of the most important effects of combining multiple conditions using one of the approaches described above is that this intrinsically normalizes scores across conditions. In other words, if a genomic window, is characterized by a bias towards higher or lower read counts (e.g., due to more or less chromatin accessibility or to biased restriction) such bias will be corrected by the normalization itself, provided that it has a multiplicative (i.e., linear) effect on the number of reads, or on the average number of reads, or on the read variance, .

Recognition site domain
Given that a GPSeq experiment comprises multiple conditions and that typically only a fraction of all the restriction sites that have been cut are sequenced in the same run, the sets of recognition sites identified in different conditions only partially overlap. In other words, even if cut, a recognition site might yield sequencing reads only in one condition, but not in another one. We can envision four different ways in which this can be taken into consideration, by defining a different domain of recognition sites to be used in the above centrality formulations.
Let us first define ( ) as the set of recognition sites, that are associated with at least one aligned read in a given condition: The universal approach considers all the genomic recognition sites, independently of whether any read was mapped to them. In this case, the number of sites in a given genomic window is: While this approach has the advantage of using the same domain for all the conditions, it also tends to include a large number of empty recognition sites (i.e., sites with no associated mapped read), which skews the centrality calculation.
The union approach considers all the genomic recognition sites that have at least one mapped read associated to them in at least one condition: As for the universal approach, also in this case the same domain is used for all conditions, while decreasing the number of empty recognition sites included in the calculation.
The intersection approach considers all the genomic recognition sites that have at least one mapped read associated to them in all the conditions: ( , ) = |{ | ∈ ⋀ ( , ) ≠ 0, = 1,2, … , }| This domain can be readily applied to all conditions, but it has the tendency to discard a large fraction of recognition sites depending on the initial overlap between the different conditions.
Finally, the separate approach uses a separate recognition site domain for each condition, defined as the set of recognition sites with at least one mapped read (see (18)): While this approach defines a different recognition site domain for each condition, it also discards all empty sites, avoiding any skewing effect that might be caused by the sub-sampling steps.

Definition of GPSeq score
To identify the best centrality score, we first used the above formulations to compute the scores for 1 Mb and 100 kb windows centered on the 68 DNA FISH probes shown in Supplementary   Fig. 1a, and then compared the scores with the median distance from the nuclear lamina measured with the same probes. Importantly, we rescaled the scores to be able to compare them across different experiments (see "Score rescaling" section below). As shown in Supplementary Fig. 1b, the rescaled probability-based formulation outperformed all the others, both at 1 Mb and at 100 kb resolution, yielding higher correlations at almost all resolutions tested. We then compared different recognition site domain approaches described above. The separate and intersection approaches yielded a higher correlation with DNA FISH in comparison to the universal and union approaches (Supplementary Note 1 Fig. 2).
Furthermore, we calculated the genome-wide rescaled score, for 1 Mb (overlapping with 100 kb steps) and 100 kb windows, using different approaches. Importantly, when calculating the score in a genome-wide fashion, we discarded the windows overlapping with a manually annotate BED file of peri-centromeric and peri-telomeric regions (Supplementary  . 3). Based on these results, we defined the rescaled centrality score calculated using the separate approach as the GPSeq score and used it throughout the manuscript. We note that the GPSeq score is free from biases affecting individual genomic windows, as such biases are inherently corrected when the windows are compared across different conditions during the GPSeq score calculation.

GPSeq score rescaling
To rescale the GPSeq score in order to make it comparable across different experiments and resolutions, we first normalize the scores to have their values ranging between 0 to 1 (excluding outliers, defined as data points with a distance from the closest quartile (1st or 3rd) equal to or greater than 1.5 times the interquartile range). Then, to avoid negative values, we raise 2 to the power of the score, shifting non-outliers scores to the [1,2] interval, using the default values --score-outliers and --score-outlier-limit of the gpseqc_estimate script described in the

Supplementary Note 2. Considerations on sequencing depth and genomic resolution
Here, we aim at recapitulating the interconnected effects that the GPSeq experiment design, sequencing depth, and genomic resolution have on the GPSeq score. We first discuss the crucial steps during the design and execution of a GPSeq experiment. We then consider the effects of sequencing depth and resolution on the distribution of GPSeq reads and their statistics and more generally on the GPSeq score.

Choosing the restriction enzyme
The first crucial step when designing a GPSeq experiment is choosing the right restriction enzyme. In this study, we have used two different restriction enzymes (the 4-base cutter MboI and the 6-base cutter HindIII), which have a different frequency and genomic distribution of recognition sites (RS) (Supplementary Note 2 Fig. 1a, b). Using a 4-base cutter implies that a larger number of genomic loci can potentially yield fragments to be sequenced. Thus, using a 4-base cutter requires a higher sequencing depth compared to using a 6-base cutter. This also implies that the density of RS will be higher in any genomic window, potentially allowing to achieve a higher resolution in radiality maps. In general, when aiming at a resolution higher than 100 kb (e.g., 10 kb), we recommend using a 4-base cutter, which requires deeper and more costly sequencing. When a lower resolution is still acceptable (e.g., 250 kb), we recommend using a 6-base cutter, which allows for a more shallow and cost-effective sequencing. The actual distribution of the RS also plays an important role. First and foremost, it is essential for the RS of the chosen restriction enzyme to be distributed as homogeneously as possible along the genome, with the least possible number of genomic windows devoid of RS ("RS deserts").
An RS desert would not be able to generate any fragment to be sequenced, rendering impossible the calculation of the GPSeq score for that region. The number of RS available in any genomic window is already accounted for by the GPSeq score formulation, as described in

Choosing the number of restriction times
In this study, we have used a specific set of digestion times, which differ in between the two enzymes that we used. Our choice of the digestion times was directed by the empirical observation of how the restriction enzyme diffuses throughout the nucleus, as revealed by YFISH. In the case of MboI, we observed a clear fluorescent band at the nuclear periphery already after only one minute of incubation in the presence of the enzyme, whereas in the case of HindIII this required a longer time (10 min). During the course of GPSeq development, we noticed that the YFISH signal pattern no longer changed beyond a certain duration of incubation with a given restriction enzyme (2 hours for HindIII; 30 min for MboI). In fact, when we compared the GPSeq score calculated using all the digestion times up to 2 hours with the GPSeq score obtained using all the times up to 6 hours, the two scores did not differ in any significant manner even when substituting the 2-hour sample with the 6-hour sample. We conclude that, as long as the longest digestion yields a YFISH signal homogenously spread throughout the whole nuclear volume, the precise timing of it is not critical. We strongly recommend to always empirically determine the set of digestion times for a given enzyme, using YFISH prior to any new GPSeq experiment. This is particularly important when applying GPSeq to a new cell line that has not been assessed by YFISH before. The number of incubation times in between the two extremes depends on how many clearly distinguishable YFISH signal patterns one is able to observe for a given cell type (the patterns should vary from a very thin fluorescence layer confined at the nuclear periphery to a homogenous staining throughout the nucleus). However, we note that, adding more times in between the shortest and longest incubations, does not substantially improve the correlation between radiality estimates by GPSeq and DNA FISH measurements (data not shown). Therefore, as a rule of thumb we suggest using 4-5 incubation times in total. We note however that this might need to be adjusted when aiming to achieve higher resolutions than in this study.

Amplification steps during library preparation
The GPSeq protocol includes two amplification steps: in vitro transcription (IVT) and polymerase chain reaction (PCR) (see step-by-step protocol at Protocol Exchange, DOI: 10.21203/rs.3.pex-570/v1). IVT allows for linear amplification of the genomic fragments of interest, while PCR is used to attach Illumina sequencing adapters and indexes to the IVTamplified fragments. At the same time, PCR exponentially amplifies DNA sequences in a potentially biased manner. In other words, PCR might lead to the over-representation of certain fragments based on their sequence. We minimize this issue by performing a low number of PCR cycles during the GPSeq library preparation (typically, 9 PCR cycles for 50-300 ng of IVT input, see Supplementary Table 2). To maximize the complexity of the final sequencing libraries, we also recommend splitting the PCR reactions into multiple tubes (e.g., 4-5 25 l PCR reactions for 300 ng of IVT input in the case of HindIII or 50 ng IVT input in the case of MboI).

Effects of sequencing depth on the GPSeq score
Once the proper restriction enzyme has been chosen and GPSeq libraries have been prepared, they need to be sequenced. The choice of the sequencing platform is crucial as it defines the maximum achievable sequencing depth. In this study, we sequenced all our GPSeq libraries on the NextSeq 550 sequencing platform from Illumina, which in single-end mode typically yields 400-500 million reads (75 nt). To investigate the effect of different sequencing depths on GPSeq, we simulated different sequencing depths by iteratively removing 5% of the original reads (Supplementary Methods). The depth of sequencing will first and foremost affect the number of reads and their overall statistics. We found that the number of reads and the number of cut RS per 1 Mb genomic window scaled linearly with the sequencing depth, both in the case of a 6-base cutter (Exp.1) and 4-base cutter (Exp.3) (Supplementary Note 2 Fig. 2a, b). At the same time, the read distribution across the RS (in terms of median and standard deviation) remained largely unaffected, until reaching a depth so shallow to hinder any type of analysis (e.g., a depth lower than 20% of the original number of reads, see Supplementary Note 2 Fig.   2a, b).
The effect of low sequencing depths on the reads will, in turn, affect the GPSeq score calculation. We first investigated how a shallower depth affects the correlation between the distance from the nuclear lamina measured by DNA FISH and the GPSeq score calculated for 1 Mb windows centered around the DNA FISH probes. The correlation remained stable until the depth was lowered to 25% and 50% of the original for HindIII and MboI, respectively (Supplementary Note 2 Fig. 2c, d). At sequencing depths above the aforementioned values, the correlation remained stable, indicating that the achieved depth was well over saturation at a resolution of 1 Mb.
Next, we focused on the effect of lower sequencing depths on the GPSeq score, calculated in a genome-wide fashion using genomic windows of 1 Mb (overlapping with a step of 100 kb) or 100 kb. To this aim, we investigated how a reference GPSeq experiment at decreasing sequencing depth correlates to the four other experiments, individually and at full depth. The results showed a consistent decrease in the inter-experiment correlation when decreasing the sequencing depth. In the case of HindIII experiments, the correlation appeared to be stable at depths above 50% and 75% of the original at 1 Mb and 100 kb resolution, respectively (Supplementary Note 2 Fig. 2e, f). On the other hand, MboI experiments did not show a stable correlation at any depth at 100 kb resolution and were correlated only above 75% depth at 1 Mb resolution (Supplementary Note 2 Fig. 2g, h). This is due to the fact that MboI libraries have a higher molecular complexity requiring a higher sequencing depth. This was further confirmed by inspecting the read count profiles corresponding to individual digestion times, which showed that the inter-experiment reproducibility was lower for MboI compared to HindIII, especially in the case of shorter digestion times (Supplementary Note 2 Fig. 2i, j).
These results underscore the need for a higher sequencing depth when using higher-frequency cutters or when aiming for higher genomic resolutions.

Effects of genomic resolution on the GPSeq score
As discussed above, depending on the desired genomic resolution, it is crucial to choose an appropriate enzyme and sequencing platform to achieve the necessary sequencing depth. In the case of HindIII and MboI, the distribution of read counts along the RS (in terms of average and standard deviation) is largely unaffected by the genomic resolution (Supplementary Note 2 Fig. 2k, l). We have investigated the effect of higher genomic resolutions on the correlation between the GPSeq score and the distance from the nuclear lamina measured by DNA FISH (Fig. 2f). In the case of HindIII (Exp.1-2), the correlation drastically decreased for resolutions higher than 100 kb. On the other hand, the correlation decreased less steeply in the case of MboI (Exp.3), with a Pearson's correlation coefficient above 0.75 even at 10 kb resolution in the case of one experiment (Exp.4) with higher sequencing depth.

Conclusions
The desired genomic resolution, choice of enzyme, and achievable depth (which depends on the available sequencing platform) are critical aspects that should be considered whenever designing a new GPSeq experiment. The higher the desired resolution, the higher the cutting frequency of the restriction enzyme should be, thus requiring a higher sequencing depth. The accuracy of the GPSeq scoredetermined by comparing it to direct measurements of radial distances by 3D DNA FISHdepends on the frequency along the genome of the RS of the restriction enzyme in use. On the other hand, the precision of the GPSeq scoreassessed by evaluating its reproducibility across multiple experimentsdepends on the coverage of the RS of the restriction enzyme in use. This means that, given a sufficient sequencing depth, a more frequent cutter will allow achieving a more accurate radiality estimate at the same genomic resolution. As a rule of thumb, a 6-base cutter such as HindIII is sufficient to obtain reproducible radiality maps at 100 kb resolution, using 5-6 different restriction times and