Conservative non-coding element identification method and system based on multi-species genome comparison

Abstract

The invention discloses a conservative non-coding element identification method and system based on multi-species genome alignment, and the method comprises the following steps: establishing an index database based on reference genome sequences of multiple species, completing whole genome alignment, and further processing to obtain a high-credibility chain alignment result; the multi-species chain type comparison results are integrated into multi-sequence comparison data in a unified format; on this basis, a neutral evolution model is constructed based on quadruple degenerate sites, and candidate conservative regions are predicted through conservative scoring; in combination with genome annotation information, a length threshold is set, and a coding region and a UTR region are rejected, so that a high-confidence non-coding conservative element is obtained; and finally, displaying a cross-species conservative distribution diagram of the CNE by utilizing a visual tool. According to the method, the CNE with a potential regulation function can be accurately, efficiently and automatically identified, and technical support is provided for regulation system analysis, functional gene mining and molecular breeding of various organisms.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of bioinformatics, and in particular to a method and system for identifying conserved non-coding elements based on multi-species genome alignment. BACKGROUND

[0002] In eukaryotic genomes, a large number of non-coding regions in addition to coding regions also have important biological functions, especially the highly conserved non-coding elements (CNEs) that maintain high sequence conservation among multiple species, which are widely believed to be involved in key life activities such as gene expression regulation, chromatin structure maintenance, development process control, and usually distributed in promoters, introns, UTRs and intergenic regions, with potential cis-regulatory functions. Studies have shown that these conserved regions often enrich regulatory elements such as transcription factor binding sites (TFBS), enhancers or gene silencing regulatory elements, and play a key role in animal and plant evolution, development and metabolic pathways.

[0003] At present, the method of genome alignment analysis mainly depends on traditional alignment algorithms such as BLAST, BWA, etc. These methods can achieve the alignment and analysis of genome sequences to some extent. However, when facing complex multi-species genomes, the existing methods often have some problems. First, the amount of genomic data is huge, and the existing alignment methods are low in efficiency when processing large-scale data. Second, in the preliminary alignment data, due to the differences in genomes between species, there are often errors in connection and mismatch, which affects the accuracy of the analysis results. Finally, how to extract and analyze non-coding conserved regions in multi-species genomes still lacks a unified and efficient solution.

[0004] Therefore, it is urgent to propose a non-coding conserved element identification method with strong reproducibility and wide applicability, which can eliminate the false judgment caused by the difference in background mutation rate, compatible with fragmented and complete genomes, directly output annotated conserved atlas, and facilitate researchers to locate functional sites. SUMMARY

[0005] Therefore, the present application provides a method and system for identifying conserved non-coding elements based on multi-species genome alignment, to solve the technical problems of high false positive rate, poor adaptability and difficult to intuitively interpret the results when identifying non-coding conserved elements in multi-species fragmented genomes in the prior art.

[0006] To achieve the above technical purpose, the present application adopts the following technical solutions: In a first aspect, the present application provides a method for identifying conserved non-coding elements based on multi-species genome alignment, comprising: establishing an index database based on reference genome sequences of multiple species; aligning the genome to be analyzed with reference genome sequences in an index database to generate initial alignment data; correcting redundancy and abnormal connections in the initial alignment data through chain structure optimization to obtain optimized alignment data; integrating multiple alignment fragments in the optimized alignment data to generate multi-species alignment results under the same data structure; constructing an evolutionary neutral model according to quadruply degenerate sites of the multi-species alignment results, and calculating a conservation score based on the model; screening a feature region with a conservation higher than a background variation rate according to a preset score threshold and the conservation score; filtering the feature region in length and function to obtain a non-coding conserved region; aligning the non-coding conserved region with multiple species using a visualization tool to generate a species conservation distribution map.

[0007] Further, the correcting redundancy and abnormal connections in the initial alignment data through chain structure optimization to obtain optimized alignment data comprises: aligning continuous sequence fragments of the genome to be analyzed with the index database for whole genome alignment to generate primary alignment results; detecting and marking cross-chromosome jumps, direction flips and redundant matches in the primary alignment results according to reference genome chromosome coordinate order; splitting or reordering the fragments using global collinearity information to output the optimized alignment results after chain structure correction.

[0008] Further, the integrating multiple alignment fragments in the optimized alignment data comprises: aligning alignment fragments of different species in the optimized alignment data according to reference genome positions; merging overlapping alignment fragments to ensure complete alignment information; when there is a missing region in any species, marking according to alignment information of other species.

[0009] Further, the multi-species alignment results adopt MAF, BED or GTF format, and the format at least includes alignment position, species information and alignment score.

[0010] Further, the constructing an evolutionary neutral model according to quadruply degenerate sites of the multi-species alignment results comprises: locating and extracting quadruply degenerate sites common to all species and located in protein coding sequence regions in the integrated alignment results; calculating an observed difference proportion for each pair of species according to the quadruply degenerate site set, and converting the difference proportion into a corrected nucleotide substitution number through a nucleotide substitution model; fitting the corrected number of nucleotide substitutions to a species phylogenetic distance, to obtain a neutral substitution rate per unit time; dividing the genome into several local windows, calculating local substitution rates in different genomic windows, and forming a local correction matrix; generating an evolutionary neutral model according to the neutral substitution rate per unit time and the local neutral substitution rate; wherein the local neutral substitution rate is used for local background correction, and the neutral substitution rate per unit time is used for global normalization.

[0011] Further, the nucleotide substitution model is a Jukes-Cantor single-parameter model, and its expression is: In the formula, d represents the corrected number of nucleotide substitutions, and p represents the observed difference ratio.

[0012] Further, the calculation of the conservation score based on the model comprises: setting a sliding window in the integrated alignment result; calculating an observed difference rate in each window and comparing it with a model expected variation rate; summarizing the conservation measures of all windows to form a conservation score in the whole genome range.

[0013] Further, the length and function filtering of the feature region comprises: length screening, used for removing fragments below a set lower limit; significance screening, used for removing regions whose conservation measures do not reach a preset significance threshold; coding region removal, used for removing regions overlapping with known protein coding regions; proximity filtering or repeated sampling test, used for improving the result stability according to research requirements.

[0014] In a second aspect, the present application provides a conserved non-coding element identification system based on multi-species genome alignment, comprising: a database establishment module, used for establishing an index database based on reference genome sequences of multiple species; an alignment module, used for aligning a to-be-analyzed genome with reference genome sequences in the index database to generate initial alignment data; and correcting erroneous connections in the initial alignment data through chain structure to obtain optimized alignment data; an integration module, used for aligning a to-be-analyzed genome with reference genome sequences in the index database to generate initial alignment data; correcting redundancy and abnormal connections in the initial alignment data through chain structure optimization to obtain optimized alignment data; and integrating multiple alignment fragments in the optimized alignment data to generate multi-species alignment results under the same data structure; a conservativeness scoring module configured to construct an evolutionary neutral model according to the quadruply degenerate sites of the multi-species alignment result, and calculate a conservativeness score based on the model; a feature region determining module configured to filter a feature region with a conservativeness higher than a background mutation rate according to a preset score threshold and the conservativeness score; a filtering module configured to perform length and function filtering on the feature region to obtain a non-coding conserved region; an output module configured to perform multi-species alignment on the non-coding conserved region by using a visualization tool to generate a species conservativeness distribution map.

[0015] In a third aspect, the present application further provides a computer readable storage medium having a computer program stored therein, wherein the computer program, when executed by a processor, implements the multi-species genome alignment based conserved non-coding element identification method described above.

[0016] Compared with the prior art, the multi-species genome alignment based conserved non-coding element identification method provided by the present application has the following advantages: 1. High integration and practicality: from index establishment, multi-species alignment, and chain structure correction to integrated alignment, model construction, filtering, and visualization, a conservativeness distribution map is output in one key, which saves the manual splicing link and is suitable for a wide range of species combinations and genome quality.

[0017] 2. More accurate recognition: a global-local double-layer neutral model is established based on quadruply degenerate sites, the mutation background is corrected segment by segment, and known coding regions and short fragments are automatically filtered in combination with annotation information, thereby improving the non-coding specificity and functional reliability of CNE.

[0018] 3. Visualization of results: after known coding regions and short fragments are automatically removed, only high-confidence non-coding conserved elements are retained; the generated visualization atlas synchronously displays adjacent genes and regulation annotations, and researchers can directly locate functional sites.

[0019] In summary, the present application can accurately, efficiently, and automatically identify CNEs with regulatory potential, thereby providing technical support for regulatory system analysis, functional gene mining, and molecular breeding of various organisms. BRIEF DESCRIPTION OF DRAWINGS

[0020] Figure 1 A flowchart of the multi-species genome alignment based conserved non-coding element identification method provided by the present application is shown in the figure. Figure 2 A structural diagram of the multi-species genome alignment based conserved non-coding element identification system provided by the present application is shown in the figure. DETAILED DESCRIPTION

[0021] Preferred embodiments of the present application will be described in detail below with reference to the drawings, wherein the drawings constitute a part of this application, and illustrate embodiments of the present application together with the principles of the present application, but are not intended to limit the scope of the present application.

[0022] Referring to Figure 1 The embodiment provides a method for identifying a conserved non-coding element based on multi-species genome alignment, comprising the following steps: Step S101: establishing an index database based on reference genome sequences of multiple species; Step S102: aligning a to-be-analyzed genome with reference genome sequences in the index database to generate initial alignment data; and correcting redundancy and abnormal connections in the initial alignment data through chain structure optimization to obtain optimized alignment data; Step S103: integrating multiple alignment fragments in the optimized alignment data to generate multi-species alignment results under the same data structure; Step S104: constructing an evolutionary neutral model according to quadruply degenerate sites of the multi-species alignment results, and calculating a conservation score based on the model; Step S105: screening a feature region with a conservation higher than a background mutation rate according to a preset score threshold and the conservation score; Step S106: performing length and function filtering on the feature region to obtain a non-coding conserved region; Step S107: performing multi-species alignment on the non-coding conserved region by using a visualization tool to generate a species conservation distribution map.

[0023] The method provided in the embodiment establishes an index based on multi-species reference genomes, obtains high-confidence alignment through chain structure correction, and integrates to a unified data structure; constructs a double-layer evolutionary neutral model relying on quadruply degenerate sites, corrects the background mutation rate window by window, and outputs a conservation score; locks a high-confidence non-coding conserved region through threshold screening, coding region elimination and robustness filtering; and finally generates a cross-species visualization map. The double-layer model can significantly reduce false positives and improve the accuracy of CNE recognition. The method is compatible with fragmented or complete genomes, and is suitable for multiple species and multiple sequencing qualities. The method outputs a conserved map with functional annotations, and accelerates downstream experiments and achievement publication.

[0024] As a preferred embodiment, in step S102, the redundancy and abnormal connections in the initial alignment data are corrected through chain structure optimization to obtain the optimized alignment data, comprising: Performing whole-genome alignment on continuous sequence fragments of the to-be-analyzed genome and the index database to generate a primary alignment result; According to the reference genome chromosome coordinate order, detecting and marking the cross-chromosome jumps, direction flips and redundant matches in the primary alignment result; The global collinearity information is used to split or reorder the contigs, and the output is the optimized alignment result with corrected chain structure.

[0025] Specifically, the preliminary alignment data after alignment is usually a pair of sequence alignment results, which shows the similarities and differences between the analyzed genome and the reference genome. These alignment results are usually stored in a certain format (such as SAM, BAM or BLAT output format), which contains information such as alignment position, matching degree, insertion / deletion, etc. In addition to the basic correspondence, the alignment result may also include alignment score, alignment quality, alignment coverage, etc. These information helps to evaluate the accuracy and reliability of the alignment. Through chain structure correction, the connection between the sequence fragments (contigs) in the analysis alignment result can be corrected.

[0026] As a preferred embodiment, in step S103, the integration of multiple alignment fragments in the optimized alignment data includes: Aligning the alignment fragments of different species in the optimized alignment data according to the reference genome position; Merging overlapping alignment fragments to ensure complete alignment information; When there is a missing region in any species, marking according to the alignment information of other species.

[0027] As a preferred embodiment, the multi-species alignment result adopts MAF, BED or GTF format, and the format at least includes alignment position, species information and alignment score.

[0028] Since in the coding sequence (CDS), the third base of some codons (i.e. quadruple degenerate site) can be mutated without affecting the encoding of amino acids. These sites are usually not under selective pressure, so they can be used as representatives of background mutation rate. As a preferred embodiment, in step S104, the construction of the evolutionary neutral model according to the quadruple degenerate sites of the multi-species alignment result includes: In the integrated alignment result, locate and extract all quadruple degenerate sites that are common to all species and located in the protein coding sequence region; According to the set of quadruple degenerate sites, calculate the observed difference proportion for each pair of species, and convert the difference proportion to the corrected nucleotide substitution number through the nucleotide substitution model; Fit the corrected nucleotide substitution number with the species phylogenetic distance to obtain the neutral substitution rate per unit time; Divide the genome into several local windows, and calculate the local substitution rate in different genome windows to form a local correction matrix; The evolutionary neutral model is generated according to the unit time neutral substitution rate and the local neutral substitution rate; wherein the local neutral substitution rate is used for local background correction, and the unit time neutral substitution rate is used for global normalization.

[0029] In some embodiments, the method further comprises performing Monte Carlo resampling or Bayesian updating on the evolutionary neutral model to determine parameter confidence intervals and output a final evolutionary neutral model.

[0030] As a preferred embodiment, the nucleotide substitution model is a Jukes-Cantor single parameter model, whose expression is: In the formula, d represents the corrected number of nucleotide substitutions, and p represents the observed difference ratio.

[0031] As a preferred embodiment, in step S104, the conservation score is calculated based on the model, including: Setting a sliding window in the integrated alignment result; Calculating the observed difference rate in each window and comparing it with the model expected variation rate; Summarizing the conservation measures of all windows to form a whole genome range conservation score.

[0032] As a specific embodiment, in step S105, by setting a threshold for the conservation score, regions with conservation higher than the background variation rate are identified. These regions are considered as candidate functional regions that are "evolutionarily invariant" and may have important biological functions.

[0033] As a preferred embodiment, in step S106, the filtering of the feature regions at least includes length screening, significance screening, coding region removal, gene proximity filtering, and robust resampling verification; The length screening is used to retain continuous regions with a length not less than a preset lower limit; The significance screening is used to remove regions whose conservation measures do not reach a preset significance threshold; The coding region removal is used to remove regions with an overlap ratio with known protein coding regions exceeding a preset upper limit; The gene proximity filtering is used to remove regions with a distance to the nearest transcription start site less than a preset threshold; The robust resampling verification is used to perform resampling inspection on the merged regions, and only regions that meet the significance requirements in a preset proportion of times are retained.

[0034] In one specific embodiment, in step S107, the conserved non-coding sequences of all species at the corresponding positions are extracted, and multi-species visual alignment is performed to provide a clear conservation distribution map. The extracted CNE sequences are subjected to multi-species alignment using a visualization tool (such as UCSC Genome Browser, IGV or MultiPipMaker) to generate a clear conservation distribution map. These tools can show sequence similarity, conservation score and structural features between different species.

[0035] The above scheme is described below through a specific embodiment.

[0036] Taking human as a representative species, monkeys, mice, pigs, lungfish and the like are selected as comparative species, and multi-species whole genome multiple alignment and CNE identification are carried out.

[0037] First, high-quality multi-species WGA alignment files are obtained, and accurate neutral models and conservation models are constructed.

[0038] A total of 40,000-50,000 CNEs with a length of ≥30 bp are obtained, and the average length is 90-120 bp.

[0039] After removing the exon region, 35,000-45,000 non-coding CNEs are retained, mainly distributed in the regulatory region (such as the promoter, intron, and the vicinity of the remote enhancer), and the retention rate is about 90%. If a more stringent threshold (such as ≥100 bp, posteriori ≥0.9) is required, the total amount can be reduced to 15,000-20,000.

[0040] Part of the region is visualized using VISTA, and the significantly conserved region corresponds to the known regulatory element (such as the HD-ZIP and bZIP transcription factor binding site).

[0041] As shown in Figure 2 The embodiment of the present application also provides a conserved non-coding element identification system based on multi-species genome alignment, which comprises: A database establishment module 201 is configured to establish an index database based on reference genome sequences of multiple species; An alignment module 202 is configured to align a to-be-analyzed genome with reference genome sequences in the index database to generate initial alignment data, and correct erroneous connections in the initial alignment data through a chain structure to obtain optimized alignment data; An integration module 203 is configured to merge alignment data of multiple gene fragments in the optimized alignment data, integrate alignment fragments of multiple species, and form an integrated alignment result under the same data structure; The conservativeness scoring module 204 is configured to construct an evolutionary neutral model according to the quadruple degenerate sites of the integrated alignment result, analyze the integrated alignment result based on the evolutionary neutral model, and output a conservativeness score; The feature region determining module 205 is configured to determine a feature region with a conservativeness higher than a background variation rate according to a preset scoring threshold and the conservativeness score. The filtering module 206 is configured to filter the feature region to obtain a filtered non-coding conservative region. The output module 207 is configured to perform multi-species alignment on the non-coding conservative region by using a visualization tool to generate a species conservativeness distribution map.

[0042] It should be noted that the modules of the method can be independently deployed or integrated into an automatic pipeline system, and are suitable for scenarios such as scientific research institutions, high-throughput functional genomics platforms, and biological breeding companies.

[0043] The above describes only the preferred specific embodiments of the present application, but the protection scope of the present application is not limited thereto, and any person skilled in the art can easily think of changes or replacements within the technical range disclosed by the present application, which should be covered within the protection scope of the present application.

Claims

1. A method for identifying conserved non-coding elements based on multi-species genome alignment, characterized by, The method comprises the following steps: establishing an index database based on reference genome sequences of multiple species; aligning a genome to be analyzed with reference genome sequences in the index database to generate initial alignment data; correcting redundancies and abnormal connections in the initial alignment data through chain structure optimization to obtain optimized alignment data; integrating multiple alignment fragments in the optimized alignment data to generate multi-species alignment results under the same data structure; constructing an evolutionary neutral model according to four-fold degenerate sites of the multi-species alignment results, and calculating a conservation score based on the model; screening a feature region with a conservation higher than a background variation rate according to a preset score threshold and the conservation score; filtering the feature region in length and function to obtain a non-coding conservative region; performing multi-species alignment on the non-coding conservative region by using a visualization tool to generate a species conservation distribution map.

2. The method for identifying a conserved non-coding element based on multi-species genome alignment according to claim 1, wherein, The method of correcting redundancies and abnormal connections in the initial alignment data through chain structure optimization to obtain optimized alignment data comprises the following steps: performing whole genome alignment of continuous sequence fragments of the genome to be analyzed with the index database to generate primary alignment results; detecting and marking cross-chromosome jumps, direction flips and redundant matches in the primary alignment results according to the order of reference genome chromosome coordinates; splitting or reordering the sequence fragments by using global collinearity information to output the optimized alignment results after chain structure correction.

3. The method for identifying conserved non-coding elements based on multi-species genome alignment according to claim 1, wherein, The method of integrating multiple alignment fragments in the optimized alignment data comprises the following steps: aligning alignment fragments of different species in the optimized alignment data according to reference genome positions; merging overlapping alignment fragments to ensure complete alignment information; when there is a missing region in any species, marking according to alignment information of other species.

4. The method for identifying a conserved non-coding element based on multi-species genome alignment according to claim 1, wherein, The multi-species alignment results adopt MAF, BED or GTF format, and the format at least includes alignment position, species information and alignment score.

5. The method for identifying conserved non-coding elements based on multi-species genome alignment according to claim 1, wherein, The method of constructing an evolutionary neutral model according to four-fold degenerate sites of the multi-species alignment results comprises the following steps: locating and extracting four-fold degenerate sites common to all species and located in protein coding sequence regions in the integrated alignment results; calculating an observed difference proportion for each pair of species according to the four-fold degenerate site set, and converting the difference proportion into a corrected nucleotide substitution number by using a nucleotide substitution model; fitting the corrected nucleotide substitution number with a species phylogenetic distance to obtain a neutral substitution rate per unit time; dividing the genome into a plurality of local windows, and calculating a local substitution rate in different genome windows to form a local correction matrix; generating an evolutionary neutral model according to the neutral substitution rate per unit time and the local neutral substitution rate; wherein the local neutral substitution rate is used for local background correction, and the neutral substitution rate per unit time is used for global normalization.

6. The method for identifying conserved non-coding elements based on multi-species genome alignment according to claim 5, wherein, The nucleotide substitution model is a Jukes-Cantor single-parameter model, and its expression is as follows: In the formula, d represents the corrected nucleotide substitution number, and p represents the observed difference proportion.

7. The method of claim 1, wherein, The method of calculating a conservation score based on the model comprises the following steps: setting a sliding window in the integrated alignment results; calculating an observed difference rate in each window, and comparing it with a model expected variation rate; The conservation scores of all windows are aggregated to form a genome-wide conservation score.

8. The method of claim 1, wherein, The length and function filtering of the feature regions comprises: Length screening for removing fragments below a set lower limit; Significance screening for removing regions whose conservation scores do not reach a preset significance threshold; Coding region removal for removing regions overlapping with known protein coding regions; Adjacent filtering or repeated sampling test for improving the stability of results according to research needs.

9. A system for identifying conserved non-coding elements based on multi-species genome alignment, characterized by, The method comprises the following steps: A database establishment module is configured to establish an index database based on reference genome sequences of multiple species; An alignment module is configured to align a to-be-analyzed genome with the reference genome sequences in the index database to generate initial alignment data, and correct erroneous connections in the initial alignment data through chain structure to obtain optimized alignment data; An integration module is configured to align a to-be-analyzed genome with the reference genome sequences in the index database to generate initial alignment data, correct redundancy and abnormal connections in the initial alignment data through chain structure optimization to obtain optimized alignment data, and integrate multiple alignment fragments in the optimized alignment data to generate multi-species alignment results in the same data structure; A conservation scoring module is configured to construct an evolutionary neutral model according to quadruple degenerate sites of multi-species alignment results, and calculate a conservation score based on the model; A feature region determination module is configured to screen feature regions with conservation higher than a background variation rate according to a preset scoring threshold and the conservation score; A filtering module is configured to perform length and function filtering on the feature regions to obtain non-coding conserved regions; An output module is configured to perform multi-species alignment on the non-coding conserved regions by using a visualization tool to generate a species conservation distribution map.

10. A computer-readable storage medium, characterized in that, The computer program stored in the computer-readable storage medium is executed by the processor to implement the method for identifying conserved non-coding elements based on multi-species genome alignment according to any one of claims 1-8.

Images