Methods for single-cell hi-c regulatory scale chromatin band analysis

Abstract

The present application relates to biological information data processing technology, and is a method for single-cell Hi-C regulatory scale chromatin band analysis, comprising the following steps: preprocessing single-cell Hi-C data to generate pseudo-bulk Hi-C data; normalizing the pseudo-bulk Hi-C data to extract the Hi-C contact matrix of each chromosome; identifying and detecting pseudo-bulk bands from the Hi-C contact matrix; projecting the pseudo-bulk bands onto the original single-cell Hi-C data to obtain single-cell bands; and quantitatively analyzing the single-cell bands in the original single-cell Hi-C data. The present application can stably identify regulatory scale chromatin bands with clear endpoints and directionality, and output statistically significant band sets.

Description

Technical Field

[0001] This invention relates to bioinformatics data processing technology, and in particular to a method for analyzing chromatin bands at the Hi-C regulatory scale in single cells. Background Technology

[0002] Genomic DNA is not in a simple linear state in the cell nucleus. Instead, it exists as chromatin with a specific high-order spatial conformation after being highly folded and condensed, and is stored in the cell nucleus as a carrier of genetic information. Genomics is the collective characterization and quantification of all genes in an organism, studying the structure, function, evolution, and location of the genome, and analyzing the interrelationships between them and their impact on the organism.

[0003] Chromatin structural bands are directional, linear, high-interaction-frequency submatrices in Hi-C (High-throughput Chromosome Conformation Capture) data, extending along continuous genomic regions and typically indicating unilateral loop extrusion processes. These bands are closely associated with extended regulatory interactions that connect promoters to distal cis-regulatory elements (such as enhancers) and are widely distributed across different genomic loci. Variations in band frequency, intensity, or directionality are often accompanied by changes in chromatin folding patterns and gene regulatory states, including developmental, neurogenesis, and tumorigenesis processes. Bands constitute an independent level in the three-dimensional organization of the genome and complement structures such as chromatin compartments, topologically associated domains (TADs), and chromatin loops. Single-cell Hi-C data offers unprecedented opportunities to characterize band-related structural variability across individual cells. However, the extreme sparsity of single-cell Hi-C data poses a major challenge to robust detection of bands at regulatory resolution.

[0004] In the detection of chromatin structural bands, various computational methods based on bulk Hi-C data have been proposed. For example, some studies have proposed identifying chromatin band anchor points by constructing signals based on edge differences and comparing the contact frequencies on both sides of the band boundaries. Furthermore, recent studies have attempted to alleviate the extreme sparsity of single-cell Hi-C data by aggregating it into a pseudo-bulk representation, and then directly applying existing bulk band detection methods for analysis. Simultaneously, enhancement methods for single-cell Hi-C data have been proposed, often processing the contact matrix at relatively coarse resolutions such as 40-50kb, and then using the enhanced data for downstream structural feature analysis, including band detection and characterization. Subsequently, some studies have explored the differences in band signals between cells by visualizing pseudo-bulk bands in selected cells. Meanwhile, by combining ultra-high resolution chromatin conformation capture data at the single-cell level with high-resolution three-dimensional genome reconstruction methods centered on gene loci of interest, researchers have been able to characterize the multi-enhancer band structure at specific marker genes in depth, revealing the spatial relationship between band structure and regulatory elements.

[0005] Despite the progress made by existing methods in strip detection, the following limitations still exist:

[0006] 1. Insufficient adaptability to single-cell Hi-C data and high dependence on data coverage: Existing band detection technologies are mainly designed for batch Hi-C data. In actual analysis of single-cell Hi-C data, the single-cell Hi-C data is usually aggregated into a pseudo-batch representation and then the batch band detection method is directly applied. However, the pseudo-batch data still has a high proportion of missing contacts, which violates the basic assumption of such technologies for high-coverage data, resulting in a significant reduction in the number of bands that can be detected under single-cell data conditions.

[0007] 2. Difficulty in supporting high-resolution analysis of control scale: For techniques that detect bands after single-cell enhancement, existing single-cell Hi-C data enhancement methods mostly operate at relatively coarse resolutions such as 40-50kb. The enhancement effect on contact maps is limited under higher resolution conditions, thus limiting the ability to perform control scale band analysis on interpolated single-cell data.

[0008] 3. Limited analysis scope and insufficient scalability: Some methods rely heavily on specific cell or local site analysis, or use high-resolution three-dimensional structure reconstruction with high computational cost, making it difficult to achieve efficient whole-genome band analysis on large-scale single-cell datasets.

[0009] These shortcomings highlight the urgent need for a specialized framework that can enable efficient and scalable whole-genome banding analysis across single cells at a controlled resolution. Summary of the Invention

[0010] To address the problems existing in the prior art, this invention proposes a method for chromatin band analysis at the regulatory scale of single-cell Hi-C. It employs a two-stage statistical framework without interpolation to detect and quantify chromatin structural bands from sparse single-cell Hi-C data at a regulatory resolution.

[0011] In an embodiment of the present invention, a method for chromatin band analysis at the Hi-C regulatory scale in single cells is provided, comprising the following steps:

[0012] S1. Preprocess the single-cell Hi-C data to generate pseudo-batch Hi-C data; normalize the pseudo-batch Hi-C data to extract the Hi-C contact matrix for each chromosome.

[0013] S2. Identify and detect pseudo-batch stripes from the extracted Hi-C contact matrix;

[0014] S3. Project the pseudo-batch bands onto the original single-cell Hi-C data to obtain single-cell bands; perform quantitative analysis on the single-cell bands in the original single-cell Hi-C data;

[0015] Step S2 includes:

[0016] S21. The Hi-C contact matrix is ​​slid-divided along the main diagonal through a preset window to extract a series of window sub-matrices;

[0017] S22. Perform eigenvalue decomposition on each window submatrix, sort the eigenvalues ​​obtained by eigenvalue decomposition in descending order according to their absolute values, and obtain a set of descending sorted eigenvectors; retain the first r eigenvectors in the descending sorted eigenvectors.

[0018] S23. For each retained feature vector, a change point detection mechanism is introduced to identify the positions where significant changes occur along the genome coordinates, so as to determine the endpoint positions of the candidate bands and obtain a list of endpoints of the candidate bands.

[0019] S24. Based on the list of endpoints of the candidate strip, select adjacent endpoint pairs and determine the width range of the candidate strip at the anchor point based on the adjacent endpoint pairs. Then, combined with the strip direction, search for the third endpoint outside the adjacent endpoint pairs on one side of the candidate strip as the far endpoint. Determine the length of the candidate strip based on the far endpoint. Thus, determine the directional candidate strip rectangle based on the adjacent endpoint pairs and the far endpoint to construct the candidate strip.

[0020] S25. Perform statistical tests on the candidate bands to obtain significant bands; multiple significant bands form a significant band set, which is used as a pseudo-batch band.

[0021] Unlike previous approaches that directly applied existing batch stripe detection algorithms to pseudo-batch data constructed from single-cell Hi-C data, this invention designs a structure-aware stripe detection strategy to address the contact sparsity and noise characteristics prevalent in low-coverage pseudo-batch data. This strategy uses the pseudo-batch data to locate the spatial position and extension direction of candidate stripes, rather than simply using the pseudo-batch data as a substitute input for the batch data.

[0022] Compared with the prior art, the technical effects achieved by the present invention specifically include:

[0023] 1. Stable and Reliable Chromatin Band Detection under Low-Coverage Pseudo-Batch Hi-C Data Conditions: This invention addresses the limited coverage and smooth transition between structural signals and noise in pseudo-batch Hi-C data by constructing a robust band detection mechanism. Through spectral selection guided by random matrix theory to suppress dominant noise components, and combined with change point detection and direction-aware candidate band assembly strategies, the mechanism rigorously compares and filters the band interiors against the matching background in a statistical sense. This allows for the stable identification of chromatin bands with clear endpoints and directions, even when the assumption of high-coverage data is difficult to uphold in traditional batch band detection techniques, and outputs a statistically significant set of bands.

[0024] 2. Significantly improves the structural consistency and biological relevance of band detection results: Compared with existing band detection technologies, the bands identified by this invention present clearer and more continuous band-like signals in the aggregate contact map (i.e., aggregate contact matrix), and show a stable and clear correspondence with chromatin structure-related features and regulatory markers (such as structural protein CTCF and histone modification marker H3K4me3) in subsequent analysis, thereby improving the reliability and usability of band analysis results in structural resolution and biological interpretation.

[0025] 3. Single-cell band quantification without interpolation or enhancement alleviates the problem of ultra-sparse data and supports large-scale analysis: In the second stage, this invention directly calculates single-cell band scores on the unbalanced, uninterpolated, and unsmoothed raw single-cell Hi-C contact counts to obtain comparable single-cell band intensity estimates. The statistical framework design in this stage maintains intercellular variability while avoiding the potential bias and high computational overhead caused by interpolation enhancement or high-resolution 3D structure reconstruction, thus enabling this invention to be used for whole-genome band quantification of large-scale single-cell data.

[0026] 4. Supporting downstream cross-cell analysis and revealing intra-band heterogeneity: The single-cell band score matrix output by this invention can be used for downstream analysis, including identifying cell type-specific bands, cell embedding, and subpopulation segmentation; and can be jointly analyzed with co-analyzed single-cell RNA-seq (Single-Cell RNA Sequencing) data to improve clustering and subtype resolution. Furthermore, this invention allows for flexible specification of the distal regions of selected bands and can divide the same band into multiple non-overlapping distal segments for segmented quantification, thereby resolving inter-cell heterogeneity within multiple enhancer bands and revealing regulatory connectivity changes that are difficult to capture using traditional pseudo-batch or "single score per band per cell" techniques. Attached Figure Description

[0027] Figure 1 This is a flowchart of a method for analyzing chromatin bands at the Hi-C regulatory scale in single cells, as described in this embodiment of the invention. Detailed Implementation

[0028] The technical solution of the present invention will be further described in detail below with reference to the embodiments and accompanying drawings, but the implementation of the present invention is not limited thereto.

[0029] Example

[0030] This embodiment proposes a method for chromatin band analysis at the regulatory scale of single-cell Hi-C. The core idea is to provide a two-stage statistical framework that does not require interpolation, for detecting and quantifying chromatin structural bands from sparse single-cell Hi-C data at regulatory resolution.

[0031] like Figure 1 As shown, the single-cell Hi-C regulatory scale chromatin band analysis method in this embodiment specifically includes the following steps:

[0032] S1. Preprocess the single-cell Hi-C data to generate pseudo-batch Hi-C data; normalize the pseudo-batch Hi-C data to extract the Hi-C contact matrix for each chromosome.

[0033] In this embodiment, the data source for extracting the Hi-C contact matrix is ​​typically a public database, such as the mouse dataset in the GEO (Gene Expression Omnibus) database. This dataset contains single-cell Hi-C data for various tissue types. The band detection method of this embodiment can be applied to other relevant single-cell Hi-C datasets to test its applicability and generalization ability in different cell types.

[0034] In this embodiment, 7469 files with the .pairs.gz extension were downloaded during the research process. These files contained Hi-C data from different single cells. The Hi-C data from different single cells were integrated using a conventional software package to generate pseudo-batch Hi-C data files in a preset format. The Hi-C contact matrices were then normalized using iterative correction and eigenvector decomposition algorithms within the software package. Subsequently, the element values ​​of the Hi-C contact matrices corresponding to each chromosome in each pseudo-batch Hi-C data were read, ultimately obtaining 141911 contact matrices. The extracted contact matrix data was used for subsequent analysis.

[0035] S2. Identify and detect pseudo-batch stripes from the Hi-C contact matrix extracted in step S1.

[0036] This step is the first stage in a two-stage statistical framework, where significant bands are detected in the Hi-C contact matrix of the pseudo-batch data and used as pseudo-batch bands.

[0037] S21. The Hi-C contact matrix is ​​slidably divided through a preset window to extract a series of window sub-matrices.

[0038] This step involves windowing the contact matrix. A series of square window sub-matrices are extracted along the main diagonal of the Hi-C contact matrix. Specifically, windows of a preset size are slidably divided along the main diagonal of the Hi-C contact matrix to extract multiple window sub-matrices; for example, each window contains... Each matrix unit corresponds to a 10 kb resolution. The genomic region is defined; the window slides with a fixed step size, which can be set as needed; in a preferred embodiment, the fixed step size is 50 matrix units. Each window submatrix serves as the basic processing unit for subsequent band significance detection.

[0039] S22. Perform eigenvalue decomposition on each extracted window submatrix, sort the eigenvalues ​​obtained by eigenvalue decomposition in descending order according to their absolute values, and obtain a set of descending sorted eigenvectors; retain the first r eigenvectors in the descending sorted eigenvectors.

[0040] This step involves eigenvector selection based on spectral decomposition and random matrix theory. Specifically, eigenvalue decomposition is performed on each window submatrix to obtain a set of eigenvectors sorted by the absolute values ​​of their eigenvalues; that is, the window submatrix can be decomposed and represented as follows:

[0041] ;

[0042] in, Represents the window submatrix, Indicates the first 1 eigenvalue, ; , The superscript T represents the corresponding eigenvector, and the superscript T indicates matrix transpose. Indicates the eigenvalues The diagonal matrix formed has the following main diagonal elements as follows: to All off-diagonal elements are zero.

[0043] Different eigenvectors correspond to spatial structural information at different scales and patterns. Eigenvectors corresponding to eigenvalues ​​with large absolute values ​​typically reflect organized chromatin structures such as stripes and TADs, while eigenvectors corresponding to eigenvalues ​​with small absolute values ​​are mainly dominated by random noise. To distinguish between chromatin structure signals and random noise, this embodiment introduces random matrix theory as the basis for spectrum selection. Under the assumption of only random noise, the eigenvalue distribution follows the Marchenko-Pastur law, whose theoretical upper bound characterizes the maximum range of noise-dominated eigenvalues. Eigenvalues ​​exceeding this upper bound are considered to have potential structural significance.

[0044] However, considering that the transition between structured signals and random noise in actual pseudo-batch Hi-C data is usually continuous rather than abrupt, directly using the theoretical upper bound as a hard threshold may lead to unstable spectral selection. Therefore, this embodiment further combines the relative change trend of the eigenvalue spectrum, sorts the eigenvalues ​​in descending order of absolute value, and identifies the first local plateau interval in the eigenvalue sequence to determine the number of eigenvectors that need to be retained.

[0045] Specifically, the eigenvalues ​​are first sorted in descending order according to their absolute values:

[0046] ;

[0047] The number of eigenvectors retained in the descending sort is then determined using the following inequality constraints. :

[0048] ;

[0049] in , In this embodiment, m is set to 1.1. The top r eigenvectors retained from the descendingly sorted eigenvectors constitute the low-dimensional representation of the input window submatrix, used to characterize the main spatial structure information contained therein. This approach avoids over-reliance on theoretical thresholds while adaptively selecting the most representative structural eigenvectors in the current window submatrix.

[0050] After feature spectrum screening guided by the above random matrix theory, only the set of low-dimensional structural feature vectors used for subsequent strip endpoint detection is retained, thereby reducing noise interference while retaining the spatial structure information related to the strip.

[0051] S23. For each retained feature vector, a change point detection mechanism is introduced to identify the locations where significant changes occur along the genome coordinates, in order to determine the endpoint positions of the candidate bands and obtain a list of endpoints of the candidate bands.

[0052] This step performs change point detection to obtain candidate band endpoints. The change point detection mechanism determines the location of change points along the genome direction of the feature vector by minimizing a piecewise cost function with a penalty term, obtaining a set of change points, and outputting an ordered set of change point indices as a list of candidate band endpoints.

[0053] Specifically, for each retained feature vector The PELT (Pruned Exact Linear Time) change point detection algorithm is applied to identify locations where significant changes occur along genomic coordinates. The PELT change point detection algorithm determines the set of change points by minimizing a piecewise cost function with a penalty term; its objective function is expressed as:

[0054] ;

[0055] in, Represents the set of points of change. Indicates the first The index of the position of each point of change in the feature vector. This indicates the number of detected change points. ; Represents the piecewise cost function. Indicates the first The index of the position of each point of change in the feature vector. Indicates the first The index of the position of each point of change in the feature vector; Indicates from position index To location index The eigenvector subsequence, i.e. the th The _th change point and its adjacent _th Sub-intervals between points of change; Subinterval The piecewise cost function.

[0056] In other words, the objective function in this embodiment consists of the piecewise cost function of each sub-interval between adjacent change points. And a penalty term proportional to the number of change points. composition; This represents the penalty parameter, used to adjust the number of variation points to prevent over-segmentation due to noise. In this embodiment, it is set to... .

[0057] In this embodiment, the piecewise cost function A radial basis function (RBF)-based approach is adopted. The piecewise cost function based on RBF measures the overall consistency of feature vectors within a given interval. It is defined as a function of the interval length and the similarity between feature vectors at different positions within the interval, where parameters are introduced... The degree to which the difference in controllable eigenvalues ​​affects the cost; the specific definition formula is as follows:

[0058] ;

[0059] in, This represents the piecewise cost function based on radial basis functions; Indicates from position index To location index eigenvector subsequences, Represents the subsequence of feature vectors The first in 1 eigenvector Represents the subsequence of feature vectors The first in 1 eigenvector; Indicates the interval length. Representing the eigenvector With feature vectors The squared Euclidean distance between them. When the change in the feature vector within the interval is small and the numerical distribution is relatively consistent, the cost of the corresponding interval is low; conversely, when the feature vector changes significantly within the interval, the cost increases significantly, thus prompting the algorithm to introduce change points at the corresponding locations.

[0060] Through the above optimization process, the PELT change point detection algorithm can adaptively determine the location of change points along the genome direction of the feature vector while ensuring computational efficiency, and finally output a set of ordered change point indices as a list of candidate band endpoints.

[0061] S24. Based on the list of endpoints of the candidate strip, select adjacent endpoint pairs and determine the width range of the candidate strip at the anchor point based on the adjacent endpoint pairs. Then, combined with the strip direction, search for the third endpoint outside the adjacent endpoint pairs on one side of the candidate strip as the far endpoint. Determine the length of the candidate strip based on the far endpoint. Thus, determine the directional candidate strip rectangle based on the adjacent endpoint pairs and the far endpoint to construct the candidate strip.

[0062] Specifically, for each window submatrix, the ordered candidate strip endpoints are: Then any pair of adjacent endpoints A width range for candidate bands was defined. Based on the extension direction of the candidate bands relative to the anchor point in the Hi-C contact matrix, candidate bands were divided into 5′-bands or 3′-bands: 5′-bands refer to bands extending upstream from the anchor point in the genome, corresponding to structures extending along the column direction in the upper triangular region of the Hi-C contact matrix; 3′-bands refer to bands extending downstream from the anchor point in the genome, corresponding to structures extending along the row direction in the upper triangular region of the Hi-C contact matrix. For each candidate band's width range, based on its band type, the distal endpoint was searched along the corresponding extension direction. Under the premise of satisfying preset statistical criteria, the farthest position was determined as the distal endpoint of the band, thus obtaining the length of the candidate band.

[0063] To mitigate the distance attenuation effect and highlight locally enriched signals, each window submatrix is ​​converted into an observed / expected (O / E) matrix. This is for the window submatrix extracted after sliding segmentation of the Hi-C contact matrix. any position in Define its corresponding distance offset as At a given distance offset Under the condition, expected value Defined as a window submatrix Mid-distance offset The average value of the corresponding diagonal elements:

[0064] ;

[0065] in Represents the window submatrix The row index that satisfies the boundary constraints. Representation matrix In the line, number The element at the column; This represents the distance offset under matrix boundary constraints. All corresponding valid row indexes The set of , specifically defined as:

[0066] ;

[0067] in Represents the window submatrix The number of matrix units in a single dimension.

[0068] Through the above methods, the expected value This characterizes the background contact level under the same distance offset condition. Furthermore, dividing the observed value by the corresponding expected value yields the O / E matrix, which represents the ratio of observed to expected values.

[0069] ;

[0070] in This represents the observations after ICE (Iterative Correction and Eigenvector Decomposition) standardization, and is a window submatrix. Middle position The element value.

[0071] After the above processing, the influence of distance attenuation on the contact matrix can be effectively eliminated, thereby providing more stable and comparable input data for subsequent strip structure detection.

[0072] In this step, for each width interval, the strip direction is first determined based on the extension direction of the candidate strip relative to the anchor point in the Hi-C contact matrix, and the far endpoint, i.e., the far termination point, is searched along one side of the corresponding strip direction. The search strategy is to gradually approach the endpoint of the width interval from the farthest end. For the width interval and the far endpoint... Determined candidate strip regions Calculate its local multiple change value :

[0073] ;

[0074] in The average O / E ratio within the candidate band region. For the entire window submatrix The average O / E ratio. When the local multiple change value first reaches a significant local maximum value exceeding the set threshold, the corresponding position is selected as the far end point of the strip, i.e., the far end point.

[0075] Finally, candidate strips that are adjacent or overlap in spatial location are merged, and further filtered based on preset constraints of maximum strip width and minimum strip length. Candidate strips that do not meet the above constraints will be eliminated.

[0076] S25. Perform statistical tests on the constructed candidate bands to obtain significant bands; multiple significant bands form a significant band set, which is used as pseudo-batch bands.

[0077] All candidate bands constructed in step S24 require statistical testing for local enrichment significance. For each candidate band, three matching background regions are constructed to evaluate the statistical significance of local enrichment. By constructing three background regions, a rigorous and direction-aware null hypothesis model is defined, thereby improving the sensitivity of band detection. The band regions are denoted as... The background area to the left of the strip area is denoted as The background area to the right of the strip area is denoted as The left background area With the background area on the right The height, width, and anchor points are all consistent with the strip area. Alignment, left background area and the right background area This is collectively referred to as the background region in the width direction of the strip region. To determine whether a candidate strip truly terminates at a certain position in its extension direction, a local segment corresponding to the end of the candidate strip is extracted near the detected distal change position, denoted as... And construct a local fragment outside the location of the local fragment. The distant background region that matches in size and orientation is denoted as .

[0078] Based on the above region definition, the strip region and its background region in the width direction (i.e., and , and ), and the local segments at the end of the strip and the distant background region (i.e. and One-sided paired statistical tests were performed to evaluate the local enrichment significance of candidate bands in the width and extension directions. When a band region showed significant enrichment relative to the corresponding background region (background region or far-end background region in the width direction), the corresponding candidate band was determined to be a significant band structure.

[0079] Specifically, when comparing candidate bands with the background regions on their left and right sides in the width direction, for each candidate band, the band regions in the matrix O / E are compared along its width direction. The corresponding elements within the range are averaged to obtain a one-dimensional average contact intensity sequence arranged along the length of the candidate strip, denoted as . Similarly, in the left background area and the right background area Within, the corresponding one-dimensional average contact intensity sequences were obtained using the same method. and To mitigate the impact of extreme values ​​on the statistical test results and improve the stability of the comparison, the above average contact intensity sequence... , and Logarithmic transformation was performed on all samples before statistical analysis. Subsequently, corresponding positions in the three average contact intensity sequences were paired one-to-one along the length direction, and a one-sided paired statistical test was performed based on the paired samples to assess whether the candidate strip region exhibited a significantly enhanced contact intensity relative to its background region. In this embodiment, paired... The test, namely the paired-samples t-test, performs the aforementioned one-sided paired statistical test and obtains the corresponding results. value.

[0080] Similarly, along the band extension direction, local segments at the ends of candidate bands... and its corresponding far-end background area Perform one-sided pairing The test also yielded the corresponding results. Value. Specifically, the corresponding elements in the matrix O / E are averaged along the strip length direction to obtain a one-dimensional average contact intensity sequence arranged along the width direction, denoted as . and And after logarithmic transformation, comparisons are made based on position-by-position pairing.

[0081] For all candidate bands obtained Values ​​were adjusted using the Benjamini–Hochberg (BH) method with multiple hypothesis testing to control for the false discovery rate (FDR). The significance of the bands was determined based on the values ​​adjusted for multiple hypothesis testing. The value is determined. In one implementation, for the two statistical tests in the width direction, the FDR-corrected value is used. A value less than 0.001 was used as the significance threshold; for statistical tests in the length direction, the value corrected by FDR was used. A value less than 0.05 is used as the significance threshold.

[0082] After performing the hypothesis tests described above, a further filtering step based on fold change is executed to constrain the local enrichment intensity of the strips in the width direction. For each candidate strip, the fold change is calculated based on the average intensity value of the strip region and its corresponding background region, and is defined as:

[0083] ;

[0084] in, , and They are respectively , and The mean; This represents the multiple change value of the left background area. This represents the multiple change value of the background area on the right.

[0085] In one implementation, the minimum threshold for the multiple change is set to 1.1. Candidate bars that simultaneously meet the significance criteria and the multiple change threshold are retained as significant bars. Multiple significant bars eventually form a significant bar set, and the rectangular coordinate information of the corresponding bars, as well as their length, width, and other structural parameters, are output. The significant bar set is the pseudo-batch bar set.

[0086] This embodiment employs Aggregated Stripe Analysis (ASA) to systematically evaluate the stripes. ASA aligns and aggregates multiple stripes according to their anchor point positions and strip directions, generating an aggregated contact matrix with aligned anchor points. This allows for the characterization and quantification of the overall spatial structural features of the stripes within a unified reference framework. In the visualization results of ASA, the stripes identified in this embodiment exhibit continuous and consistent strip-shaped contact enrichment trajectories along the anchor point towards the distal end. The stripe boundaries are clear, the structure is intact, and they demonstrate good spatial consistency and directional stability. This result indicates that the stripes identified in this embodiment are more likely to correspond to repeatable and stable chromatin structure signals, rather than being caused by local noise or scattered contacts.

[0087] Furthermore, in terms of biological relevance, the bands identified in this embodiment exhibit a consistent correspondence with chromatin structures and regulatory markers: signal enrichment related to structural protein binding and active promoters can be observed at the band anchor points, such as CTCF (CCCTC-binding factor) and H3K4me3 (trimethylation of lysine 4 at histone H3); and after aligning multiple bands according to the anchor point positions, the relevant contact signals form obvious enrichment peaks near the anchor points, gradually weakening along the direction of the band towards the distal end, showing a directional attenuation characteristic extending from the anchor point to the distal end. This spatial distribution pattern conforms to the typical structural characteristics of bands as anchored regulatory structural units.

[0088] Therefore, this embodiment demonstrates the reliability and usability of band detection results in structural analysis and biological interpretation through the consistency between the ASA aggregation results and the structural / regulatory markers.

[0089] S3. Project the pseudo-batch bands obtained in step S2 onto the original single-cell Hi-C data to obtain single-cell bands; perform quantitative analysis on the single-cell bands in the original single-cell Hi-C data.

[0090] This step is the second stage in the two-stage statistical framework. The pseudo-batch bands obtained in step S2 are projected onto the original single-cell Hi-C data to directly quantify the variability of single-cell bands among cells at a resolution of 10 kb.

[0091] This embodiment achieves quantitative analysis of single-cell bands by calculating single-cell band scores. Specifically, the single-cell band score is used to measure band-related contact enrichment under ultra-sparse conditions by comparing the contact intensity between the band's interior and the matching background regions on its left and right sides.

[0092] The single-cell band score can be calculated under two complementary settings: one is a genome-wide analysis, which calculates the overall score for each single-cell band in each cell for all single-cell bands; the other is an intra-band tiling analysis, which calculates multiple segment-level scores in each cell for a selected single-cell band to analyze its internal heterogeneity. Multiple segment-level scores are used to characterize the internal structural distribution characteristics of the band within a single cell. By analyzing the spatial variation patterns of the segment-level score vectors, such as score differences between different segments, the location and continuity of high-scoring local segments, the spatial non-uniformity of contact enrichment within single-cell bands can be characterized. One of the single-cell band score calculation methods can be selected according to actual needs.

[0093] S31. Calculate the overall score of single-cell bands.

[0094] For a pseudo-batch stripe and a cell By comparing cells The observed contact number within a specific sub-region of the central stripe is compared with the contact number in the matching left and right background regions to define the overall score of the single-cell stripe. Let... This represents the overall stripe region corresponding to the pseudo-batch stripe. This represents the latter half of the subregion of the band extending along the direction of the single-cell band; denoted as... For cells In the latter half of the strip subregion The sum of the number of contacts observed in the middle, To be similar in size and shape to the latter half of the sub-region of the strip Matching left background area, To be similar in size and shape to the latter half of the sub-region of the strip For a matching right-side background region, the overall score of the single-cell band is defined as the ratio with pseudo-count stabilization. :

[0095] ;

[0096] in, For cells In the left background area The number of contacts observed in the middle, For cells In the background area on the right The number of contacts observed; the addition of a constant 1 is used to avoid a denominator of 0 and to stabilize the ratio under extremely sparse conditions, so that the overall score of the single-cell strip can reflect the enrichment of the strip region relative to the local background region.

[0097] It should be noted that the overall score of a single-cell band differs from the fold change statistic used in the first stage to screen for pseudo-batch candidate bands, i.e., the local fold change value. The first-stage fold change was calculated on the Hi-C contact matrix of a pseudo-batch with low coverage but low sparsity, using the O / E average of the entire pseudo-batch strip relative to... The window submatrix O / E average value is compared; while the second stage operates on the original single-cell observation and uses the matched left and right background regions for comparison.

[0098] The overall score of a single-cell band is the same as the whole-genome band score. This embodiment sets the whole-genome band score as follows: by default, the distal 50% of each single-cell band is used as the latter half of the band sub-region. This setting reduces the confounding effect of strong near-diagonal contacts: near the diagonal, short-range interactions are more concentrated and tend to be similarly elevated in the band and the background regions on both sides, thus reducing the ratio. Approaching 1 reduces sensitivity to band-specific enrichment. Confounding effects are particularly pronounced in most single-cell bands spanning multiple TADs. When calculating the overall score, focusing on the distal half of the band better captures enrichment from the anchor point to the distal end, while reducing near-diagonal bias.

[0099] Based on the calculated overall scores of single-cell bands, this embodiment can construct a single-cell band score matrix. Rows in the matrix correspond to different bands, and columns correspond to different single cells. Each element in the matrix represents the overall score of the corresponding band in its corresponding single cell, used to quantify the presence intensity and structural characteristics of the band at the single-cell level. The resulting single-cell band score matrix supports downstream analysis, including identifying cell type marker bands, cell embedding, and integration with co-tested single-cell RNA-seq data to achieve more refined clustering and subtype resolution.

[0100] S32. Calculate the segment-level scores of single-cell bands.

[0101] First, an internal tiling score is set for the selected single-cell band. By calculating scores for multiple non-overlapping distal segments of the same single-cell band, the heterogeneity within the band can be analyzed. This setting is particularly suitable for promoter-anchored multi-enhancer bands, where different distal segments correspond to different enhancers. For example, in EBF1 (Early B Cell Factor 1) promoter-anchored multi-enhancer bands, the main body of the single-cell band is divided into six consecutive distal segments, and each distal segment is scored separately. The number of segments and the division method are only examples and do not constitute a limitation, thereby revealing the heterogeneity of enhancer-promoter connectivity states and the progressive truncation pattern.

[0102] In this study, multiple segment-level scores are not summed or merged, but are preserved as a multidimensional feature vector describing the internal structural distribution of the same band within a single cell. By arranging the multiple segment-level scores corresponding to the same band in a single cell according to their spatial order in the genome, a segment-level score vector reflecting the contact enrichment distribution of the band along the distal direction can be obtained.

[0103] In this embodiment, the spatial variation pattern of the segment-level score vector is analyzed to characterize the non-uniformity of contact enrichment within the strip at the single-cell level. Specifically, the distribution difference of segment-level scores in the distal direction in different cells can be used to characterize the difference in the coverage of the strip in different cells; the relative high and low relationship of segment-level scores between adjacent segments can be used to reflect the gradient change of contact intensity within the strip; high-score or low-score regions appearing in local segments can be used to identify functional sub-regions within the strip; and the spatial continuity or discontinuity of high-score segments is used to characterize the continuous or discontinuous configuration state of the strip in a single cell.

[0104] Based on the aforementioned segment-level score vectors, individual cells can be further grouped or classified using any suitable method to distinguish cell populations with different internal structural distribution characteristics of bands. For cell populations with similar segment-level score vector patterns, their corresponding single-cell contact matrices can be aggregated to obtain a population-level band structure representation, thereby verifying the differences in band configurations corresponding to different segment-level score distribution patterns.

[0105] Through the above method, this embodiment achieves the analysis of the structural heterogeneity within a single-cell strip without interpolating the single-cell contact matrix or reconstructing the explicit three-dimensional structure. It is particularly suitable for the case of multiple enhancer strips anchored by promoters, where different distal segments can correspond to different enhancer regions, thereby revealing the differences and change patterns of enhancer-promoter connectivity states in different cells.

[0106] This invention employs a two-stage statistical framework without interpolation, providing a statistically sound and scalable framework for whole-genome band detection and single-cell analysis of single-cell Hi-C data at adjustable resolution. In the first stage, a spectral feature selection and change point detection process guided by random matrix theory is combined, and three statistical tests are used to quantitatively assess the separation between bands and background, thereby defining bands in the pseudo-batch data and obtaining a greater number of higher-quality bands with stronger biological relevance. In the second stage, the original single-cell Hi-C data is directly utilized. By comparing the band's internal region with the matching background region, a single-cell band score is calculated for each band defined in the pseudo-batch data, thus alleviating the extreme sparsity of single-cell Hi-C data and enabling quantitative analysis of intercellular variation without interpolation or 3D reconstruction. These single-cell band spectra support integrated analysis with co-analyzed single-cell RNA-seq data, resulting in more refined resolution of cell subtypes. In the example of EBF1 promoter anchoring bands, single-cell band scores successfully resolved enhancer connectivity differences between different cells in a 5kb resolution scMicro-C (single-cell Micro-C, a single-cell nucleosome resolution chromatin conformation capture technique) dataset.

[0107] The above embodiments are preferred embodiments of the present invention, but the embodiments of the present invention are not limited to the above embodiments. Any changes, modifications, substitutions, combinations, or simplifications made without departing from the spirit and principle of the present invention shall be considered equivalent substitutions and shall be included within the protection scope of the present invention.

Claims

1. A method for single-cell Hi-C regulatory scale chromatin banding analysis, characterized in that, The method comprises the following steps: S1, preprocessing single-cell Hi-C data to generate pseudo-bulk Hi-C data; S2, identifying and detecting pseudo-bulk bands from the extracted Hi-C contact matrix; S3, projecting the pseudo-bulk bands onto the original single-cell Hi-C data to obtain single-cell bands; S4, quantitatively analyzing the single-cell bands in the original single-cell Hi-C data; Step S2 comprises: S21, slidingly dividing the Hi-C contact matrix along the main diagonal direction through a preset window to extract a series of window sub-matrices; S22, performing eigenvalue decomposition on each window sub-matrix, and obtaining a set of descendingly sorted eigenvectors by descendingly sorting the eigenvalues obtained by the eigenvalue decomposition according to absolute values; the first r eigenvectors in the descendingly sorted eigenvectors are reserved; S23, for each reserved eigenvector, introducing a change point detection mechanism to identify positions where significant changes occur along the genomic coordinates, to determine the endpoint positions of the candidate bands, and obtaining an endpoint list of the candidate bands; S24, selecting adjacent endpoint pairs according to the endpoint list of the candidate bands, and determining the width interval of the candidate bands at the anchor points according to the adjacent endpoint pairs, and then combining a third endpoint outside the adjacent endpoint pairs on one side of the candidate bands as a distal endpoint according to the direction of the band, to determine the length of the candidate bands, and thus determining the directional candidate band rectangle according to the adjacent endpoint pairs and the distal endpoint to construct the candidate bands; S25, performing statistical test on the candidate bands to obtain significant bands; a plurality of significant bands form a significant band set, and the significant band set is taken as the pseudo-bulk bands. In step S22, the window sub-matrix is decomposed and expressed as:

2. The method for single-cell Hi-C regulatory scale chromatin banding analysis of claim 1, wherein, The eigenvalues are descendingly sorted according to the absolute values of the eigenvalues: ； wherein, denotes a window sub-matrix, denotes the th eigenvalue, ; , denotes the corresponding eigenvector, and the superscript T denotes the matrix transpose; denotes a diagonal matrix composed of eigenvalues , the main diagonal elements of which are to in turn, and the off-diagonal elements are all zero; The first r eigenvectors finally reserved in the descendingly sorted eigenvectors constitute a low-dimensional representation of the input window sub-matrix, which is used to represent the main spatial structure information contained therein. ； The number of eigenvectors to be reserved in the descendingly ordered eigenvectors is determined by the following inequality constraint condition : ； wherein , ; The change point detection mechanism of step S23 determines the change point positions of the eigenvectors along the genomic direction by minimizing a piecewise cost function with a penalty term, obtains a set of ordered change point indexes as the endpoint list of the candidate bands.

3. The method for single-cell Hi-C regulatory scale chromatin banding analysis of claim 1, wherein, The objective function of the piecewise cost function with a penalty term is:

4. The method for single-cell Hi-C regulatory scale chromatin banding analysis of claim 3, wherein, The definition formula of the piecewise cost function based on the radial basis function is: ； wherein, denotes a set of change points, denotes a position index of the change point in the feature vector, denotes the number of detected change points, ; denotes a segmentation cost function, denotes a position index of the change point in the feature vector, denotes a position index of the change point in the feature vector; denotes a sub-interval between the change point and its adjacent change point; denotes a segmentation cost function of the sub-interval ; denotes a penalty term; denotes a penalty parameter for adjusting the number of change points.

5. The method for single-cell Hi-C regulatory scale chromatin banding analysis of claim 4, wherein, The piecewise cost function adopts a piecewise cost function based on a radial basis function, which is defined as a function of interval length and similarity between feature vectors at each position in the interval, and the parameter controls the degree of influence of the difference between the control feature values on the cost.

6. The method for single-cell Hi-C regulatory scale chromatin banding analysis of claim 5, wherein, According to the extension direction of the candidate bands relative to the anchor points in the Hi-C contact matrix, the candidate bands are divided into 5'-bands or 3'-bands: the 5'-band refers to a band extending from the anchor point to the upstream direction of the genome, which corresponds to a structure extending along the column direction in the upper triangular region of the Hi-C contact matrix; the 3'-band refers to a band extending from the anchor point to the downstream direction of the genome, which corresponds to a structure extending along the row direction in the upper triangular region of the Hi-C contact matrix; ； wherein, represents a piecewise cost function based on a radial basis function; represents a sub-sequence of feature vectors from position index to position index represents the i-th feature vector in the sub-sequence of feature vectors represents the i-th feature vector in the sub-sequence of feature vectors represents an interval length, represents the squared Euclidean distance between the feature vector and the feature vector .​​​​​ 7. The method for single-cell Hi-C regulatory-scaled chromatin banding analysis of claim 1, wherein, In step S24, for each window sub-matrix, the ordered candidate slice endpoints are then any adjacent pair of endpoints defines a width interval of a candidate slice; For the width interval of each candidate band, according to the type of the band, the distal endpoint is searched along the corresponding extension direction, and under the premise of meeting the preset statistical criterion, the most distant position is determined as the distal endpoint of the band, so as to obtain the length of the candidate band; ​ The candidate bands that are adjacent or overlapped in spatial position are merged, and the candidate bands that do not satisfy the constraint conditions are removed based on the preset maximum band width and minimum band length constraint conditions.

8. The method for single-cell Hi-C regulatory-scaled chromatin banding analysis of claim 7, wherein, In step S24, each window sub-matrix is converted into an observation-to- expectation ratio (O / E) matrix; for a window sub-matrix at any position , the corresponding distance offset is defined as , the expectation is defined as the average of the diagonal elements in the window sub-matrix corresponding to the distance offset , given the distance offset . Defining window sub-matrices Mid position The element values are the observed values; dividing the observed values by the corresponding expected values gives the observed-to-expected ratio, O / E matrix.

9. The method for single-cell Hi-C regulatory-scaled chromatin banding analysis of claim 1, wherein, In step S25, for each candidate strip, three background regions matching the strip are constructed: the strip region is denoted as , the left background region of the strip region is denoted as , and the right background region of the strip region is denoted as , where the left background region is aligned with the strip region in height, width and anchor point, and the right background region is also aligned with the strip region in height, width and anchor point. The left background region and the right background region are collectively referred to as the background regions of the strip region in the width direction; a local segment corresponding to the end of the candidate strip is extracted near the detected far-end change position, denoted as , and a far-end background region matching the local segment in size and direction is constructed outside the position of the local segment, denoted as ; Respectively for the band regions and the distal background regions and the local fragments at the band ends Unpaired t-test was used to assess the local enrichment significance of the candidate bands in the width direction and the extension direction. When the band region is significantly enriched relative to the corresponding background region, the corresponding candidate band is determined as a significant band structure.

10. The method for single-cell Hi-C regulatory-scaled chromatin banding analysis of claim 1, wherein, Step S3 realizes quantitative analysis of the single-cell band by calculating a single-cell band score; the calculation method of the single-cell band score is: calculating a single-cell band overall score, or calculating a plurality of section-level scores of the single-cell band; S31, calculating a single-cell band overall score, comprising: For a pseudo-bulk band and a cell , the overall score of a single cell band is defined by comparing the observed number of contacts in a specific sub-region of the band region with the number of contacts in the matching left and right background regions of the cell ​ Let denote the overall band region corresponding to the pseudo-batch band, denote the band posterior sub-region along the direction of the single-cell band extension; let denote the cell contact count observed in the band posterior sub-region , denote the left background region matching the band posterior sub-region in size and shape, denote the right background region matching the band posterior sub-region in size and shape, then the overall score of the single-cell band is defined as the ratio of the band pseudo-count stabilized : ； in, For cells In the left background area The number of contacts observed in the middle, For cells In the background area on the right The number of contacts observed in the middle; According to the calculated single-cell band overall score, a single-cell band score matrix is constructed; the rows of the single-cell band score matrix correspond to different bands, and the columns correspond to different single cells; each element in the matrix represents the band overall score of the corresponding band in the corresponding single cell, which is used to quantify the existence intensity and structural characteristics of the band at the single-cell level; S32, calculating a plurality of section-level scores of the single-cell band, comprising: Setting the internal tiling score of the selected single-cell band, calculating the score of each of a plurality of non-overlapping distal sub-sections of the same single-cell band to analyze the heterogeneity inside the band; The plurality of section-level scores of the same band in a single cell are arranged in the spatial order in the genome to obtain a section-level score vector reflecting the contact enrichment distribution of the band along the distal direction.

Images