DNase high-throughput sequencing data preference correction method based on Bayesian method and application

Abstract

The invention relates to the technical field of biological information correction, in particular to a DNase high-throughput sequencing data preference correction method based on a Bayesian method and application, and the method comprises the following steps: firstly, constructing a DNase-seq library, processing data, and obtaining sequencing coverage data; secondly, a Bayesian model is constructed, observation signals are modeled to be in Poisson distribution, real signals are in Gamma distribution, and preference signals are in logarithmic normal distribution; then, training the Bayesian model by using a Markov chain Monte Carlo (MCMC) method to obtain posteriori distribution of the real signal and the preference signal, and calculating a posteriori mean value; and finally, the observation signal is corrected by using the posterior mean value, and the real signal is recovered. According to the method, through combining prior distribution and observation data, accurate estimation of real signals is provided, and preference noise in the data is effectively reduced. Taking corn stem tissue as an example, the coefficient of variation of the corrected coverage data is significantly reduced, the volatility is reduced by about 86.40%, and the false positive risk in the data is effectively reduced.

Description

Technical Field

[0001] The present invention relates to the technical field of biological information correction, and in particular to a DNase high-throughput sequencing data bias correction method based on a Bayesian method and an application thereof. Background Art

[0002] DNase-seq (DNase I hypersensitive sites sequencing) is a high-throughput sequencing technology widely used to study open regions of the genome. This technology uses DNase I enzyme to cut open chromatin regions and perform high-throughput sequencing on these regions to identify regulatory elements in the genome. However, DNase-seq data are often affected by sequencing biases, such as sequencing depth, enzyme cutting efficiency, etc. These biases can significantly affect the quality and reliability of DNase-seq data. For example, the existence of biases may cause some real regulatory elements to be misidentified or omitted, misleading researchers' conclusions. Therefore, how to effectively correct these biases and restore the real biological signals in the data has become an important challenge in DNase-seq data analysis.

[0003] Traditional bias correction methods are usually based on empirical models and rely on pre-set assumptions and parameters. These methods often fail to make full use of all the information in the data and lack reasonable handling of uncertainty. In addition, empirical models often find it difficult to provide accurate and reliable correction results when faced with complex bias patterns. With the development of computing power and statistical methods, the Bayesian method provides a new approach to bias correction with its flexibility and superiority in handling uncertainty. The Bayesian method establishes a statistical model, estimates the model parameters using observed data, and infers through the posterior distribution, thereby achieving bias correction and uncertainty quantification. The Bayesian method can not only make full use of the information in the data, but also provide uncertainty assessment of the correction results, thereby improving the reliability and accuracy of the analysis. Therefore, the DNase-seq data bias correction method based on the Bayesian method has important theoretical significance and application value, and provides a new approach to data bias correction. Summary of the invention

[0004] The purpose of the present invention is to solve the shortcomings in the prior art, and to propose a DNase high-throughput sequencing data bias correction method and application based on the Bayesian method, by building a Bayesian model, decomposing the observed data into real signals and bias signals, and correcting the bias by the Bayesian inference method, thereby restoring the real signal. The method can effectively correct sequencing bias and improve the accuracy and reliability of DNase-seq data.

[0005] In order to achieve the above object, the present invention adopts the following technical solutions:

[0006] A method for correcting DNase high-throughput sequencing data bias based on a Bayesian approach, the correction method comprising the following steps:

[0007] Step 1: DNase-seq library construction and data processing:

[0008] Take corn stem tissue as an example to carry out DNase-seq experiment, obtain the sequencing coverage data of DNase-seq experiment, and record the sequencing depth of each site in the genome as the initial observation data y i ;

[0009] Step 2: Build a Bayesian model:

[0010] Assume that the observed signal y i is the real signal λ i and the preference signal θ i As a result of the combined effect, the observed signal is modeled as a Poisson distribution, whose mean is obtained by multiplying the true signal and the preference signal:

[0011] y i ~Poisson(λ i ·θ i )

[0012] Among them, y i is the observed signal, λ i is the real signal, θ i Signal preference;

[0013] True signal λ i Gamma distribution:

[0014] λ i ~Gamma(α,β)

[0015] Among them, λ i is the real signal, α and β are the shape parameter and scale parameter of the Gamma distribution;

[0016] Preference signal θ i Follows a lognormal distribution:

[0017] θ i ~LogNormal(μ θ ,σ θ )

[0018] Among them, θ i is a preference signal, μ θ and σ θ are the mean and standard deviation of the lognormal distribution;

[0019] Step 3: Model training:

[0020] The Bayesian model is trained using the Markov Chain Monte Carlo (MCMC) method to obtain the parameter λ i and θ i The posterior distribution of ;

[0021] P(λ i ,θ i ∣y i )∝P(yi∣λ i ,θ i )·P(λ i )·P(θ i )

[0022] Among them, P(λ i ,θ i ∣y i ) is given the observed data y i In the case of i and the preference signal θ i The joint posterior distribution of i ∣λ i ,θ i ) is the likelihood function, which means that given the real signal and the preference signal, the observed signal y i The probability of i ) and P(θ i ) are the prior distributions of the true signal and the preference signal respectively;

[0023] Based on the obtained posterior distribution, the true signal λ at each position is calculated i and preference θ i The posterior mean of ;

[0024] For the real signal λ i The posterior mean estimate of is:

[0025]

[0026] in, is i The posterior mean of i (j) represents the λ obtained by the hth sampling i value, N is the number of sampling;

[0027] For the preference θ i The posterior mean estimate of is:

[0028]

[0029] in, is θ i The posterior mean of i (j) represents the θ obtained by the jth sampling i value, N is the number of sampling;

[0030] Step 4: Preference Correction:

[0031] Using the estimated preference parameters Correct the observed signal to obtain the corrected real signal:

[0032]

[0033] Among them, corrected_coverage i represents the value of the observed signal after correction of the preference parameter, y i To observe the signal, is the posterior mean of the preference signal;

[0034] Save the corrected data as a new file for subsequent analysis.

[0035] Preferably, in step one, DNase-seq library construction and data processing include: DNase-seq library construction; using trim_galore software for data quality control; using Bowtie2 software to align DNase-seq data to the maize reference genome; and using bedtools software to generate a genome coverage file.

[0036] Preferably, in step 2, the prior parameters α and β in the prior distribution of the real signal are set to 2, and the prior parameter μ in the prior distribution of the preference signal is set to θ and σ θ Set to 0 and 1 respectively.

[0037] The present invention also provides an application of the correction method in plants, wherein the plant is corn, and after the DNase-seq data of corn stem tissue is corrected, the data bias can be reduced.

[0038] By adopting the above technical solution: the present invention provides an accurate estimation of the real signal by combining the prior distribution with the observed data, effectively reducing the bias noise in the data. Taking corn stem tissue as an example, the coefficient of variation of the coverage data after correction is significantly reduced, and the volatility is reduced by about 86.40%. This correction method effectively reduces the risk of false positives in the data and provides reliable data guarantee for genome function research and application.

[0039] Compared with the prior art, the present invention has the following beneficial effects:

[0040] 1. The present invention adopts the Bayesian method, and by introducing the prior distribution, it effectively separates the real signal from the preference signal and realizes accurate signal correction. The Bayesian method comprehensively utilizes the prior knowledge and observation data, and provides a reasonable estimate of the real signal, thereby providing a solid data guarantee for subsequent data analysis and application.

[0041] 2. After correcting the DNase-seq data of corn stem tissue, the Bayesian model used in the present invention significantly reduced the coefficient of variation of the coverage data from 2.2037 before correction to 0.2996 after correction, reducing the volatility by about 86.40%. This shows that the signal distribution after correction is smoother and more stable, effectively reducing the risk of false positives introduced by preferential signals, and providing a more reliable basis for subsequent genome function research. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] Figure 1 It is a schematic diagram of the variation of the original coverage and the corrected coverage with the genomic position of the present invention;

[0043] Figure 2 It is a schematic diagram of the scattered distribution of the original coverage and the corrected coverage of the present invention;

[0044] Figure 3 Schematic diagram of the coefficient of variation of the original coverage and the corrected coverage of the present invention. DETAILED DESCRIPTION

[0045] The technical solutions in the embodiments of the present invention will be described clearly and completely below in conjunction with the accompanying drawings, so that those skilled in the art can better understand the advantages and features of the present invention, thereby making a clearer definition of the protection scope of the present invention. The embodiments described in the present invention are only a part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments in the present invention, all other embodiments obtained by ordinary technicians in the field without making creative work are within the scope of protection of the present invention.

[0046] A method for correcting DNase high-throughput sequencing data bias based on a Bayesian approach, the correction method comprising the following steps:

[0047] (1) DNase-seq library construction and data processing

[0048] DNase-seq experiments were performed using maize (Zea mays) stem tissue as a sample. First, fresh maize stem tissue was collected, immediately frozen in liquid nitrogen to maintain its original state, and then ground into a uniform powder. After grinding, about 1 gram of powder sample was weighed and suspended in an equal volume of pre-cooled nuclear isolation buffer (NIB; 10mM Tris-HCl, 80mM KCl, 10mM EDTA, 1mM spermidine, 1mM spermine, 0.15% mercaptoethanol, 0.5M sucrose, pH 9.5) to isolate cell nuclei. Next, the extracted nuclei were resuspended in nuclear digestion buffer (NDB; 10mM Tris-HCl, 10mM NaCl, 3mM MgCl2, pH 7.5), and the nuclear DNA was digested using gradient concentrations of DNase I enzyme at 37°C for 10 minutes.

[0049] The degree of DNase I digestion was assessed by pulsed field gel electrophoresis (PFGE). High molecular weight DNA after DNase I digestion was extracted and end-repaired with T4 DNA polymerase. Subsequently, biotin-labeled adapter I (5′Bio ACAGGT TCA GAG TTC TAC AGT CCG AC and 5′P-GTC GGA CTG TAG AAC TCT GAA C) was attached to the DNA molecules. After MmeI digestion, magnetic beads were used to enrich the DNase I digested DNA ends. Adapter II (5′P-TCG TATGCC GTC TTC TGC TTG and 5′-CAA GCA GAA GAC GGC ATA CGA NN) was attached to the MmeI treated DNA ends. PCR amplification was performed using adapter-specific primers (5′-CAA GCA GAA GAC GGC ATA CGA and 5′-AAT GAT ACG GCGACC ACC GAC AGG TTC AGA GTT CTA CAG TCC GA), and 90 bp DNA fragments were separated by polyacrylamide gel electrophoresis (PAGE). Finally, sequencing was performed on the Illumina HiSeq 2500 platform.

[0050] A total of 86,837,684 reads were measured, ensuring the coverage and depth of the sequencing data, which was suitable for subsequent analysis. The obtained sequencing data were first quality controlled by trim_galore software to remove low-quality reads and adapter sequences to ensure the overall reliability of the data. Then, Bowtie2 software was used to align the processed reads to the reference genome of maize (version: B73 RefGen_v4). The alignment parameters were set to: sensitive mode (--sensitive) to improve the alignment accuracy, and the maximum acceptable alignment mismatch number was 0 (-N 0) to reduce the mismatch caused by mismatch. Finally, the bedtools software was used to generate the genome coverage file, and the generated coverage.bedgraph file recorded the sequencing depth information of each site in detail, providing a reliable data basis for the subsequent Bayesian bias correction method.

[0051] (2) Building a Bayesian model

[0052] A Bayesian model is constructed to explain the coverage data observed in DNase-seq experiments. The core idea of ​​this model is to regard the observed signal as the result of the interaction between the real signal and the bias signal. The specific process is as follows:

[0053] First, assume that the coverage data y for each site in the DNase-seq experiment is i is driven by two main factors: the true signal λ i and the preference signal θ i Observation signal y i It follows a Poisson distribution with a mean of λ i and θ i This means that the observed coverage y i Not only does it reflect the true activity of the genomic region (λ i ), and is also affected by various preference factors during the experiment (such as sequencing preference, nuclease cleavage preference, etc.) (θ i ). Therefore, y i It can be expressed as:

[0054] y i ~Poisson(λ i ·θ i )

[0055] Among them, y i is the observed signal, λ i is the real signal, θ i Signal preference;

[0056] Next, for the real signal λ i Modeling. iGamma distribution is subject to Gamma distribution. Gamma distribution can flexibly model positive data of normal distribution, and is especially suitable for describing biological signals with positive values ​​and small fluctuations. The shape parameter α and scale parameter β of Gamma distribution determine the concentration and volatility of the signal. The details are as follows:

[0057] λ i ~Gamma(α,β)

[0058] Among them, λ i is the real signal. The shape parameter α is set to 2, indicating that the expected real signal is relatively stable and has small fluctuations. The scale parameter β is also set to 2, so that the mean of the Gamma distribution is more concentrated, which is in line with biological expectations, that is, the real signal should not have too large deviations.

[0059] In order to capture the preference signal θ in the experiment i The heterogeneity of θ i Obey the lognormal distribution. The lognormal distribution can effectively describe the wide distribution of positive data, and is particularly suitable for modeling preference signals that come from multiple factors and have wide fluctuations. The details are as follows:

[0060] θ i ~LogNormal(μ θ ,σ θ )

[0061] Among them, θ i is a preference signal. The logarithmic mean parameter μ θ is set to 0, which means there is no obvious preference. The standard deviation parameter σ θ Set to 1, this large standard deviation allows the bias signal to have a wider range of fluctuations between sites to cover a variety of possible sources of experimental bias.

[0062] (3) Model training

[0063] The constructed Bayesian model is trained by the Markov Chain Monte Carlo (MCMC) method to obtain the true signal λ i and the preference signal θ i To ensure the convergence of the model and the reliability of the sampling results, the present invention adopts the following training process:

[0064] First, No-U-Turn Sampler (NUTS) is selected as the sampling algorithm. NUTS is a variant of Hamiltonian Monte Carlo (HMC), which can automatically adjust the step size to improve the efficiency and accuracy of sampling. NUTS is particularly effective in processing high-dimensional parameter spaces, and can ensure that the present invention obtains high-quality posterior samples in complex Bayesian models.

[0065] Before starting formal sampling, perform 1000 initial iterations. These initial iterations are used to help the MCMC chain find high probability areas for model parameters and ensure that the MCMC chain has converged. Eliminate the influence of the initial values ​​of the parameters on the sampling results and ensure the reliability of subsequent sampling. During the formal sampling process, set the number of sampling to 10,000. A larger number of samples helps improve the estimation accuracy of the posterior distribution and ensure the accuracy of the real signal λ. i and the preference signal θ i The estimate of has a smaller variance.

[0066] Based on the obtained MCMC sampling results, calculate the true signal λ at each site i and the preference signal θ i These posterior means provide us with the best point estimate of the signal. The specific calculation method is as follows:

[0067] For the real signal λ i The posterior mean estimate of is:

[0068]

[0069] in, is i The posterior mean of i (j) represents the λ obtained by the jth sampling i value, N is the number of sampling;

[0070] For the preference θ i The posterior mean estimate of is:

[0071]

[0072] in, is θ i The posterior mean of i (j) represents the θ obtained by the jth sampling i value, N is the number of samplings.

[0073] (4) Correction of coverage data

[0074] Based on the calculated posterior mean, bias correction is performed on the sequencing data to improve the accuracy and reliability of the data. The specific steps are as follows:

[0075] Using the posterior mean of the preference signal obtained in the previous step, the coverage data y of each site is i Correction formula is:

[0076]

[0077] Among them, corrected_coverage i represents the value of the observed signal after correction of the preference parameter, y i To observe the signal, is the posterior mean of the bias signal. The purpose of this correction step is to eliminate the fluctuations caused by the bias signal to obtain coverage data that is closer to the actual genomic activity. The corrected coverage data will be saved in a new file for subsequent bioinformatics analysis.

[0078] (5) Compare the data before and after correction

[0079] Plot the curves of original coverage and corrected coverage versus genomic position ( Figure 1 ), showing that the corrected data is smoother and less volatile, reducing the fluctuations caused by the bias signal. Draw a scatter plot of the original coverage and the corrected coverage ( Figure 2 ), analyze the consistency of the data before and after correction. The closer the points are distributed along the diagonal line (y=x), the more consistent the coverage before and after correction. The coverage after correction is more uniform, which reflects the effectiveness of the correction method. In addition, the coefficient of variation ( Figure 3 ), which is used to evaluate the relative volatility of the data. The coefficient of variation before correction was high, indicating that the original data had large fluctuations. After correction, the coefficient of variation was significantly reduced, indicating the effectiveness of the correction method in reducing data volatility and improving data stability.

[0080] In summary, the present invention significantly improves the accuracy and reliability of DNase-seq sequencing data through the bias correction technology based on the Bayesian method, has high versatility and practical value, and is of great significance in promoting bioinformatics research, improving data analysis quality, and promoting scientific discovery.

[0081] The description and practice disclosed in the present invention are easy to think and understand for ordinary technicians in the technical field, and several improvements and modifications can be made without departing from the principles of the present invention. Therefore, modifications or improvements made without departing from the spirit of the present invention should also be regarded as the protection scope of the present invention.

Claims

1. A Bayesian-based bias correction method for DNase high-throughput sequencing data. It is characterized in that The calibration method comprises the following steps: Step 1: DNase-seq library construction and data processing: Obtain the sequencing coverage data of the DNase-seq experiment and record the sequencing depth of each site in the genome as the initial observation data y i ; Step 2: Build a Bayesian model: Assume that the observed signal y i is the real signal λ i and the preference signal θ i As a result of the combined effect, the observed signal is modeled as a Poisson distribution, whose mean is obtained by multiplying the true signal and the preference signal: y i ~Poisson(λ i ·i i ) Among them, y i is the observed signal, λ i is the real signal, θ i Signal preference; True signal λ i Gamma distribution: l i ~Gamma(α,β) Among them, λ i is the real signal, α and β are the shape parameter and scale parameter of the Gamma distribution; Preference signal θ i Follows a lognormal distribution: i i ~LogNormal(μ θ ,s θ ) Among them, θ i is a preference signal, μ θ and σ θ are the mean and standard deviation of the lognormal distribution; Step 3: Model training: The Bayesian model is trained using the Markov Chain Monte Carlo method to obtain the parameter λ i and θ i The posterior distribution of ; P(λ i ,i i ∣y i )∝P(y i ∣λ i ,i i )·P(λ i )·P(θ i ) Among them, P(λ i ,θ i ∣y i ) is given the observed data y i In the case of i and the preference signal θ i The joint posterior distribution of i ∣λ i ,θ i ) is the likelihood function, which means that given the real signal and the preference signal, the observed signal y i The probability of i ) and P(θ i ) are the prior distributions of the true signal and the preference signal respectively; Based on the obtained posterior distribution, the true signal λ at each position is calculated i and preference θ i The posterior mean of ; For the real signal λ i The posterior mean estimate of is: in, is i The posterior mean of i (j) represents the λ obtained by the jth sampling i value, N is the number of sampling; For the preference θ i The posterior mean estimate of is: in, is θ i The posterior mean of i (j) represents the θ obtained by the jth sampling i value, N is the number of sampling; Step 4: Preference Correction: Using the estimated preference parameters Correct the observed signal to obtain the corrected real signal: Among them, corrected_coveragei represents the value of the observed signal after correction of the preference parameter, y i To observe the signal, is the posterior mean of the preference signal; Save the corrected data as a new file for subsequent analysis.

2. A method for correcting DNase high-throughput sequencing data bias based on a Bayesian approach according to claim 1, It is characterized in that In step one, DNase-seq library construction and data processing include: DNase-seq library construction; data quality control using trim_galore software; DNase-seq data alignment to the maize reference genome using Bowtie2 software; and genome coverage file generation using bedtools software.

3. A method for correcting DNase high-throughput sequencing data bias based on a Bayesian approach according to claim 1, It is characterized in that In step 2, the prior parameters α and β in the prior distribution of the true signal are set to 2, and the prior parameter μ in the prior distribution of the preference signal is set to θ and σ θ Set to 0 and 1 respectively.

4. Use of the correction method according to any one of claims 1 to 3 in plants.

5. The use according to claim 4, It is characterized in that The plant is corn, and after the DNase-seq data of corn stem tissue is corrected, the data bias can be reduced.

Images