An easy-to-use and improved cross-species Hi-C contact map comparison framework

Abstract

This invention relates to the field of 3D genomics analysis technology, specifically to an easy-to-use and improved cross-species Hi-C contact map comparison framework, including modules for data processing, correction, optimization, diagnosis, interpretation, and full-process management. It utilizes standardized Cooler / .mcool format input and output, and achieves efficient cross-species Hi-C data comparison through intelligent preprocessing, multiple error correction, dynamic index optimization, intelligent visualization diagnosis, multi-omics association, and full-process monitoring. Its core lies in integrating technologies such as collinearity confidence grading, parameter self-optimization, and interactive diagnosis, addressing pain points in existing processes such as inconsistent formats, large errors, poor usability, and difficulty in reproducing results. This application significantly improves the accuracy and reproducibility of cross-species comparisons, lowers the operational threshold, expands analytical dimensions, and ensures analytical stability and efficiency, providing strong support for research on 3D genome evolutionary constraints and lineage-specific changes.

Description

Technical Field

[0001] This invention relates to the field of 3D genomics analysis technology, and more specifically, to an easy-to-use and improved cross-species Hi-C contact map comparison framework. Background Technology

[0002] Eukaryotic genomes possess a hierarchical three-dimensional organization, and their structure is crucial for gene expression regulation, DNA replication, and genome stability. Cross-species Hi-C contact map comparisons are a core tool for revealing evolutionary constraints and lineage-specific changes in 3D genomes. However, existing cross-species Hi-C analysis workflows suffer from significant drawbacks: they rely on multiple disparate scripts, lack standardized input / output formats, and lack supporting diagnostic tools, leading to unreliable and difficult-to-reproduce results. Traditional methods often depend on TAD or chromatin loop boundary calls and homology comparisons, which are easily affected by analysis tools and parameter settings, resulting in inconsistent boundary calls, strong resolution dependence, and amplified collinearity mapping errors. Although tools such as C-InterSecture have introduced relevant statistical measures, they still require handling complex configuration files, and I / O and visualization support are limited. These problems severely restrict the efficiency, accuracy, and reproducibility of cross-species 3D genome comparisons, necessitating an improved, user-friendly framework with standardized input / output, integrated error correction, and intelligent diagnostic capabilities to meet the needs of 3D genomics comparative research. Summary of the Invention

[0003] In view of this, the present invention addresses the shortcomings of the prior art by proposing an easy-to-use and improved cross-species Hi-C contact map comparison framework, aiming to solve at least one of the problems mentioned in the background art.

[0004] This invention provides an easy-to-use and improved cross-species Hi-C contact map comparison framework, comprising: a data processing module configured to receive input Hi-C data files, and to perform data preprocessing through resolution consistency verification, bin coverage grading processing, and interchromosomal interaction adaptation; wherein resolution consistency verification includes dynamic selection of downsampling algorithms, bin coverage grading processing includes edge region weight adjustment and low coverage region filtering logic, and interchromosomal interaction adaptation includes proportion judgment and classification labeling, and finally outputs a standardized contact data table;

[0005] A correction module, whose input is connected to the output of the data processing module, is configured to receive the standardized contact data table and complete the mapping of homologous segments through collinear chain confidence grading, chain fragmentation merging processing, and real-time verification of mapping results. The collinear chain confidence grading is based on three-level classification of identity and coverage and assigned weights. The chain fragmentation merging processing includes interval and direction judgment. The real-time verification of mapping results includes coordinate validity and rejection ratio control. Finally, the module outputs mapped contact pair data with confidence levels.

[0006] An optimization module, whose input is connected to the output of the correction module, is configured to receive the mapped contact pair data with confidence levels, perform comparative index calculation and parameter optimization based on logical judgment driven by research objectives and data characteristics, and output differential contact area and quantitative index data.

[0007] The diagnostic module has its input connected to the output of the optimization module. The diagnostic module is configured to receive the differential contact area and quantitative index data, and output multiple types of visualization results and interactive diagnostic information through result feature-driven visualization type recommendation and interactive anomaly localization logic.

[0008] The interpretation module has a first input terminal connected to the output terminal of the diagnostic module, and a second input terminal of the interpretation module receives external multi-omics data. It generates results through multi-omics association analysis and evolutionary analysis. The multi-omics association analysis includes the linkage between contact differences and gene expression. Finally, it outputs association results and a hierarchical interpretation report.

[0009] The processing module has its monitoring terminals connected to the running nodes of the aforementioned modules. The processing module is configured to achieve full-process management through node parameter verification, anomaly backtracking and solution recommendation, and historical result reuse judgment. The anomaly backtracking and solution recommendation includes tool error handling, and its control signal output terminals are connected to the parameter adjustment interfaces of the aforementioned modules to dynamically adjust the running parameters of each module.

[0010] In some embodiments, the data processing module is configured to receive an input Hi-C data file and perform data preprocessing through resolution consistency verification, bin coverage grading, and inter-chromosome interaction adaptation. Resolution consistency verification includes dynamic selection of downsampling algorithms; bin coverage grading includes edge region weight adjustment and low-coverage region filtering logic; and inter-chromosome interaction adaptation includes proportion judgment and classification labeling. When finally outputting a standardized contact data table, the module includes:

[0011] The dynamic selection of the resolution consistency check includes:

[0012] When the resolution difference of the input samples is less than or equal to a preset multiple, and the coefficient of variation of the contact count distribution is less than a preset first coefficient threshold, mean aggregation is used for downsampling.

[0013] When the resolution difference of the input samples is less than or equal to a preset multiple, and the coefficient of variation of the contact count distribution is greater than or equal to the first coefficient threshold, median aggregation is used for downsampling.

[0014] When the resolution difference is greater than the preset multiple, and the Pearson correlation coefficient of the coverage of adjacent bins is greater than or equal to the preset second coefficient threshold, linear interpolation is used to complete the count of high-resolution bins.

[0015] In some embodiments, the data processing module is configured to receive an input Hi-C data file and perform data preprocessing through resolution consistency verification, bin coverage grading, and inter-chromosome interaction adaptation. The resolution consistency verification includes dynamic selection of downsampling algorithms; the bin coverage grading includes edge region weight adjustment and low-coverage region filtering logic; and the inter-chromosome interaction adaptation includes proportion judgment and classification labeling. When finally outputting a standardized contact data table, the module further includes:

[0016] The edge region weight adjustment in the bin coverage hierarchical processing includes:

[0017] When the coverage of a bin is within a preset coverage threshold and is located in the middle of a chromosome arm, K-nearest neighbor interpolation is used to complete the missing count, and the K value is adjusted according to the average coverage of the surrounding bins.

[0018] When the average coverage is greater than the preset coverage range, K is 3;

[0019] When the average coverage is lower than the preset coverage range, K is 5;

[0020] When it is in the centromere or telomere region, it is marked as an edge bin and the weight is set to the coverage percentage.

[0021] In some embodiments, the input of the correction module is connected to the output of the data processing module. The correction module is configured to receive the standardized contact data table and complete the mapping of homologous fragments through collinear chain confidence grading, chain fragmentation merging, and real-time verification of mapping results. The collinear chain confidence grading is based on a three-level classification of identity and coverage, with assigned weights. The chain fragmentation merging includes interval and direction judgment. The real-time verification of mapping results includes coordinate validity and rejection ratio control. When finally outputting mapped contact pair data with confidence levels, the following steps are included:

[0022] The three-level classification and weight assignment of the collinear chain confidence level include:

[0023] When the identity is greater than or equal to the preset first setting value, and the coverage is greater than or equal to the preset first coverage setting value, it is determined to be a high confidence chain and assigned the first preset weight;

[0024] When the identity is greater than or equal to the preset second setting value and less than the first setting value, or the coverage is greater than or equal to the preset second coverage setting value and less than the first coverage setting value, it is determined to be a middle confidence chain and assigned a second preset weight.

[0025] When the identity is less than the second set value and the coverage is less than the second coverage set value, it is determined to be a low-confidence chain and assigned a third preset weight.

[0026] In some embodiments, the input of the correction module is connected to the output of the data processing module. The correction module is configured to receive the standardized contact data table and complete the mapping of homologous segments through collinear chain confidence grading, chain fragmentation merging, and real-time verification of mapping results. The collinear chain confidence grading is based on a three-level classification of identity and coverage, with assigned weights. The chain fragmentation merging includes interval and direction judgment. The real-time verification of mapping results includes coordinate validity and rejection ratio control. When finally outputting mapped contact pair data with confidence levels, the module further includes:

[0027] In the correction module, the interval and direction determination of the chain fragmentation merging process includes: for short chains with a length less than a preset third setting value;

[0028] When the interval between the chain and the adjacent chain is less than the preset fourth setting value and the comparison direction is consistent, they are merged into a long chain and the original coordinates are recorded.

[0029] If the interval with an adjacent chain is greater than or equal to the fourth setting value or the comparison direction is opposite, it is retained as a fragmented chain and is only allowed to be associated with a high-confidence chain for mapping.

[0030] In some embodiments, the input of the optimization module is connected to the output of the correction module. The optimization module is configured to receive the mapped contact pair data with confidence levels, perform comparative index calculation and parameter optimization based on research target adaptation and data feature-driven logical judgment, and output differential contact area and quantitative index data, including:

[0031] In the optimization module, the preferred selection of comparison indicators based on the research objective includes: when the research objective is a conservative analysis, the correlation coefficient between PBAD and Pearson is calculated first.

[0032] When the research objective is to locate differential structures, log ratio and strip score should be calculated first.

[0033] In some embodiments, the input of the optimization module is connected to the output of the correction module. The optimization module is configured to receive the mapped contact pair data with confidence levels, perform comparative index calculation and parameter optimization based on research target adaptation and data feature-driven logical judgment, and output differential contact area and quantitative index data, including:

[0034] The data feature-driven parameter self-optimization includes: when the noise level is such that the standard deviation of the contact count in the background area is greater than or equal to the preset fifth setting value, the logarithmic ratio is Gaussian smoothed, and the smoothing kernel size is adjusted according to the noise distribution, with a 3x3 kernel for concentrated noise and a 5x5 kernel for dispersed noise;

[0035] The bandwidth of the band score is adjusted with the genome resolution: 40kb resolution corresponds to 200kb bandwidth, 20kb resolution corresponds to 100kb bandwidth, and 10kb resolution corresponds to 50kb bandwidth.

[0036] In some embodiments, the input of the diagnostic module is connected to the output of the optimization module. The diagnostic module is configured to receive the differential contact area and quantitative index data, and output multiple types of visualization results and interactive diagnostic information through result feature-driven visualization type recommendation and interactive anomaly localization logic, including:

[0037] The recommended visualization types based on differential features include:

[0038] When the overlap rate of extreme percentile contact pairs between two species is less than the preset sixth setting value, it is recommended to generate an extreme value clustering diagram.

[0039] When the contact intensity difference of a certain distance bin is greater than or equal to the preset seventh setting value, it is recommended to generate a magnified map of the local deviation and mark the chromosome regions with significant differences.

[0040] In some embodiments, the first input terminal of the interpretation module is connected to the output terminal of the diagnostic module, and the second input terminal of the interpretation module receives external multi-omics data, generating results through multi-omics association analysis and evolutionary analysis; wherein the multi-omics association analysis includes the linkage between contact differences and gene expression, and the final output of association results and hierarchical interpretation reports includes:

[0041] The interpretation module includes multi-omics association analysis, which involves associating differential contact regions with TSS regions of differentially expressed genes in RNA-seq. When the number of associated genes is greater than or equal to a preset eighth setting value, GO enrichment analysis is performed and a contact differential-gene expression association heatmap is generated.

[0042] The differential contact region is compared with the H3K27acChIP-seq reinforcing sub-region. When the overlap rate is greater than or equal to the preset ninth setting value, the correlation between contact strength and reinforcing sub-activity is calculated. When the correlation coefficient is greater than or equal to the preset correlation coefficient threshold, it is determined to be a significant correlation.

[0043] In some embodiments, the monitoring terminal of the processing module is connected to the running nodes of each of the above modules respectively. The processing module is configured to achieve full-process management through node parameter verification, anomaly backtracking and solution recommendation, and historical result reuse judgment. The anomaly backtracking and solution recommendation includes tool error handling, and its control signal output terminal is connected to the parameter adjustment interface of each of the above modules respectively. When dynamically adjusting the running parameters of each module, it includes:

[0044] In the processing module, the anomaly backtracking and solution recommendation include: when the asm2link stage fails due to insufficient memory, adjusting parameters to compare in batches;

[0045] When the contact-stat stage overflows with memory, the graphs generated by chromosomes are merged.

[0046] When the tool fails to be invoked, the corresponding conda installation command or environment variable configuration method will be displayed.

[0047] Compared with existing technologies, the advantages of this invention are as follows: Standardized data input and output and intelligent preprocessing unify Cooler / .mcool format adaptation, avoiding multi-script dependencies and format confusion, significantly improving the usability and reproducibility of the analysis process; integration of multiple error correction mechanisms such as collinearity confidence grading, fragmentation merging, and mapping verification effectively reduces the impact of issues such as inconsistent boundary calls and collinearity mapping errors, improving the accuracy of cross-species comparisons; dynamic optimization of indicator selection and parameter configuration based on research objectives and data characteristics, coupled with intelligent visualization recommendations and interactive diagnostics, enables precise localization of differential structures and intuitive interpretation of results; integration of multi-omics association and evolutionary analysis linkage functions expands the analysis dimensions, providing more comprehensive support for 3D genome structure and function association research; full-process monitoring and anomaly handling modules ensure analysis stability, and the historical result reuse function reduces redundant calculations and improves research efficiency.

[0048] The above general description and the following detailed description are exemplary and explanatory only, and are not intended to limit this disclosure.

[0049] Other features and aspects of this disclosure will become clearer from the following detailed description of exemplary embodiments with reference to the accompanying drawings. Attached Figure Description

[0050] To more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the specific embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0051] Figure 1 Functional block diagram of an easy-to-use and improved cross-species Hi-C contact diagram comparison framework provided in embodiments of the present invention. Detailed Implementation

[0052] Exemplary embodiments of the present disclosure will now be described in more detail with reference to the accompanying drawings. While exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure may be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided to enable a more thorough understanding of the present disclosure and to fully convey the scope of the disclosure to those skilled in the art. It should be noted that, unless otherwise specified, embodiments and features in the embodiments of the present invention can be combined with each other. The present invention will now be described in detail with reference to the accompanying drawings and embodiments.

[0053] See Figure 1 As shown, an easy-to-use and improved cross-species Hi-C contact map comparison framework according to an embodiment of this application includes:

[0054] The data processing module is configured to receive input Hi-C data files and perform data preprocessing through resolution consistency verification, bin coverage classification processing, and inter-chromosome interaction adaptation. The resolution consistency verification includes dynamic selection of downsampling algorithms, the bin coverage classification processing includes edge region weight adjustment and low coverage region filtering logic, and the inter-chromosome interaction adaptation includes proportion judgment and classification labeling. Finally, a standardized contact data table is output.

[0055] A correction module, whose input is connected to the output of the data processing module, is configured to receive the standardized contact data table and complete the mapping of homologous segments through collinear chain confidence grading, chain fragmentation merging processing, and real-time verification of mapping results. The collinear chain confidence grading is based on three-level classification of identity and coverage and assigned weights. The chain fragmentation merging processing includes interval and direction judgment. The real-time verification of mapping results includes coordinate validity and rejection ratio control. Finally, the module outputs mapped contact pair data with confidence levels.

[0056] An optimization module, whose input is connected to the output of the correction module, is configured to receive the mapped contact pair data with confidence levels, perform comparative index calculation and parameter optimization based on logical judgment driven by research objectives and data characteristics, and output differential contact area and quantitative index data.

[0057] The diagnostic module has its input connected to the output of the optimization module. The diagnostic module is configured to receive the differential contact area and quantitative index data, and output multiple types of visualization results and interactive diagnostic information through result feature-driven visualization type recommendation and interactive anomaly localization logic.

[0058] The interpretation module has a first input terminal connected to the output terminal of the diagnostic module, and a second input terminal of the interpretation module receives external multi-omics data. It generates results through multi-omics association analysis and evolutionary analysis. The multi-omics association analysis includes the linkage between contact differences and gene expression. Finally, it outputs association results and a hierarchical interpretation report.

[0059] The processing module has its monitoring terminals connected to the running nodes of the aforementioned modules. The processing module is configured to achieve full-process management through node parameter verification, anomaly backtracking and solution recommendation, and historical result reuse judgment. The anomaly backtracking and solution recommendation includes tool error handling, and its control signal output terminals are connected to the parameter adjustment interfaces of the aforementioned modules to dynamically adjust the running parameters of each module.

[0060] It should be understood that the cross-species Hi-C contact map comparison framework of this application achieves standardized, high-precision, and intelligent analysis of cross-species Hi-C data through a six-level linkage mechanism of "preprocessing-correction-optimization-diagnosis-interpretation-full-process control". The six modules are arranged in a serial + parallel monitoring architecture: the data processing module completes the cleaning and standardization of raw data, providing high-quality input for subsequent analysis; the correction module solves the error problem of cross-species homologous fragment mapping, ensuring the accuracy of mapping; the optimization module dynamically adapts the analysis strategy according to research needs, improving the relevance of the results; the diagnosis module transforms quantitative data into intuitive visualization results, assisting in anomaly localization; the interpretation module integrates multi-omics data, expanding the depth of analysis and application value; and the processing module monitors the running status of each node throughout the process, ensuring that the process is stable, efficient, and reproducible.

[0061] Data processing module: Its core function is to eliminate format differences and quality defects in the original Hi-C data. It unifies the analytical scale of different samples through resolution consistency verification, avoiding comparison biases caused by resolution mismatch; it repairs the counting information of low-quality bins and filters invalid data through bin coverage grading, ensuring the reliability of each bin; and it classifies inter-chromosome / intra-chromosome interaction data through inter-chromosome interaction adaptation, providing a clearly categorized and standardized contact data table for subsequent analysis.

[0062] Correction Module: Focuses on error correction for cross-species homologous fragment mapping. Collinearity chain confidence grading is used to screen high-quality homologous chains, avoiding mapping errors introduced by low-reliability chains; chain fragmentation and merging processes address mapping breakage caused by short chain dispersion; mapping results are verified in real time to remove invalid mapping pairs, controlling errors within acceptable ranges, and finally outputting mapping contact pairs with confidence levels, providing a high-quality data foundation for subsequent index calculations.

[0063] Optimization Module: Enables personalized adaptation of analysis strategies. It selects core comparison indicators based on research objectives, avoiding blind selection of indicators; it optimizes parameter configuration based on data characteristics (such as noise level and sample size) to address the result bias caused by fixed parameters under different data characteristics, accurately identifying areas of difference and quantifying the degree of difference.

[0064] Diagnostic Module: Building a bridge between "data, visualization, and diagnosis." It recommends suitable visualization types based on the characteristics of the differences, allowing users to quickly focus on the core discrepancies; through interactive anomaly location functions, it supports users in exploring the details of the discrepancy areas, providing intuitive clues for subsequent interpretation.

[0065] Interpretation Module: This module transforms "data differences" into "biological significance." It establishes associations between 3D genome structural differences and functional phenotypes such as gene expression and enhancer activity through multi-omics association analysis; it reveals the evolutionary patterns of structural differences through evolutionary analysis, ultimately outputting hierarchical interpretation reports to meet the needs of users with diverse professional backgrounds.

[0066] Processing Module: The framework's "central control system." Node parameter verification ensures the consistency and integrity of data transmission between modules; anomaly backtracking and solution recommendations enable rapid response to operational failures, lowering the operational threshold; historical result reuse reduces redundant calculations and improves analysis efficiency, while dynamically adjusting the parameters of each module ensures the process is always in optimal operating condition.

[0067] In some specific embodiments, the data processing module is configured to receive input Hi-C data files and perform data preprocessing through resolution consistency verification, bin coverage grading, and interchromosomal interaction adaptation. Resolution consistency verification includes dynamic selection of downsampling algorithms; bin coverage grading includes edge region weight adjustment and low-coverage region filtering logic; and interchromosomal interaction adaptation includes proportion judgment and classification labeling. When finally outputting a standardized contact data table, the module includes:

[0068] The dynamic selection of the resolution consistency check includes:

[0069] When the resolution difference of the input samples is less than or equal to a preset multiple, and the coefficient of variation of the contact count distribution is less than a preset first coefficient threshold, mean aggregation is used for downsampling.

[0070] When the resolution difference of the input samples is less than or equal to a preset multiple, and the coefficient of variation of the contact count distribution is greater than or equal to the first coefficient threshold, median aggregation is used for downsampling.

[0071] When the resolution difference is greater than the preset multiple, and the Pearson correlation coefficient of the coverage of adjacent bins is greater than or equal to the preset second coefficient threshold, linear interpolation is used to complete the count of high-resolution bins.

[0072] It should be understood that resolution consistency verification is a core step in data preprocessing. Its purpose is to unify Hi-C samples with different resolutions to the same analytical scale, avoiding distortion in contact intensity contrast caused by resolution differences. This embodiment uses a dual logic of "resolution difference judgment + data feature adaptation" to dynamically select the optimal processing scheme, ensuring resolution uniformity while preserving the effective information of the original data to the greatest extent.

[0073] The preset multiplier is 2; the preset first coefficient threshold (coefficient of variation threshold) is 0.3; and the preset second coefficient threshold (Pearson correlation coefficient threshold) is 0.6.

[0074] When the input sample resolution difference is ≤2 times (e.g., sample A is 20kb and sample B is 40kb), a downsampling strategy is used to unify the data to a lower resolution. In this case, the data dispersion is determined by the coefficient of variation (COP) of the contact count distribution: a COP < 0.3 indicates a uniform data distribution with no obvious extreme values, and mean aggregation (averaging the counts from multiple high-resolution bins and merging them into a single low-resolution bin) is used, balancing efficiency and accuracy; a COP ≥ 0.3 indicates the presence of extreme high / low counts, and median aggregation is used to avoid interference from extreme values ​​in the merged counts, ensuring data accuracy.

[0075] When the resolution difference is greater than 2 times (e.g., sample A is 10kb and sample B is 100kb), if the Pearson correlation coefficient of the coverage of adjacent bins of the low-resolution sample is ≥0.6, it indicates that the contact pattern continuity of adjacent bins is good. Linear interpolation (fitting the expected count of the high-resolution bin based on the count of adjacent bins) is used to complete the count of the high-resolution bin, thereby improving the resolution and ensuring the precision of the comparison. If the correlation coefficient is <0.6, it indicates poor continuity and the interpolation result is unreliable. The user will be prompted to supplement data of the same resolution to avoid errors introduced by invalid interpolation.

[0076] In some specific embodiments, the data processing module is configured to receive input Hi-C data files and perform data preprocessing through resolution consistency verification, bin coverage grading, and inter-chromosome interaction adaptation. Resolution consistency verification includes dynamic selection of downsampling algorithms; bin coverage grading includes edge region weight adjustment and low-coverage region filtering logic; inter-chromosome interaction adaptation includes proportion judgment and classification labeling; and when finally outputting a standardized contact data table, it also includes:

[0077] The edge region weight adjustment in the bin coverage hierarchical processing includes:

[0078] When the coverage of a bin is within a preset coverage threshold and is located in the middle of a chromosome arm, K-nearest neighbor interpolation is used to complete the missing count, and the K value is adjusted according to the average coverage of the surrounding bins.

[0079] When the average coverage is greater than the preset coverage range, K is 3;

[0080] When the average coverage is lower than the preset coverage range, K is 5;

[0081] When it is in the centromere or telomere region, it is marked as an edge bin and the weight is set to the coverage percentage.

[0082] It should be understood that bin coverage directly reflects the richness of effective interaction information for each genome bin in Hi-C data. Bins with low coverage will lead to distorted contact counts, affecting the accuracy of subsequent analyses. This embodiment uses a logic of "coverage grading + regional localization + dynamic interpolation" to differentiate bins of different qualities, thereby both restoring effective data and filtering out invalid data, while reducing interference from edge regions.

[0083] The preset coverage threshold range is 50%-80%; the preset mean threshold (average coverage threshold of surrounding bins) is 60%.

[0084] When the coverage of a bin is between 50% and 80% (medium coverage), its position is first determined by genomic coordinates: if it is located in the middle of a chromosome arm (far from the core functional region of the centromere and telomere), K-nearest neighbor interpolation is used to complete the missing count—when the mean coverage of surrounding bins is >60% (high mean, strong data reliability), K is set to 3 (only referencing 3 adjacent high-quality bins to avoid interference from low-quality data at a distance); when the mean is <60% (low mean, weak data reliability), K is set to 5 (referencing more adjacent bins to improve interpolation stability).

[0085] If a medium coverage bin is located in a centromere or telomere region (located through the species genome annotation file), it is marked as an "edge bin": These regions have special genomic structures (such as centromeres rich in repetitive sequences), and the reliability of Hi-C data is naturally lower. Their weight is set to the coverage percentage (e.g., a weight of 0.6 for 60% coverage) to reduce their impact on the results in subsequent index calculations and avoid misleading conclusions from edge region data.

[0086] When the coverage of a bin is less than 50% (low coverage), if the proportion of low-coverage bins on the chromosome is greater than or equal to 10%, it is determined to be a globally low-quality chromosome region, prompting the user to check the original sequencing data (such as sequencing depth and alignment rate) to prevent low-quality data from entering the subsequent process.

[0087] In some specific embodiments, the input of the correction module is connected to the output of the data processing module. The correction module is configured to receive the standardized contact data table and complete the mapping of homologous fragments through collinear chain confidence grading, chain fragmentation merging, and real-time verification of mapping results. The collinear chain confidence grading is based on a three-level classification of identity and coverage, with assigned weights. The chain fragmentation merging includes interval and direction judgment. The real-time verification of mapping results includes coordinate validity and rejection ratio control. When finally outputting mapped contact pair data with confidence levels, the following steps are included:

[0088] The three-level classification and weight assignment of the collinear chain confidence level include:

[0089] When the identity is greater than or equal to the preset first setting value, and the coverage is greater than or equal to the preset first coverage setting value, it is determined to be a high confidence chain and assigned the first preset weight;

[0090] When the identity is greater than or equal to the preset second setting value and less than the first setting value, or the coverage is greater than or equal to the preset second coverage setting value and less than the first coverage setting value, it is determined to be a middle confidence chain and assigned a second preset weight.

[0091] When the identity is less than the second set value and the coverage is less than the second coverage set value, it is determined to be a low-confidence chain and assigned a third preset weight.

[0092] It should be understood that collinear chains are the foundation of cross-species homologous fragment mapping, and their reliability directly determines the accuracy of the mapping results. This embodiment uses a dual metric of "sequence identity" and "alignment integrity" to classify collinear chains into three levels of confidence and assign weights, thereby screening high-quality chains, labeling medium-quality chains, and restricting low-quality chains, thus reducing mapping errors from the source.

[0093] The first preset value (identity threshold) is 90%; the first preset coverage value is 85%; the second preset value (identity threshold) is 70%; the second preset coverage value is 70%; the first preset weight is 1.0; the second preset weight is 0.7; and the third preset weight is 0.3.

[0094] High confidence chain: identity ≥ 90% and coverage ≥ 85%, indicating that the fragment sequences of the two species are highly consistent and the alignment is complete, which is a reliable homologous fragment correspondence. It is assigned a weight of 1.0 and is used first in subsequent mapping as the core mapping basis.

[0095] Medium confidence chain: identity ≥ 70% and < 90%, or coverage ≥ 70% and < 85%, indicating that the sequence consistency or alignment integrity is somewhat insufficient, but still has homology characteristics. It is assigned a weight of 0.7 and can be used as an auxiliary mapping basis. It needs to be used in conjunction with the high confidence chain to reduce errors.

[0096] Low-confidence chains: identity < 70% and coverage < 70%, indicating poor sequence consistency, incomplete alignment, and low homology confidence. They are assigned a weight of 0.3 and are only used as a supplement when there are no high / medium confidence chains. Low-weight mapping pairs are strictly filtered in subsequent index calculations (e.g., only mapping pairs with a weight ≥ 0.5 are retained).

[0097] In some specific embodiments, the input of the correction module is connected to the output of the data processing module. The correction module is configured to receive the standardized contact data table and complete the mapping of homologous segments through collinear chain confidence grading, chain fragmentation merging, and real-time verification of mapping results. The collinear chain confidence grading is based on a three-level classification of identity and coverage, with assigned weights. The chain fragmentation merging includes interval and direction judgment. The real-time verification of mapping results includes coordinate validity and rejection ratio control. When finally outputting mapped contact pair data with confidence levels, the module further includes:

[0098] In the correction module, the interval and direction determination of the chain fragmentation merging process includes: for short chains with a length less than a preset third setting value;

[0099] When the interval between the chain and the adjacent chain is less than the preset fourth setting value and the comparison direction is consistent, they are merged into a long chain and the original coordinates are recorded.

[0100] If the interval with an adjacent chain is greater than or equal to the fourth setting value or the comparison direction is opposite, it is retained as a fragmented chain and is only allowed to be associated with a high-confidence chain for mapping.

[0101] It should be understood that fragmented short strands (collinear strands of excessively short length) are easily generated during genome alignment. These short strands may lead to breaks or misalignments in the mapping of homologous fragments, affecting the accuracy of the correspondence between cross-species contact pairs. This embodiment uses the logic of "length judgment + interval and direction verification" to merge or restrict the use of short strands, optimize the collinear strand structure, and improve the continuity and reliability of mapping.

[0102] The preset third setting value (short chain length threshold) is 10kb; the preset fourth setting value (chain spacing threshold) is 5kb.

[0103] For short chains with a length of <10kb, first determine the distance between them and adjacent collinear chains and the alignment direction: if the distance is <5kb and the alignment direction is consistent (e.g., both are forward alignments), it means that the two short chains belong to the same continuous homologous region. They are merged into one long chain and the original coordinates are recorded (for subsequent backtracking verification) to avoid mapping breakage caused by fragmentation.

[0104] If the short chain is ≥5kb apart from the adjacent chain, or the alignment direction is opposite (e.g., one in the forward direction and one in the reverse direction), it means that the two do not belong to the continuous homogeneous region and are retained as fragmented chains. However, their use is restricted—they are only allowed to be associated with high-confidence chains for mapping (e.g., the fragmented chain must overlap with the high-confidence chain by ≥50% to participate in the mapping), to avoid mapping misalignment errors caused by using fragmented chains alone.

[0105] If the fragmented chain ratio in the .mark file is ≥20%, the MCScanX tool will be automatically invoked to re-identify collinear blocks and optimize the chain structure. The optimized chain will only be used when the fragmented chain ratio is reduced by ≥10% after optimization, to ensure the overall reliability of the chain structure.

[0106] In some specific embodiments, the input of the optimization module is connected to the output of the correction module. The optimization module is configured to receive the mapped contact pair data with confidence levels, and to perform comparative index calculation and parameter optimization based on logical judgment driven by research objectives and data features. When outputting the difference contact area and quantitative index data, the optimization module includes:

[0107] In the optimization module, the preferred selection of comparison indicators based on the research objective includes: when the research objective is a conservative analysis, the correlation coefficient between PBAD and Pearson is calculated first.

[0108] When the research objective is to locate differential structures, log ratio and strip score should be calculated first.

[0109] It should be understood that the research objectives of cross-species Hi-C comparisons differ (such as focusing on conservation and differential structure), and the suitability of different comparison indicators varies. This embodiment dynamically selects core comparison indicators based on research objectives to ensure that the indicators are highly matched with research needs, improve the relevance and practicality of the results, and avoid the omission of key information due to a "one-size-fits-all" approach to indicator selection.

[0110] The core requirement is to identify regions with highly similar 3D genome structures between two species. Priority is given to calculating PBAD (Percentile Background Adjusted Difference) and Pearson correlation coefficient: PBAD quantifies the difference in contact intensity after removing background noise; the closer the value is to 0, the stronger the conservation. The Pearson correlation coefficient reflects the linear correlation between the contact intensity of the two species; the closer the value is to 1, the higher the conservation. Combining these two methods allows for precise screening of conserved regions. Simultaneously, the significance threshold is adjusted according to the sample size (p < 0.01 when sample size ≥ 3, p < 0.05 when sample size = 1) to ensure statistical reliability.

[0111] The core requirement is to identify species-specific 3D genome structures (such as unique TADs and loop structures). Priority is given to calculating the logarithmic ratio (log2((obs+)). ) / (ctrl+ The log-ratio visually represents the multiple difference in the contact intensity between two species (positive values ​​indicate OBS species enrichment, and negative values ​​indicate Ctrl species enrichment), facilitating rapid location of differential regions. The band score focuses on the unilateral enrichment differences on both sides of the diagonal, accurately capturing local structural differences (such as band enrichment caused by chromatin loops unique to a particular species). The combination of the two enables precise location and quantification of differential structures.

[0112] In some specific embodiments, the input of the optimization module is connected to the output of the correction module. The optimization module is configured to receive the mapped contact pair data with confidence levels, and to perform comparative index calculation and parameter optimization based on logical judgment driven by research objectives and data features. When outputting the difference contact area and quantitative index data, the optimization module includes:

[0113] The data feature-driven parameter self-optimization includes: when the noise level is such that the standard deviation of the contact count in the background area is greater than or equal to the preset fifth setting value, the logarithmic ratio is Gaussian smoothed, and the smoothing kernel size is adjusted according to the noise distribution, with a 3x3 kernel for concentrated noise and a 5x5 kernel for dispersed noise;

[0114] The bandwidth of the band score is adjusted with the genome resolution: 40kb resolution corresponds to 200kb bandwidth, 20kb resolution corresponds to 100kb bandwidth, and 10kb resolution corresponds to 50kb bandwidth.

[0115] It should be understood that Hi-C data vary in noise levels, resolution, and other characteristics, and fixing parameters can lead to distorted index calculation results. This embodiment dynamically optimizes parameter configuration based on data characteristics, and improves the accuracy and reliability of index calculation through noise reduction processing, bandwidth adaptation, and other means, ensuring that the results can truly reflect the differences in 3D genome structure across species.

[0116] The default fifth setting value (background area contact count standard deviation threshold) is 1.0.

[0117] Noise level adaptation (Gaussian smoothing optimization): When the standard deviation of contact count in the background region (non-TAD, non-ring structure non-functional region) is ≥1.0, it is judged as high noise data, and the logarithmic ratio is Gaussian smoothed. When the noise is concentrated (such as noise mainly distributed in a specific distance bin), a 3x3 kernel is used (small kernel smoothing to avoid excessive blurring of effective differences); when the noise is dispersed (such as noise spread across multiple distance bins), a 5x5 kernel is used (large kernel smoothing to ensure noise reduction effect). The dynamic kernel size balances noise reduction and effective information preservation.

[0118] Resolution adaptation (band score bandwidth optimization): The bandwidth of the band score directly affects the ability to capture local differences and needs to be matched with the genome resolution: 40kb resolution corresponds to 200kb bandwidth (low resolution requires a larger bandwidth to cover a sufficient area to avoid missing differences); 20kb resolution corresponds to 100kb bandwidth; 10kb resolution corresponds to 50kb bandwidth (high resolution requires a smaller bandwidth to ensure the precision of difference localization). By dynamically adjusting the bandwidth, accurate capture of differential structures under different resolution data can be achieved.

[0119] In some specific embodiments, the input terminal of the diagnostic module is connected to the output terminal of the optimization module. The diagnostic module is configured to receive the differential contact area and quantitative index data, and output multiple types of visualization results and interactive diagnostic information through result feature-driven visualization type recommendation and interactive anomaly localization logic, including:

[0120] The recommended visualization types based on differential features include:

[0121] When the overlap rate of extreme percentile contact pairs between two species is less than the preset sixth setting value, it is recommended to generate an extreme value clustering diagram.

[0122] When the contact intensity difference of a certain distance bin is greater than or equal to the preset seventh setting value, it is recommended to generate a magnified map of the local deviation and mark the chromosome regions with significant differences.

[0123] It should be understood that visualization is key to interpreting cross-species Hi-C comparison results, and different differences require appropriate visualization types to be clearly displayed. This embodiment intelligently recommends visualization types based on data difference characteristics, allowing users to quickly focus on core differences. Simultaneously, by marking areas of significant difference, it improves the efficiency and accuracy of interpreting visualization results. The preset sixth setting value (extreme percentile contact pair overlap rate threshold) is 40%; the preset seventh setting value (contact intensity difference threshold) is 30%.

[0124] Visualization of extreme value differences: When the overlap rate of extreme percentile contact pairs (0-10 quantile, 90-100 quantile, corresponding to extremely strong and extremely weak contact respectively) is <40%, it indicates that the extreme contact patterns of the two species are significantly different. It is recommended to generate an extreme value clustering map: use K-means clustering (K value is automatically determined to be 3-5 through the elbow rule) to divide extreme value contact pairs into three categories: "Species A specific", "Species B specific", and "Shared", and label them with different colors. At the same time, calculate the proportion and statistical significance (p<0.05) of each category to intuitively display the species-specific extreme contact characteristics.

[0125] Distance binning difference visualization: When the contact intensity difference of a certain distance bin (e.g., 1-2Mb) is ≥30%, it indicates that there is a significant difference in the chromatin interaction pattern within that distance range. It is recommended to generate a local distance deviation magnification map: extract the contact pairs of the bin separately, sort them by chromosome coordinates, and display the contact intensity distribution of the two species. Use a t-test (p<0.05) to mark the chromosome regions with significant differences to help users accurately locate the specific location of the difference.

[0126] In some specific embodiments, the first input terminal of the interpretation module is connected to the output terminal of the diagnostic module, and the second input terminal of the interpretation module receives external multi-omics data, generating results through multi-omics association analysis and evolutionary analysis; wherein the multi-omics association analysis includes the linkage between contact differences and gene expression, and the final output of association results and hierarchical interpretation reports includes:

[0127] The interpretation module includes multi-omics association analysis, which involves associating differential contact regions with TSS regions of differentially expressed genes in RNA-seq. When the number of associated genes is greater than or equal to a preset eighth setting value, GO enrichment analysis is performed and a contact differential-gene expression association heatmap is generated.

[0128] The differential contact region is compared with the H3K27acChIP-seq reinforcing sub-region. When the overlap rate is greater than or equal to the preset ninth setting value, the correlation between contact strength and reinforcing sub-activity is calculated. When the correlation coefficient is greater than or equal to the preset correlation coefficient threshold, it is determined to be a significant correlation.

[0129] It should be understood that the biological significance of cross-species 3D genome structural differences needs to be interpreted in conjunction with multi-omics data. This embodiment reveals the association between structural differences and functional phenotypes by linking 3D genome contact differences with functional data such as gene expression and enhancer activity. At the same time, enrichment analysis and correlation tests are used to improve the biological explanatory power of the results, thus achieving a leap from "structural differences" to "functional significance".

[0130] The preset eighth setting (threshold for the number of associated genes) is 20; the preset ninth setting (threshold for overlap rate) is 30%; and the preset correlation coefficient threshold is 0.4.

[0131] Linkage between contact differences and gene expression: Differential contact regions (|PBAD|>0.5) were associated with the TSS regions (transcription start sites ±10kb, core regulatory regions) of differentially expressed genes in RNA-seq. When the number of associated genes was ≥20, it indicated a systematic association between differential contact regions and changes in gene expression. GO enrichment analysis was automatically performed (using the clusterProfiler tool) to screen GO terms related to "chromatin organization" and "transcriptional regulation," generating a contact difference-gene expression association heatmap that visually displays the changes in gene expression levels corresponding to differential contact regions, revealing the regulatory role of structural differences in gene expression.

[0132] Contact differences and enhancer activity linkage: The differential contact regions were compared with the H3K27ac ChIP-seq enhancer regions (H3K27ac is a specific marker of active enhancers). When the overlap rate was ≥30%, it was indicated that the differential contact regions may play a role by affecting enhancer activity. The Pearson correlation coefficient between the contact intensity and the H3K27ac signal intensity of each overlapping region was calculated. A correlation coefficient ≥0.4 was considered a significant association, and the significantly associated enhancer-promoter interactions were marked, providing clues to the functional mechanism of structural differences.

[0133] In some specific embodiments, the monitoring terminal of the processing module is connected to the running nodes of the aforementioned modules respectively. The processing module is configured to achieve full-process management through node parameter verification, anomaly backtracking and solution recommendation, and historical result reuse judgment. The anomaly backtracking and solution recommendation includes tool error handling, and its control signal output terminal is connected to the parameter adjustment interface of each of the aforementioned modules. When dynamically adjusting the running parameters of each module, it includes:

[0134] In the processing module, the anomaly backtracking and solution recommendation include: when the asm2link stage fails due to insufficient memory, adjusting parameters to compare in batches;

[0135] When the contact-stat stage overflows with memory, the graphs generated by chromosomes are merged.

[0136] When the tool fails to be invoked, the corresponding conda installation command or environment variable configuration method will be displayed.

[0137] It should be understood that the stability and reproducibility of the entire process are key to cross-species Hi-C analysis. By monitoring the operating status of each module in real time, intelligent backtracking and solution recommendations are made for abnormal situations. At the same time, historical results are reused to reduce redundant calculations, ensuring that the analysis process is efficient, stable and reproducible, and lowering the operating threshold for users.

[0138] Asm2link Stage Anomaly Handling: The asm2link stage is the core step in collinear strand generation. It relies on the minimap2 tool for genome alignment. If it fails due to insufficient memory, the minimap2 --split-prefix parameter is automatically adjusted to align the genome in batches (e.g., by chromosome) to avoid excessive memory consumption in a single alignment and ensure smooth strand generation.

[0139] Contact-stat stage anomaly handling: The contact-stat stage is responsible for generating diagnostic maps. When the number of contact pairs exceeds 1 million, memory overflow is likely to occur. Automatic chromosome-based chart generation (each chromosome generates a heatmap and scatter plot separately) is used to automatically generate a whole genome chart after generation. This solves the memory shortage problem without affecting the integrity of the results.

[0140] Tool call exception handling: The framework depends on bioinformatics tools such as bedtools and minimap2. If the call fails due to the tool not being installed or the environment variable not being configured, the system will automatically suggest the corresponding conda installation command (such as conda install -c bioconda bedtools) or environment variable configuration method (such as manually adding the tool installation path to the PATH environment variable) to help users quickly solve the problem and lower the technical threshold.

[0141] Historical results reuse: When a user resubmits an analysis task, the input parameters (species, resolution, collinear chain generation parameters, etc.) of the new task are compared with those of the historical task. If the parameters are completely consistent and the historical task was successfully completed, the historical intermediate results (such as .mark files and .liftContacts tables) are automatically reused, and only subsequent modules (such as metrics and visualizations) are rerun, reducing redundant calculations and improving analysis efficiency.

[0142] Research Motivation: Cross-species Hi-C comparisons can reveal evolutionary constraints and lineage-specific variations in 3D genome structure. However, existing workflows are fragile: reliance on multiple scripts, inconsistent input / output formats, and a lack of diagnostic tools make the results unreliable and difficult to reproduce.

[0143] Results: We propose Coolsecture, a Python 3-based command-line toolkit for simplifying cross-species Hi-C analysis. This framework standardizes input and output using the Cooler format (.cool / .mcool) as its core, streamlining the lift-over process for cross-species contact relationships and introducing an uncertainty-aware percentile scoring system along with corresponding diagnostic visualizations. The toolkit generates ready-to-use Cooler output files and provides compact graphical results, including percentile heatmaps, distance heatmaps, ratio scatter plots, and two-panel local difference plots.

[0144] Architecture and Installation

[0145] Coolsecture is a modular Python package defined via a pyproject.toml file. It provides a single CLI (command-line interface) entry point, Coolsecture, which dispatches commands to various subcommands, each implemented within its own module. It requires Python 3.8 or later and depends on Cooler for processing .cool / .mcool files, NumPy for vectorized numerical computations, and Matplotlib for plotting. Optional dependencies include minimap2 and paftools (for generating genome alignments and collinear strands), and Snakemake (for workflow automation).

[0146] Standardized input / output

[0147] Unlike C-InterSecture, which relies on HiC-Pro (Servant et al. 2015) or Juicer (Durand et al. 2016) output formats, Coolsecture directly reads the contact matrix from the Cooler file. The prepare subcommand accepts a .cool file or a multi-resolution .mcool::resolutions / RES path and outputs two files: (1) a rich contact table containing information for each pixel, including the chromosome of the interacting bin, start and end coordinates, global bin ID, percentile ranking, strict and loose confidence labels, coverage of each bin, and distance bins; (2) a summary table by distance bins, including sample size, quantiles, and median of the 100 percentiles. The distance bins are defined by the genomic distance (bp); interactions between chromosomes are classified into bin-1 (only included when --inter is set). If bin weights exist in the Cooler file, the interaction count is normalized by dividing by the product of the two bin weights before assigning percentiles; otherwise, the raw count is used.

[0148] Previous C-InterSecture processes categorized contacts into high, medium, and low based solely on percentile ranking, without quantifying uncertainty. Our approach calculates percentile ranking within each distance bin using the empirical cumulative distribution function (CDF). For a pixel with count c, we first calculate its Poisson confidence interval. Then, the upper and lower bounds of the interval are mapped to the reference distribution, respectively, to obtain two rankings: , When there are ties, a mid-rank is used. These labels enable subsequent analysis steps to distinguish between robust and critical signals.

[0149] Lift-over and Chain File Conversion

[0150] We implemented a three-step process for cross-species lift-over: (1) asm2link: encapsulates the genome cross-alignment of minimap2 (Li2021) and compresses the PAF alignment into six collinear strands (chromosomes, start position, and end position of the query and target species), while preserving the alignment direction. (2) link2mark: converts the strands generated by asm2link or other collinearity tools into .mark files, supporting optional block merging and reverse complementation. (3) liftcontacts: using the .mark file and the genome FASTA index, we perform actual lift-over on the prepared contacts, generating a .liftContacts table that stores percentile and distance information between the two species.

[0151] The `lift2cool` subcommand then converts the generated `.liftContacts` table into the standard Cooler format. Specifically, it outputs two cooler matrices: `Observed.cool` stores the contact counts of the original species, and `Control.cool` stores the corresponding counts of the mapped species. Both are preserved in the source genome coordinate system and represented numerically in the 0–100 percentile, thus providing a unified data structure for downstream visualization and quantitative comparisons.

[0152] Indicator Calculation and Diagnostic Visualization

[0153] The metric subcommand outputs window-based statistical results in the form of bedGraph tracks and summary charts: (1) PBAD (percentile-based background-adjusted difference), which inherits the calculation method of C-InterSecture; (2) logarithmic ratio: the formula is and in Add small constants to avoid division by zero; (3) Strip score: calculate the one-sided enrichment on both sides of the diagonal, defined as the difference (or normalized difference) between the two diagonal wings within a fixed bandwidth; (3) Correlation coefficient: calculate Pearson's r or Spearman's ρ between matching pixels. The output BedGraph can be visualized using IGV and similar tools to help locate structurally conservative / differential chromatin regions.

[0154] The `contact-stat` subcommand generates three diagnostic plots: (1) a percentile heatmap—comparing percentile scores between species; (2) a distance heatmap—showing the relationship between contact intensity and genomic distance; and (3) a ratio scatter plot—showing a scatter plot of the ratio of contact distance to contact frequency, highlighting species-specific contractions or expansions. These plots can be exported as PDF, PNG, or SVG formats and include randomized controls to assess significance. The `plot-cross` subcommand generates a two-panel heatmap that highlights corresponding portions of two species in the same region, providing a quick visual check of lift-over accuracy.

[0155] Command line interface and workflow integration

[0156] All functions are exposed via CLI, with the calling format: Coolsecture<subcommand>. This avoids running scattered scripts or manually editing INI files, and help information is uniformly processed, providing concise parameter descriptions. The entire workflow can be automated using the provided Snakemake file, sequentially calling prepare, link2mark, liftcontacts, metric, contact-stat, and lift2cool, and using predictable prefixed naming for the output results.

[0157] To demonstrate the functionality of Coolsecture, we utilized publicly available sequencing data from previous studies of *Arabidopsis suecica* (Asu) (Burns et al. 2021; Jiang et al. 2021) and *Arabidopsis thaliana* (Ath) (Feng et al. 2014; Wang et al. 2022) to generate cooler-format contact matrices for each. Each matrix was processed at 40 kb resolution using the prepare subcommand with default parameters. Collinear strands between species pairs were generated using asm2link, subsequently converted to .mark files by the link2mark subcommand, and used for execution of the liftcontacts subcommand.

[0158] The `metric` subcommand displays the windowed PBAD score distribution for Asu–Ath. The peak of the observed distribution is concentrated around 0.15, while the peak of the random distribution is skewed towards a slightly higher range. Furthermore, the observed distribution exhibits a more pronounced tail in the high-score range, while the tail of the random distribution decays more rapidly. This peak misalignment and tail morphology difference suggests a potential systematic shift in PBAD between the observed and random data, implying a spatial clustering deviation beyond the distance matching expectations.

[0159] The diagnostic results generated by contact-stat confirmed the robustness of percentile scores. Percentile heatmaps showed significant enrichment at the edges of the axes. This indicates that in Asu-Ath comparisons, extremely strong or weak contacts do not necessarily correspond to the same intensity in another species, but rather exhibit a systematic shift. Distance heatmaps showed enrichment along the diagonal, suggesting that the distance-attenuation relationship between the two species is generally comparable. However, some shifts exist, with Asu's contact span being shorter or longer than Ath's in some regions, suggesting differences in chromosome folding scales between species. Ratio scatter plots showed a butterfly distribution, with most points clustered near the origin, indicating that the distance-intensity relationship between the two species is generally similar; however, asymmetric enrichment exists in the negative region of the horizontal axis and the diagonal quadrants, suggesting that Ath has shorter distances in some cases, and that distance and intensity differences are coupled in some regions. These patterns collectively indicate that chromosome folding is conserved overall, but measurable, site-specific deviations also exist.

[0160] To visually illustrate local, distance-matching deviations, a plot-cross comparison of Asu and Ath is used, taking the sT1:0–10 Mb region of Asu as an example. In the observation-control plot, the lower triangle displays the Asu percentile (observed), while the upper triangle displays the Ath percentile (control) on the same genomic axis. We observed that local enrichments in the observation panel were weakened or absent in the control, indicating that these enrichments represent relatively non-conserved structures across species. The observation-difference plot summarizes this contrast, displaying percentile differences in the upper triangle heatmap while keeping the lower triangle unchanged, thus presenting a clearer comparison. The combined use of these two complementary views provides a diagnostic test for local structures, demonstrating Coolsecture's ability to visualize distance-matching deviations in user-specified genomic regions.

[0161] We developed Coolsecture, a modern framework for cross-species Hi-C comparisons. It improves upon earlier workflows in three key areas: (1) Standardized input and output: based on the Cooler format, it avoids reliance on specific upstream tools and generates .cool files compatible with other Hi-C tools; (2) Uncertainty-aware percentile scoring—introducing strict and lenient rankings to distinguish robust signals from marginal signals in downstream analysis; and (3) Integrated diagnostic visualization—providing percentile heatmaps, distance heatmaps, ratio scatter plots, and region cross plots, enabling users to quickly assess the plausibility of results. Testing on Asu–Ath, we demonstrated that Coolsecture generates multi-level statistical and visualization results in cross-species comparisons, revealing both conservation and variation in chromosome folding from whole-genome distribution to local regions. The results show that the tool not only captures overall conserved patterns but also pinpoints site-specific deviations beyond distance matching expectations. Overall, Coolsecture provides a more robust and reproducible approach to cross-species Hi-C analysis. It can meet the needs of beginners to get started quickly, and also provide advanced users with a scalable workflow, thereby promoting comparative research in three-dimensional genomics.

[0162] Those skilled in the art will understand that embodiments of this application can be provided as methods, systems, or computer program goods. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, this application can take the form of a computer program goods embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0163] This application is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program goods according to embodiments of this application. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart... Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0164] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0165] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0166] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit it. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the specific implementation of the present invention. Any modifications or equivalent substitutions that do not depart from the spirit and scope of the present invention should be covered within the scope of protection of the claims of the present invention.

Claims

1. An easy-to-use and improved cross-species Hi-C contact map comparison framework, characterized in that, include: The data processing module is configured to receive input Hi-C data files and perform data preprocessing through resolution consistency verification, bin coverage grading, and interchromosomal interaction adaptation. The resolution consistency check includes dynamic selection of downsampling algorithm, the bin coverage classification processing includes edge region weight adjustment and low coverage region filtering logic, the chromosome interaction adaptation includes proportion judgment and classification labeling, and finally outputs a standardized contact data table. A correction module, whose input is connected to the output of the data processing module, is configured to receive the standardized contact data table and complete the mapping of homologous segments through collinear chain confidence grading, chain fragmentation merging processing, and real-time verification of mapping results. The collinear chain confidence grading is based on three-level classification of identity and coverage and assigned weights. The chain fragmentation merging processing includes interval and direction judgment. The real-time verification of mapping results includes coordinate validity and rejection ratio control. Finally, the module outputs mapped contact pair data with confidence levels. An optimization module, whose input is connected to the output of the correction module, is configured to receive the mapped contact pair data with confidence levels, perform comparative index calculation and parameter optimization based on logical judgment driven by research objectives and data characteristics, and output differential contact area and quantitative index data. The diagnostic module has its input connected to the output of the optimization module. The diagnostic module is configured to receive the differential contact area and quantitative index data, and output multiple types of visualization results and interactive diagnostic information through result feature-driven visualization type recommendation and interactive anomaly localization logic. The interpretation module has a first input terminal connected to the output terminal of the diagnostic module, and a second input terminal of the interpretation module receives external multi-omics data and generates results through multi-omics association analysis and evolutionary analysis. The multi-omics association analysis includes the linkage between exposure differences and gene expression, and finally outputs the association results and hierarchical interpretation reports. The processing module has its monitoring terminals connected to the running nodes of the above modules respectively. The processing module is configured to achieve full-process management through node parameter verification, anomaly backtracking and solution recommendation, and historical result reuse judgment. The anomaly backtracking and solution recommendation includes tool error handling, and its control signal output terminals are connected to the parameter adjustment interfaces of the above modules respectively, for dynamically adjusting the running parameters of each module. The three-level classification and weight assignment of the collinear chain confidence level include: When the identity is greater than or equal to the preset first setting value, and the coverage is greater than or equal to the preset first coverage setting value, it is determined to be a high confidence chain and assigned the first preset weight; When the identity is greater than or equal to the preset second setting value and less than the first setting value, or the coverage is greater than or equal to the preset second coverage setting value and less than the first coverage setting value, it is determined to be a middle confidence chain and assigned a second preset weight. When the identity is less than the second set value and the coverage is less than the second coverage set value, it is determined to be a low-confidence chain and assigned a third preset weight; In the optimization module, the preferred selection of comparison indicators based on the research objective includes: when the research objective is a conservative analysis, the correlation coefficient between PBAD and Pearson is calculated first. When the research objective is to locate differential structures, log ratio and strip score should be calculated first. The data feature-driven parameter self-optimization includes: when the noise level is such that the standard deviation of the contact count in the background area is greater than or equal to the preset fifth setting value, the logarithmic ratio is Gaussian smoothed, and the smoothing kernel size is adjusted according to the noise distribution, with a 3x3 kernel for concentrated noise and a 5x5 kernel for dispersed noise; The bandwidth of the band score is adjusted with the genome resolution: 40kb resolution corresponds to 200kb bandwidth, 20kb resolution corresponds to 100kb bandwidth, and 10kb resolution corresponds to 50kb bandwidth.

2. The easy-to-use and improved cross-species Hi-C contact map comparison framework according to claim 1, characterized in that, The data processing module is configured to receive input Hi-C data files and perform data preprocessing through resolution consistency verification, bin coverage grading, and inter-chromosome interaction adaptation. Resolution consistency verification includes dynamic selection of downsampling algorithms; bin coverage grading includes edge region weight adjustment and low-coverage region filtering logic; and inter-chromosome interaction adaptation includes proportion judgment and classification labeling. The final output of the standardized contact data table includes: The dynamic selection of the resolution consistency check includes: When the resolution difference of the input samples is less than or equal to a preset multiple, and the coefficient of variation of the contact count distribution is less than a preset first coefficient threshold, mean aggregation is used for downsampling. When the resolution difference of the input samples is less than or equal to a preset multiple, and the coefficient of variation of the contact count distribution is greater than or equal to the first coefficient threshold, median aggregation is used for downsampling. When the resolution difference is greater than the preset multiple, and the Pearson correlation coefficient of the coverage of adjacent bins is greater than or equal to the preset second coefficient threshold, linear interpolation is used to complete the count of high-resolution bins.

3. The easy-to-use and improved cross-species Hi-C contact map comparison framework according to claim 2, characterized in that, The data processing module is configured to receive input Hi-C data files and perform data preprocessing through resolution consistency verification, bin coverage grading, and inter-chromosome interaction adaptation. Resolution consistency verification includes dynamic selection of downsampling algorithms; bin coverage grading includes edge region weight adjustment and low-coverage region filtering logic; inter-chromosome interaction adaptation includes proportion judgment and classification labeling; and when finally outputting a standardized contact data table, it also includes: The edge region weight adjustment in the bin coverage hierarchical processing includes: When the coverage of a bin is within the preset coverage threshold range and is located in the middle of a chromosome arm, K-nearest neighbor interpolation is used to complete the missing count, and the K value is adjusted according to the average coverage of the surrounding bins. When the average coverage is greater than the preset coverage range, K is 3; When the average coverage is lower than the preset coverage range, K is 5; When it is in the centromere or telomere region, it is marked as an edge bin and the weight is set to the coverage percentage.

4. The easy-to-use and improved cross-species Hi-C contact map comparison framework according to claim 3, characterized in that, The input of the correction module is connected to the output of the data processing module. The correction module is configured to receive the standardized contact data table and complete the mapping of homologous segments through collinear chain confidence grading, chain fragmentation merging, and real-time verification of mapping results. The collinear chain confidence grading is based on a three-level classification of identity and coverage, with assigned weights. The chain fragmentation merging includes interval and direction determination. The real-time verification of mapping results includes coordinate validity and rejection ratio control. When finally outputting mapped contact pair data with confidence levels, the module also includes: In the correction module, the interval and direction determination of the chain fragmentation merging process includes: for short chains with a length less than a preset third setting value; When the interval between the chain and the adjacent chain is less than the preset fourth setting value and the comparison direction is consistent, they are merged into a long chain and the original coordinates are recorded. If the interval with an adjacent chain is greater than or equal to the fourth setting value or the comparison direction is opposite, it is retained as a fragmented chain and is only allowed to be associated with a high-confidence chain for mapping.

5. The easy-to-use and improved cross-species Hi-C contact map comparison framework according to claim 4, characterized in that, The input terminal of the diagnostic module is connected to the output terminal of the optimization module. The diagnostic module is configured to receive the differential contact area and quantitative index data, and output multiple types of visualization results and interactive diagnostic information through result feature-driven visualization type recommendation and interactive anomaly localization logic, including: Visualization types based on differential features include: When the overlap rate of extreme percentile contact pairs between two species is less than the preset sixth setting value, it is recommended to generate an extreme value clustering diagram. When the contact intensity difference of a certain distance bin is greater than or equal to the preset seventh setting value, it is recommended to generate a magnified map of the local deviation and mark the chromosome regions with significant differences.

6. The easy-to-use and improved cross-species Hi-C contact map comparison framework according to claim 5, characterized in that, The first input terminal of the interpretation module is connected to the output terminal of the diagnostic module, and the second input terminal of the interpretation module receives external multi-omics data and generates results through multi-omics association analysis and evolutionary analysis. The multi-omics association analysis includes the linkage between exposure differences and gene expression. The final output of the association results and hierarchical interpretation report includes: The interpretation module includes multi-omics association analysis, which involves associating differential contact regions with TSS regions of differentially expressed genes in RNA-seq. When the number of associated genes is greater than or equal to a preset eighth setting value, GO enrichment analysis is performed and a contact differential-gene expression association heatmap is generated. The differential contact region is compared with the H3K27acChIP-seq reinforcing sub-region. When the overlap rate is greater than or equal to the preset ninth setting value, the correlation between contact strength and reinforcing sub-activity is calculated. When the correlation coefficient is greater than or equal to the preset correlation coefficient threshold, it is determined to be a significant correlation.

7. The easy-to-use and improved cross-species Hi-C contact map comparison framework according to claim 6, characterized in that, The monitoring terminals of the processing modules are respectively connected to the running nodes of the aforementioned modules. The processing modules are configured to achieve full-process management through node parameter verification, anomaly backtracking and solution recommendation, and historical result reuse judgment. The anomaly backtracking and solution recommendation includes tool error handling, and its control signal output terminals are respectively connected to the parameter adjustment interfaces of the aforementioned modules. When dynamically adjusting the running parameters of each module, this includes: In the processing module, the anomaly backtracking and solution recommendation include: when the asm2link stage fails due to insufficient memory, adjusting parameters to compare in batches; When the contact-stat stage overflows with memory, the graphs generated by chromosomes are merged. When the tool fails to be invoked, the corresponding conda installation command or environment variable configuration method will be displayed.

Images