Methods, apparatuses, and media for predicting non-coding mutation function based on single-cell multi-omics

Abstract

The application discloses a method, device and medium for predicting non-coding mutation functions based on single-cell multi-omics, and belongs to the technical field of single-cell multi-omics. The application provides a unified deep learning framework, which can directly start from DNA sequences and jointly predict chromatin accessibility and gene expression levels at single-cell resolution. The application proposes two new quantitative scoring systems, regulatory effect scores and expression effect scores, which can respectively improve the accuracy and interpretability of non-coding variation function annotation based on the prediction effect of chromatin accessibility and gene expression. The application can be applied to the research of diseases such as Alzheimer's disease, and can analyze the pathogenic cascade process of transcription factor binding change, cis-regulatory remodeling and key gene-specific transcription disorder involved in candidate mutations.

Description

Technical Field

[0001] This application relates to the field of single-cell multi-omics technology, and in particular to methods, devices and media for predicting non-coding mutation functions based on single-cell multi-omics. Background Technology

[0002] Non-coding regions in the human genome are widely considered to play a central role in gene expression regulation and disease pathogenesis. Genome-wide association studies (GWAS) have revealed the crucial role of non-coding variants in regulating complex phenotypes. However, GWAS mutation sites are often limited by linkage disequilibrium interference, making it difficult to accurately distinguish between causal and related variants and to precisely identify key risk genes.

[0003] In human brain diseases, numerous non-coding genetic variations have been found to be associated with disease susceptibility. These variations may regulate cis-regulatory elements such as enhancers and promoters, as well as their associated long-range regulatory networks. Abnormal regulation caused by non-coding mutations can disrupt the function of specific brain cell types, thereby inducing a variety of neurological diseases. However, traditional batch-level omics methods struggle to resolve the highly heterogeneous cellular composition of the human brain and cannot effectively reveal the complex interactions between polygenic effects and genetic heterogeneity, thus limiting a deeper understanding of the regulatory mechanisms of individual mutations.

[0004] The rapid development of single-cell sequencing technology has made it feasible to establish higher-resolution blueprints of human cells, promoting the development of a new generation of models. Some transfer learning models, such as Huatuo and ExPectoSC, can infer gene expression at the cell type level based on pre-trained model parameters; other models, such as scBasset, NvwaCE, and Nvwa, directly predict chromatin accessibility or gene expression from DNA sequences. These tools provide some technical support for the analysis of functional non-coding variations. However, there are still two key limitations: (1) the lack of a unified framework that can simultaneously achieve joint prediction of single-cell chromatin accessibility and gene expression; and (2) the lack of the ability to identify brain tissue-specific regulatory logic, making it difficult for existing models to accurately decode the cell type-specific gene regulatory syntax rules in the human nervous system. Summary of the Invention

[0005] To solve at least one of the above-mentioned technical problems, the technical solution adopted in this application is as follows.

[0006] The first aspect of this application provides a method for constructing a deep learning model for predicting the function of DNA sequences based on single-cell multi-omics data, wherein the single-cell multi-omics data includes an epigenome dataset and a single-cell chromatin accessibility dataset, and the method includes the step of constructing a regulatory module for predicting the chromatin accessibility probability of sequence fragments: Pre-training is performed based on the epigenome dataset, and parameters are further optimized based on the single-cell chromatin accessibility dataset to obtain general epigenomic features of sequence fragments. Single-cell specific features are extracted based on the single-cell chromatin accessibility dataset. The general epigenomic features and the single-cell specific features are mapped.

[0007] In some embodiments of this application, the epigenetic dataset is derived from public epigenomic databases, including but not limited to the ENCODE database, the Roadmap Epigenomics database, and the Cistrome database. In some embodiments of this application, the epigenetic dataset contains 21,907 epigenomic features, including 9,471 transcription factor binding sites, 10,064 histone modification signal regions, and 2,372 chromatin accessibility peak regions.

[0008] In some embodiments of this application, the apparent dataset is prepared as follows: (1) Sequence Segmentation and Tag Generation: The human reference genome GRCh38 was divided into non-overlapping 200 bp sequence fragments. The inventors compared each sequence fragment with the peak regions of the aforementioned public epigenome dataset. For any epigenomic feature, if the overlap between the sequence fragment and the peak region corresponding to that feature exceeded 50%, the tag position corresponding to that feature was marked as 1; otherwise, it was marked as 0. Thus, each sequence fragment was assigned a 21907-dimensional binary tag vector, where each dimension corresponds to a specific epigenetic regulatory feature.

[0009] (2) Sequence Expansion and Encoding: Using the midpoint of each sequence fragment as the reference coordinate, the sequence is expanded 1000 bp upstream and downstream to obtain a DNA sequence fragment of length 2000 bp. Subsequently, the inventors used one-hot encoding to digitally represent each base in the sequence, where the bases "A", "C", "G", and "T" are encoded as [1,0,0,0], [0,1,0,0], [0,0,1,0], and [0,0,0,1], respectively. This method can transform the DNA sequence into a multidimensional tensor suitable for deep learning input.

[0010] (3) Dataset partitioning: To ensure the independence of model evaluation and generalization performance, a chromosome-level data partitioning strategy is adopted. Specifically, sequence fragments on chromosomes 7 and 8 are used as the test set, fragments on chromosome 9 are used as the validation set, and the remaining chromosome fragments are used as the training set. This partitioning method can avoid the information leakage problem caused by overlapping sequence positions and ensure the independence between model training and test data.

[0011] In some embodiments of this application, the single-cell chromatin accessibility dataset is derived from public databases, including but not limited to the GEO database and the ATAC database.

[0012] In some embodiments of this application, the preprocessing steps for the single-cell chromatin accessibility dataset are as follows: (1) Sequence extraction and encoding: Taking the center point of each chromatin accessibility peak as a reference, extend 1000 bp upstream and downstream to generate a 2000 bp genome sequence fragment. Then, perform one-hot encoding on each fragment to convert the DNA sequence into a four-channel numerical matrix. (2) Tag generation: Based on the cell-peak binary matrix, assign a corresponding open state tag to each sequence fragment, where "1" indicates that the peak is in an open state in this cell type and "0" indicates that it is not open; (3) Data partitioning: The dataset partitioning strategy is consistent with that of the pre-training stage, that is, the sequence fragments of chromosomes 7 and 8 are used as the test set, the chromosome 9 fragment is used as the validation set, and the remaining chromosome fragments are used as the training set.

[0013] In some embodiments of this application, the single-cell multi-omics data further includes a single-cell transcriptome dataset, and the method further includes the step of constructing an expression module for predicting the expression level of genes corresponding to sequence fragments: Based on the aforementioned regulatory module, the probability of chromatin accessibility in the region near the gene transcription start site is predicted, and the expression level of the corresponding gene for the sequence fragment is further predicted based on the single-cell transcriptome dataset.

[0014] In some embodiments of this application, the single-cell transcriptome dataset is derived from public databases, including but not limited to the GEO database.

[0015] In some embodiments of this application, the preprocessing procedure for the single-cell transcriptome dataset is as follows: (1) Expression matrix normalization: The original scRNA-seq expression matrix was normalized using the CPM method to eliminate the influence of sequencing depth differences on gene expression values; (2) Sparsity repair and imputation: The MAGIC algorithm is used to imput data in the normalized matrix. The potential true expression pattern is restored by a diffusion model based on the expression similarity of neighboring cells. (3) Data partitioning: The dataset is partitioned according to chromosome stratification. Genes on chromosome 8 are set as the test set, and genes on the other chromosomes are used as the training set.

[0016] In some embodiments of this application, the single-cell chromatin accessibility dataset and the single-cell transcriptome dataset need to be aligned: first, the single-cell chromatin accessibility dataset and the single-cell transcriptome dataset are normalized and dimensionality reduced to extract representative low-dimensional embedding features; then, anchor point identification is performed in the latent representation space of the two datasets to calculate the mutual information and feature similarity between cells to determine cross-modal correspondence.

[0017] The second aspect of this application provides a first model constructed using the method described in the first aspect of this application, wherein the first model includes only a control module.

[0018] A third aspect of this application provides a second model constructed using the method described in the first aspect of this application, the second model including an expression module.

[0019] The fourth aspect of this application provides a method for predicting the functional effects of non-coding mutation sites, including the step of predicting the effect of non-coding mutations on chromatin accessibility using the first model described in the second aspect of this application. Specifically, the first model is used to predict the chromatin accessibility probabilities of wild-type sequences and mutant sequences containing the non-coding mutation sites, respectively. Based on the difference in the chromatin accessibility probabilities of the two sequences, a non-coding mutation regulatory effect score is determined, and the effect of non-coding mutations on chromatin accessibility is further determined.

[0020] In some embodiments of this application, the method further includes the step of predicting the effect of non-coding mutations on gene expression using the second model described in the third aspect of this application. Specifically, the second model is used to predict the expression levels of the corresponding genes of wild-type sequences and mutant sequences containing the non-coding mutation sites, respectively, and the effect of the non-coding mutation on gene expression is determined by the difference in the expression levels of the corresponding genes.

[0021] The fifth aspect of this application provides a method for identifying cell type-specific regulatory mutations. For different mutation sites, the method described in the fourth aspect of this application is used to determine their regulatory effect scores in different cell types, and cell type-specific regulatory mutations are obtained based on statistical methods.

[0022] In some embodiments of this application, it further includes: The perturbation of transcription factor binding sites by cell type-specific regulatory mutations was assessed, and the transcription factor binding sites perturbed by cell type-specific regulatory mutations were enriched based on the software analysis results. Based on the aforementioned regulatory module, single-cell level prediction of transcription factor activity is achieved, and combined with the expression of transcription factors in the corresponding single-cell transcriptome data, cell type-specific positive regulatory transcription factors are inferred. Based on the above results, it is inferred that cell type-specific regulatory mutations mediate gene regulatory programs at transcription factor binding sites.

[0023] The sixth aspect of this application provides a method for constructing a cell type-specific "mutation-cis-regulatory element-gene" regulatory circuit, comprising the following steps: By integrating the cell type-specific regulatory mutations obtained in the fifth aspect of this application with eQTL data, mutation linkage disequilibrium, cis-regulatory element interactions, and chromatin interaction data for verification, a complete "mutation-cis-regulatory element-gene" regulatory circuit is established to elucidate the specific mechanism by which cell type-specific mutation sites regulate their target genes.

[0024] A seventh aspect of this application provides a computer device, comprising: a memory for storing a computer program; and a processor for executing the computer program to implement the steps of the methods described in the first, fourth, fifth, and sixth aspects of this application.

[0025] The eighth aspect of this application provides a computer-readable storage medium storing a computer program that, when executed by a processor, implements the steps of any of the methods described in the first, fourth, fifth, and sixth aspects of this application.

[0026] The eighth aspect of this application provides the use of a gene mutation detection reagent in the preparation of a kit for diagnosing Alzheimer's disease, wherein the gene mutation is A→G of rs10821626.

[0027] Compared with the prior art, this application has the following advantages: This application provides a unified deep learning framework that can directly predict chromatin accessibility and gene expression levels at single-cell resolution, starting from DNA sequences.

[0028] This application proposes two novel quantitative scoring systems: regulatory effect score and expression effect score. These systems can be based on the predictive effects of chromatin accessibility and gene expression, respectively, which significantly improve the accuracy and interpretability of non-coding variant functional annotation.

[0029] This application provides a new paradigm for understanding non-coding variants that reshape cell type-specific gene regulatory networks. By integrating systematic analysis of single-cell transcription factor activity, transcription factor expression, and mutation regulatory effects, it reveals cell type-specific regulatory mechanisms.

[0030] This application can be used to construct an integrated analysis process of "variation-cis-regulatory elements-genes", effectively identify cell type-specific regulatory circuits, and open up a mechanistic cognitive pathway between non-coding variations and regulatory dysfunction.

[0031] This application can be applied to research on diseases such as Alzheimer's disease, to elucidate the pathogenic cascade processes involving altered transcription factor binding, cis-regulatory remodeling, and key gene-specific transcriptional dysregulation related to candidate mutations.

[0032] It should be understood that the description in this section is not intended to identify key or essential features of the embodiments of this application, nor is it intended to limit the scope of this application. Other features of this application will become readily apparent from the following description. Attached Figure Description

[0033] The above and other objects, features, and advantages of exemplary embodiments of this application will become readily apparent from the following detailed description taken in conjunction with the accompanying drawings. Several embodiments of this application are illustrated in the drawings by way of example and not limitation, in which: Figure 1 UMAP projection maps of single-cell chromatin accessibility and transcriptome atlases of the amygdala, hippocampus, cerebellum, midbrain, basal ganglia, pons, cortex and thalamus in Example 1 of this application are shown.

[0034] Figure 2 The diagram shows the framework of the deep learning model in Embodiment 2 of this application. The model consists of a regulation module and an expression module, which can achieve joint prediction of single-cell chromatin accessibility and gene expression level.

[0035] Figure 3 The comparison of the prediction performance of the control module and the expression module in the deep learning model of Embodiment 2 of this application is shown.

[0036] Figure 4 A flowchart illustrating the de novo prediction of non-coding mutation regulation and expression effects using a model in Embodiment 3 of this application is shown.

[0037] Figure 5 The diagram shows a scatter plot illustrating the correlation between the PCC-based assessment scores of the regulatory effect and expression effect and the experimental observation results in Embodiment 3 of this application.

[0038] Figure 6 The flowchart shown illustrates the process of screening cell type-specific regulatory mutation sites based on the top four scores of the regulatory effect score in Example 4 of this application.

[0039] Figure 7The z-score heatmap based on AUCell analysis in Embodiment 4 of this application is shown, quantifying the enrichment of inferred regulatory variants in cell type-specific chromatin accessibility regions. The x-axis represents cell types from eight brain regions, and the y-axis represents chromatin accessibility regions for specific cell types.

[0040] Figure 8 The following are representative transcription factors showing the enrichment effect of microglia-specific regulatory mutation sites, obtained by MotifBreakR and SNP2TFBS calculation methods respectively in Example 4 of this application.

[0041] Figure 9 The flowchart illustrating the prediction of single-cell level transcription factor activity based on the regulatory module in Embodiment 4 of this application is shown.

[0042] Figure 10 The scatter plot shows the differential gene expression of transcription factors (x-axis) and the inferred activity score of the regulatory module (y-axis) in Example 4 of this application.

[0043] Figure 11 The effect of the mutant rs10203826 on the binding site of transcription factor IKZF1 in Example 4 of this application is shown, and the single-cell level regulatory effect score of this mutant and the corresponding transcription factor regulatory activity inferred by the regulatory module are displayed using UMAP projection diagram.

[0044] Figure 12 The flowchart of the analysis of the "variation-cis-regulatory element-gene" regulatory circuit in Example 5 of this application is shown, combined with various verification methods such as regulatory module prediction, eQTL intersection, cis-regulatory element interaction, chromatin interaction, and transcription factor binding effect.

[0045] Figure 13 This paper presents some candidate microglia-specific regulatory mutations that can establish a "mutation-cis-regulatory element-gene" regulatory circuit in Example 5 of this application, and uses a heatmap to show the associated genes and supporting evidence.

[0046] Figure 14 This document illustrates a specific example of how the microglia-specific candidate mutation (rs1293303) affects the associated gene (CTSB) in Example 5 of this application.

[0047] Figure 15 The analysis process for screening AD-related pathogenic mutations in Embodiment 6 of this application is shown.

[0048] Figure 16 The comparison of significance, regulatory effect, and expression effect scores of candidate pathogenic mutations and all mutations in the GWAS dataset within a 1Mb region in Embodiment 6 of this application is shown.

[0049] Figure 17 The representative AD pathogenic mutations screened in Example 6 of this application are shown, and the associated genes and supporting evidence are displayed using a heatmap.

[0050] Figure 18 This document illustrates a specific example of how the candidate AD pathogenic mutation (rs10821626) affects the regulatory mechanism of the downstream AD risk gene (CCDC6) in Example 6 of this application.

[0051] Figure 19 The example in Example 6 of this application shows a genomic region 100 bp upstream and downstream of the candidate AD pathogenic mutation (rs10821626) as the center, and the nucleotide contribution score calculated based on simulated saturation mutations presents a microglia-specific sequence pattern.

[0052] Figure 20 The comparison of ATAC-seq read counts at the rs10821626 site in microglia derived from human pluripotent stem cells in Example 6 of this application is shown. The accessibility of the reference allele (A) is significantly higher than that of the risk allele (G), and the statistical significance is assessed by a one-sided t-test.

[0053] Figure 21 The diagram shows the qPCR results of CCDC6 expression after base editing (A>G) at the rs10821626 site in Example 6 of this application, with GAPDH (left) and ACTB (right) as internal reference genes. Compared with the control group, the introduction of the risk allele (G) reduced CCDC6 expression. Detailed Implementation

[0054] Unless otherwise stated, implied from the context, or as is customary in the art, all parts and percentages in this application are based on weight, and all testing and characterization methods used are concurrent with the filing date of this application. Where applicable, any patent, patent application, or disclosure relating to this application is incorporated herein by reference in its entirety, and its equivalent patent families are also incorporated herein by reference, particularly the definitions of relevant terms in the art disclosed in such documents. If any definition of a specific term disclosed in the prior art is inconsistent with any definition provided in this application, the definition provided in this application shall prevail.

[0055] The numerical ranges used in this application are approximate values ​​and therefore may include values ​​outside the range unless otherwise stated. The numerical range includes all values ​​from the lower limit to the upper limit, increasing by one unit, provided that there is an interval of at least two units between any lower and any higher value. For ranges containing values ​​less than 1 or fractions greater than 1 (e.g., 1.1, 1.5, etc.), one unit is appropriately considered as 0.0001, 0.001, 0.01, or 0.1. For ranges containing single digits less than 10 (e.g., 1 to 5), one unit is generally considered as 0.1. These are merely specific examples of what is intended to be expressed, and all possible combinations of values ​​between the listed lowest and highest values ​​are considered to be clearly described in this application.

[0056] The terms “comprising,” “including,” “having,” and their derivatives do not exclude the presence of any other components, steps, or processes, regardless of whether such other components, steps, or processes are disclosed in this application. To eliminate any doubt, unless expressly stated otherwise, all compositions using the terms “comprising,” “including,” or “having” in this application may contain any additional additives, excipients, or compounds. Conversely, except for those necessary for operational performance, the term “substantially constitutes…” excludes any other components, steps, or processes described below with respect to that term. The term “consisting of…” does not include any components, steps, or processes not specifically described or listed. Unless expressly stated otherwise, the term “or” refers to the individual members listed or any combination thereof.

[0057] To make the technical problems, technical solutions and beneficial effects solved by this application clearer, the following detailed description is provided in conjunction with embodiments.

[0058] The following examples are used to illustrate preferred embodiments of this application. Those skilled in the art will understand that the techniques disclosed in the examples represent technologies discovered by the inventors that can be used to implement this application, and therefore can be considered preferred embodiments of this application. However, those skilled in the art should understand from this specification that many modifications can be made to the specific embodiments disclosed herein, still yielding the same or similar results, without departing from the spirit or scope of this application.

[0059] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this application pertains, and all materials cited herein and referenced by them are incorporated herein by reference.

[0060] Those skilled in the art will recognize, or can learn through routine experimentation, many equivalents of the specific embodiments of the invention described herein. These equivalents will be included in the claims.

[0061] Unless otherwise specified, the experimental methods used in the following examples are conventional methods. Unless otherwise specified, the instruments and equipment used in the following examples are all conventional laboratory instruments and equipment; unless otherwise specified, the experimental materials used in the following examples were all purchased from conventional biochemical reagent stores.

[0062] Example 1: Data Processing of Human Brain Single-Cell Omics Atlas 1. Collection of scATAC-seq and scRNA-seq datasets The inventors collected single-cell ATAC sequencing (scATAC-seq) and single-cell RNA sequencing (scRNA-seq) datasets from the GEO database and publicly available normal adult brain tissue, as shown in Table 1: Table 1: Statistical analysis of human brain scRNA-seq and scATAC-seq data

[0063] The scATAC-seq dataset integrates approximately 60 samples, covering over 1.2 million single cells. The samples originate from multiple key brain regions, including the cerebral cortex, hippocampus, thalamus, pons, basal ganglia, cerebellum, and amygdala. This data possesses high representativeness and organizational diversity, comprehensively reflecting the chromatin accessibility of different brain regions.

[0064] The scRNA-seq dataset comprises approximately 230 samples, totaling over 3.7 million single cells. These samples cover key brain regions including the cerebral cortex, hippocampus, thalamus, pons, basal ganglia, cerebellum, and amygdala. This dataset exhibits high representativeness and organizational diversity, enabling a comprehensive characterization of gene expression profiles across different brain regions.

[0065] 2. Preprocessing of human brain scATAC-seq data To ensure the reliability and comparability of the scATAC-seq data, the inventors performed rigorous quality control and preprocessing on the raw data. The specific steps are as follows: (1) Dual threshold filtering strategy: adopt dual index screening criteria. First, remove low-quality cells with fewer than 1,000 unique fragments to exclude samples with insufficient sequencing depth or high noise signal. Second, calculate the transcription start site (TSS) enrichment score of each cell and retain only high-quality cells with a TSS enrichment score greater than 7.

[0066] (2) Peak Recall and Redundancy Elimination: Using cell type as the unit of analysis, MACS2 software was used to perform chromatin accessibility peak recall on the screened high-quality data. To unify the peak sets among different cell types and eliminate repetitive signals, a fixed-width peak integration strategy in the ArchR workflow was adopted. Specifically, the center position (peak summit) of each peak was used as the anchor point, and extended 250 bp to both sides (total length 500 bp) to construct a standardized peak region of fixed length. If the distance between the peak centers of different cell types was less than 250 bp, they were considered the same peak and merged. After merging, further screening was performed based on signal intensity and repeatability, retaining only the peak regions that were repeatedly detected in multiple samples and whose significance level met the threshold (q-value < 0.05). Finally, a unified, non-redundant, and high-confidence peak set was obtained for the subsequent construction of the cell-peak matrix.

[0067] (3) Constructing the cell-peak matrix: The cell-peak binary matrix is ​​generated using the `addPeakMatrix` and `getMatrixFromProject` functions in the ArchR package. To ensure the representativeness of peaks in the population, only peak rows that are open in at least 25 cells are retained. The matrix is ​​in binary encoding, where "1" indicates that the peak is open in the corresponding cell, and "0" indicates that the peak is closed in that cell.

[0068] After the above preprocessing, the resulting scATAC-seq dataset can accurately reflect the chromatin accessibility characteristics of different cell types, laying a high-quality data foundation for subsequent cross-omics integration and deep learning modeling of epigenetic regulation at the single-cell level.

[0069] 3. Preprocessing of scRNA-seq datasets To ensure data quality and comparability across studies, the inventors employed a multi-dimensional quality control (QC) and filtering strategy to preprocess the scRNA-seq dataset. Specifically, this included: (1) Transcript filtering: Low-quality cells with fewer than 200 detected transcripts are removed to eliminate samples with insufficient sequencing depth or low transcriptional activity. (2) Mitochondrial gene ratio control: Filter cells with a mitochondrial gene expression ratio exceeding 10% to exclude cells that may be in apoptosis or stress state; (3) Double cell detection: Based on the doublet score, possible mixed signals of two cells are removed, and only high-confidence single cells with doublet scores less than 0.3 are retained; (4) Cell type annotation and screening: retain high-quality cell type annotations defined manually or by algorithm in the original study, and delete cells whose cell identity could not be clearly determined; at the same time, to ensure statistical reliability, remove rare cell types with fewer than 100 cells.

[0070] After the above processing, the resulting scRNA-seq dataset exhibits a high signal-to-noise ratio, complete cell type annotation, and good data consistency. This high-quality single-cell transcriptome data can provide a reliable foundation for subsequent integration analysis with scATAC-seq data, model training, and cross-omics feature learning.

[0071] 4. Cell type alignment using scRNA-seq and scATAC-seq The inventors utilized the FindTransferAnchors and IntegrateData modules of the Seurat software to achieve cross-omics integration and cell type alignment of single-cell transcriptome (scRNA-seq) and chromatin accessibility (scATAC-seq) data. This method is based on the principle of Canonical Correlation Analysis (CCA), establishing a shared low-dimensional representation space across different omics datasets by identifying cross-modal feature anchors. The specific steps are as follows: (1) Data preprocessing and feature extraction: Standardization and high-variability gene screening were performed on scRNA-seq and scATAC-seq data, respectively. The scRNA-seq data were normalized using the LogNormalize method and high-variability genes were extracted. The gene activity matrix of the scATAC-seq data was calculated, and normalization and dimensionality reduction were performed using TF-IDF (Term Frequency-Inverse Document Frequency) and SVD (Singular Value Decomposition) to obtain representative latent feature vectors.

[0072] (2) Cross-modal feature mapping and CCA calculation: Using Seurat's FindTransferAnchors function, with the CCA algorithm as the core, canonical correlation analysis is performed on the low-dimensional feature vectors of the two sets of data to identify the feature combination with the greatest linear correlation between the transcriptome and the accessibility group.

[0073] (3) Anchor selection and joint embedding construction: The set of highly relevant anchor pairs is used to construct a shared potential space across omics. The anchor weight normalization and batch effect correction are performed by the IntegrateData function to realize the joint representation of different omics data in a unified coordinate system.

[0074] After processing with this integration strategy, scATAC-seq data can be transferred to obtain high-quality cell type information from scRNA-seq, achieving consistent annotation across omics and providing an accurate cell type comparison and classification basis for subsequent model training and regulatory mechanism analysis.

[0075] 5. Establishment of a multi-omics database for single-cell brain tissue Building upon this foundation, the inventors further developed a comprehensive human brain single-cell multi-omics database for systematically integrating, storing, and displaying the single-cell regulatory and expression characteristics of different brain regions. The platform uses the Bootstrap framework to construct the front-end interface, possessing excellent adaptability and interactivity; the back-end is driven by PHP, R, and a MySQL database, supporting efficient data querying, visualization, and dynamic analysis. A browsing module is included to display the integrated human brain single-cell multi-omics dataset, including sample size, cell type distribution, and statistical information for each brain region. Users can quickly retrieve dataset summaries by brain region, cell type, or data type.

[0076] Figure 1 This presentation showcases integrated single-cell multi-omics data from eight major brain regions: the amygdala, hippocampus, cerebellum, midbrain, basal ganglia, pons, cortex, and thalamus. The single-cell multi-omics data for each brain region are derived from high-quality, rigorously quality-controlled samples, comprehensively reflecting the cell type composition and epigenetic regulatory characteristics of different brain regions.

[0077] 6. Creation of the pre-training dataset for the regulation module The inventors pre-trained the regulatory module using high-quality epigenomic data from the ENCODE database (815 epigenomics), the Roadmap Epigenomics database (1187 epigenomics), and the Cistrome database (19905 epigenomics). The dataset contains 21907 epigenomic features, including 9471 transcription factor binding sites, 10064 histone modification signal regions, and 2372 chromatin accessibility peak regions. By integrating multi-source epigenomic data, comprehensive learning and modeling of regulatory signals from different sources can be achieved. The specific production steps are as follows: (1) Sequence Segmentation and Tag Generation: The human reference genome GRCh38 was divided into non-overlapping 200 bp sequence fragments. The inventors compared each sequence fragment with the peak regions of the aforementioned public epigenome dataset. For any epigenomic feature, if the overlap between the sequence fragment and the peak region corresponding to that feature exceeded 50%, the tag position corresponding to that feature was marked as 1; otherwise, it was marked as 0. Thus, each sequence fragment was assigned a 21907-dimensional binary tag vector, where each dimension corresponds to a specific epigenetic regulatory feature.

[0078] (2) Sequence Expansion and Encoding: Using the midpoint of each sequence fragment as a reference coordinate, the sequence is expanded 1000 bp upstream and downstream to obtain a DNA sequence fragment of length 2000 bp. Subsequently, the inventors used one-hot encoding to digitally represent each base in the sequence, where the bases "A", "C", "G", and "T" are encoded as [1,0,0,0], [0,1,0,0], [0,0,1,0], and [0,0,0,1], respectively. This method can transform the DNA sequence into a multidimensional vector suitable for deep learning input.

[0079] (3) Dataset partitioning: To ensure the independence of model evaluation and generalization performance, a chromosome-level data partitioning strategy is adopted. Specifically, sequence fragments on chromosomes 7 and 8 are used as the test set, fragments on chromosome 9 are used as the validation set, and the remaining chromosome fragments are used as the training set. This partitioning method can avoid the information leakage problem caused by overlapping sequence positions and ensure the independence between model training and test data.

[0080] The pre-training dataset constructed through the above steps can comprehensively cover the multidimensional regulatory features of the human genome, providing rich prior knowledge of regulatory semantics for regulatory modules.

[0081] 7. Creation of the single-cell chromatin accessibility dataset for formal training of the regulatory module. Based on the cell-peak binary matrix obtained from the aforementioned "collection and preprocessing of human brain scATAC-seq data," the inventors constructed a single-cell chromatin accessibility dataset for the formal training of the regulatory module. The specific process is as follows: (1) Sequence extraction and encoding: Taking the center point of each chromatin accessibility peak as a reference, extend 1000 bp upstream and downstream to generate a 2000 bp genomic sequence fragment. Then, one-hot encoding is performed on each fragment to convert the DNA sequence into a four-channel numerical matrix as input to the deep learning model.

[0082] (2) Tag generation: Based on the cell-peak binary matrix, each sequence fragment is assigned a corresponding open state label. "1" indicates that the peak is open in this cell type, and "0" indicates that it is not open. Finally, a cell type-specific two-dimensional label matrix is ​​obtained, which is used to supervise the model to learn cell-specific chromatin accessibility patterns.

[0083] (3) Data partitioning: The dataset partitioning strategy is consistent with that of the pre-training stage, that is, the sequence fragments of chromosomes 7 and 8 are used as the test set, the chromosome 9 fragment is used as the validation set, and the remaining chromosome fragments are used as the training set. This partitioning method ensures the independence between the model training data and the test data.

[0084] 8. Creation of the Expression Module Dataset To train the expression module, the inventors used a high-quality single-cell transcriptome matrix obtained from the collection and preprocessing of human brain scRNA-seq data to construct an expression prediction training dataset. The specific process is as follows: (1) Expression matrix normalization: The original scRNA-seq expression matrix was normalized using the CPM (Counts Per Million) method to eliminate the influence of sequencing depth differences on gene expression values.

[0085] (2) Sparsity Repair and Imputation: Due to amplification bias and low capture efficiency in single-cell RNA sequencing, a large number of zero values ​​exist in the expression matrix, especially in low-expression genes. The inventors used the MAGIC algorithm to imput data in the normalized matrix and restored the potential true expression patterns through a diffusion model based on the expression similarity of neighboring cells. This step significantly solved the sparsity problem of the expression matrix and provided more stable input features for model training.

[0086] (3) Data partitioning: The dataset is partitioned according to chromosome stratification, with genes on chromosome 8 serving as the test set and genes on the remaining chromosomes serving as the training set. This partitioning method is also to ensure the independence between the model training and test data.

[0087] Example 2: Joint prediction of chromatin accessibility and gene expression levels based on a unified deep learning framework This embodiment, based on the high-quality model dataset obtained from the "Human Brain Multi-omics Comprehensive Database" constructed in Embodiment 1, develops a unified deep learning framework consisting of a regulation module and an expression module, used to predict chromatin accessibility and gene expression levels at single-cell resolution from DNA sequences, respectively. The specific framework is as follows: Figure 2 As shown.

[0088] 1. Control Module The regulation module adopts a dual-branch hybrid architecture, consisting of a DeepSEA pre-training branch and a ResNet18 end-to-end learning branch, which are integrated through a feature fusion module. This dual-branch structure allows the model to simultaneously learn universal regulatory rules shared across samples and cell-specific local accessibility features, achieving the fusion expression of multi-level regulatory information.

[0089] (1) DeepSEA pre-trained branch structure: This branch utilizes the DeepSEA structure pre-trained on a public epigenome dataset to capture general epigenetic regulatory signals. The model contains three one-dimensional convolutional layers (with a kernel size of 7 and channel numbers of 320, 480, and 960 respectively) and two fully connected layers, ultimately outputting 21,907-dimensional epigenomic features to characterize the feature expression of DNA sequences at the general regulatory level.

[0090] (2) ResNet18 end-to-end learning branch structure: This branch is used to learn specific regulatory patterns at the single-cell level from scratch. The model structure includes an input layer, a residual network layer, and an output layer. The input layer is a one-dimensional convolutional layer (with a kernel size of 4 and 48 channels). The residual network layer consists of 8 progressive residual blocks with channel numbers of 96, 96, 128, 128, 256, 256, 512, and 512, respectively. Each residual block contains two convolutional layers and a ReLU activation function, used to progressively extract nonlinear regulatory features from the sequence. The output layer is a two-layer fully connected layer used to generate single-cell-specific feature representations of the input sequence fragments.

[0091] (3) Feature fusion module: The epigenomic features output by the DeepSEA pre-training branch are fused with the single-cell specific features generated by the ResNet18 end-to-end learning branch. The fusion result is nonlinearly mapped through three fully connected layers, and finally outputs the chromatin accessibility probability prediction value of each input sequence fragment at the single-cell level.

[0092] 2. Expression Module Architecture The expression module consists of three main parts: epigenomic feature extraction, spatial transformation dimensionality reduction, and ridge regression, as detailed below: (1) Epigenomic Feature Extraction: Based on the regulatory module trained on scATAC-seq data, for each of the 24,339 genes annotated by Gencode (Release 19, GRCh37.p13), a region 20,000 bp upstream and downstream of the transcription start site (TSS) was defined, and a sliding window scan was performed within this region. Each window was 2,000 bp long with a step size of 200 bp, resulting in 200 overlapping windows within the TSS ±20 kb range of each gene. These sequences were input into the regulatory module to obtain the corresponding predicted values ​​of chromatin accessibility features.

[0093] (2) Spatial Transformation Dimensionality Reduction: To reduce the dimensionality of the predicted chromatin accessibility features, this embodiment employs a weighted dimensionality reduction strategy. Specifically, position-dependent weights are assigned to 200 sliding windows within the TSS ± ​​20 kb region to reflect the biological hypothesis that the influence of distal regulatory elements on gene expression gradually weakens. A weight matrix is ​​constructed. It contains ten exponentially decaying functions; for each decay constant Calculate the weights of the upstream and downstream windows separately; the weight calculation formula is:

[0094] in Display window The distance to TSS. The resulting accessibility prediction matrix. (in (For the number of single cells), dimension reduction projection is performed using matrix multiplication: finally, a low-dimensional matrix is ​​obtained. This matrix preserves location-dependent regulatory information in a compact form and serves as input to the subsequent linear regression module.

[0095] (3) Ridge regression model: A ridge regression model is trained independently for each cell type, using the transformed chromatin accessibility feature matrix. Using the input as input, predict the expression level of the corresponding gene. Ridge regression suppresses feature collinearity through L2 regularization, improving the robustness and generalization ability of the model.

[0096] 3. Model Training The training scheme for the regulation module consists of two phases: pre-training and joint training. In the pre-training phase, the DeepSEA model is trained using public epigenetic datasets (from ENCODE, Roadmap Epigenomics, and Cistrome databases) to learn general epigenomic features. In the joint training phase, the DeepSEA pre-training structure and the ResNet18 end-to-end learning structure are jointly trained using a processed single-cell chromatin accessibility dataset. The DeepSEA pre-training structure employs deep transfer learning, inheriting pre-trained weights and continuously optimizing parameters through backpropagation; the ResNet18 module has completely randomized initial parameters to learn single-cell-level chromatin accessibility patterns from scratch. Both phases are trained on NVIDIA 1080Ti GPUs, using binary cross-entropy as the loss function, employing an early stopping strategy and a cosine annealing learning rate scheduler, and performing gradient comparison to select the optimal model hyperparameters, such as batch size and learning rate. The expression module uses a ridge regression linear model, and the L2 regularization coefficient is adjusted using grid search during the training phase. And select the parameter that performs best on the validation set ( ).

[0097] 4. Model Evaluation The performance evaluation of the regulatory module is mainly based on the following two aspects: First, evaluation at the cellular level. For each single cell, the area under the receiver operating characteristic curve (AUROC) is calculated for the peak values ​​reserved in the test set. This measures the model's ability to distinguish between open and closed regions at the cellular level. Second, evaluation at the peak level: For each peak value in the test set, the AUROC at the peak level is calculated, reflecting the model's accuracy in predicting the openness of different regulatory elements.

[0098] The expression module uses Pearson coefficient (PCC) and AUROC to evaluate predictive performance. PCC measures the linear correlation between the model's predicted gene expression levels and actual experimental data; AUROC evaluates the model's ability to distinguish between high and low gene expression. Figure 3 As shown, the prediction performance of both the regulation module and the expression module of this deep learning framework is superior to existing models.

[0099] Example 3: De novo prediction of regulatory and expression effects by non-coding mutations This embodiment provides a mutation function evaluation method based on a unified deep learning framework. This method can simultaneously predict the regulatory effect score and expression effect score of mutations at single-base resolution. The specific process is as follows: Figure 4 As shown.

[0100] 1. Predicting the regulatory effects of mutations A regulatory module trained on single-cell chromatin accessibility data can analyze the regulatory effects of single-base mutations at the single-cell level. The specific steps are as follows: A DNA sequence centered on the mutation site and covering 1000 bp upstream and downstream is constructed to simulate the chromatin environment of the reference allele and the alternative allele, respectively. This sequence is then input into the regulatory module, which outputs predicted chromatin accessibility values ​​for both allele sequences. The difference in chromatin accessibility caused by a single-base mutation is defined as the regulatory effect score of the mutation. Its mathematical expression is:

[0101] in, and These represent the predicted accessibility of reference and alternative allele sequences, respectively.

[0102] The larger the absolute value of the regulatory effect score, the more significant the effect of the mutation on the open state of chromatin; a positive value indicates that the mutation may enhance chromatin accessibility, while a negative value indicates that it has an inhibitory effect.

[0103] 2. Predicting the expression effects of mutations For each single-base mutation, the regulatory module generates allele-specific accessibility features, and predicts the corresponding gene expression changes using a ridge regression model independently modeled for each cell type. The inventors score the mutation's effect on expression. Defined as the absolute value of the difference between the expression levels inferred from the reference allele sequence and the alternative allele sequence, its mathematical expression is:

[0104] in, and These represent the predicted expression levels of the reference and alternative allele sequences, respectively.

[0105] Expression effect score This reflects the extent to which mutations affect the expression level of target genes at the transcriptional level. A higher score indicates a more significant change in expression caused by the mutation, potentially leading to stronger biological functional consequences.

[0106] 3. Evaluate the accuracy of the model in predicting mutation regulation and expression effects. To verify the reliability and biological consistency of the model's predictions, the inventors systematically evaluated the model using publicly available chromatin accessibility QTL (caQTL) and expression QTL (eQTL) datasets.

[0107] (1) Evaluation data sources: The caQTL dataset comes from the research results of Zeng et al., Genetic regulation of cell type-specific chromatin accessibility shapes brain diseaseetiology (Science, 2024); the eQTL dataset comes from the GTEx database (website: https: / / www.gtexportal.org / home / downloads / adult-gtex / qtl).

[0108] (2) Evaluation indicators and methods: The correlation between the predicted regulatory effect score and the expression effect score and the experimental observation value was calculated using the Pearson Correlation Coefficient (PCC).

[0109] (3) Results Analysis: The evaluation results showed that the caQTL site regulatory effect scores predicted by the regulatory module were significantly positively correlated with the experimentally measured changes in chromatin accessibility; the eQTL expression effect scores predicted by the expression module also showed high consistency with the measured expression changes (e.g., Figure 5 (As shown).

[0110] Example 4: Regulatory program for recognizing cell type-specific non-coding variants This embodiment, based on the unified deep learning framework established in Embodiment 2, combines publicly available single-cell multi-omics data to achieve the identification and analysis of cell type-specific non-coding mutations and their regulatory mechanisms. The overall process is as follows: Figure 6 As shown.

[0111] 1. Screening for cell type-specific regulatory mutation sites The inventors sorted the mutation regulatory effect scores predicted by the regulatory module from high to low, and selected mutations in the upper quartile distribution range as cell type-specific candidate regulatory mutation sites.

[0112] To verify the rationality of the cell type specificity of the selected mutations, the inventors introduced an external validation dataset—the human brain single-cell chromatin accessibility map (covering eight brain regions and more than 30 cell subtypes) from Li et al. (Science, 2023). The AUCell algorithm was used to validate the distribution of the predicted cell type-specific regulatory mutations' accessibility across different cell types.

[0113] The verification results are as follows Figure 7As shown: the predicted cell type-specific regulatory mutations exhibit significant open chromatin enrichment signals in their corresponding cell types; the hierarchical clustering structure among different cell subtypes conforms to their biological lineage.

[0114] 2. Analysis of the perturbation effect of mutations on transcription factor binding To clarify the potential impact of mutations on transcription factor (TF) binding sites, the inventors used two independent bioinformatics tools, MotifBreakR and SNP2TFBS, to conduct quantitative analysis on cell type-specific mutations in microglia.

[0115] (1) MotifBreakR analysis: The upstream and downstream sequences of each mutation site were scanned, and the relative entropy difference (ΔEntropy) between the reference allele and the alternative allele was calculated based on the TF binding motif in the JASPAR2022 database to measure the change in TF binding affinity caused by the mutation.

[0116] (2) SNP2TFBS analysis: By comparing the mutation location with a pre-calculated set of TF binding sites, the change in TF binding affinity before and after the mutation is quantified, and the TF binding perturbation score is obtained.

[0117] The results are as follows Figure 8 As shown, both MotifBreakR and SNP2TFBS analyses revealed significantly enriched perturbed TF types in the candidate mutation set, including IRF and IKZF family transcription factors.

[0118] 4. Prediction of transcription factor activity at single-cell resolution To further characterize the dynamics of transcriptional regulation under the influence of mutations, the inventors proposed a single-cell level transcription factor activity scoring method using a trained regulatory module. The specific process is as follows (e.g.) Figure 9 (as shown) (1) Background sequence generation: 2000 peaks are randomly sampled from the scATAC-seq peak set and the sequence order is randomly shuffled to construct a background sequence library; (2) Inserting transcription factor binding sequences: The standard binding sequences of the TFs to be studied were extracted from the JASPAR database and inserted into the center of the background sequence; (3) Activity inference: The regulatory module is used to predict the inserted sequence and the background sequence respectively, calculate the difference in chromatin accessibility between the two, and obtain the TF activity score at single-cell resolution by Z-score standardization.

[0119] This method overcomes the limitations of traditional motif enrichment analysis, enabling quantitative prediction of transcription factor regulatory activity at the single-cell level, thereby revealing the dynamic characteristics of cell type-specific transcriptional regulatory networks.

[0120] 5. Identification of positively regulated transcription factors The inventors further combined the model-inferred single-cell TF activity fractions with corresponding single-cell transcriptome data to identify positively regulatory transcription factors that simultaneously exhibit high expression and high regulatory activity. The screening criteria are as follows: (1) The differential expression level (log2FC) of transcription factors in the single-cell expression matrix is ​​greater than 0.25; (2) The activity fraction inferred by the regulatory module is located in the upper quartile range.

[0121] The filtering results are as follows Figure 10 As shown, transcription factors such as IRF8 and IKZF1 exhibit both high expression and high activity in microglia, indicating their central role in cell type-specific regulation. This strategy can systematically identify key TTFs with regulatory drivers in specific cell types, providing a basis for constructing cell-specific regulatory networks.

[0122] 5. Identification of regulatory programs for non-coding mutations To reveal the regulatory patterns of cell type-specific non-coding mutations, the inventors combined information about the perturbation of TF binding by mutations with model-inferred single-cell TF activity and transcriptome expression data to form a comprehensive regulatory inference framework. Taking the mutation rs10203826 as an example (e.g.) Figure 11 As shown in the figure): In microglia, this mutation exhibits a significant regulatory effect score; MotifBreakR analysis results show that the mutation disrupts the binding site of the IKZF1 transcription factor; IKZF1 has a specific high activity score and high expression level in microglia.

[0123] The combined results indicate that this mutation alters the expression of downstream target genes by perturbing the IKZF1 binding site, thereby producing a cell type-specific transcriptional regulatory effect. It also suggests that the screened cell type-specific mutations specifically affect the binding sites of transcription factors that are highly active in their respective cell types.

[0124] Example 5: Construction of a microglial cell-specific "mutation-cis-regulatory element-gene" regulatory circuit This embodiment aims to combine deep learning model predictions with multi-omics evidence to systematically analyze how microglia-specific non-coding mutations affect target gene expression through cis-regulatory elements, thereby constructing a "mutation-cis-regulatory element-gene" regulatory loop. The overall analysis process is as follows: Figure 12 As shown.

[0125] 1. Integrate eQTL data to determine the association between mutations and gene expression. The inventors first integrated eQTL (expression Quantitative Trait Loci) datasets at the tissue and cell type levels. By integrating eQTL data, the functional impact of mutations can be verified at the gene expression level, so as to establish the functional link between mutations and gene expression.

[0126] (1) Data sources: including tissue-level eQTLs and cell-type-level eQTLs from the GTEx project (from Bryois et al., Nature Neuroscience, 2022).

[0127] (2) Intersection analysis: The cell type-specific regulatory mutations obtained by deep learning model screening are matched with eQTL sites by coordinate intersection; (3) Target gene identification: If the mutation site overlaps with the eQTL signal, the target gene corresponding to the eQTL is included in the mutation-associated gene set.

[0128] 2. Chain Disequilibrium Analysis To further capture potential functional mutations associated with and regulating mutation signals, the inventors used PLINK software to calculate the linkage disequilibrium coefficient (LD) between the mutation and neighboring sites. By introducing LD information, allelic variations that have a co-genetic effect with the mutation can be captured.

[0129] (1) Parameter settings: Calculate the LD value of each regulatory mutation and the SNP within a 1 Mb range; (2) Screening criteria: Mutations with LD>0.8 are retained – eQTL pairings are considered to have shared genetic signals; (3) Expanding the candidate mutation set: Including eQTL target genes that share signals into the mutation-associated gene set, thereby expanding the "mutation-gene" regulatory network.

[0130] 3. Interaction Analysis of Cis-Regulator Elements (1) Enhancer-promoter co-accessibility analysis Based on the cell-peak matrix obtained from the scATAC-seq data, co-accessibility analysis of chromatin accessibility signals was performed using Cicero software to identify potential three-dimensional interactions between enhancers and promoters. Specifically, chromatin open peaks were divided into window intervals according to genomic location, and the degree of co-accessibility between each peak at the single-cell level was calculated. The co-accessibility score calculated by the Cicero algorithm was used to quantitatively describe the spatial connection strength between two peaks.

[0131] The inventors screened enhancer-promoter pairings with accessibility scores greater than 0.1 as significant interaction events, and further located the corresponding target genes based on these pairings. Subsequently, by comparing known mutation sites with the enhancer-promoter regions, it was determined whether the mutation fell within or at the boundary of a key regulatory element. If the mutation is located on a highly co-accessible enhancer, it may weaken the regulatory connection between the enhancer and promoter by disrupting transcription factor binding motifs or altering chromatin conformation, thereby affecting the expression of the target gene.

[0132] (2) ABC model analysis The inventors further integrated scATAC-seq data with publicly available human brain Hi-C spatial conformation data to construct a three-dimensional chromatin regulatory model based on the Activity-by-Contact (ABC) principle, which is used to identify potential cis-regulatory elements and gene connections.

[0133] Specifically, enhancer activity is first calculated based on chromatin accessibility and histone modification signals (such as H3K27ac) in enhancer regions from scATAC-seq data. Then, Hi-C data is used to estimate the spatial contact probability between enhancers and promoters of various genes. The ABC score of each enhancer-gene pair is obtained by multiplying the two scores, which is used to quantitatively characterize its regulatory potential.

[0134] The inventors screened enhancer-gene connections with ABC scores greater than 0.015 as high-confidence cis-regulatory relationships. For each high-scoring enhancer region, they further located mutation sites and inferred the possible target genes and regulatory mechanisms of the mutations.

[0135] 4. Validation of Plac-seq experimental data The inventors downloaded Plac-seq (Proximity ligation-assisted ChIP-seq) experimental data from the UCSC database, accessible at: https: / / genome.ucsc.edu / s / nottalexi / glassLab_BrainCellTypes_hg19. This dataset, derived from human brain tissue samples, directly reflects the spatial folding state of chromatin and the physical interactions between different regulatory regions, providing structural evidence for the functional annotation of non-coding mutations.

[0136] To ensure coordinate system consistency across different data sources, the inventors used the LiftOver genome coordinate transformation tool, along with the hg19ToHg38.over.chain.gz chain file, to accurately convert Plac-seq data from the hg19 version of the genome to the hg38 version. The converted data maintains consistent genomic coordinates with mutation sites and scATAC-seq accessibility data, thus enabling unified analysis across platforms.

[0137] Building upon this foundation, the inventors performed peak identification and interval annotation on the Plac-seq interaction signals to locate significant chromatin interaction regions. Combined with mutation site information, they screened for events where mutations fell within these interaction regions. If the mutation is located in an enhancer region and this region significantly interacts with a gene promoter via the Plac-seq signal, it indicates that the mutation may regulate the target gene by altering the stability or frequency of the enhancer-promoter physical connection.

[0138] 5. Analysis of the impact of mutations on transcription factor binding The inventors used the MotifBreakR algorithm to evaluate the perturbation effect of mutations on transcription factor binding sites. This perturbation can be categorized into two types: (1) Destructive mutations: disrupt the original transcription factor binding sequence and reduce binding affinity; (2) Degenerative mutation: forming new transcription factor binding sites and enhancing the binding potential of specific transcription factors.

[0139] The inventors ultimately screened 186 high-confidence microglia-specific regulatory mutations that met multi-dimensional support (eQTL co-occurrence, LD sharing, and spatial interaction verification), such as Figure 13 As shown, taking the mutant rs1293303 (G>C) as an example: this mutation disrupts the binding site of the ELF3 transcription factor, and the regulatory effect score in microglia is significantly lower than that in other cell types, indicating that this mutation significantly downregulates chromatin accessibility; at the same time, this mutation is a microglia-specific eQTL, and the associated gene is CTSB; and the regulatory effect of this mutation is supported by Cicero co-accessibility analysis, ABC model prediction, and Plac-seq spatial interaction data. Figure 14 ).

[0140] Example 6: Screening for Alzheimer's Disease-Related Causal Mutations and Elucidation of Microglial-Specific Pathogenic Mechanisms This embodiment, based on the aforementioned deep learning framework and multi-omics integrated analysis, aims to systematically identify potential causal mutations associated with Alzheimer's Disease (AD) and reveal their specific regulatory mechanisms in microglia. The overall process is as follows: Figure 15 As shown.

[0141] 1. Screening for potential pathogenic mutations related to Alzheimer's disease (AD). (1) Collection and filtering of GWAS genetic mutation sites The inventors collected GWAS results from multiple large-scale AD genetic studies (Jansen et al., 2019; Kunckle et al., 2019; Wightman et al., 2021; Bellenguez et al., 2022), and after integration, obtained approximately over 10,000 relevant mutations. These mutations were then analyzed according to a significance threshold of P-value < 5 × 10⁻⁶. -8 High-confidence mutations were screened and functionally annotated using the Ensembl Variant Effect Predictor (VEP) tool, retaining only mutations located in non-coding regions to avoid interference from protein-coding variations in the results.

[0142] (2) Screening for potential pathogenic mutations Using each mutation site as the center, a 2,000 bp DNA sequence is extracted as input and fed into both the regulatory and expression modules to predict the regulatory and expression effect scores of the mutation. The inventors then weighted and combined these two effect scores to screen out a set of potential AD-related pathogenic mutations, such as... Figure 16 As shown, the regulatory and expression effect scores of these candidate mutations were significantly higher than those of GWAS mutations in the surrounding 1Mb region, indicating that this strategy can accurately locate AD-related mutation sites.

[0143] 2. Analysis of gene regulatory programs for potentially pathogenic mutations Based on the "mutation-cis-regulatory element-gene" regulatory loop analysis process constructed in Example 5, the inventors further analyzed the downstream target genes and regulatory mechanisms of AD-related candidate mutations. For example... Figure 17 As shown, a total of 10 representative AD candidate mutation-gene pairs were identified.

[0144] The mutant rs10821626 is located in a microglia-specific open chromatin region. MotifBreakR analysis revealed that this mutation disrupts the binding motif of the transcription factor SPIC. Further Cicero co-accessibility analysis showed that the SPIC binding site has a significant spatial interaction with the promoter of the AD risk gene CCDC6 (e.g., Figure 18 As shown in the figure, this mutation may reduce the transcriptional activation of CCDC6 by disrupting SPIC binding.

[0145] This step, combining multidimensional validation of mutation sites, regulatory elements, and target genes, reveals the molecular regulatory pathways of AD-related mutations, providing clear targets for subsequent functional validation.

[0146] 3. Simulated saturation mutation analysis Simulated saturation mutation is an interpretability analysis method based on deep learning models, used to quantify the importance of each nucleotide in a DNA sequence to the model's prediction results. The specific steps are as follows: (1) Reference sequence construction and mutation simulation A 2,000 bp reference sequence was extracted centered on the target mutation site. All possible single-base substitutions were performed at each nucleotide position, generating 3 × 2,000 mutant sequences. All mutant sequences were input into the regulatory module to obtain a single-cell resolution chromatin accessibility prediction matrix.

[0147] (2) Score standardization and locus contribution calculation The prediction results for each cell were Z-score normalized, and the site contribution score of each reference nucleotide was calculated to assess the strength of the influence of the base on local chromatin accessibility.

[0148] (3) Site-specific analysis For the mutation rs10821626, saturation mutation simulation was performed within ±50 bp upstream and downstream of it.

[0149] The results showed that the G allele disrupted the SPIC binding motif, thereby reducing chromatin accessibility at that site. Figure 19 As shown, the contribution score of the bases in the SPIC binding motif in microglia is significantly higher than that in the other five brain cell types, indicating that this motif is specific and functionally critical in microglia.

[0150] 4. Allele-specific open chromatin analysis The inventors downloaded ATAC-seq alignment data (BAM format) of microglial cells differentiated from human induced pluripotent stem cells from the GEO database (GSE164315, GSE206479, GSE226690, GSE245524). After alignment using the hg38 version of the reference genome, allele imbalance analysis was performed using the WASP2 workflow.

[0151] To ensure statistical reliability, only samples with ≥5 valid alignments at the target mutation site (e.g., rs10821626) were retained. A one-sided t-test was used to assess the bias of the two allele signals. Results are as follows: Figure 20 As shown, the A→G mutation at rs10821626 significantly reduced chromatin accessibility, further supporting the repressive regulatory effect of the G allele.

[0152] 5. Base editing experiments to verify AD-related mutations Base editing experiments are a crucial step in functionally validating model predictions, used to verify whether single-base mutations at specific sites will induce changes in gene expression or epigenetic regulation in real biological systems. This experiment, based on the CRISPR / Cas9 system, uses a ribonucleoprotein complex (RNP) to introduce the Cas9 protein and specific sgRNA into human cells, thereby achieving precise single-base substitutions at specific sites in the genome.

[0153] Compared with traditional double-strand break repair methods, base editing methods have advantages such as high editing efficiency, low off-target effects, and low cytotoxicity, making them suitable for verifying the functional effects of non-coding mutations. In this embodiment, the inventors introduced a target mutation (such as the A→G substitution in rs10821626) into 293FT cells and, combined with subsequent genome sequencing and quantitative PCR (qPCR) analysis, were able to directly observe the transcriptional regulatory changes induced by the mutation, thereby experimentally validating the prediction results of the deep learning model. The specific steps are as follows: (1) Cell culture 293FT cells (Thermo Fisher Scientific, R70007) were selected and cultured routinely at 37 °C and 5% CO2 in DMEM medium (Gibco, 10566) containing 10% fetal bovine serum (Gibco, 16000-044) and 1% penicillin and streptomycin (Invitrogen, 15140122).

[0154] (2) CRISPR-Cas9-mediated mutation editing Before transfection, 0.8 × 10 5 Up to 1.0×10 5 Two 293FT cells were seeded in 24-well plates and cultured for 24 hours. Subsequently, 36 pmol of Cas9 protein (Genscript, Z03702), 36 pmol of sgRNA (GAUAUGAGCUUUGAGAAAAA), and 0.9 nmol of Oligo ssDNA template were co-transfected using Lipofectamine™ 2000 (Thermo Fisher Scientific, 11668019).

[0155] (3) DNA extraction and sequencing 72 hours after transfection, genomic DNA was extracted using QuickExtract™ DNA Extraction Solution (Epicentre, QE09050) and sequenced to confirm the successful editing of the A→G mutation.

[0156] (4) Target gene expression analysis The edited cells were seeded in 60 mm culture dishes (density 1×10⁶). 5 Total RNA was extracted from each well ( / well), treated with DNase, and then subjected to reverse transcription PCR to obtain cDNA. The expression level of the target gene CCDC6 was detected by qPCR using the Vazyme Q711-02AA commercial kit. Figure 21 As shown, the A→G mutation led to a significant decrease in CCDC6 expression, and the results were consistent when using different internal reference genes (GAPDH and ACTB).

[0157] Furthermore, it should be understood that after reading the foregoing teachings of this application, those skilled in the art can make various alterations or modifications to this application, and these equivalent forms also fall within the scope defined by the appended claims.

Claims

1. A method for constructing a deep learning model for predicting the function of non-coding mutations based on single-cell multi-omics data, characterized in that, The single-cell multi-omics data includes an epigenome dataset and a single-cell chromatin accessibility dataset. The method includes the step of constructing a regulatory module for predicting the chromatin accessibility probability of sequence fragments. The regulatory module consists of a DeepSEA pre-trained branch and a ResNet18 end-to-end learning branch, and the two are integrated through a feature fusion module. DeepSEA pre-training branch structure: pre-training is performed based on the epigenome dataset, and parameters are further optimized based on the single-cell chromatin accessibility dataset to obtain the general epigenomic features of sequence fragments; ResNet18 end-to-end learning branch structure: Extract single-cell specific features based on the single-cell chromatin accessibility dataset; Feature fusion module: This module fuses the general epigenomic features output from the DeepSEA pre-training branch and the single-cell-specific features generated by the ResNet18 end-to-end learning branch. The fusion result is then non-linearly mapped through three fully connected layers, ultimately outputting the predicted chromatin accessibility probability value for each input sequence fragment at the single-cell level. The single-cell multi-omics data also includes a single-cell transcriptome dataset, and the method further includes the step of constructing an expression module for predicting the expression level of genes corresponding to sequence fragments: Based on the aforementioned regulatory module, the probability of chromatin accessibility in the region near the gene transcription start site is predicted, and the expression level of the corresponding gene for the sequence fragment is further predicted based on the single-cell transcriptome dataset.

2. A method for predicting the functional effect of non-coding mutations, characterized in that, The deep learning model constructed using the method described in claim 1 predicts the impact of non-coding mutations on chromatin accessibility and the impact of non-coding mutations on gene expression. Specifically, The regulatory module is used to predict the chromatin accessibility probability of wild-type and mutant sequences containing the non-coding mutation site, respectively. Based on the difference in chromatin accessibility probability between the two, the non-coding mutation regulatory effect score is determined, and the impact of non-coding mutation on chromatin accessibility is further determined. The expression module is used to predict the expression levels of the corresponding genes of wild-type and mutant sequences containing the non-coding mutation site, respectively. The difference in expression levels between the two corresponding genes is used to determine the effect of the non-coding mutation on gene expression.

3. A method for identifying cell type-specific regulatory mutations, characterized in that, For different mutation sites, the regulatory effect score of each mutation in different cell types is determined using the method described in claim 2, and cell type-specific regulatory mutations are obtained based on statistical methods.

4. A computer device, characterized in that, include: Memory, used to store computer programs; A processor for executing the computer program to implement the steps of the method as described in any one of claims 1-3.

5. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a computer program that, when executed by a processor, implements the steps of the method as described in any one of claims 1-3.

Images