Genomic word framework analysis for genomic methylation data

Abstract

Genome word framework (GWF) analysis of DNA methylation relates to methylation motif analysis and digital signal processing. GWF is a DNA sequence fragment covering a differential methylation position (DMP). This analysis allows for the identification of DNA sequence methylation motifs found in genes with potential epigenetic regulatory functions, such as genes induced by environmental changes or diseases. Analytical heuristics can be implemented and used to identify DNA sequences of methylation motifs, with higher order DNA base interdependencies with respect to methylated cytosine, as well as statistically non-random base distributions. The discovery lays a foundation for further model prediction. For example, such model predictions may be used to identify and treat patients suffering from autism, cancer, and other diseases that benefit from early diagnosis.

Description

[0001] CROSS-REFERENCE TO RELATED APPLICATIONS

[0002] This application claims priority under 35 U.S.C. §119 to U.S. Provisional Patent Application Serial No. 63 / 323,690, filed on March 25, 2022. The entirety of this provisional patent application is incorporated herein by reference, including but not limited to the specification, claims, and abstract, and any figures, tables, appendices, or diagrams therein.

[0003] STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH OR DEVELOPMENT

[0004] This invention was made with government support under Grant No. GM134056 awarded by the National Institutes of Health. The government has certain rights in this invention. Technical Field

[0005] The present disclosure relates generally to improvements in computer science that have application to any industry that can benefit from the study of genes, phenotypes and / or DNA / RNA.More particularly, but not exclusively, the present disclosure relates to genomic word frame analysis of genomic methylation data. Background Art

[0006] The background description provided herein provides context for the present disclosure. The work of the presently named inventors and aspects of the description that may not have originally qualified as prior art at the time of filing are neither explicitly nor implicitly admitted as prior art.

[0007] DNA sequences carry not only information about how to build proteins, but also regulatory information about the survival and reproduction of organisms, including but not limited to epigenetic information that controls the behavior of chromatin.

[0008] The existence of many diseases and many phenotypes is unclear because the phenotypic penetrance and expression of a gene vary depending on the combination of modifying alleles present in one genetic background versus another.

[0009] Therefore, there is a need in the art for a device to address the possibility of superimposing more than one structured language in the nucleotide sequence of a genome. Summary of the invention

[0010] The following objectives, features, advantages, aspects and / or embodiments are not exhaustive and do not limit the entire disclosure. No single embodiment is required to provide each and every objective, feature or advantage. Any object, feature, advantage, aspect and / or embodiment disclosed herein may be integrated with each other in whole or in part.

[0011] The primary objectives, features and / or advantages of the present disclosure are to improve or overcome deficiencies in the art.

[0012] Another object, feature and / or advantage of the present disclosure is to provide an extension for a general programming language or a statistical programming language. In an embodiment, the extension includes an algorithm for analyzing methylation signals on DNA sequence fragments. In an embodiment, these DNA sequences are characterized by (i) methylation information and (ii) physicochemical information around each methylated cytosine. In an embodiment, the algorithm includes one or more functions that can estimate a distance matrix for a selected set of regions of the DNA sequence; perform a hierarchical clustering analysis on the selected set of regions; group the selected set of regions into a specified number of clusters; and align multiple DNA sequences from these clusters into methylation motifs.

[0013] In some embodiments, the extension is written in the R statistical language.

[0014] Another object, feature and / or advantage of the present disclosure is to reduce the amount of computation. For example, differentially methylated genes (GLM) identified by methylation analysis can be integrated into a gene network to identify network hubs via protein-protein interaction network analysis and weighted correlation network analysis. There is no precedent for weighted correlation network analysis based on prior knowledge of DMP and DMG so far.

[0015] Another object, feature and / or advantage of the present disclosure is to provide a computerized heuristic. In some embodiments, the computerized heuristic includes high-order DNA base interdependencies for methylated cytosines; and a statistically non-random base distribution.

[0016] According to some aspects of the present disclosure, the heuristic method may include (1) a statistic (sum, mean, density, etc.) of information divergence (ID) estimated for each gene that carries at least one DMP on it (on the gene body or on the promoter region); (2) principal component analysis (PCA), in which the top k components carrying 1% or more of the total sample variance are considered in downstream analysis; (3) calculation of a correlation matrix of pairwise gene correlations carrying vectors represented as PCs; (4) analysis of the correlation matrix of the network; and (5) the contribution of each gene to phenotypic discrimination, which is evaluated based on the fraction of the cumulative variance that accounts for the total sample variance carried by the gene.

[0017] According to some additional aspects of the present disclosure, the ID is selected from the group consisting of: Hellinger divergence / distance, J divergence, total variation distance, etc.

[0018] According to some additional aspects of the present disclosure, PCA can be applied with the pcaLDA function described in the '986 patent. Thus, genes are represented as k-dimensional vectors of PCs, where the square of each coordinate carries the vector contribution (in terms of variance) to the treatment versus control group distinction. The correlation matrix can be mathematically equivalent (in terms of information) to a weighted correlation network (WCN).

[0019] According to some additional aspects of the present disclosure, the WCN is analyzed as done for the network, which can be a PPI network. The new knowledge retrieved from the WCN is derived from the raw methylation data, and it does not rely on our beliefs or biological knowledge about the genes represented in the network. The results from the WCN and PPI networks were compared with the identification of consistent relationships and epigenetic gene contributions to the phenotype.

[0020] According to some additional aspects of the present disclosure, the heuristic method further includes calculating the magnitude as the Euclidean norm of the gene represented as a vector of k PCs.

[0021] Another object, feature and / or advantage of the present disclosure is to selectively construct motif libraries. For example, methylation motifs in all DMGs can be identified, which provide raw materials for constructing motif libraries. These libraries can then serve as the required basic data sets for predictive models used in building plant science and biomedical research.

[0022] The genomic word frameworks and genomic methylation data disclosed herein can be used for a variety of applications. For example, such GWF-based model predictions can be used to identify and treat patients with autism, cancer, and other diseases that would benefit from early diagnosis. The models can also help provide further understanding to discover the causes of (1) currently unknown phenotypes and (2) multifactorial diseases that appear to be caused by genetic and environmental factors (such as diabetes and alcoholism).

[0023] The visual presentation of the genomic word frame analysis can be automatically and intuitively configured to quickly convey meaning to the interpreter. Thus, at least one embodiment disclosed herein can include a unique aesthetic appearance. The decorative aspects included in such an embodiment can help people further understand the potential relationship between genomic methylation data and applications in the physical world (e.g., phenotypes).

[0024] Methods that facilitate the use, production, assembly, maintenance, and repair of DNA methylation motif libraries can be implemented to achieve some or all of the foregoing objectives.

[0025] Another object, feature and / or advantage of the present disclosure is to provide a method for analyzing methylation signals on DNA sequence fragments. In some embodiments, the method includes performing a hierarchical clustering analysis on regions of the DNA sequence; hierarchically grouping a set of selected regions into a specified number of clusters; aligning potential DNA sequence motifs from the clusters; and applying digital signal processing to the encoded methylation and physicochemical signals.

[0026] Library creation and maintenance can be further incorporated into an automated, heuristic analysis process that continually refines and improves the DNA base interdependencies regarding methylated cytosines until they achieve a statistically non-random base distribution.

[0027] These and / or other objects, features, advantages, aspects and / or embodiments will become apparent to those skilled in the art after reading the following brief and detailed description of the drawings. In addition, the present disclosure encompasses aspects and / or embodiments that are not explicitly disclosed but can be understood by reading the present disclosure, including at least: (a) combinations of the disclosed aspects and / or embodiments and / or (b) reasonable modifications not shown or described. BRIEF DESCRIPTION OF THE DRAWINGS

[0028] Several embodiments in which the present disclosure may be practiced are illustrated and described in detail, wherein like reference numerals represent like components throughout the several views. The drawings are presented for illustrative purposes and are not to scale unless otherwise specified.

[0029] Figure 1 The msh1 system consisting of four msh1-derived epigenetic states is shown. In Arabidopsis, four distinct plant states are derived from MSH1 knockdown or knockout. States 1 and 2 are directly derived from msh1 disruption, resulting in highly stress-responsive phenotypes. State 1 is variable under short days and includes a low-frequency “perennial-like” phenotype16. States 3 and 4 involve the interaction of the msh1-modified genome with the native (wild-type) genome through grafting or hybridization, resulting in growth vigor phenotypes.

[0030] Figures 2A to 2HThe characteristics of epi strains derived by hybridization of msh1 T-DNA mutants with isogenic wild type are shown.2A shows the phenotype of different epi strain F3 populations when days after planting (days after planting, DAP) is 34. These strains are derived from WT x msh1 hybridization, wherein Epi 24 and Epi 8 are from a parental hybridization, and Epi10 and Epi19 are from a secondary parental hybridization.2B, 2C and 2D show the total leaf area (34DAP), days to bolting, seed weight (mg) of the four populations and WT controls.2B to 2D show the bars representing mean ± SE.The Mann-Whitney U-test (Mann-Whitney U-test) with a two-sided alternative hypothesis is used to test the significance of the mean value difference between each Epi F3 population and the WT control.2E shows the root phenotype (33DAP) of four different Epi F3 populations grown in sand. 2F shows the total leaf area (33DAP), dry leaf weight (mg), and dry root weight (mg) of the four populations grown on sand and the WT control. 2F to 2H show that the bars represent the mean ± SE. The Mann-Whitney U-test with a two-sided alternative hypothesis was used to test the significance of the mean difference between each Epi F3 population and the WT control. Significance code: *p<0.05, **p<0.01, ***p<0.001, ns-not significant.

[0031] Figures 3A to 3H Tomato plants for methylome sequencing are shown. 3A shows Rutgers wild-type plants at 4 weeks. 3B to 3D show comparisons of Rutgers msh1 memory plants (state 2, on the right side of each figure) with WT at 4 weeks. 3E to 3G show comparisons of Rutgers msh1 RNAi transgenic plants (state 1, on the right side of each figure) with WT. 3H shows epi lines with higher vegetative growth (3 plants on the left), represented as high epi lines (Epi-lineHigh, Epi-H) in subsequent analysis, and epi lines with lower vegetative growth (3 plants on the right), represented as low epi lines (Epi-line Low, Epi-L). Young leaves of 4-week-old plants were collected for whole-genome bisulfite sequencing.

[0032] Figures 4A to 4CShown the reversion phenotype in Arabidopsis and soybean.4A shows the plant growth phenotype of three F3 epi strain populations, Epi 10 and Epi 6 derived by hybridization with msh1 T-DNA mutant and Epi 9 derived by hybridization with msh1 memory strain.The dotted circle represents the revertant of inference.Col-0 wild type and msh1 memory are shown as control.4B shows soybean (cv.Thorne) F3 epi strain of field planting, shows the evidence (arrow) of spontaneous reversion.Harvested in October after severe cold, revertant shows the cold tolerance similar to msh1 memory.4C shows the principal component analysis-linear discriminant analysis of methylation group data from three epi strain non-revertants and three revertant full sibling F3 progeny compared to three independently derived msh1 memory plants.

[0033] Figure 5 The relative DMP frequencies between different Arabidopsis epi lines are shown. DMPs were assigned to genic regions, TE-associated regions, and other genomic regions. The centroid of the wild-type sample was used as a reference. The relative DMP frequency in each genomic feature was estimated as the number of DMPs divided by the number of total genomic cytosine positions in each genomic feature. A DMP was defined as “high” if the difference in cytosine methylation levels was greater than 0, and as “low” if the difference in methylation levels was less than 0.

[0034] Figures 6A to 6C The total high and low DMP counts in the epi lines are shown. 6A shows a comparison of Arabidopsis epi lines with WT. 6B shows a comparison of tomato low epi lines (BE) with high epi lines (GE). 6C shows a comparison of soybean epi-F4 with epi-F6. Each bar graph represents a single plant in a given population. The cytosine background (CG, CHG, CHH) of each plant is shown separately. If the difference in cytosine methylation levels is greater than 0, DMP is defined as "high"; and if the difference in methylation levels is less than 0, DMP is defined as "low". The control centroid is used as a reference.

[0035] Figures 7A to 7C Figure 7A shows the hierarchical clustering results of gene methylation group data from three different msh1 states in Arabidopsis: msh1 null mutant (state 1), Col-0 / Col-0 msh1Grafted progeny (state 3), Epi F3 population (state 4), and related Col-0 controls. 7B shows the hierarchical clustering results of the individual plant (P) dataset from three different msh1 states in tomato: Rutgers MSH1-RNAi+ (state 1), Rutgers / Rutgers msh1-RNAi Grafted (HEG) progeny (state 3) and Epi F3 populations (Epi-L and Epi-H) (state 4) and associated Rutgers controls. 7C shows the hierarchical clustering results of TEDMP data from the same Arabidopsis dataset presented in 7A. Individuals are represented as vectors of the sum of Hellinger divergences (HD) at DMP positions within gene regions (7A, 7B) or TE regions (7C). Hierarchical clusters were constructed using the Ward agglomerative method and the centroids of the corresponding wild-type samples were used to calculate Hellinger divergences (HD). The HD formula is reported elsewhere.

[0036] Figures 8A to 8E Shown are GO pathways significantly enriched by DMG analysis in epi lines.8A shows the enriched GO pathways shared by Epi 8 and Epi 24 in Arabidopsis.8B shows the enriched GO pathways shared by Epi 10 and Epi 19 in Arabidopsis.8C shows the enriched GO pathways identified by the F1 hybrid DMG from C24 x Ler hybridization by our analysis.The raw methylation data of F1 hybrids have been reported previously.8D shows the enriched GO pathways derived from the DMG analysis of low epi lines and high epi lines in tomato.8E shows the enriched GO pathways based on the DMG analysis of Epi-F4 and Epi-F6 comparisons in soybean.For 8A to 8C, the bar graphs represent the number of DMGs, and the dotted lines represent the -log10 (p value) of the enrichment test.For 8D to 8E, the bar graphs represent the number of DMGs, and the dotted lines represent the enrichment multiples of each enriched pathway. The enriched GO pathways were identified using the Network Enrichment Analysis Test (NEAT). The p-values ​​were calculated using the one-sided hypergeometric distribution test implemented in NEAT. Thus a complete list of enriched pathways can be provided.

[0037] Figures 9A to 9C Shown are the GO pathways of significant enrichment based on DEG analysis of epi lines in Arabidopsis.9A shows the enrichment GO pathways of DEGs shared from Epi 8 and Epi 10 leaf tissues.9B shows the enrichment GO pathways of DEGs shared from Epi 8 and Epi 24 flower stem tissues, and 9C shows root tissues.Bar graphs represent the number of DEGs, and dotted lines represent the enrichment multiples of each enrichment pathway.One-sided Fisher's exact test is used to calculate the FDR implemented in DAVID GO.Therefore a complete list of enrichment pathways can be provided.

[0038] Figures 10A to 10B Protein-protein interaction (PPI) networks from the hubs of a subset of network-associated DMGs and DEGs in Arabidopsis Epi 8 and Epi 24 full-sib lines are shown. Epi 8 and Epi 24 represent daughter lines derived from the same cross, with Epi 8 exhibiting enhanced vegetative growth rates and Epi 24 significantly enhancing seed yield. The main subnetworks of the hubs were obtained by applying machine learning K-means clustering to a collection of 3647 (Epi 8, 10A) and 3523 (Epi24, 10B) network-associated DMGs and DEGs identified in Arabidopsis epi lines compared to WT. The analysis generated 153 (10A) and 346 (10B) hub genes, forming closely related subnetworks. The GO network enrichment analysis of the string application in Cytoscape was used to identify gene functional pathways enriched within the network. The size of each node is proportional to its node degree value, while the label font size is proportional to its betweenness centrality.

[0039] Figures 11A to 11B Protein-protein interaction (PPI) networks derived from hubs of a subset of network-associated DMGs and DEGs in soybean epi lines and tomato grafted offspring (HEG) lines are shown. 11A shows the main subnetworks of the hubs obtained by applying a machine learning K-means clustering algorithm to a set of 1264 network-associated DMGs and DEGs identified in a comparison of soybean Epi-F6 and Epi-F4 (state 4). The analysis generated 421 hub genes that formed closely related subnetworks. GO network enrichment analysis using the string application in Cytoscape was used to identify gene functional pathways enriched within the networks. 109 genes related to gene expression regulation were selected for display. 11B ... comparing the network functions of the 1264 network-associated DMGs and DEGs in tomato Rutgers / Rutgers msh1-RNAi The main subnetwork of hubs obtained by applying the machine learning K-means clustering algorithm to 3173 network-related DMGs and DEGs identified in comparison with the Rutgers / Rutgers graft (HEG) offspring (state 3). The analysis generated 305 hub genes, forming subnetworks. GO network enrichment analysis using the string application in Cytoscape was used to identify gene functional pathways enriched within the network. 126 genes related to gene expression regulation were selected for display. For both 11A and 11B, the size of each node is proportional to its node degree value, while the label font size is proportional to its betweenness centrality. The datasets of soybean DEGs and tomato HEGs (DMGs and DEGs) have been reported previously.

[0040] Figures 12A to 12DThe relationship between 871 core hub genes and 4 different msh1-derived states and a core network with biological significance are shown. 12A shows a Venn diagram of Arabidopsis DMGs identified from four different msh1-derived states (Col-0 genetic background): msh1 mutants (state 1), msh1 memory (state 2), grafted progeny (HEG, state 3) and epi strains (Epi 24, state 4). 12B shows an overview of the PPI network and individual 871 core hub genes. 12C shows the hierarchical clustering of individual plant data sets from four different msh1-derived states based on the sum of Bayesian methylation level differences of DMPs from 871 core genes from 12A. Hierarchical clustering was constructed using the Ward agglomeration method. Bayesian methylation level differences were calculated as described previously. 12D shows the main subnetwork of the hub obtained by applying the machine learning K-means clustering algorithm to the set of 871 core genes from 12A. GO network enrichment analysis using the string application in Cytoscape was used to identify enriched gene functional pathways within the network. 67 genes involved in the enriched network were identified. This 67-gene set provides RdDM candidate target genes for further study. For both 12B and 12D, the size of each node is proportional to its node degree value, and the label font size is proportional to its betweenness centrality.

[0041] Figures 13A to 13DShown the research of gene methylation remodeling (repatterning) in candidate RdDM target gene, it distinguishes four kinds of msh1 derived states in Arabidopsis.13A shows two inferred cluster motifs identified based on differential gene methylation in four kinds of msh1 derived states.Hierarchical clustering is carried out to the set of 14-bp region comprising DMP, these regions are from seven (TOR, SYD, NRPB1, NRPD1, UBP26, SUVR5, UPF1) in 67 core hub loci, then DNA multiple sequence alignment is carried out to each cluster, thereby allows identification methylation motif.13B shows the difference of methylation level of gene body DMP in motif cluster 11 in inferred RdDM target gene UBP26.The variation of DMP motif methylation remodeling is shown by chromosome and position.For each plant measured under each msh1 state, the methylation variation detected separately is shown as the point with cross hatching, wherein positive indication DMP, and negative indication does not detect methylation variation.Each line represents single plant data set. 13C shows sample DMPs within motif cluster 11 in UBP26 and UPF1 that show dcl2 / dcl3 / dcl4 sensitivity in grafted rootstock (state 1) and grafted progeny (state 3) from a comparative mutant rootstock experiment. Note that dcl2 / dcl3 / dcl4 effects can only be determined for msh1 and grafted progeny data. 13D shows sample putative cluster motifs identified based on differential gene methylation in the four msh1-derived states, following analysis of the 67 core central loci and multiple DNA sequence alignment of methylation motifs for each cluster.

[0042] Figures 14A to 14D Exemplary cluster motifs containing DMPs on the gene AT1G50030 (TOR) in Arabidopsis thaliana are shown. 14A shows the cluster motif for cluster 1. 14B shows the cluster motif for cluster 9. 14C shows the cluster motif for cluster 12. 14D shows the cluster motif for cluster 17. Individual consensus nucleotides are evident at different distances from the target cytosine, which is nucleotide 7 within each motif (on either the '+' strand or the '-' strand).

[0043] Figures 15A to 15B The power spectral density of the control (15A) and treatment (msh1 mutant, dwarf phenotype) (15B) in Arabidopsis is shown. The dots help to visually identify the peaks where spectral differences between the control and the treatment can be found.

[0044] Figures 16A to 16BSpectrograms of AT1G50030 (TOR) in Arabidopsis are shown. In both 16A and 16B, the upper figure is the control and the lower figure is the treatment (dwarf phenotype). 16A is the spectrogram of the first region and 16B is the spectrogram of the third region. The methylation signal decomposes the power spectrum energy around the periodicity at about 1 / 3 frequency (white dotted line).

[0045] Figures 17A to 17C Wavelet power spectra (WPS) of regions AT1G50030.1 (17A), AT1G50030.2 (17B) and AT1G50030.3 (17C) in the control group (left graph in each figure) and the treated group (dwarf phenotype) (right graph in each figure) are shown. The differences in WPS observed between the control group and the treated group correspond to changes in bases.

[0046] Fig.18 Correlation plots based on 17A WPS are shown. Control vs. control are shown in the upper graph, and control vs. treatment (dwarf phenotype) are shown in the lower graph. In addition to the expected correlation break at region 55-60, the correlation between the energy spectra at regions 30-35 and 55-60 was also lost in the control vs. treatment.

[0047] Fig.19 Subnetworks of 81 DMG hubs are shown. The PPI network was built on a collection of 751 DMGs identified by principal component analysis and analyzed with the STRING Cytoscape application. Circles without cross-hatching highlight DMGs involved in biological processes, including nervous system, nervous system development, synapses, neuronal projections, ion transport, CNS diseases, axon guidance, neurogenesis, ion / cation binding, ion / cation transport across membranes, voltage-gated channels, and TRP channels.

[0048] Fig. 20 Shown in Fig.19 Network enrichment analysis of DMG hubs identified in the primary network.

[0049] One of ordinary skill in the art need not view the nearly infinite number of different arrangements of features described in the following detailed description in isolated figures to facilitate an understanding of the present disclosure. DETAILED DESCRIPTION

[0050] The present disclosure is not limited to what is described herein. Mechanical, electrical, chemical, procedural and / or other changes may be made without departing from the spirit and scope of the present disclosure. Unless otherwise specified, the features shown or described are not essential to allow the basic operation of the present disclosure.

[0051] The genomic word framework (GWF) analysis of DNA methylation relates to methylation motif analysis and digital signal processing. GWF is a DNA sequence fragment covering differential methylation positions (DMPs). However, GWF should not be confused with the concept of DNA sequence motifs. Word framework (WF) can include one or more motifs. In other words, "WF sentences" are also GWFs. DMPs can be identified by methylation analysis utilizing an extension program of a statistical programming language, such as the R package described by the inventor of the present application in U.S. Patent No. 10,913,986. This analysis allows identification of DNA sequence methylation motifs found in genes with potential epigenetic regulatory functions, including genes induced by environmental changes or diseases.

[0052] One potential embodiment of the analysis heuristics described herein has been implemented in an R package named GenomicWordFramework. More specifically, GenomicWordFramework is a utility package for identifying potential genomic word framework (GWF) regions for hypothetical languages. GWF facilitates further analysis of DNA sequences carrying genomic signals by applying digital signal processing (DSP) and machine learning (ML) tools in other R packages. GenomicWordFramework includes several functions to complete the reading and conversion of data into a suitable form so that different statistical methods can be applied for data analysis, such as clustering algorithms and statistical tests.

[0053] The GenomicWordFramework can pre-identify DMPs using R packages dedicated to methylation analysis, such as those described in U.S. Pat. No. 10,913,986. Thus, using the methylation analysis and GWF analysis described herein can form a pipeline that implements signal detection and machine learning methods that allow filtering signals from noise at a high rate.

[0054] The discriminatory power of the GenomicWordFramework was first tested in a small dataset of the Arabidopsis msh1 mutant system (a model organism) and subsequently applied to a published study analyzing DNA methylation in placental tissue from typically developing children and children with autism. The GenomicWordFramework analysis is a proven concept that is compatible with a virtually unlimited number of differentially methylated network hubs. The differentially methylated network hubs previously identified by methylation and network analysis in autism research are related to biological processes including: nervous system, nervous system development, synapses, neuronal projections, central nervous system diseases, axon guidance, neurogenesis, ion / cation binding, ion / cation transmembrane transport, and voltage-gated channels.

[0055] The results show that the GWF-based heuristic can identify DNA sequences with methylation motifs, high-order DNA base interdependencies with respect to methylated cytosines, and statistically non-random base distributions. These findings lay the foundation for further model predictions for patients, not only for autism, but also for other diseases such as cancer that would benefit from early diagnosis.

[0056] In other words, GWF analysis enables the identification of sets / clusters of synonymous methylation term frameworks within genes that undergo targeted methylation changes and participate in gene networks involving biological processes relevant to the system under study. GWF further lays the foundation for the creation of DNA methylation motif libraries to be used for patient diagnosis and prognosis. Therefore, GWF analysis greatly improves the utility and value of identifying DMPs and differentially methylated genes (DMGs) through methylation analysis.

[0057] GenomicWordFramework R package and GWF identification

[0058] GenomicWordFramework is an extension written in a static programming language. More specifically, GenomicWordFramework is an R package designed to analyze methylation signals on DNA sequence fragments that are characterized by not only methylation information but also the physicochemical information around each methylated cytosine (and adenine in the case of animals and bacteria). The package's functions convert methylation data into a suitable format for further DSP analysis outside of the R package.

[0059] GWFs are identified in the methylome in two possible ways: (1) applying the algorithm described by Sanchez et al. in "Information Thermodynamics of Cytosine DNA Methylation" (published on March 10, 2016), which is incorporated herein by reference in its entirety; and (2) directly extracting DNA sequence fragments covering a specific number of DNA bases upstream and downstream of the DMP. The GWFs obtained using the first method can be further analyzed using digital signal processing (DSP) tools. The GWFs obtained using the second method can be clustered into alignment motif groups and further tested to evaluate the deviation of each in the multiple sequence alignment (MSA) from the random Monte Carlo simulated MSA. Statistically significant motif regions are expanded for further application of DSP analysis.

[0060] Binary encoding of methylation signals

[0061] All the information needed to search for DNA motifs, DNA sequences, and modified nucleotide bases can be stored in a single binary string, using a triplet representation of each base, taking into account the number of hydrogen bonds in Watson-Crick base pairs, the chemical type of the DNA base (pyrimidine and purine), and the base modification state. By convention, the DNA sequence is called the positive strand. For example, the following DNA base encoding can be used to study methylation signals:

[0062] Methylation status: 1 (methylated) and 0 (unmethylated)

[0063] Strong and weak base pair interactions: 1 (three hydrogen bonds, S type) and 0 (two hydrogen bonds, W type)

[0064] Chemical types: 1 (pyrimidine, Y class) and 0 (purine, R class).

[0065] Free radical content: 1 (base containing amino group (NH_3, M type)) and 0 (base containing keto group (C=O) group (K type))

[0066] Each base will form a binary string:

[0067] Signal mC on the plus strand, a methylated pyrimidine with three hydrogen bonds: 111.

[0068] Signal mC on the negative strand, methylated pyrimidine with three hydrogen bonds ('-' base G on the positive strand): 110.

[0069] C, pyrimidine with three hydrogen bonds, M class: 011

[0070] G, purine with three hydrogen bonds, K class: 010

[0071] T, pyrimidine with two hydrogen bonds, K class: 001

[0072] Signal mA on the positive strand, methylated purine with two hydrogen bonds: 100.

[0073] Signal mA on the negative strand, methylated purine with two hydrogen bonds ('+' base T on the positive strand): 101

[0074] A, purine with two hydrogen bonds, M class: 000

[0075] ·U, same as 'T': 001.

[0076] Each base will form a complex number (a+ib), which is a cyclic group integrated by the eighth root of unity: Where k = 0,…,7.

[0077] ·Signal mC on the positive strand, a methylated pyrimidine with three hydrogen bonds: 1.

[0078] • Signal mC on the negative strand, methylated pyrimidine with three hydrogen bonds ('+' base G on the positive strand): -1.

[0079] C, pyrimidine with three hydrogen bonds, S type:

[0080] G, purine with three hydrogen bonds, S type:

[0081] T, pyrimidine with two hydrogen bonds, W class:

[0082] • Signal mA on the positive strand, methylated purine with two hydrogen bonds: i.

[0083] Signal mA on the negative strand, methylated purine with two hydrogen bonds ('+' base T on the positive strand): -i

[0084] A, purine with two hydrogen bonds, class W:

[0085] U, same as 'T':

[0086] There may be more than one complex encoding. For DSP analysis, an encoding based on a group structure (here the Abelian group) may be more desirable. This is a group structure defined for the set of DNA bases that include methylated cytosine and adenine members of the alphabet, as described further below. Complex signals may be difficult to analyze in the GWF R package, and the complex encoding and signals can be exported to other languages ​​later, such as C++, Python or MatLab.

[0087] Specific details for implementing this encoding are provided in the documentation for the function, which encodes the 0 and 1 binary signals previously detected from the DNA sequence into a user-defined numerical code. Given two objects, one carrying the signal and the other carrying the DNA sequence, such a function can perform the encoding set by the user. Binary encoding of the methylation signal allows the incorporation of physicochemical information into the DSP analysis of DNA sequence motifs.

[0088] Motif clustering

[0089] Small regions spanning DMPs, typically 7-30 bp on at least three samples, are identified and considered as DNA methylation motif candidates in 67 genes. A distance matrix can be estimated for a selected region set using a function in the R package, which calculates the matrix distances between the aligned sequences from each multiple sequence alignment (MSA). Next, a hierarchical clustering analysis of the selected region set (using the distance matrix estimated before) can be completed by utilizing a function of a selected information divergence matrix to group the selected regions into 100 clusters. It should be understood that Hellinger divergence is only one of the possible information divergences that can be estimated and applied here. For example, J divergence is more suitable for applications that are intended to extract new knowledge in terms of information thermodynamics of epigenomic phenomena. J divergence is a symmetric version of relative entropy.

[0090] The unweighted pairwise method of arithmetic mean ("UPGMA") was used as the agglomerative algorithm. In one embodiment, clusters with less than ten (10) regions were discarded. DNA multiple sequence alignment on each sequence cluster can be performed using the MUltiple Sequence Comparison by Log-Expectation (MUSCLE) algorithm implemented in the R package to analyze and understand genomic data generated by molecular biology wet lab experiments.

[0091] Further partitioning of the motif set may be undertaken for downstream analysis. A variety of clustering algorithms may be applied. For example, a distance-based k-medoid clustering algorithm may be used to apply and implement fast k-medoid clustering: such as simple and fast k-medoid clustering, sorted k-medoid clustering, and k-medoid clustering that increases the number of clusters. In one embodiment, the algorithm is an algorithm included in the kmed R package.

[0092] The clustering results can be plotted in a labeled bar graph or a principal component analysis biplot. The final clustering division depends on the applied clustering algorithm and its corresponding parameter settings, including the type of metric applied to calculate the distance matrix required by the clustering algorithm. Methylation motifs are objective DNA sequence features, and the applied clustering algorithm is only a supporting tool to achieve motif identification.

[0093] Motif score

[0094] Motif score s of aligned sequences j and k jk It can be defined in an intuitive way as the base-2 logarithm of the number of matching bases found in the alignment. Formally:

[0095]

[0096] in For each base position i on sequences j and k, the maximum motif score is: Max{s jk} = log2N. Next, the motif score in MSA is defined as:

[0097]

[0098] For an MSA with M sequences of length N, the number of pairwise comparisons is: Therefore, for a fixed value of motif size, a perfect MSA of a DNA sequence motif will have a maximum score: Max{S} = log2N. That is, in this modeling, the maximum amount of information carried by MSA is: I max =log2NI max =log2N, and the amount of information (uncertainty change) carried by MSA is given by the following expression:

[0099] I=-(S-log2N)

[0100] The same results are obtained if instead of matching, the letter frequencies in MSA and Shannon entropy before (perfect alignment) and after are estimated. The alignment information is then calculated as:

[0101] I=-(H 之后 -H 之前 )

[0102] Monte Carlos Testing of MSA Randomness

[0103] The Monte Carlo test (MCT) of how a given DNA multiple sequence alignment differs from a randomly generated MSA was implemented in the R function contained in the GenomicWordFramework R package. To perform the MCT, assume that:

[0104] A count vector n summarizing the observed DNA base frequencies in N columns k (k=A,C,G,T) will be based on the parameter p k The multinomial distribution of .

[0105] In the Bayesian framework, since the Dirichlet distribution is the conjugate prior distribution of the multinomial distribution, we assume one of two ways to generate the parameter vector:

[0106] o If the number of genes or genomic regions from which motifs are identified is small: parameter vector P = (pA,

[0107] pC, pG, pT) are derived from the Dirichlet distribution D_αP, where the parameter vector α = (α1, α2, α3,

[0108] α4).

[0109] o If the number of genes or genomic regions from which motifs are identified is large: parameter vector P = (pA,

[0110] pC, pG, pT) are drawn from a Dirichlet mixture distribution, where each Dirichlet component D_(α_iQ_i) has parameters α_i = (α_i1, α_i2, α_i3, α_i4) and Q = (qA, qC, qG, qT). Formally:

[0111]

[0112] where ∑ i π i =1.

[0113] DNA bases are generated independently according to the distribution pk.

[0114] Dirichlet distribution estimation

[0115] The parameter vector of a specific Dirichlet distribution is α = (α1,…,α k ) is estimated based on the entire set of identified DNA motif candidates. Given a DNA methylation motif matrix with N columns corresponding to the entire set of identified motif candidates, the frequency of each DNA base in each column is calculated, resulting in a four-dimensional count vector for each column in the dataset, where each vector coordinate carries the absolute frequency of one of the four DNA bases in a given alignment column. The resulting N×4 count matrix is ​​the raw count data used in the estimation of the parameters of the Dirichlet distribution, and the functions in the R package are applied to estimate the multivariate generalization of the continuous multivariate probability distribution family: the Beta distribution. Next, random DNA MSA sequences are generated based on the estimated Dirichlet distribution and the probability density function (PDF) or cumulative density function (CDF) in the R package. That is, the random DNA MSA is generated by sampling from the estimated Dirichlet distribution.

[0116] Monte Carlo p-value

[0117] For MSAs of fixed length N, log2N is a constant, so it is sufficient to proceed as follows for the purposes of MCT: Consider the score statistic given in Equation 2 and evaluate how statistically different the observed alignment is from the Monte Carlo simulated alignment. The Monte Carlo p-value is estimated as:

[0118]

[0119] Where S0 represents the alignment score of the MSA to be tested, S iis the alignment score of the ith Monte Carlo simulation MSA, and

[0120] It is worth noting that the raw frequencies observed from a small motif matrix are usually not very close to the distribution of DNA bases across all motifs that the model is supposed to represent. However, for a typical analysis of 50 genes or more, the motif matrix required for Dirichlet distribution model estimation will typically carry thousands of motifs.

[0121] Digital signal processing of binary-coded methylation signals

[0122] The binary coded methylation signal is the raw data for digital signal processing (DSP) tools. DSP tools have a large number of possible applications. The GenomicWordFramework R package can include some applications of wavelet spectrograms via wavelet transforms, as well as traditional Fourier power spectra and spectrograms.

[0123] Multiplicative groups of the expanded DNA alphabet with methylated bases

[0124] make is an ordered set of the following: DNA bases plus methylated adenine (A m ) and cytosine (C m ) and methylated adenine in the negative strand and cytosine Let's use the multiplication operation '×' Define a multiplication group Among them C m It is a multiplication unit, and the unmethylated complementary bases of DNA are also algebraically complementary, that is: C×G=C m and A×T=C m If base C is a group generator, that is, generator condition C 8 =C m And the condition C×G=C m Means C 7 =G, then the final algebraic biophysical constraint holds. Since C 3 ×C 5 =C 8 , set C 3 =T means C 5 =A. To preserve the order of the bases in the set, we let C×C=C 2 =A m ,

[0125] The group defined above is an Abelian cyclic group, isomorphic to the cyclic group of the 8th root of unity integral: where k = 0,…,7. Although we can complete symbolic algebraic operations on , but for specific applications in computational biology and bioinformatics, the set It is more convenient to operate on the cyclic group defined above. The elements of this group are bijective mappings Write a sequence set as: where i is the imaginary unit defined in the complex number set.

[0126] Complex encoding of methylation signals using the GWF R package

[0127] Library (GenomicWordFramework)

[0128] We will use methylation signals from three gene regions in Arabidopsis

[0129] data(at_signal,at_gene_seq,package="GenomicWordFramework")

[0130] Each base will form a complex number (a+ib), which is a cyclic group integrated by the eighth root of unity. Where k = 0,…,7.

[0131] ·Signal mC on the positive strand, a methylated pyrimidine with three hydrogen bonds: 1.

[0132] • Signal mC on the negative strand, methylated pyrimidine with three hydrogen bonds ('+' base G on the positive strand): -1.

[0133] C, pyrimidine with three hydrogen bonds, S type:

[0134] G, purine with three hydrogen bonds, S type:

[0135] T, pyrimidine with two hydrogen bonds, W class:

[0136] • Signal mA on the positive strand, methylated purine with two hydrogen bonds: i.

[0137] Signal mA on the negative strand, methylated purine with two hydrogen bonds ('+' base T on the positive strand): -i

[0138] A, purine with two hydrogen bonds, class W:

[0139] U, same as 'T':

[0140] In the GWF R package, the encoding of methylation signals is done by the function signalEncoding. In the current case, we use the group generator: To easily set the encoding:

[0141]

[0142]

[0143] The signal to be exported can be retrieved using the getEncoding function:

[0144]

[0145] The entire matrix of the encoded signal can be retrieved:

[0146]

[0147]

[0148] The matrix carrying the signal can be exported to Python or MatLab.

[0149]

[0150] The entire exportable numeric matrix (showing two rows and only the first 20 columns):

[0151] Example

[0152] The following non-limiting numbered examples also form part of the present disclosure:

[0153] 1. An extension program for a general programming language or a statistical programming language, the extension program comprising: an algorithm for analyzing methylation signals on a DNA sequence fragment, the DNA sequence being characterized by:

[0154] i. methylation information; and

[0155] ii. The physicochemical information around each methylated cytosine;

[0156] Wherein the algorithm comprises one or more functions, the one or more functions may:

[0157] estimating a distance matrix for a set of selected regions of the DNA sequence;

[0158] performing a hierarchical cluster analysis on the selected region set;

[0159] grouping the selected set of regions into a specified number of clusters; and

[0160] Multiple DNA sequences from the clusters are aligned into methylation motifs.

[0161] 2. The extended program according to embodiment 1, wherein the extended program is written in R statistical language.

[0162] 3. An extension program according to any one of embodiments 1 to 2, further comprising a digital signal processing (DSP) tool usable in another programming language.

[0163] 4. The extension program according to embodiment 3, wherein the other programming language is C++, Python or MatLab.

[0164] 5. The expansion program according to any one of embodiments 1 to 4, further comprising:

[0165] Encoding the physical and chemical properties of DNA sequences and DNA bases into numerical or complex signals further includes DSP.

[0166] 6. An extension procedure according to embodiment 5, wherein the encoding is based on a group structure.

[0167] 7. The extension procedure according to embodiment 6, wherein the group structure is an Abelian group.

[0168] 8. The expansion procedure according to any one of embodiments 1 to 7, wherein the set of selected regions are grouped into groups of at least 100 clusters.

[0169] 9. An extension procedure according to any one of embodiments 1 to 8, wherein clusters with less than ten regions are discarded.

[0170] 10. The extended procedure according to any one of embodiments 1 to 9, wherein the plurality of DNA sequences are aligned using multiple sequence comparison performed by a logarithmic expectation algorithm.

[0171] 11. The expansion program according to embodiment 5 further comprises:

[0172] Export the numerical or complex signal to C++, Python or MatLab.

[0173] 12. The extended procedure according to any one of embodiments 1 to 12, wherein the methylation motifs are further grouped for downstream analysis.

[0174] 13. An extended procedure according to embodiment 12, wherein the methylation motifs are further grouped using a clustering algorithm.

[0175] 14. An extended program according to embodiment 13, wherein the clustering algorithm is a distance-based k-center point clustering algorithm.

[0176] 15. The expansion program according to any one of embodiments 1 to 14, further comprising:

[0177] Encode the methylation and physicochemical signals of DNA bases into a binary, numerical or complex signal.

[0178] 16. The extension procedure of embodiment 15, wherein the encoding is based on a group structure.

[0179] 17. The extension procedure according to embodiment 16, wherein the group structure is an Abelian group.

[0180] 18. The expansion program according to any one of embodiments 15 to 17, further comprising:

[0181] The encoded methylation and physicochemical signals were subjected to power spectrum analysis.

[0182] 19. An extended procedure according to embodiment 18, wherein the power spectrum analysis is wavelet power spectrum analysis (WPS).

[0183] 20. The expansion program according to embodiment 19 further comprises:

[0184] A correlogram is calculated based on the WPS.

[0185] 21. A computerized heuristic method comprising:

[0186] Higher-order DNA base interdependencies regarding methylated cytosine; and

[0187] Statistically non-random distribution of bases.

[0188] 22. A computerized heuristic method according to embodiment 21, wherein said high-order DNA interdependencies and said base distributions are obtained from analyzing regions of DNA sequences by performing at least one hierarchical clustering and aligning multiple DNA sequences from the clusters into methylation motifs.

[0189] 23. A method for analyzing methylation signals on a DNA sequence fragment, the method comprising:

[0190] performing hierarchical cluster analysis on regions of the DNA sequence;

[0191] Hierarchically group the selected set of regions into a specified number of clusters;

[0192] aligning potential DNA sequence motifs from the clusters; and

[0193] Apply digital signal processing to the encoded methylation and physicochemical signals.

[0194] 24. The method of embodiment 23, wherein the set of selected regions are grouped into groups of at least 100 clusters.

[0195] 25. The method of any one of embodiments 23 to 24, wherein clusters with less than ten regions are discarded.

[0196] 26. The method according to any one of embodiments 23 to 25, wherein the encoded methylation and physicochemical signals are encoded based on group structure.

[0197] 27. The method of embodiment 26, wherein the group structure is an Abelian group.

[0198] 28. The method according to any one of embodiments 23 to 27, wherein the potential DNA sequences are aligned using a multiple sequence comparison performed by a logarithmic expectation algorithm.

[0199] 29. The method according to any one of embodiments 23 to 28, further comprising:

[0200] The encoded methylation and physicochemical signals were subjected to power spectrum analysis.

[0201] 30. The method according to embodiment 29, wherein the DSP is applied by using the Genome Word Framework (GWF) R package.

[0202] 31. The method of embodiment 30, wherein the GWFR package includes one or more clustering algorithms.

[0203] 32. The method of any one of embodiments 30 to 31, further comprising exporting the GWFR package for analysis with other DSP tools in different programming languages.

[0204] 33. The method of embodiment 32, wherein the different programming language is C++, Python, or MatLab.

[0205] 34. The method according to any one of embodiments 23 to 33, further comprising:

[0206] Methylation analysis is used to identify differentially methylated positions (DMPs) that constitute part of the methylation signal.

[0207] 35. The method according to any one of embodiments 23 to 34, further comprising:

[0208] Each of the multiple sequence alignments (MSAs) was evaluated for deviation from a random Monte Carlo simulated MSA.

[0209] 36. The method according to any one of embodiments 23 to 35, further comprising:

[0210] Analyze genes with epigenetic regulatory functions to diagnose or treat diseases.

[0211] 37. The method according to embodiment 36, wherein the disease is autism or cancer.

[0212] 38. A computerized heuristic method for analyzing genetic data, comprising:

[0213] (1) Information Diffusion (ID) statistics estimated for each gene carrying at least one differentially methylated position (DMP) on the gene or in the promoter region;

[0214] (2) principal component analysis (PCA), in which the top k components carrying 1% or more of the entire sample variance are considered in downstream analyses;

[0215] (3) calculating the correlation matrix of pairwise gene correlations carrying vectors represented as PCs;

[0216] (4) Analyze the network's correlation matrix; and

[0217] (5) The contribution of each gene to phenotypic discrimination, which is assessed based on the fraction of the cumulative variance accounted for by the total sample variance carried by the gene.

[0218] 39. The computerized heuristic method of embodiment 38, wherein the ID is selected from the group consisting of: Hellinger divergence / distance, J divergence, and total variance distance.

[0219] 40. A computerized heuristic according to embodiment 39, wherein the ID is the J-divergence.

[0220] 41. A computerized heuristic method according to any one of embodiments 38 to 40, wherein the PCA is applied with a function such that genes are represented as a k-dimensional vector of PCs, and further wherein the square of each coordinate carries a vector contribution in terms of variance that distinguishes the treatment from the control group.

[0221] 42. A computerized heuristic method according to embodiments 38 to 41, wherein the PCA is applied together with the pcaLDA function.

[0222] 43. The computerized heuristic method according to any one of embodiments 38 to 42, wherein the correlation matrix comprises a weighted correlation network (WCN), wherein the WCN and the network are analyzed, and the network is a protein-protein interaction (PPI) network.

[0223] 44. The computerized heuristic method of embodiment 43, further comprising comparing the results from the WCN and the PPI network with the identified consistent relationships and epigenetic gene contributions to the phenotype.

[0224] 45. The computerized heuristic method according to any one of embodiments 38 to 44, further comprising calculating the magnitude of the Euclidean norm of the gene represented as a vector of k PCs.

[0225] 46. ​​The computerized heuristic method of any one of embodiments 38 to 45, wherein the network involves one or more biological processes.

[0226] 47. A method for reducing the amount of calculation, the method comprising:

[0227] Use methylation analysis to identify differentially methylated genes (DMGs) on DNA sequence fragments;

[0228] integrating the DMG into a gene network; and

[0229] Identify network hubs.

[0230] 48. The method according to embodiment 47, wherein the network hubs are identified via protein-protein interaction (PPI) network analysis and weighted correlation network (WCN) analysis.

[0231] 49. The method of embodiment 48, further comprising comparing the results from the WCN and the PPI network to identify consistent relationships and epigenetic gene contributions to phenotypes.

[0232] 50. The method according to any one of embodiments 47 to 49, wherein the network hubs are involved in one or more biological processes.

[0233] Example 1: Arabidopsis dataset

[0234] In order to perform specific applications on raw genomic data, the concepts, algorithms, and formulas must be implemented in some computing language. It should be understood that although a variety of computing languages ​​can be used, this example chooses the R statistical language and relies on an R package called GenomicWordFramework. The steps and results obtained in this example come from applying the heuristics to a specific (and small) experimental data set. Results for larger data sets (in humans) are presented in a later example. library(GenomicWordFramework)

[0235] library(MethylIT.utils)

[0236] library(ggplot2)

[0237] library(ggseqlogo)

[0238] library(rtracklayer)

[0239] The goal of this package is to derive objects that can be used for further applications of DSP and ML tools available in other R packages. However, some applications of basic DSP tools are also provided.

[0240] Signal analysis of Arabidopsis experimental dataset

[0241] Examples of empirical methylation signal data are illustrated using the dataset included in the package. The experimental dataset carries the methylation levels of the Arabidopsis Columbia-0 ecotype (Col-0) and the msh1 mutant (dwarf phenotype). The methylation data was derived from a previous methylation analysis application and is included as a dataset in the package. The DMP dataset can be loaded from the package:

[0242]

[0243]

[0244] Retrieving genomic signals from DNA sequences

[0245] Here, the coordinates of the Arabidopsis annotated genes are given in the GTF file downloaded from ftp: / / ftp.ensemblgenomes.org / pub / plants / release-49 / gtf / arabidopsis_thaliana / .

[0246]

[0247] Use the coordinates of the gene to obtain the corresponding DNA sequence carrying the signal. For simplicity, the DNA sequence is provided in the package.

[0248]

[0249] A sequence for this instance is provided with the package.

[0250]

[0251] Next, the binary signal on the gene region was retrieved using the function getSignalAtRegions. Methylation signals were provided for three genes in five different Col-0 samples and five msh1 mutant (dwarf) samples. Only DMPs were included.

[0252] In different samples, methylation can be located at different positions, so there is a set complete.region = TRUE to retrieve information from the entire region, and zero_signal = TRUE to allow positions with no signal.

[0253]

[0254] This dataset is included in the following R package:

[0255] data(at_signal)

[0256] Example Results

[0257] Additional results supporting our heuristic for identifying methylation motifs in the Arabidopsis methylome are provided below:

[0258] The msh1 epi lines contain distinct non-genetic states based on phenotype. Manipulation of the msh1 mutant resulted in four distinct phenotypes, where states 1 and 2 are characterized by slow growth, delayed flowering, and a persistent stress response, while states 3 and 4 produce enhanced growth vigor and higher seed set (seeded) than the wild type (WT) ( Figure 1 ). State 4 was obtained by making Col-0msh1 mutant (state 1) or memory strain (state 2) x Col-0 WT orthogonal or reciprocal crosses and generating epi-F2 and epi-F3 families in Arabidopsis. Similar results were obtained regardless of the direction of hybridization. Compared with state 3 grafted progeny, the progeny populations in the F2 generation showed a more variable distribution of growth enhancement phenotypes, and each epi strain population could be classified as enhanced vegetative growth, higher seed setting rate, or both (Figure 2).

[0259] Four F3 populations followed hybrid-derived epi lines. Epi 8 and Epi 24 are sibling lines from a single WT x msh1 hybrid event, while Epi 10 and Epi 19 are sibling lines from a secondary WT x msh1 hybrid event. All four F3 epi lines showed consistent phenotypes within each population, but there were significant differences between the four populations (Figure 2). Epi line populations showed increases in aboveground vegetative growth and belowground root development ( Figure 2A , 2E ). Three populations, Epi 8, 10 and 19, had significantly higher total leaf areas than the WT control ( Figure 2B , 2F ), and all four populations had higher dry leaf weights than WT ( Figure 2G The four groups also showed higher dry root weights than WT ( Figure 2H ).

[0260] Enhanced reproductive growth was observed in three epi lines, Epi 8, 10, and 19, with Epi 8 showing early flowering ( Figure 2C). Importantly, full-sib populations Epi 8 and Epi 24 from the same primary cross displayed different phenotypes. Among the epi lines, Epi 24 had the highest seed weight and the lowest leaf area; whereas Epi 8 had the highest leaf area and the lowest seed weight. These observations indicate that within closely related populations, growth vigor may vary due to vegetative versus reproductive allocation. A similar phenomenon may occur in tomato epi lines (Figure 3), where populations of sibling epi lines also differ in terms of vegetative growth enhancement ( Figure 3H ).

[0261] By the fifth or sixth generation (S5, S6), the epi strain phenotype reverted to the wild type. In these successive generations, occasional reversions to a phenotype similar to the memory (state 2) were observed ( Figure 4A , 4B ), characterized by reduced growth rate, altered leaf morphology, and increased stress response. This reversion has not been observed in grafted progeny. The frequency of reversion varies between populations and species (Table 1).

[0262] Table 1. Reversion frequency within epi strains

[0263]

[0264]

[0265] *n, i, and e refer to mild (normal), intermediate, and extreme phenotypes of the msh1 memory strain used in crosses to generate the epi strain population. (TDNA) refers to the msh1 TDNA insertion mutant used as a control in crosses.

[0266] Plant features of the four msh1-derived states in Arabidopsis were recapitulated in tomato (Fig. 3) and soybean lines used in the study. Cross-species similarities included epi-line reversion (Fig. 4; Table 1). What is now considered epi-line reversion was previously reported as instability in msh1-altered sorghum. The observation of reversion across species implies that the epi-line state 4 condition is closely related to the memory state 2 condition. Reversion occurred only in epi-line populations derived from crosses with an msh1 memory line (state 2) parent, but not when an msh1 null mutant (state 1) was the parent (Table 1). One explanation for this observation is that msh1 signaling is weaker in the memory line than in the null mutant, resulting in less stable epi-line results.

[0267] By genome-wide methylation analysis, msh1 states 1 to 4 comprise discrete epigenetic stages. Significant changes in DNA methylation were detected in four Arabidopsis epi lines (F3), with gene-related changes primarily occurring in the CG environment ( Figure 5 ). Similar DNA methylation changes occurred in Arabidopsis, tomato, and soybean epi lines, showing variable hypermethylation and hypomethylation in CHH environments within and between epi lines of the three species (Figure 6). In a comparison of soybean Epi-F4 and Epi-F6, the growth vigor of the Epi-F4 line was enhanced, but declined to the wild-type phenotype by the F6 generation, allowing the identification of potential growth vigor-related DNA methylation changes ( Figure 6C ).

[0268] To estimate the association between differential methylation of genes and changes in plant phenotype, high-resolution methylome analysis was used. This procedure combines signal detection and machine learning to distinguish high-probability methylation changes associated with treatments within gene regions. Hierarchical clustering of methylome data from individual plants of all four states of Arabidopsis used changes in methylation levels at differentially methylated positions (DMPs) in gene regions (calculated as Hellinger divergence). The results showed clustering of individual plants (biological replicates) from the same population, with plants from different states divided into six branches that were different from the WT control group (Figure 7). Although derived from the same cross, Epi 8 and Epi 24 were divided into two subgroups ( Fig. 7A ), which seems to be consistent with their different phenotypes (Figure 2). The data show that DNA methylation remodeling of Epi 24 gene is more closely related to grafted offspring (HEG; state 3) than to its full sibling Epi 8 ( Fig. 7A Epi 10 and Epi 19 also formed two closely related clusters consistent with their phenotypic similarities. In contrast, hierarchical clustering of transposable element-associated DMPs in the four epi lines yielded results consistent with the lineages, such that Epi 8 and Epi 24 now clustered together, as did Epi 10 and Epi 19 in a single subcluster, and both clusters separated from the grafted offspring of state 3 (HEG) ( Figure 7C ).

[0269] Similar results were found in the hierarchical clustering analysis of tomato gene methylation data. The tomato msh1 mutant (state 1), HEG (state 3), and epi line data sets (state 4) formed three branches, consistent with the unique gene methylation effects (Figure 4b). Two epi lines with different plant growth vigor phenotypes, Epi High (H; biomass increased compared to WT) and EpiLow (L; biomass equivalent to WT) were divided into two subgroups within the same main cluster ( Figure 7B ). These observations suggest that gene-related methylation remodeling can be discerned to reflect plant phenotypes in lines predicted to be genetically uniform.

[0270] Methylation data support that epi line revertants are more closely associated with msh1 memory state 2. Arabidopsis epi line (state 4) revertant samples were included in the methylome analysis to test for evidence of methylation remodeling in revertant and non-revertant full-sib samples. Principal component and linear discriminant analysis (PCA-LDA) of these datasets using Hellinger divergence produced distinct clusters of non-revertants and revertants from the same progeny population. Epigenomic features that distinguish revertant and non-revertant full-sib progeny presumably arise spontaneously. The revertant methylome dataset showed a closer relationship with the data for msh1 memory state 2 ( Figure 4C ). These msh1 memory plants were derived from an independent population. The relationship between the methylome signatures of the revertants and msh1 memory is consistent with the similarities they exhibit in plant phenotypes (Figure 4). These results support previous indications that msh1 memory plays a role in methylome signatures. 18 The results showed that epigenetic state 4 (enhanced growth activity) was unique in the control group and confirmed the low-frequency conversion of epigenetic state 4 (enhanced growth activity) to state 2 (a persistent stress response phenotype).

[0271] Epi strain methylome datasets were altered in growth-related gene networks. Differentially methylated genes (DMGs) were identified for each Epi strain population by applying generalized linear regression analysis (GLM) to test the significance of differences between group DMP counts (WT and Epi strains) in genes. This analysis identified 3204, 2860, 3208, and 2797 DMGs in Epi 8, Epi 10, Epi 19, and Epi 24, respectively. To investigate the consistency of DMG functional relationships in each population, a gene network-based enrichment analysis was performed. DMGs from different epi strains shared enrichment of specific functional networks. For example, full siblings Epi 8 and Epi 24 shared common enrichment networks for responses to auxin, responses to red or far-red light, responses to nutrient levels, photoperiodism, detection of abiotic stimuli, responses to stress, and catabolic processes ( Fig. 8A ), while Epi 10 and Epi 19 shared enriched networks for responses to cyclopentenone, responses to auxin, cellular responses to auxin stimulation, and auxin-activated signaling pathways ( Figure 8B ). This result reflects the non-stochastic quality of the DMP dataset, meaning that loci targeted by the methylation machinery and their respective networks can be identified.

[0272] Methylome data from previously reported Arabidopsis F1 heterotic crosses of ecotypes C24 and Ler were analyzed using the methylome analysis method applied to the msh1 dataset. The enriched networks emerging from the F1 hybrid (C24 x Ler and Ler x C24) data showed significant agreement with those identified in the msh1-derived epi-line data, focusing on the following pathways: response to red or far-red light, response to auxin, growth regulation, cellular response to auxin stimulation, and auxin-activated signaling pathways ( Figure 8C ). DMG overlap between the hybrid data and epi lines ranged from 47% overlap with Epi 10 to 54% overlap with Epi 8. This clear alignment of methylome-identified gene pathways between heterotic hybrids and msh1 epi lines supports the connection of phenotype to methylome remodeling and links enhanced vigor to changes in auxin signaling. These observations are also consistent with the profile of DMG-derived enriched networks in tomato and soybean epi line datasets, which identified prominent auxin-related networks ( Figures 8D to 8E ). Soybean research focuses on phenotypic changes, with Epi-F4 lines showing significantly enhanced growth vigor, which then weakened back to wild-type growth in the Epi-F6 generation. Therefore, methylome signatures were analyzed to distinguish Epi-F4 from Epi-F6 datasets as a direct comparison to distinguish methylome signatures associated with vigor within a single lineage.

[0273] Examination of gene expression changes in Arabidopsis epi line populations involved sampling three different tissues for RNAseq analysis. Differential gene expression in leaf tissues of Epi 8 and Epi 10 was compared with wild type by similar sampling as performed in the methylome study, and Epi 8 and Epi 24 were additionally compared by analyzing scape and root tissues. From leaf tissue, 1884 differentially expressed genes (DEGs) were identified in Epi 8 and 992 in Epi 10 relative to wild type. Epi 8 and Epi 24 analysis revealed 1991 DEGs in Epi 8 and 1650 in Epi 24 relative to wild type from scape, and 1133 DEGs in Epi 8 and 1111 in Epi 24 from roots.

[0274] Network enrichment analysis of the derived DEG datasets revealed shared pathways altered by msh1 effects in leaf, flower stem, and root tissues. The largest enrichments involved abiotic and biotic stress responses, with circadian rhythm and phytohormone response-related networks also significantly enriched in these three tissues. Transcriptional regulation was particularly prominent in flower stems where MSH1 accumulated (Figure 9). These results appear to be consistent with the epi line DMG results (Figure 8), making the integration of the DEG dataset with the DMG dataset help amplify the enrichment of networks associated with msh1 relative to either dataset alone.

[0275] Differential methylation and expression analysis identified central gene hubs in the epi strain state. To investigate the interaction between the DMG dataset and the DEG dataset in the epi strain, the DMG and DEG data were input into Cytoscape to construct a protein-protein interaction (PPI) graph, followed by K-means clustering analysis to identify putative core networks carrying central gene hubs. The K-means clustering machine learning algorithm uses betweenness centrality, closeness centrality, average shortest path length, clustering coefficient, degree, and eccentricity as parameters, allowing the identification of clusters containing the most concentrated nodes (proteins) in the PPI network.

[0276] In Arabidopsis Epi 8, a total of 3647 unique loci from DMGs and DEGs were used for analysis to generate a PPI network formed by 430 genes. Functional enrichment analysis of these putative hub genes was performed using the STRING database functional enrichment tool, revealing a PPI network of 153 hub genes and associated functional networks ( Fig. 10A ). These core networks included developmental processes, responses to hormones (particularly auxin), responses to cold, chromosome organization, mRNA processing, spliceosome, and ribosome biogenesis. Compared to Epi 24, a larger PPI network of 346 genes was formed using a total of 3523 unique loci from DMGs and DEGs, indicating that some of these same processes were represented ( Fig. 10B ). However, the Epi 24 analysis also revealed more prominent changes in the expression of development-related genes and, most notably, strong and overlapping effects in the expression of genes involved in ribosome biogenesis ( Fig. 10B ). Ribosome biogenesis is integral to the growth response. Since Epi 8 and Epi 24 are full siblings from the same cross, these different characteristics are expected to be related to their different phenotypes, with Epi 8 being enhanced in vegetative growth and Epi 24 showing significantly increased reproductive vigor (Figure 2).

[0277] The analysis of soybean epi line data included assigning Arabidopsis orthologs to each identified soybean gene using BLASTP and obtaining the soybean epi line PPI network of 109 core hub genes and their functional networks ( Fig.11A The resulting enriched networks included developmental processes, ribosomes, chromatin organization, mRNA processing, and response to light stimulation. These results suggest that there is cross-species overlap in the identified pathways ( Fig.11A ), so that the effects of msh1 lead to parallel phenotypic results in different species through similar networks. The annotated gene overlap rate for core hub genes between soybean and Arabidopsis was 16% (18 genes).

[0278] To assess the relationship between the PPI network output of the epi lines and the grafted HEGs, a core hub PPI network was constructed for the tomato HEG (state 3) phenotype. This network contains ribosome biogenesis, developmental processes, and chromatin organization ( Fig. 11B ). Core hub proteins shared by epi lines and HEG datasets of soybean, Arabidopsis, and tomato have been identified, including TOR, MOM, RFC1, UBQ family, NRPD family, and HAC family proteins. Although several overlapping hubs are conserved across species, detailed cross-species comparisons are limited by differences in gene annotations.

[0279] Comparison of Epi states in Arabidopsis reveals conserved MSH1 epigenomic targets in biologically significant gene networks. To investigate the relationship between gene methylation remodeling between the four distinct MSH1-derived states, DMG overlap was assessed. Fig. 12A 871 DMGs were shown to be shared between the four msh1 states, including 17.6% of msh1 (state 1), 39.7% of memory (state 2), 31.6% of grafted progeny (HEG, state 3), and 31.1% of epi strain (state 4, Epi 24) DMGs. Overlapping established a conserved msh1 "core" DMG dataset ( Fig. 12A Each state was also classified by state-specific DMGs, accounting for 33.8% in state 1, 15.4% in state 2, 10.2% in state 3, and 15.2% in state 4 ( Fig. 12A Of the 871 msh1-related “core” DMGs, 531 (61%) could be placed in known networks.

[0280] Components of the RdDM pathway have been shown to be required for the induction of msh1 state 1, the transition from state 1 to state 2, and the generation of state 3 after grafting. The methylome dataset compares msh1 with the dcl2 / dcl3 / dcl4 / msh1 quadruple mutant (state 1) or Col-0 / Col-0msh1 Col-0 / Col-0 dcl1 / dcl3 / dcl4 / msh1 The grafted progeny of the graft (state 3) helped to classify DMGs as RdDM (sRNA) dependent. These subtractive datasets confirmed that 674 (77%) of the 871 core DMG loci were predicted to be DCL2, DCL3, DCL4 dependent by obtaining overlap between the 871 DMG core dataset and the msh1 and dcl2 / dcl3 / dcl4 / msh1 datasets. Fig. 12B The PPI core hub network of DMGs is shown, with those predicted to be DCL2, DCL3, and DCL4-dependent with cross-hatching. These results provide evidence for a relationship between the epigenetic effects directed by msh1 and the effects of the RdDM pathway on target methylation remodeling.

[0281] Although 871 core DMGs are shared between the four msh1 states, methylation changes found within these 871 loci also help distinguish the four states. Fig. 12C Results of hierarchical clustering of individual plants from four epigenetic states using only DMP methylation information on 871 msh1 core genes, which separated into four distinct clades, are shown. Based on the PPI network of 871 DMGs, K-means clustering was performed to identify putative central hubs and their functional enrichment analysis. Fig.12D The resulting network of 67 candidate hubs is shown, which are involved in responses to stress, developmental (growth) processes, gene expression, spliceosomes, histone modifications, and chromosome organization networks. More than 80% of the 67 hubs are predicted to be RdDM targets, several of which play a central role in plant development and environmental stress responses, such as the cell growth regulator rapamycin target TOR, the SWI2 / SNF2-like meristem recognition gene SPLAYED (SYD), its close homolog BRAHMA (BRM), and the chromatin remodeling factors PICKLE (PKL) and PICKLE-RELATED 1 (PKR1 / CHR4).

[0282] Transposable elements and sRNAs are associated with candidate RdDM target loci in 871 msh1 DMGs. The RdDM pathway is known to actively target transposable element (TE) sequences, which has motivated the investigation of TE associations with msh1-responsive loci. Looking at the 871 DMGs shared between msh1 states, we detected associations of these loci with TE and sRNA (20-24nt) clusters above the genome-wide level (61% DMGs within 2kb TEs vs. 47% DMGs genome-wide; 77% DMGs within 2kb sRNA clusters vs. 46% DMGs genome-wide). This enrichment for TE and sRNA cluster proximity was further increased when dcl2, dcl3, dcl4-sensitive DMGs (65% DMGs within 2kb TEs; 81% DMGs within 2kb sRNA clusters) were subsets. However, when focusing on the 67 hub DMGs derived by k-means clustering, only sRNA cluster enrichment was seen (49% of DMGs within 2 kb TEs; 72% of DMGs within 2 kb sRNA clusters). Fig.12D Genes associated with only TEs (7), genes associated with only sRNAs (22), and genes associated with both TEs and sRNAs (26) were shown. Thus, while TEs may affect the methylation status of neighboring genes and serve as RdDM targets, sRNA cluster associations can define RdDM targets for msh1 effects regardless of TE proximity.

[0283] Some possible evidence for an association between TE families and DMG proximity was found. Comparison of the observed and expected numbers of TE family members revealed a significant overrepresentation of members of the L1 and Gypsy families, and an underrepresentation of members of the Helitron family, in 871 DMGs shared between msh1 states and in 674 DMGs sensitive to dcl2, dcl3, and dcl4. Further studies are needed to reveal any biological significance of these associations. Based on different analyses, four criteria facilitated the classification of RdDM target loci in our study ( Fig.12D ): Positive loci were (1) confirmed DMGs in all four msh1 states, (2) contained putative network hubs based on k-means clustering, had genetic evidence for (3) dcl2, dcl3, dcl4 sensitivity, and / or (4) association with sRNA clusters.

[0284] Detailed methylation analysis of selected RdDM target genes in four distinct non-genetic states revealed sequence motifs containing dcl2 / 3 / 4-sensitive DMPs. Annotation of 67 candidate hub loci supported their relevance to the phenotypic effects observed in the four msh1 states. In addition to gene networks of altered gene expression and chromatin behavior, major overlapping networks appeared to reflect the observed transitions between stress response and growth ( Fig.12D ), which supports the relationship between msh1 and memory stress phenotypes and growth vigor in HEG and epi strain states. Methylation signatures of 67 candidate hub loci that distinguish four msh1 states support a model for precise targeted remodeling ( Fig. 12C ).

[0285] Figure 13 and Tables 2 and 3 show the results of the analysis of seven selected candidate RdDM target loci: target of rapamycin (TOR), a cell growth regulator; SPLAYED (SYD), a chromatin remodeling factor component; ubiquitin-specific protease 26 (UBP26), required for heterochromatin silencing; nuclear RNA polymerase D1A (NRPD1), the largest subunit of RNA polymerase IV; RNA polymerase II large subunit involved in transcription (NRPB1); SU(VAR)3-9-related protein 5 (SUVR5), a gene involved in H3K9me2 deposition and gene silencing; and UP-FRAMESHIFT1 (UPF1), a gene involved in RNA interference and splicing.

[0286] Table 2. Summary of gene methylation data within motif cluster 11*

[0287]

[0288] *+ indicates the presence of DMP in at least two samples, - indicates the absence of DMP, ! indicates dcl2, dcl3, dcl4 sensitivity. **The numbers in brackets refer to the number of epi strains that displayed DMP in at least 2 samples.

[0289] Table 3. Summary of gene methylation data within motif cluster 18*

[0290] Gene ID TOR NRPD1 NRPD1 SUVR5 SYD SYD UBP26 NRPB1 UPF1 UPF1 State 1 +！ - - +！ - - - +！ +！ +！ State 2 - - - + - - - + - - State 3 - - - - + - +！ - +！ +！ State 4 - +(4) +(1) +(4) +(1) +(1) +(3) - - -

[0291] *+ indicates the presence of DMP in at least two samples, - indicates the absence of DMP, ! indicates dcl2, dcl3, dcl4 sensitivity. **The numbers in brackets refer to the number of epi strains that displayed DMP in at least 2 samples.

[0292] Differential methylation analysis within seven loci revealed evidence of state-specific remodeling (Tables 2 and 3). The methylation changes at each locus were associated with a recognizable sequence motif spanning approximately 14 nucleotides. Fig.13ATwo sample composite motifs are shown in Fig. 13B and 13C Examples of state-specific methylation remodeling of UBP26 and UPF1 within cluster 11 are shown in . Systematic analysis of 67 loci revealed similar state-specific changes. Although multiple DMPs could be detected within a single identified motif domain, DMPs were analyzed separately. DMPs were distinguished in the replicated datasets, with several loci shown to be dcl2 / dcl3 / dcl4-sensitive in both state 1 and state 3 ( Fig. 13C ; Tables 2 and 3). It was feasible to test the dcl2 / dcl3 / dcl4 sensitivity only in the two states containing the msh1 rootstock. Tables 2 and 3 summarize the scoring results of the cluster 11 and cluster 18 motifs in the seven selected genes.

[0293] The identified DMPs showed no obvious exonic, intronic or junctional localization patterns and contained multiple motif sites per gene. Evaluation of clustered motifs containing DMPs within 67 msh1 core hub loci revealed evidence of higher-order dependencies in many cases. Multiple sequence alignment (MSA) of a given DNA motif can reveal dependencies between two nucleotides located at different positions within the motif, reflected in the frequency of their co-occurrence. First-order dependencies refer to adjacent nucleotides and are often found in the context of CG methylation, second-order dependencies refer to nucleotides separated by two nucleotides, and higher-order dependencies refer to nucleotides separated by more than two nucleotides. These relationships are derived from the study of Markov dependencies in DNA sequences, which is the basis for the application of hidden Markov models for motif discovery.

[0294] For the identified motifs, individual consensus nucleotides are evident at different distances from the target cytosine, which is nucleotide 7 within each motif (on either the plus or minus strand). For example, the motif from cluster 65 shows an invariant T at position 14 and a consensus A at position 12, while the motif from cluster 66 shows an invariant G at position 14 and a consensus AG pair at positions 2 and 3, respectively ( Fig.13D ). In motif 76, only Gs are located at positions 6 and 9, with a consensus T at position 5, while motifs 82 and 86 show a consensus or invariant T at position 2. These observations support the hypothesis that the target methylation sites within these genes are characterized by invariant or nearly invariant sequence features.

[0295] Small regions (14 bp) containing DMPs in at least three samples were identified and considered as candidates for DNA methylation motifs among the 67 identified msh1 core hub genes. The distance matrix was estimated for the selected set of regions using the function dist.dna in the ape R package (version 5.5). Hierarchical clustering analysis of the selected set of regions (using the previously estimated distance matrix) was completed using the function hclust in the stats R package (version 4.1.1) and grouped into 100 clusters. The UPGMA method was applied as the agglomerative algorithm. Clusters with less than 10 regions were discarded. DNA multiple sequence alignment on each sequence cluster was performed using the MUSCLE algorithm implemented on the Bioconductor R package muscle (version 3.14). Fig.13A The motifs presented in are the result of a specific analysis of the genes TOR, SYD, NRPB1, NRPD1, UBP26, SUVR5, and UPF1.

[0296] Example 2: Detection of signal peaks

[0297] A rapid way to discover DNA sequence features associated with methylation is to determine the locations of methylation sites in a sample. This approach relies on evidence that methylation occurs at specific, non-random sites, a DNA sequence motif.

[0298]

[0299] These data are included in the following packages:

[0300]

[0301] The function signal_peaks provides a potential DNA sequence motif with coordinates centered on the methylation peak and covering 3-up and 4-down bp around the signal. It requires the presence of methylation signals in at least three samples (cutoff = 3L).

[0302] Clustering of potential DNA sequence motifs

[0303] Applying DNA multiple sequence alignment MUSCLE and hierarchical clustering to identify DNA sequence motif clusters ( Figures 14A to 14D ).

[0304]

[0305] The clustered motifs containing DMPs on the gene AT1G50030 (TOR) revealed evidence for higher-order dependencies. That is, in a multiple sequence alignment (MSA) of a given DNA motif, base Y at a given site k depends on base X at the previous site j if bases X and Y are both observed with high frequencies in the MSA. In particular, if k–j = 1, k is a first-order dependency commonly found in the context of CG methylation; if k–j = 2, k is a second-order dependency; and when k–j > 2, k is a higher-order dependency.

[0306] For these motifs, individual consensus nucleotides are evident at different distances from the target cytosine, which is nucleotide 7 within each motif (on either the '+' strand or the '-' strand). For example, the motif from cluster 1 ( Fig.14A ) shows only T at position 1 (maximum score is 2 positions), a consensus G at position 9, T at positions 12 and 15 (about 1.5 positions), while motif 17 ( Fig.14D ) shows only T at position 1. In motif 12 ( Fig. 14C ), a clear consensus of A was found at positions 5 and 10, and at the level of methylated C at position 7.

[0307] Get DNA sequence motifs by gene:

[0308]

[0309] Applying Digital Signal Processing (DSP)

[0310] DSP can be applied to the previously obtained signal. Next, the power spectrum of the chromosome 2 and 3 signals is calculated.

[0311] Functions can be applied that encode previously detected 0 and 1 binary signals from the DNA sequence into a user-defined numerical code. Possible encodings can be binary, real and complex numbers. GWF also supports encoding of DNA methylation sequences using complex numbers. Encoding using ordinary real numbers is also supported. The basic idea is to encode the physicochemical properties of the DNA bases. In this case, by applying different DSP tools, periodicity and correlation of the encoded signal for the superposition of the methylation signal with the physicochemical signal can be searched. Currently, DSP analysis of complex signals with R does not work well. However, the methylation signal can be encoded with GWF and then exported to, for example, MatLab or Pythom and the DSP analysis can be done there.

[0312] Given two objects, one carrying the signal and the other carrying the DNA sequence, this function will perform an encoding set by the user. This function can be used to re-encode previously detected 0 and 1 binary signals from a DNA sequence into a user-defined numeric code. In particular, information about the physicochemical properties of neighboring DNA bases can be incorporated: the number of hydrogen bonds and the base chemical type. This can be the default value used by the function:

[0313] 1. Methylation status: 1 (methylated) and 0 (unmethylated)

[0314] 2. Hydrogen bonds: 1 (three hydrogen bonds) and 0 (two hydrogen bonds)

[0315] 3. Chemical type: 1 (pyrimidine) and 0 (pyrimidine)

[0316] Therefore, each base generates a binary string:

[0317] 1. Signal mC, methylated pyrimidine with three hydrogen bonds: 111

[0318] 2.C, pyrimidine with three hydrogen bonds: 011

[0319] 3.G, purine with three hydrogen bonds: 010

[0320] 4. T, pyrimidine with two hydrogen bonds: 001

[0321] 5.A, purine with two hydrogen bonds: 000

[0322]

[0323]

[0324] Next, divide the signal into non-overlapping window intervals of 90 bits (30 bp):

[0325]

[0326]

[0327] There are 592 potential DNA sequence motifs:

[0328]

[0329] Overlap can be searched between the motifs with a 30-bp wide sliding window.

[0330]

[0331] That is, 20 of the 592 genomic word frames were completely located within a 30-bp region on gene AT1G50030. However, in the dwarf sample 'dw2', this region was covered by only 3 significant motifs.

[0332]

[0333]

[0334] These three motifs are embedded at the beginning of the observed signal region.

[0335]

[0336] Sequences from clusters 12 and 9 ( Fig. 14B , 14C ) differ only in the base shift (position 3).

[0337]

[0338] Power spectrum analysis of experimental data sets

[0339] The power spectrum of a binary signal can be obtained from a SignalMatrix-class object using the function plot_power_spectral (Figure 15). The dots help to visually identify the control (wt) ( Fig.15A ) and treatment (msh1 mutant, dwarf phenotype) ( Fig. 15B ) have peaks with spectral differences between them.

[0340]

[0341] The peak at 1 / 3 (0.33) indicates that the observed region is a protein coding region ( Fig.15A ). Methylation may or may not affect the 1 / 3 periodicity. In protein coding regions, peaks other than the peak at 1 / 3 are generated in the present case by periodicity derived from the physicochemical properties of the DNA bases encoded in the binary vector representing the base sequence.

[0342] Spectrum

[0343] A spectrogram is a visual representation of the spectrum of the signal's frequencies varying over time. In our case, the axis "time" will be represented by bits or individual DNA base positions. A spectrogram can be obtained from any region of a coding signal using the function "spectrogram" in the R package phonTools.

[0344] To prepare the methylation signal as a digital signal:

[0345]

[0346] Next, the values ​​from the control ( Fig.16A , above) and processing ( Fig.16A , below) of the same gene:

[0347]

[0348] The methylation signal is at about 1 / 3 frequency ( Fig.16A The dotted line in the figure) and the power spectrum energy around the decomposed periodicity at box 16. In the third region we find:

[0349]

[0350] In this case, the methylation effect is located around bins 6-10, and the energy power is shifted down by 1 / 3 frequency ( Fig. 16B ). In essence, the effect of methylation on signal power and periodicity depends on the DNA sequence context. This is an expected result that has been ignored by traditional methylation analysis so far.

[0351] Wavelet Power Spectrum

[0352] The wavelet coefficients yield information about the correlation between the wavelet (at a specific scale) and the data array (at a specific location). Larger positive amplitudes mean higher positive correlations, while larger negative amplitudes mean higher negative correlations.

[0353] Wavelet power spectra provide a useful way to determine the distribution of energy within a data array. By looking for regions within the high-power wavelet power spectrum (WPS) in the plot, we can determine which features of the signal are important and which can be ignored. Here, the term "energy" is not arbitrary, but is borrowed from its application in artificially constructed communication systems. The energy level represented in the WPS is proportional to the energy consumed to transmit a binary signal of the same size to a given receiver through an artificially constructed communication grid.

[0354] Furthermore, with the advancement of information theory and its application to biomolecular methods, it is well known that in order to complete a single methylation change, each methyltransferase / demethylase must consume a minimum amount of energy to process the information associated with the change. This energy is determined by Landauer's principle, according to which a molecular machine must consume a minimum amount of energy ε = k at each step of a genetic logic operation (including proofreading). B T ln2 (about 3x 10-21 joules per bit at room temperature).

[0355] Wavelet power spectrum analysis of experimental data sets

[0356] The wavelet power spectrum analysis of the previously estimated at_signal_diff dataset (Figure 17) is done in the following way:

[0357]

[0358] Some methylation motifs have the same methylation status in the same region in both Col-0 and msh1 dwarf groups. This is the case for the following region: AT1G50030.3 ( Fig. 17C ) and AT1G50030.6. The changes observed in the wavelet power (WPS) spectrum correspond to the bases where the signal is located.

[0359] In AT1G50030.1( Fig.17A ), WPS indicates that the change occurred at position 60, which means that the base at position 20 = 60 / 3. The table below confirms our estimate. The methylated cytosine is located at base 20 of the Col-0 negative strand. The reported base in the positive strand is G, which is encoded by the binary string 110 when the complementary base C is methylated. The sample base position in msh1Dwarf-2 (G:010) is not methylated.

[0360]

[0361] In the same way, the correlation graph can be calculated based on WPS ( Fig.18 ):

[0362]

[0363]

[0364] In addition to the expected disruption of the correlation in the region 55-60, the correlation between the energy spectra in the regions 30-35 and 55-60 is lost ( Fig.18 ).

[0365] Example 3: Human Dataset

[0366] Completed methylation signal analysis of 81 genes related to the development of the human nervous system.

[0367] Fig.19is a subnetwork of 81 DMG hubs. The PPI network was built on a collection of 751 DMGs identified by principal component analysis and analyzed by the STRING Cytoscape application. Circles without cross-hatching highlight DMGs involved in biological processes, including nervous system, nervous system development, synapses, neuronal projections, ion transport, CNS diseases, axon guidance, neurogenesis, ion / cation binding, ion / cation transport across membranes, voltage-gated channels, and TRP channels.

[0368] The entire DMG set was derived from 751 DMGs, which were selected based on their contribution to the classification of patients into two groups, namely "typical" and "autistic". Specifically, these 751 DMGs contributed more than 1% of the total variance to the main principal components of PCA. These DMGs were analyzed using the STRING Cytoscape application, and the K-means clustering method was applied to identify the main sub-networks of hubs using the network centrality index as a variable. All DMGs were obtained through methylation analysis.

[0369] Fig. 20 Is from Fig.19 Network enrichment analysis identified in the main network of the DMG hub of the Autism-related gene network. The identification of methylation motifs in the 81 DMG members of the main subnetwork derived from the hub of the Autism-related gene network can be achieved by the methods described herein.

[0370] From the foregoing it can be seen that the present disclosure achieves at least all stated objectives.

[0371] Glossary

[0372] Unless otherwise defined, all technical and scientific terms used above have the same meaning as commonly understood by one of ordinary skill in the art to which embodiments of the present disclosure belong.

[0373] The terms "a," "an," and "the" include both singular and plural referents.

[0374] The term "or" is synonymous with "and / or" and refers to any one member or combination of members of a particular list.

[0375] The terms "invention," "present invention," "disclosure" or "present disclosure" are not intended to refer to any single embodiment of the particular invention, but rather encompass all possible embodiments described in the specification and claims.

[0376] As used herein, the term "about" refers to slight variations in the numerical amount of any quantifiable variable. Unintentional errors may occur, for example, due to use of typical measuring techniques or equipment, or due to differences in the manufacture, origin, or purity of the components.

[0377] The term "substantially" refers to a large or significant degree. Thus, given the appropriate context, "substantially" may refer to a plurality, a majority, and / or a majority of the quantifiable variable.

[0378] The term "generally" encompasses "about" and "substantially".

[0379] The term "configuration" describes a structure that is capable of performing a task or adopting a specific configuration. The term "configuration" can be used interchangeably with other similar phrases such as build, arrange, adjust, manufacture, etc.

[0380] Terms characterizing order, position and / or orientation are not limiting and are referred to merely in terms of the views presented.

[0381] "Phenotype" refers to the set of observable characteristics of an individual resulting from the interaction of his or her genotype with the environment.

[0382] "Epigenetic" refers to non-genetic influences that arise from gene expression.

[0383] An “R package” is an extension to the R statistical programming language. R packages contain code, data, and documentation in a standardized, integrated format that R users can typically install from centralized software repositories such as CRAN.

[0384] "Expression" is the degree to which an individual with a particular genotype expresses a phenotype.

[0385] The "scope" of the present invention is defined by the appended claims and the full scope of equivalents authorized by such claims. The scope of the present invention is further identified as including any possible modifications to any aspect and / or embodiment disclosed herein, which modifications will result in other embodiments, combinations, sub-combinations, etc. that are obvious to those skilled in the art.

Claims

1. An extension program for a general programming language or a statistical programming language, the extension program comprising: An algorithm for analyzing methylation signals on a DNA sequence fragment characterized by: iii. Methylation information; and iv. The physicochemical information around each methylated cytosine; Wherein the algorithm comprises one or more functions, the one or more functions may: estimating a distance matrix for a set of selected regions of the DNA sequence; performing a hierarchical cluster analysis on the selected region set; grouping the set of selected regions into a specified number of clusters; as well as Multiple DNA sequences from the clusters are aligned into methylation motifs.

2. The extension program according to claim 1, wherein the extension program is written in R statistical language.

3. The extension of claim 2, further comprising a digital signal processing (DSP) tool available in another programming language.

4. The extension program according to claim 3, wherein the other programming language is C++, Python or MatLab.

5. The expansion program according to claim 3, further comprising: Encoding the physical and chemical properties of DNA sequences and DNA bases into numerical or complex signals further includes DSP. The expansion program according to claim 5 , wherein the encoding is based on a group structure.

7. The expansion program according to claim 6, wherein the group structure is an Abelian group.

8. The expansion program of claim 1, wherein the set of selected regions are grouped into clusters of at least 100 clusters.

9. The expansion program of claim 1, wherein clusters having less than ten regions are discarded.

10. The extension program of claim 1, wherein the plurality of DNA sequences are aligned using multiple sequence comparisons performed using a logarithmic expectation algorithm.

11. The expansion program according to claim 5, further comprising: Export the numerical or complex signal to C++, Python or MatLab.

12. The expansion program of claim 1, wherein the methylation motifs are further grouped for downstream analysis.

13. The extended program according to claim 12, wherein the methylation motifs are further grouped using a clustering algorithm. The expansion program according to claim 13 , wherein the clustering algorithm is a distance-based k-centers clustering algorithm.

15. The expansion program according to claim 1, further comprising: Encode the methylation and physicochemical signals of DNA bases into a binary, numerical or complex signal. The expansion program according to claim 15 , wherein the encoding is based on a group structure.

17. The expansion program according to claim 16, wherein the group structure is an Abelian group.

18. The expansion program according to claim 15, further comprising: Power spectrum analysis of the encoded methylation and physicochemical signals.

19. The expansion program according to claim 18, wherein the power spectrum analysis is wavelet power spectrum analysis (WPS).

20. The expansion program according to claim 19, further comprising: A correlogram is calculated based on the WPS.

21. A computerized heuristic method comprising: Higher-order DNA base interdependencies regarding methylated cytosine; as well as Statistically non-random distribution of bases.

22. The computerized heuristic method of claim 21, wherein the higher order DNA interdependencies and the base distributions are obtained from analysis of regions of DNA sequences by performing at least one hierarchical clustering and aligning a plurality of DNA sequences from the clusters into methylation motifs.

23. A method for analyzing methylation signals on a DNA sequence fragment, the method comprising: performing hierarchical cluster analysis on regions of the DNA sequence; Hierarchically group the selected set of regions into a specified number of clusters; aligning potential DNA sequence motifs from the clusters; as well as Apply digital signal processing to the encoded methylation and physicochemical signals.

24. The method of claim 23, wherein the set of selected regions are grouped into groups of at least 100 clusters.

25. The method of claim 23, wherein clusters having fewer than ten regions are discarded.

26. The method of claim 23, wherein the encoded methylation and physicochemical signals are encoded based on population structure.

27. The method of claim 26, wherein the group structure is an Abelian group.

28. The method of claim 23, wherein the potential DNA sequences are aligned using multiple sequence comparisons performed using a logarithmic expectation algorithm.

29. The method of claim 23, further comprising: The encoded methylation and physicochemical signals were subjected to power spectrum analysis.

30. The method of claim 29, wherein the DSP is applied by using the Genome Word Framework (GWF) R package.

31. The method of claim 30, wherein the GWFR package includes one or more clustering algorithms.

32. The method of claim 30, further comprising exporting the GWFR package for analysis with other DSP tools in different programming languages.

33. The method of claim 32, wherein the different programming language is C++, Python or MatLab.

34. The method of claim 23, further comprising: Methylation analysis is used to identify differentially methylated positions (DMPs) that constitute part of the methylation signal.

35. The method of claim 23, further comprising: Each of the multiple sequence alignments (MSAs) was evaluated for deviation from a random Monte Carlo simulated MSA.

36. The method of claim 23, further comprising: Analyze genes with epigenetic regulatory functions to diagnose or treat diseases.

37. The method of claim 36, wherein the disease is autism or cancer.

38. A computerized heuristic method for analyzing genetic data, comprising: (1) Information Diffusion (ID) statistics estimated for each gene carrying at least one differentially methylated position (DMP) on the gene or in the promoter region; (2) principal component analysis (PCA), in which the top k components carrying 1% or more of the entire sample variance are considered in downstream analyses; (3) calculating the correlation matrix of pairwise gene correlations carrying vectors represented as PCs; (4) Analyze the network’s correlation matrix; as well as (5) The contribution of each gene to phenotypic discrimination, which is assessed based on the fraction of the cumulative variance accounted for by the total sample variance carried by the gene.

39. The computerized heuristic of claim 38, wherein the ID is selected from the group consisting of: Hellinger divergence / distance, J divergence, and total variance distance.

40. The computerized heuristic of claim 39, wherein the ID is J-divergence.

41. The computerized heuristic of claim 38, wherein the PCA is applied with a function such that genes are represented as a k-dimensional vector of PCs, and further wherein the square of each coordinate carries a vector contribution in terms of variance in distinguishing treatment from control groups.

42. The computerized heuristic of claim 38, wherein the PCA is applied with a pcaLDA function.

43. The computerized heuristic method of claim 38, wherein the correlation matrix comprises a weighted correlation network (WCN), wherein the WCN and the network are analyzed, and the network is a protein-protein interaction (PPI) network.

44. The computerized heuristic method of claim 43, further comprising comparing results from the WCN and the PPI networks with identified consistent relationships and epigenetic gene contributions to the phenotype.

45. The computerized heuristic method of claim 38, further comprising computing a magnitude of the Euclidean norm of a gene represented as a vector of k PCs.

46. ​​The computerized heuristic method of claim 38, wherein the network involves one or more biological processes.

47. A method for reducing the amount of calculation, the method comprising: Use methylation analysis to identify differentially methylated genes (DMGs) on DNA sequence fragments; integrating the DMG into a gene network; as well as Identify network hubs.

48. The method of claim 47, wherein the network hubs are identified via protein-protein interaction (PPI) network analysis and weighted correlation network (WCN) analysis.

49. The method of claim 48, further comprising comparing results from the WCN and the PPI network with identified consistent relationships and epigenetic gene contributions to phenotype.

50. The method of claim 47, wherein the network hubs are involved in one or more biological processes.

Images