Group epigenetic analysis method and system based on consensus peak recognition

Abstract

The invention provides a group epigenetic analysis method and system based on consensus peak recognition, and belongs to the technical field of bioinformatics, and the method comprises the steps: S1, recognizing a consensus peak from a multi-sample peak file based on a position offset threshold and a minimum sample proportion; s2, calculating the coverage depth of a consensus peak based on the segment / read segment to obtain an original peak intensity matrix; s3, eliminating GC preference through GC content inter-partition local correction; s4, performing cross-sample normalization by adopting methods such as DESeq2, TMM and the like; and S5, automatically removing a hidden batch effect by utilizing potential variable analysis. The system correspondingly comprises five functional modules. According to the method, the problems of inconsistent peak coordinates, incomparable signals and large technical deviation interference in multi-sample epigenetic data are solved through an integrated process, full-automatic generation from original data to a high-quality and high-comparability peak intensity matrix is realized, and the analysis efficiency and accuracy of researches on large-scale groups ChIP-seq, ATAC-seq and the like are remarkably improved.

Description

Technical Field

[0001] This invention relates to the field of bioinformatics technology, and in particular to a population epigenetic analysis method and system based on consensus peak identification. Background Technology

[0002] With the development of high-throughput sequencing technology, population epigenetics studies based on technologies such as ChIP-seq and ATAC-seq are increasing, aiming to elucidate the variation patterns of chromatin states in a population and their association with traits. The core of such research lies in constructing a cross-sample comparable and accurate "peak." The "sample" matrix. However, existing technical processes face a series of challenges when processing multi-sample data.

[0003] First, due to factors such as experimental noise and alignment differences, the start and end coordinates of independently identified enriched peaks on the genome often shift, making it difficult to directly establish peak correspondences across samples. Simple coordinate merging strategies cannot accurately identify true common peaks.

[0004] Secondly, the quantification of peak intensity (such as coverage) is systematically interfered with by various technical factors, including sequencing depth, library type (single-end SE or paired-end PE), GC content bias (GC-bias), and the batch effect. Uncorrected peak intensity matrices cannot reflect true biological differences, seriously affecting the accuracy of subsequent differential analysis, clustering, and epigenetic quantitative trait locus (QTL) mapping.

[0005] Existing tools often focus on single steps in the data processing workflow (such as peak detection, normalization, or batch correction), lacking an integrated, automated solution that covers everything from defining common peaks, accurate signal quantization, and systematic bias correction to batch effect removal. Researchers typically need to manually connect multiple independent tools, resulting in a cumbersome, error-prone process that makes it difficult to ensure consistency of parameters and logic between steps, failing to meet the data quality and analytical efficiency requirements of large-scale population studies.

[0006] Therefore, there is an urgent need to develop an integrated analysis method and system that can systematically solve the above problems. Summary of the Invention

[0007] The present invention aims to solve at least one of the technical problems existing in the prior art, and proposes a population epigenetic analysis method and system based on consensus peak identification.

[0008] In a first aspect, embodiments of the present invention provide a population epigenetic analysis method based on consensus peak identification, comprising:

[0009] S1. Read the enriched peak region files of each sample; scan the peak regions after sorting by chromosome, judge and merge peaks belonging to the same enriched region according to the preset position offset threshold; filter out the regions that appear stably within the population based on the preset minimum sample ratio, and generate a consensus peak set and a sample list in which each consensus peak appears.

[0010] S2. For each region in the consensus peak set, read the sequence alignment file of the corresponding sample; calculate the coverage depth of the region based on the fragment or read according to the sequencing library type; calculate the intensity index of each sample in each consensus peak based on the coverage depth to form the original peak intensity matrix.

[0011] S3. Extract the genome sequence of each consensus peak region and calculate its GC content; divide the GC content into multiple intervals; generate a correction factor based on the relationship between the local statistics and global statistics of the peak intensity in each GC content interval, and correct the original peak intensity matrix to obtain the GC-corrected matrix;

[0012] S4. The GC-corrected matrix is ​​scaled using a selected normalization method to obtain the normalized peak intensity matrix.

[0013] S5. Perform data transformation on the normalized peak intensity matrix; estimate the number of latent variables with hidden batch effects; remove the influence represented by the latent variables through a regression model; perform inverse transformation and post-processing to output the final corrected peak intensity matrix.

[0014] Furthermore, in step S1, the position offset threshold ranges from 50 to 150 base pairs; the minimum sample ratio ranges from 0.1 to 0.3.

[0015] Furthermore, in step S2, the intensity index includes at least one of maximum coverage, average coverage, median coverage, and peak area; for paired-end sequencing data, the fragment coverage interval is calculated based on the template length of the paired reads; for single-end sequencing data, the coverage interval is calculated based on the read length.

[0016] Furthermore, in step S3, the GC content is uniformly divided into 50 to 200 intervals; the correction factor is determined based on the ratio of the local mean of the peak intensity in each GC interval to the global mean of all peak intensities.

[0017] Furthermore, in step S4, the normalization method includes at least one of the following: DESeq2 size factor method, TMM method, CPM method, TPM method, and upper quartile method.

[0018] Furthermore, step S5 includes: performing an asinh transformation on the data matrix; automatically estimating the number of latent variables using latent variable analysis; constructing a linear model containing a design matrix of known biological variables and the latent variables for regression correction; and performing an inverse sinh transformation and minimum value shifting on the residuals.

[0019] Furthermore, in step S1, for each finally determined consensus peak, its start and end coordinates are the average of the start and end coordinates of all merged sample peak regions.

[0020] Furthermore, the method also includes step S6: visualizing the quality assessment of the final corrected peak intensity matrix, including generating at least one of a peak intensity distribution map, a heatmap of inter-sample correlation, a principal component analysis map, and a MA map.

[0021] The technical advantages of the population epigenetic analysis method based on consensus peak identification disclosed in this invention are as follows: This method, through an innovative algorithm, achieves fully automated processing from multiple raw sample data to a high-quality, comparable peak intensity matrix. It can automatically and robustly identify a high-quality consensus peak set from hundreds of samples. Through parameterized control, users can flexibly adjust the degree of "sharedness" of the consensus peaks according to research needs (pursuing sensitivity or specificity). This method solves the problems of blurred peak boundaries and high false positive / false negative rates caused by manual merging or simple overlap methods, providing a precise structural basis for subsequent quantitative comparison of all samples in a unified coordinate system.

[0022] Secondly, the present invention also provides a population epigenetic analysis system for implementing the method, comprising:

[0023] The consensus peak identification module is used to execute step S1 to generate a consensus peak set and a sample occurrence list; the peak intensity calculation module is used to execute step S2 to generate the original peak intensity matrix.

[0024] The GC correction module is used to perform step S3 and generate the GC-corrected matrix;

[0025] The normalization module is used to perform step S4, providing multiple normalization methods and outputting a normalization matrix;

[0026] The batch effect correction module is used to perform step S5, remove implicit batch effects, and output the final matrix.

[0027] Furthermore, it also includes a visualization module for performing step S6 and generating a data quality assessment chart. Attached Figure Description

[0028] Figure 1A flowchart illustrating a population epigenetic analysis method based on consensus peak identification, provided in an embodiment of the present invention;

[0029] Figure 2 A schematic diagram of a population epigenetic analysis system framework based on consensus peak identification is provided in an embodiment of the present invention;

[0030] Figure 3 A schematic diagram illustrating the framework configuration of the various modules and optional methods provided in the embodiments of the present invention;

[0031] Figure 4 This is a comparative diagram of the GC correction module before and after processing provided in an embodiment of the present invention;

[0032] Figure 5 This is a comparative schematic diagram of uniformity correction before and after the present invention. Figure 6 This is a schematic diagram comparing the standardization effect of GC content provided in an embodiment of the present invention. Detailed Implementation

[0033] To enable those skilled in the art to better understand the technical solutions of the present invention, exemplary embodiments of the present invention are described below in conjunction with the accompanying drawings, including various details of the embodiments of the present invention to aid understanding. These should be considered merely exemplary. Therefore, those skilled in the art should recognize that various changes and modifications can be made to the embodiments described herein without departing from the scope and spirit of the present invention. Similarly, for clarity and brevity, descriptions of well-known functions and structures are omitted in the following description.

[0034] A population epigenetic analysis method and system based on consensus peak identification, referenced Figures 1 to 5 As shown, the specific steps include:

[0035] S1. Consensus Peak Identification: Read the enriched peak region files of each sample; scan the peak regions after sorting by chromosome, judge and merge peaks that may belong to the same enriched region based on the preset position offset threshold; filter out regions that appear stably within the population based on the preset minimum sample ratio, and generate a consensus peak set and a sample list of each consensus peak.

[0036] Read the narrowPeak file for each sample; scan the peaks after sorting by chromosome; determine the peaks that can be merged based on the position offset threshold (max_position_shift); output the common peaks that meet the minimum sample ratio (min_sample_fraction); generate the start and end coordinates of consensus and a sample list for each common peak.

[0037] The algorithm mainly comprises an input parsing unit, a sorting and scanning unit, a position offset judgment unit, a sample ratio screening unit, and a consensus peak generation unit, forming a logical flow that can run automatically in a computing environment. First, the input parsing unit reads the narrowPeak file of each sample from the storage system, extracts the chromosome number, start and end coordinates, and related quantitative information for each peak, and stores these peak data in a structured list according to chromosome division to ensure partition independence in subsequent processing. Then, the sorting and scanning unit sorts all peak objects under the same chromosome according to their starting coordinates and processes each peak sequentially using a linear scanning method to ensure controllable time complexity of the algorithm. During the scanning process, the position offset judgment unit detects the coordinate difference between the current peak and its neighboring peaks, and determines whether two peaks might originate from the same true enrichment region based on a preset position offset threshold `max_position_shift`. If the condition is met, they are classified into the same candidate region. Next, the sample ratio screening unit counts the number of samples containing peaks in the candidate region and compares it with the total number of samples. If the ratio exceeds a preset `min_sample_fraction`, the region is identified as a shared peak; otherwise, it is discarded. Finally, the consensus peak generation unit calculates the average value of the peak positions of each sample in the qualified region, forms the start and end coordinates of the consensus, and records a list of all samples participating in the consensus peak, generating a structured consensus peak output object. The entire module relies on the process described in the above logic diagram, enabling the program to systematically merge multi-sample peak signals into a structurally unified whole, forming the standardized input required for subsequent peak intensity calculation and normalization processing.

[0038] Implementation Process: The narrowPeak files of 533 ATAC-seq samples with open chromatin from rice were collected as input. The consensus peak identification module (Script1) of this invention was run. The module first reads all files and sorts the peak regions by chromosome. The position shift threshold (max_position_shift) was set to 100 base pairs (bp), meaning that during the scan, if the starting positions of two peak regions differ by less than 100 bp, they are considered to originate from the same region and are grouped into the same candidate cluster. After traversing all samples, each candidate cluster was screened: the minimum sample ratio (min_sample_fraction) was set to 0.2, meaning that a candidate cluster must appear in at least 533 samples. Only among approximately 107 samples (approximately 0.2) can a consensus peak be considered valid. For a cluster that passes the test, the average of the start and end coordinates of all peaks within the cluster is calculated to generate the final start and end coordinates (chrom, start, end) of that consensus peak, and a list of samples containing that peak (samples_list) is recorded.

[0039] To verify parameter sensitivity, tests were conducted with min_sample_fraction values ​​of 0.1, 0.2, and 0.3. At a fraction of 0.1, approximately 150,000 consensus peaks were identified; at 0.2, approximately 120,000 peaks were identified; and at 0.3, approximately 90,000 peaks were identified. Simultaneously, max_position_shift values ​​of 50bp and 150bp were tested. The results show that the algorithm runs stably under different parameters, and the number of consensus peaks decreases as min_sample_fraction increases, demonstrating the rigor of the selection process.

[0040] This step demonstrates that the method of this invention can automatically and robustly identify a high-quality consensus peak set from hundreds of samples. Through parameterized control, users can flexibly adjust the degree of "sharedness" of the consensus peaks according to research needs (pursuing sensitivity or specificity). This method solves the problems of blurred peak boundaries and high false positive / false negative rates caused by manual merging or simple overlap methods, providing a precise structural basis for subsequent quantitative comparison of all samples in a unified coordinate system.

[0041] S2. Precise peak intensity calculation: For each region in the consensus peak set, read the sequence alignment file of the corresponding sample; calculate the coverage depth of the region based on the fragment or read according to the sequencing library type; calculate the intensity index of each sample in each consensus peak based on the coverage depth to form the original peak intensity matrix.

[0042] For each shared peak region:

[0043] For PE libraries: Pair reads and calculate fragment coverage based on template length. For SE libraries: Take the read coverage interval; calculate multiple metrics: maximum coverage, mean, peak area, peak position, CV, etc. The output is "Common peaks". The rawpeakheight matrix of the "sample".

[0044] Its peak coverage calculation module consists of a library type identification unit, a fragment construction unit, a coverage accumulation unit, and a peak index calculation unit, which can be executed in parallel via CPU multi-threading on standard computing devices. First, the library type identification unit determines whether the library type of the sample corresponding to the common peak is PE or SE based on the metadata of the input file, and selects different processing paths accordingly. In the PE path, the fragment construction unit performs name matching on paired reads read from the BAM file, filters out secondary and supplementary alignments, calculates the actual coverage interval of the DNA fragment based on the template length of the read pair, and maps the start and end positions of the fragments to the local coordinate system of the common peak. In the SE path, the system directly uses the alignment start position and length of a single read to determine the coverage interval. Subsequently, the coverage accumulation unit constructs an integer array within the peak length range, performs an accumulation operation of incrementing the entire segment for each fragment or read-covered position to obtain the base-by-base coverage depth within the peak region. Next, the peak index calculation unit calculates a series of quantitative indicators based on the coverage array, including maximum coverage, average coverage, median coverage, peak area (sum of coverage), peak position (point of maximum coverage), and coefficient of variation (CV), and combines these with the peak coordinates and sample identifiers to form a structured record. The calculation results of all samples are summarized into a "shared peak" list in the order of common peaks. The original peak height matrix of the "sample" provides normalized input for subsequent GC correction, global normalization, and batch effect processing. This module ensures that peak intensity quantization matches different library structures through fragment-level calculation, thereby achieving accurate peak-level signal modeling for ChIP-seq / ATAC-seq data.

[0045] Implementation Process: The consensus peak set generated in S1 is processed using the peak coverage calculation module (Script2) of this invention. The module inputs are the consensus peak list and the BAM file for each sample. For each consensus peak and each sample, the module first determines that the sample's BAM file is a paired-end (PE) library. For PE data, the module reads the read pairs aligned to the peak region, filtering out secondary and supplementary alignments. Key Step: Based on the alignment start position and template length of each read pair, the module calculates the precise coverage region of the original DNA fragment it represents on the genome, rather than simply calculating the coverage of the read itself. Then, the coverage regions of all fragments are summed to obtain a fragment-based coverage depth profile for that region. Based on this profile, multiple intensity indicators for the sample in this peak are calculated: maximum coverage (max_cov), average coverage (mean_cov), and peak area (peak_area, i.e., the sum of coverage depths).

[0046] In comparison, the same peak area was calculated using a traditional read count-based method. 1000 high-confidence open regions, validated through other experiments, were selected as the gold standard. The Pearson correlation coefficients between the peak areas obtained by both methods and the signal strength estimated by the gold standard were calculated. The results show that the correlation coefficient of the traditional read count method is lower. The correlation coefficient was 0.61, while the correlation coefficient r of the fragment-level calculation method of this invention was improved to 0.80. In addition, examination of the coverage curve revealed that the curve obtained by the method of this invention is smoother at the peak center and decays more naturally at the boundary, reducing the jagged noise introduced by the read segment boundary.

[0047] This embodiment demonstrates that the fragment-level coverage calculation method used in this invention can significantly improve the accuracy of peak intensity quantification. By utilizing paired information from PE data to reconstruct the true fragments, this method is more consistent with the physical experimental nature of ATAC-seq / ChIP-seq, effectively avoiding biases introduced by single-end counting or erroneous fragment estimation, and providing more reliable raw data input for subsequent GC correction and normalization.

[0048] S3. GC Content Bias Correction: Extract the genomic sequence of each consensus peak region and calculate its GC content; divide the GC content into multiple intervals; generate a correction factor based on the relationship between the local and global statistics of peak intensity within each GC content interval, and correct the original peak intensity matrix to obtain the GC-corrected matrix. The GC correction module can effectively eliminate sequence-specific bias.

[0049] Sequences of each peak were extracted from the reference genome; GC content was calculated; GC content was divided into 100 intervals; normalization factor was calculated based on the local mean of each GC group; and GC preference correction was applied to the peak height.

[0050] Its GC correction module consists of a sequence extraction unit, a GC calculation unit, a binning and classification unit, a local mean modeling unit, and a peak correction unit. These units operate collaboratively in a sequential workflow within the computing device. First, the sequence extraction unit reads the corresponding base sequences from the reference genome index file segment by segment based on the chromosome number and start / end coordinates of the common peaks. This unit relies on the FASTA index structure and random access functionality to ensure rapid parsing of large-scale peak sets. Subsequently, the GC calculation unit traverses the sequences to count the total number of G and C values, and divides this by the total sequence length to obtain the GC ratio of the peak region. The resulting GC ratio is input as a continuous variable to the binning and classification unit. In this unit, the system maps the GC content of all peaks to 100 equally wide intervals according to preset interval division rules, clustering peaks with similar GC values ​​into the same group. Next, the local mean modeling unit calculates the mean of peak heights (calculated by the coverage module) within each GC group to obtain the local expected coverage level for that GC interval. The system then calculates the overall mean of all peak heights and constructs a normalization factor for that group using the ratio of the overall mean to the group mean. Finally, the peak correction unit multiplies each peak height by its corresponding GC group's normalization factor to form the corrected peak intensity matrix. This module utilizes the statistical dependency between GC content and sequencing coverage, and approximates and eliminates the systematic bias caused by GC preference at the algorithmic level through bucket modeling, making peak intensity independent of sequence composition. This provides structurally consistent numerical input for subsequent normalization and batch correction modules.

[0051] The sequence of each consensus peak was extracted from the rice reference genome, and the GC content was calculated. The GC content range (0% to 100%) was uniformly divided into 100 intervals (buckets). For the original peak intensity matrix before GC correction (from Example 2), the mean of all peak intensities (taking mean_cov as an example) within each GC interval was calculated as the "local observed mean" for that GC interval. Simultaneously, the global mean of all peak intensities was calculated. The correction factor for each GC interval was defined as global mean / local observed mean. The intensity of each peak was multiplied by the correction factor of its corresponding GC interval to obtain the GC-corrected intensity.

[0052] S4. Cross-sample normalization: The GC-corrected matrix is ​​scaled using a selected normalization method to obtain the normalized peak intensity matrix.

[0053] Depending on the user's selected method: use the size factor of DESeq2; or the TMM factor of edgeR; or the normalization strategy such as CPM / TPM / UQ / CQN; output the normalized peak height matrix.

[0054] Its peak height normalization module consists of a method selection unit, a factor calculation unit, a numerical transformation unit, and an output matrix generation unit. Each unit executes sequentially to ensure the uniformity of peak intensities from different sources on a numerical scale. First, the method selection unit automatically determines the normalization path based on the user-input parameter `method`, including different statistical models such as DESeq2, TMM, CPM, TPM, UQ, or CQN. Then, the factor calculation unit performs the corresponding mathematical operations according to the selected method: in the DESeq2 path, this unit calculates the geometric mean across samples for each peak and obtains the size factor for each sample using the median ratio method; in the TMM path, the system trims and averages the data based on the M and A values ​​of the reference and target samples to obtain a correction factor for library size differences; in the CPM or TPM paths, the system calculates the normalization ratio based on the total count or peak length, respectively; in the UQ path, the 75th percentile of the sample peak intensity is used as the scaling factor; and in the CQN path, a conditional quantile model based on peak length or sequence attributes is constructed. Next, the numerical transformation unit divides each sample column in the original peak height matrix by its corresponding normalization factor. If it is in TPM form, it is further scaled by a factor of millions. Finally, the output matrix generation unit re-integrates the normalized values ​​with the peak position information to form a structured peak. The sample matrix is ​​used for subsequent batch correction or differential analysis. The above modules automate the execution of various normalized mathematical models through a unified interface and logical flow, enabling peak intensity data under different experimental conditions to be mapped to a comparable numerical space, thereby meeting the needs of large-scale epigenetic data processing.

[0055] For the GC-corrected matrix, cross-sample normalization was performed using the DESeq2 size factor method and the TMM method, respectively. The DESeq2 method calculates a size factor for each sample to align the median ratio of peak intensities between samples. The TMM method calculates a scaling factor by trimming the mean. Performance evaluation: Before GC correction, the correlation coefficients between all peak intensities and their corresponding GC contents were calculated. On three test sets of different sizes (100, 300, and 533 samples), the correlation coefficients were 0.34, 0.33, and 0.32, respectively, indicating a significant GC bias. After GC correction, the correlation coefficients decreased to 0.02, 0.01, and 0.00, respectively, indicating that the GC bias was effectively eliminated. For normalization, the two samples A and B with the largest difference in sequencing depth were selected; before normalization, their median peak intensity ratio was 2.1:1. After DESeq2 normalization, the ratio becomes 1.05:1; after TMM normalization, it becomes 1.08:1, and the scale difference between samples is successfully corrected.

[0056] The normalization module robustly corrects for technical scale differences between samples. These steps are crucial for ensuring the comparability of peak intensities among samples from different experimental conditions and sequencing depths, laying a numerical foundation for subsequent detection of real biological differences.

[0057] S5. Removal of latent batch effects: Perform data transformation on the normalized peak intensity matrix; estimate the number of latent batch effect variables; remove the influence represented by the latent variables through a regression model; perform inverse transformation and post-processing to output the final corrected peak intensity matrix.

[0058] Construct the model matrix mod; perform an asinh transformation on the peak heights; automatically estimate the number of latent variables (n_sv) using Python-SVA; regress the latent variables and eliminate batch effects; perform inverse transformation and minimum normalization. Output the final batch-corrected peak height matrix.

[0059] The batch effect correction process consists of a model matrix construction unit, a data transformation unit, a latent variable estimation unit, a batch regression correction unit, and an inverse transformation normalization unit. Each unit works collaboratively with the numerical optimization process through matrix operations within the computing device. First, the model matrix construction unit generates a design matrix `mod` based on whether the input samples contain phenotypic or fixed effect information. This matrix typically includes an intercept term and necessary biological variables to ensure that the subsequent regression process can distinguish between the target signal and non-target biases. Subsequently, the data transformation unit performs an asinh transformation on the peak height matrix, maintaining approximately linearity in the low-value region and achieving stable compression in the high-value region, avoiding numerical instability caused by directly using the original non-negative discrete distribution. Based on this, the latent variable estimation unit uses the Python-SVA method to perform eigenvalue decomposition on the residual matrix obtained after removing the portion explained by `mod` from the peak matrix, and automatically determines the number of latent variables `n_sv` using permutation tests or statistical thresholds, thereby constructing a latent variable matrix characterizing unknown batch effects. The batch regression correction unit then performs joint regression on the model matrix and the latent variable matrix after the asinh transformation, calculating and removing the regression components corresponding to the latent variables, so that the matrix retains only the corrected biologically relevant components. Finally, the inverse transformation normalization unit applies an inverse sinh transformation to the corrected matrix, mapping the data back to the original dimensions, and using minimum value normalization to ensure that all values ​​are positive when negative values ​​occur. After all processing is complete, the program outputs a peak × sample matrix after batch effect elimination, combining the corrected peak height and peak location information, providing consistent input for subsequent multi-sample comparisons.

[0060] Implementation Process: Using the matrix normalized by DESeq2 as input, the batch effect correction module (Script6) is run. First, an asinh transformation is performed on the matrix to stabilize the variance. Then, the SVA (Surrogate Variable Analysis) method implemented in Python is used: given the design matrix (mod) for biological grouping of the samples, the algorithm performs singular value decomposition on the residual matrix and automatically estimates the number of latent variables (surrogate variables) n_sv through permutation tests. In this example, with a dataset of 533 samples, n_sv is automatically estimated to be 7. These 7 latent variables are used as covariates, and together with the original biological variables, a linear model is constructed. Regression is used to remove the variation represented by the latent variables from the asinh-transformed data. Finally, an inverse sinh transformation is performed, and the entire matrix is ​​shifted to its minimum value to ensure all values ​​are positive, resulting in the final corrected peak intensity matrix.

[0061] Quality Assessment and Downstream Applications: 1) PCA Analysis: Before correction, the PCA plot showed that samples were clearly clustered by experimental batch; after correction, samples were mainly clustered by biological type (e.g., different tissues), and batch effects were eliminated. 2) Sample Correlation Heatmap: After correction, samples within the same biological group showed higher correlation, while the abnormally high correlation between different batches disappeared. 3) Downstream Difference Analysis: Using the corrected matrix for peak analysis (t-test) between two groups of samples, the number of detected peaks was more reasonable compared to using the uncorrected matrix, and the enriched biological pathways were more relevant with less noise. 4) Machine Learning Prediction: Using the corrected matrix as features to predict a certain trait of the samples (e.g., drought resistance), the prediction accuracy was improved by approximately 15% compared to using the uncorrected or partially corrected matrix using 5-fold cross-validation.

[0062] This step demonstrates that the SVA batch effect correction module integrated in this invention can automatically identify and remove hidden technical variations in data, significantly improving data purity and the efficiency of downstream statistical analysis and machine learning modeling. The final output peak intensity matrix is ​​a high-quality foundational data layer that can be directly used for epigenetic analysis of various populations.

[0063] S6, Quality Control Visualization (Script7).

[0064] Peak height distribution violin plot, correlation heatmap, PCA, MA plot, help to check the normalization effect and outliers.

[0065] Its quality control visualization module consists of a distribution visualization unit, a sample correlation analysis unit, a dimensionality reduction analysis unit, and a two-sample comparison unit. Each unit executes sequentially after data normalization and batch correction, used to detect the stability, consistency, and outliers of the results in a computer environment. First, the distribution visualization unit reads the normalized peak height matrix and constructs a violin plot based on a statistical plotting program, using kernel density estimation to display the overall distribution of peak intensities for each sample. This unit determines whether there are scale differences or long-tailed anomalies by comparing the distribution shapes of different samples; its generation process includes: density estimation for each column of values, generating a symmetrical distribution graph, and mapping the coordinate axes. Subsequently, the sample correlation analysis unit constructs a correlation matrix based on the Pearson or Spearman correlation coefficients of peak intensities between samples and plots it as a heatmap. This unit uses standard matrix operations and color mapping rules to make the correlation levels visually easily distinguishable, providing a basis for identifying outlier samples or insufficiently corrected batch differences. Next, the dimensionality reduction analysis unit performs PCA, compressing the information of the samples in the multidimensional peak space into the coordinates of the first two principal components, and plots a scatter plot to show the overall clustering relationship between the samples. This is based on the eigenvalue decomposition of the covariance matrix, which visually reflects the remaining systematic bias in the graph. Finally, the two-sample comparison unit constructs an MA plot, performing scatter analysis on the log ratio (M) of the peak intensities and the mean intensity (A) between the two samples to detect whether systematic bias still exists after normalization. These visualization units, through statistical modeling and graphical mapping, enable users to systematically evaluate the normalization quality and potential anomalies within a logical graph framework, forming a complete data quality audit process.

[0066] The beneficial effects of this invention are as follows:

[0067] High-precision consensus peak identification: Through the dual mechanisms of "position offset threshold" and "sample ratio screening", it can robustly identify the true shared enrichment regions across samples, effectively overcoming the problem of inconsistent peak coordinates between samples, and laying the foundation for building a unified analysis framework.

[0068] Precise signal quantification: Based on fragment-level coverage calculations, it distinguishes between PE and SE libraries, providing peak intensity estimates that more closely approximate the distribution of real DNA fragments, thus improving quantification accuracy.

[0069] Systematic bias correction: It integrates local GC content correction, multiple robust normalization methods, and automated batch effect removal of latent variables, forming a coherent correction chain that minimizes technical noise and enables the peak intensity matrix to truly reflect biological variation.

[0070] Full-process automation and high throughput: The above steps are integrated into a complete system process, supporting parallel computing and capable of efficiently processing large-scale population data of tens to hundreds of samples, significantly improving analysis efficiency and reproducibility.

[0071] Flexible and scalable: The modular design allows users to select different parameters and normalization strategies based on data characteristics. The visualization module facilitates result evaluation and process optimization. The system has good versatility and scalability.

[0072] refer to Figure 4 By comparing the peak height distribution changes at each stage of the data processing flow, the box width represents the data dispersion, and the median line shows the central trend. This allows observation of the effect of each processing step; the box should gradually narrow. (Reference) Figure 5 As shown, the effect of data standardization on improving peak height distribution can be observed. Colors represent density distribution, and the MAD value measures the data dispersion, providing a quantitative assessment of the standardization effect. An MAD improvement of >50% is considered excellent. (Reference) Figure 6 As shown, the GC content bias and its correction effect are displayed. The color depth represents the data point density, and the red trend line shows the degree of GC bias. The changes in the correlation coefficient before and after correction can be compared.

[0073] Based on the same concept, embodiments of the present invention also provide a population epigenetic analysis system for implementing the method, comprising:

[0074] (1) Consensus peak identification module, used to execute step S1 and generate consensus peak set and sample occurrence list.

[0075] This invention presents a method for constructing a set of common peaks on 533 ChIP-seq samples, along with its characterization data, to support the technical features regarding the common peak identification algorithm and parameter ranges in the claims. The invention tests the parameter `min_sample_fraction` at three values: 0.1, 0.2, and 0.3, as set in the technical solution, to evaluate the sensitivity of the number of common peaks to parameter changes. When `min_sample_fraction` = 0.1, the system identifies approximately n1 common peaks; when `min_sample_fraction` = 0.2, it identifies approximately n2 common peaks; and when `min_sample_fraction` = 0.3, it identifies approximately n3 common peaks, where n1 ≥ n2 ≥ n3, indicating that the algorithm can flexibly adjust the conservatism of the common peak set according to the sample ratio. Under all parameter conditions, candidate peak regions are matched based on `max_position_shift` (using three values: 50bp, 100bp, and 150bp), and after coordinate offset correction, peak objects with consensus start and end coordinates are generated. The structured matrix of the common peak results includes fields such as chrom, start, end, peak_id, and samples_list, which can be used as structured input for subsequent peak coverage calculations. This embodiment demonstrates that the method of the present invention can stably identify consistent enriched regions across sample signals at a multi-sample scale, without relying on the construction of artificial peak sets, and has repeatability and scalability.

[0076] (2) Peak intensity calculation module, used to execute step S2 and generate the original peak intensity matrix.

[0077] How can the fragment-level coverage calculation method based on BAM files improve the accuracy of peak area estimation to support the technical features regarding the peak coverage calculation logic in the claims? Under PE library conditions, this invention calculates fragment coverage depth using template length and infers fragment intervals by comparing the reference_start and template_length of the readpair, ensuring that the coverage reflects the actual DNA fragment structure. Under SE library conditions, the coverage interval is calculated based on the read start position and read length. To verify the algorithm's performance, this embodiment sets three test conditions: using read-level estimation (control group), using fragment-level estimation (in this invention's method, the template length upper limit is 600bp, 800bp, and 1000bp), and a mixed-mode test. In 1000 known true enrichment regions, the correlation coefficient between the peak area estimation of this invention and the true value is higher than that of the control group. =0.61 promoted to fragment-level method =0.80, an improvement of over 30%. Peak area distribution plots and coverage curves also show that this method provides smoother, abrupt coverage curves at both the peak center and peak boundaries, improving the continuity of the peak intensity matrix. This embodiment demonstrates that fragment-level coverage modeling is a crucial technique for ensuring the accuracy of peak quantification.

[0078] (3) GC correction module, used to perform step S3 and generate GC-corrected matrix.

[0079] To verify that the normalization module can be implemented under different numerical conditions, this embodiment selects three sets of peak height matrices (with sample sizes of 100, 300, and 533 points), and calculates the Pearson correlation coefficient between peak height and GC content before and after GC correction. The results show that under uncorrected conditions, the correlations of the three sets of data are r = 0.34, 0.33, and 0.32, respectively. After processing with the GC binning correction method of this invention, the correlations are reduced to r = 0.02, 0.01, and 0.00. The histogram further shows that the correlation is weakened in almost all GC intervals. Furthermore, to verify the applicability of the DESeq2 normalization method, this embodiment selects the sample with the highest and lowest median peak height as test objects. Before normalization, the ratio of their median peak heights is approximately 2.1:1, which converges to approximately 1.05:1 after normalization, reducing the difference by approximately 50%, consistent with the normalization mechanism in the claims. The normalized peak matrix was used in PCA to demonstrate the structural relationships between samples. All three datasets showed that samples were clustered by biological category, rather than by sequencing batch, indicating that numerical standardization eliminated major non-biological variations in the matrix. This embodiment supports the stable operation of the normalization module under different sample sizes through both mathematical modeling and performance testing.

[0080] (4) Normalization module, used to execute step S4, providing a variety of normalization methods and outputting a normalization matrix.

[0081] For three datasets with different sample sizes (100, 300, and 533 samples), Python-SVA automatically estimated the number of latent variables, n_sv, to be 4, 6, and 7, respectively. After latent variable regression and minimum value normalization, the peak height distribution showed good cross-sample consistency in the violin plot, and the MAplot showed no systematic shift. The correlation heatmap also no longer showed batch clustering structure. The corrected peak × sample matrix was used for two downstream tasks: one was differential enrichment analysis, using t-tests to screen regions with significant peak intensity differences; the other was machine learning tasks (such as random forests and XGBoost) to predict sample categories or biological traits. The input data for the above tasks were all standardized matrices generated by this invention, and their stability remained consistent under different training / test split ratios. This embodiment shows that the matrix obtained by the method of this invention not only meets the data consistency requirements but can also be directly used in large-scale epigenetic analysis workflows, providing a reproducible and verifiable implementation basis for the technical solution described in the claims.

[0082] (5) Batch effect correction module, used to perform step S5, remove implicit batch effects and output the final matrix.

[0083] Example embodiments have been disclosed herein, and while specific terminology has been used, it is for illustrative purposes only and should be construed as such, and is not intended to be limiting. In some instances, it will be apparent to those skilled in the art that features, characteristics, and / or elements described in conjunction with particular embodiments may be used alone, or in combination with features, characteristics, and / or elements described in conjunction with other embodiments, unless otherwise expressly indicated. Therefore, those skilled in the art will understand that various changes in form and detail may be made without departing from the scope of the invention as set forth in the appended claims.

Claims

1. A population epigenetic analysis method based on consensus peak identification, characterized in that, include: S1. Read the enriched peak region files of each sample; scan the peak regions after sorting by chromosome, judge and merge peaks belonging to the same enriched region according to the preset position offset threshold; filter out the regions that appear stably within the population based on the preset minimum sample ratio, and generate a consensus peak set and a sample list in which each consensus peak appears. S2. For each region in the consensus peak set, read the sequence alignment file of the corresponding sample; Based on the sequencing library type, the coverage depth of the region is calculated based on fragments or reads; based on the coverage depth, the intensity index of each sample at each consensus peak is calculated to form the original peak intensity matrix; S3. Extract the genome sequence of each consensus peak region and calculate its GC content; divide the GC content into multiple intervals; generate a correction factor based on the relationship between the local statistics and global statistics of the peak intensity in each GC content interval, and correct the original peak intensity matrix to obtain the GC-corrected matrix; S4. The GC-corrected matrix is ​​scaled using a selected normalization method to obtain the normalized peak intensity matrix. S5. Perform data transformation on the normalized peak intensity matrix; estimate the number of latent variables with hidden batch effects; remove the influence represented by the latent variables through a regression model; perform inverse transformation and post-processing to output the final corrected peak intensity matrix.

2. The method according to claim 1, characterized in that, In step S1, the position offset threshold ranges from 50 to 150 base pairs; the minimum sample ratio ranges from 0.1 to 0.

3.

3. The method according to claim 1, characterized in that, In step S2, the intensity index includes at least one of maximum coverage, average coverage, median coverage, and peak area; for paired-end sequencing data, the fragment coverage interval is calculated based on the template length of the paired reads. For single-end sequencing data, the coverage area is calculated based on the read length.

4. The method according to claim 1, characterized in that, In step S3, the GC content is uniformly divided into 50 to 200 intervals; the correction factor is determined based on the ratio of the local mean of the peak intensity in each GC interval to the global mean of all peak intensities.

5. The method according to claim 1, characterized in that, In step S4, the normalization method includes at least one of the following: DESeq2 size factor method, TMM method, CPM method, TPM method, and upper quartile method.

6. The method according to claim 1, characterized in that, Step S5 includes: performing an asinh transformation on the data matrix; automatically estimating the number of latent variables using latent variable analysis; constructing a linear model containing a design matrix of known biological variables and the latent variables for regression correction; and performing an inverse sinh transformation and minimum value shift on the residuals.

7. The method according to claim 1, characterized in that, In step S1, for each finally determined consensus peak, its start and end coordinates are the average of the start and end coordinates of all merged sample peak regions.

8. The method according to claim 1, characterized in that, The method further includes step S6: visualizing the quality assessment of the final corrected peak intensity matrix, including generating at least one of a peak intensity distribution map, a heatmap of inter-sample correlation, a principal component analysis map, and a MA map.

9. A population epigenetic analysis system for implementing the method according to any one of claims 1-8, characterized in that, include: The consensus peak identification module is used to perform step S1 and generate a consensus peak set and a sample occurrence list; The peak intensity calculation module is used to execute step S2 and generate the original peak intensity matrix; The GC correction module is used to perform step S3 and generate the GC-corrected matrix; The normalization module is used to perform step S4, providing multiple normalization methods and outputting a normalization matrix; The batch effect correction module is used to perform step S5, remove implicit batch effects, and output the final matrix.

10. The system according to claim 9, characterized in that, It also includes a visualization module for performing step S6 and generating a data quality assessment chart.

Images