Classification method for predicting promoter fragments in scATAC-seq based on transcription factor motifs

Abstract

The present invention discloses a classification method for predicting promoter fragments in scATAC-seq based on transcription factor motifs, which solves the problems of cumbersome and low-efficiency methods for determining promoters through multi-omics integration in the prior art. It includes analyzing multiple omics data in the same tissue, and then screening out the top 2000 highly variable genes through RNA expression analysis. Then, using gene IDs to determine whether these genes have corresponding RNA expression levels, and finding the scATAC-seq sequences at the transcription start sites (TSS) of these genes, while distinguishing the sequences in non-promoter regions, and splitting the sequences in the promoter regions into 100-bp fragments for subsequent feature extraction. Subsequently, motif features are used to extract features from these fragments to establish a prediction model for predicting promoter regions. Finally, the prediction model is evaluated to ensure its good prediction performance and accuracy. The advantages of the present invention are: it can accurately and quickly identify promoter regions by learning multi-omics data.

Description

Technical Field

[0001] The present invention relates to the field of bioinformatics technology, and particularly relates to a classification method for predicting promoter fragments in scATAC-seq based on transcription factor motifs. Background Art

[0002] Promoters, which play a crucial role in directly regulating gene expression levels, are located in the transcriptional start region of genes and are the basis for understanding the molecular mechanisms of gene regulation. In the gene regulatory network, promoters, as key regulatory elements responsible for gene activation, do not operate in isolation. Instead, they are part of a complex regulatory system, intertwined with transcription factors (TFs), enhancers, and other elements, jointly driving target genes to achieve specific expression levels in specific cell types. With the rapid development of high-throughput sequencing technologies, especially the breakthrough of single-cell sequencing technology, it is now possible to analyze the gene expression diversity and dynamic changes at single-cell resolution with unprecedented precision. Single-cell resolution analysis can detect subtle differences in gene expression among individual cells within the same tissue or organ.

[0003] Chromatin accessibility analysis is the basis of epigenetics research, providing key insights into the regulatory mechanisms coordinating gene expression. ATAC-seq can utilize transposase to label accessible regions of chromatin and then reveal key information about cis-regulatory elements through sequencing. Single-cell RNA sequencing (scRNA-seq) is the basis of transcriptomics, enabling the mapping of gene expression patterns at single-cell resolution. Integrating scRNA-seq and scATAC-seq data can reveal the subtle connections between gene expression and chromatin states and accurately evaluate the regulatory differences driven by promoters in various cell types. This comprehensive cross-omics approach deepens the understanding of the complexity and heterogeneity of the gene expression regulatory network. However, it is cumbersome to determine promoters through multi-omics integration because it requires using a large number of R language-dependent packages and performing complex processes.

[0004] In summary, the existing methods require a large amount of manpower and material resources to download and install dependent packages and execute complex processes, resulting in low efficiency. Summary of the Invention

[0005] The object of the present invention is to address the above problems and provide a classification method for predicting promoter fragments in scATAC-seq based on transcription factor motifs.

[0006] To achieve the above object, the present invention adopts the following technical solutions: A classification method for predicting promoter fragments in scATAC-seq based on transcription factor motifs, the method comprising the following steps:

[0007] S1. Obtain multi-omics data of the same tissue and compare them to ensure that the number of cells in the two omics data is consistent;

[0008] S2. Determine the top 2000 highly variable genes through the gene expression level data of RNA;

[0009] S3. Correlate the scATAC-seq data with the transcription start sites (TSS) of genes corresponding to the RNA expression levels determined by identifiers. If there is a correspondence, it is classified as a promoter; if there is no correspondence, it is classified as a non-promoter;

[0010] S4. Split the promoter and non-promoter data into sequence lengths of 100 bp, and divide them into a training set, a test set, and an independent set;

[0011] S5. Extract features from the sequences using motifs to construct a feature matrix;

[0012] S6. Based on the extracted features, use a CNN model to construct a prediction model for promoters;

[0013] S7. Evaluate the performance of the established prediction model.

[0014] In step S1, the multi-omics data are scRNA-seq and scATAC-seq data of PBMC from the 10x Genomics database. The common cells are identified by comparing the cell barcodes in the two datasets, and the number of cells in each dataset is counted to ensure the consistency of the two datasets.

[0015] In step S2, determining the top 2000 highly variable genes specifically includes the following steps:

[0016] S21. Obtain the scRNA-seq dataset after quality control, standardization, and normalization;

[0017] S22. Calculate the standard deviation of the expression level of each gene in all cells to measure the variability of gene expression;

[0018] S23. Use the function in the Seurat software package to sort the genes according to the standard deviation of gene expression, and set parameters to select the top 2000 genes as highly variable genes.

[0019] In step S3, the classification of promoters and non-promoters specifically includes the following steps:

[0020] S31. Extract the cell identifiers in the scATAC-seq data and match them with the cell identifiers in the scRNA-seq data to ensure that the different omics data of the same cell are analyzed;

[0021] S32. For each matching cell identifier, determine the transcription start site (TSS) of the gene in the scRNA-seq data;

[0022] S33. Align each fragment in the scATAC-seq data with the gene TSS corresponding to the RNA expression level determined by the cell identifier;

[0023] S34. If a fragment in the scATAC-seq data has a corresponding relationship with the TSS of a certain gene, divide the fragment into the promoter region;

[0024] S35. If a fragment in the scATAC-seq data has no corresponding relationship with the TSS of any gene, divide the fragment into the non-promoter region.

[0025] In the above classification method for predicting promoter fragments in scATAC-seq based on transcription factor motifs, step S4 specifically includes the following steps.

[0026] S41. Split the promoter and non-promoter fragments and divide them into a training set, a test set, and an independent set;

[0027] S42. Extract sequence fragments of the promoter region and the non-promoter region from the scATAC-seq data.

[0028] S43. Uniformly divide the sequence of each promoter and non-promoter region into several 100bp subsequences.

[0029] In step S43, for sequences with a length less than 100bp, use the zero-padding strategy for filling to ensure that all sequences have the same length and format in subsequent analysis.

[0030] In the above classification method for predicting promoter fragments in scATAC-seq based on transcription factor motifs, the number of samples of promoters and non-promoters in the training set is 6000 each, and the distribution numbers of promoters and non-promoters in the test set are 2000 each, and the distribution numbers of promoters and non-promoters in the independent set are 2000 each.

[0031] In step S5, use the motif to extract features from the sequence. The specific steps include:

[0032] S51. Obtain the position frequency matrix (PFM) of the motif with a length of 12 for 600 transcription factors in the JASPAR database;

[0033] S52. Convert the sequence fragment into a One-Hot encoding matrix, where each base corresponds to a unique encoding vector;

[0034] S53. Take the position frequency matrix PFM as a sliding window and slide it position by position onto the One-Hot encoding matrix;

[0035] S54. At each position, extract the corresponding frequency value in the position frequency matrix PFM according to the position of 1 in the One-Hot encoding as the feature value of this position;

[0036] S55. Through the sliding integration operation, generate a new feature matrix for each sequence segment;

[0037] In step S6, use the CNN model to construct a promoter prediction model, which specifically includes the following steps:

[0038] S61. According to the feature matrix extracted in step S5, determine the input dimension of the prediction model as (600, 89, 1);

[0039] S62. Design the initial layer of the prediction model as a two-dimensional convolutional layer with 32 filters of size 7×7, and use batch normalization and the LeakyReLU activation function, and retain the output result;

[0040] S63. Add another convolutional layer with 64 filters of size 3×3 to the prediction model, use batch normalization and the LeakyReLU activation function, and integrate the output of the initial convolutional layer into this layer through a residual connection before applying the LeakyReLU;

[0041] S64. Use an average pooling layer of size 2×2 to reduce the dimension of the model feature map, convert the obtained feature map into a one-dimensional vector, and add a dropout layer;

[0042] S65. Add a fully connected layer with 256 units to the prediction model, and use batch normalization and the LeakyReLU activation function;

[0043] S66. The output layer of the prediction model is a neuron using the sigmoid activation function to achieve binary classification, use the Adam optimizer, and use binary cross-entropy as the loss function.

[0044] In step S7, evaluate the performance of the prediction promoter model constructed by the CNN model.

[0045] Compared with the prior art, the advantages of the present invention are as follows: it can more accurately predict the promoter region, thereby improving the reliability of gene expression regulation research; it can quickly process a large amount of data, shorten the time required for promoter prediction, improve research efficiency, reduce manual operations and experimental costs, and provide an efficient, accurate and practical tool for bioinformatics research and gene regulatory network analysis, which has important scientific and economic value for promoting biological research and clinical applications. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] Figure 1 It is the overall flow chart for predicting promoters in the present invention;

[0047] Figure 2 It is the calculated feature matrix diagram in the present invention;

[0048] Figure 3 It is the structural diagram of the prediction model in the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0049] The present invention will be further described in detail below with reference to the drawings and specific embodiments.

[0050] As Figures 1-3 shown, a classification method for predicting promoter fragments in scATAC-seq based on transcription factor motifs, the method includes the following steps:

[0051] S1. Obtain multi-omics data of the same tissue and compare them to ensure that the number of cells in the two omics data is the same;

[0052] In the embodiment of the present invention, the data is scRNA-seq and scATAC-seq data of PBMC from the 10x Genomics database. The common cells are identified by comparing the cell barcodes in the two datasets, and the number of cells in each dataset is counted to ensure the consistency of the two datasets.

[0053] S2. Determine the top 2000 highly variable genes through the gene expression level data of RNA;

[0054] In the implementation process of the present invention, an scRNA-seq dataset after quality control, standardization and normalization is obtained; the standard deviation of the expression level of each gene in all cells is calculated to measure the variability of gene expression; functions in the Seurat software package are used to sort the genes according to the standard deviation of gene expression, and parameters are set to select the top 2000 genes as highly variable genes.

[0055] S3. Correlate the scATAC-seq data with the transcription start sites (TSS) of genes corresponding to the RNA expression levels determined by the identifiers. If there is a correlation, it is classified as a promoter; if there is no correlation, it is classified as a non-promoter.

[0056] In the implementation process of the present invention, cell identifiers in the scATAC-seq data are extracted and matched with the cell identifiers in the scRNA-seq data to ensure that different omics data of the same cell are analyzed.

[0057] For each matched cell identifier, determine the transcription start site (TSS) of genes in the scRNA-seq data, and compare each fragment in the scATAC-seq data with the gene TSS corresponding to the RNA expression level determined by the cell identifier.

[0058] If a fragment in the scATAC-seq data has a corresponding relationship with the TSS of a certain gene, the fragment is classified as a promoter region; if a fragment in the scATAC-seq data has no corresponding relationship with the TSS of any gene, the fragment is classified as a non-promoter region.

[0059] S4. Split the promoter and non-promoter data into sequence lengths of 100 bp, and divide them into a training set, a test set, and an independent set.

[0060] In the implementation process of the present invention, sequence fragments of the promoter region and the non-promoter region are extracted from the scATAC-seq data, and the sequences of each promoter and non-promoter region are evenly split into multiple 100-bp subsequences.

[0061] For sequences with a length less than 100 bp, a strategy of padding with 0 is adopted to ensure that all sequences have the same length and format in subsequent analysis.

[0062] The number of samples of promoters and non-promoters in the training set is 6000 each, the distribution of promoters and non-promoters in the test set is 2000 each, and the distribution of promoters and non-promoters in the independent set is 2000 each.

[0063] S5. Extract features from the sequences using motifs and construct a feature matrix.

[0064] In the implementation process of the present invention, the position frequency matrix PFM of motifs with a length of 12 of 600 transcription factors was obtained from the JASPAR database. The sequence fragments were converted into a One-Hot encoding matrix, where each base corresponded to a unique encoding vector. The PFM matrix was used as a sliding window and slid over the One-Hot encoding matrix position by position. At each position, the corresponding frequency value in the PFM matrix was extracted according to the position of 1 in the One-Hot encoding as the eigenvalue of that position. Through the sliding integration operation, a new 600*89 feature matrix was generated for each sequence fragment.

[0065] S6. Based on the extracted features, a prediction model for promoters was constructed using a CNN model;

[0066] In the implementation process of the present invention, according to the features extracted in step S5, the input dimension of the prediction model was determined to be (600, 89, 1). The initial layer of the model was designed as a two-dimensional convolutional layer with 32 filters of size 7×7, and batch normalization and the LeakyReLU activation function were used, and the output result was retained. Another convolutional layer with 64 filters of size 3×3 was added to the prediction model, and batch normalization and the LeakyReLU activation function were used. Before applying LeakyReLU, the output of the initial convolutional layer was integrated into this layer through a residual connection. An average pooling layer of size 2×2 was used to reduce the dimension of the model feature map, the obtained feature map was converted into a one-dimensional vector, a dropout layer with a dropout rate of 0.3 was added, a fully connected layer with 256 units was added, and batch normalization and the LeakyReLU activation function were used. The output layer of the prediction model was a neuron using the sigmoid activation function to achieve binary classification. The Adam optimizer with a learning rate of 0.01 was used, and binary cross-entropy was used as the loss function.

[0067] S7. Evaluate the performance of the established prediction model.

[0068] Finally, the performance of the prediction model was evaluated through four evaluation metrics, and the specific definitions are as follows:

[0069] ,

[0070] In the implementation process of the present invention, a comparison was also made with the existing prediction model. It can be found from Table 1 that in the independent set, the evaluation parameters Sn, Sp, ACC, and MCC of the promoter prediction model all exceeded the results of the existing prediction model.

[0071] Table 1

[0072]

[0073] In summary, the principle of this embodiment is as follows: First, by analyzing multiple omics data in the same tissue, it is ensured that the number of cells under study remains consistent across various omics. Then, through RNA expression analysis, the top 2000 highly variable genes are screened out. Next, using gene IDs, it is determined whether these genes have corresponding RNA expression levels, and the scATAC-seq sequences at the transcription start sites (TSSs) of these genes, i.e., the promoter regions, are found. At the same time, the sequences in the non-promoter regions are distinguished, and the sequences in the promoter regions are segmented into 100-bp fragments for subsequent feature extraction. Subsequently, motif features are used to extract features from these fragments, a prediction model is established to predict the promoter regions, and finally, the prediction model is evaluated to ensure its good prediction performance and accuracy.

[0074] The specific embodiments described herein are merely illustrative of the spirit of the present invention. Those skilled in the art to which the present invention pertains can make various modifications or supplements to the described specific embodiments or use similar means for substitution, but will not deviate from the spirit of the present invention or exceed the scope defined by the appended claims.

Claims

1. A classification method for predicting promoter fragments in scATAC-seq based on transcription factor motifs, characterized in that This method includes the following steps: S1. Obtain multi-omics data of the same tissue and compare them to ensure that the number of cells in the two omics data is the same; S2. Determine the top 2000 highly variable genes based on the gene expression level data of RNA; S3. Correlate the scATAC-seq data with the transcription start site (TSS) of the gene corresponding to the RNA expression level determined by the identifier. If there is a correspondence, it is classified as a promoter; if there is no correspondence, it is classified as a non-promoter; S4. Split the promoter and non-promoter data into sequence lengths of 100 bp, and divide them into a training set, a test set, and an independent set; S5. Extract features from the sequences using motifs to construct a feature matrix; S6. Based on the extracted features, use a CNN model to construct a prediction model for promoters; S7. Evaluate the performance of the established prediction model.

2. The classification method for predicting promoter fragments in scATAC-seq based on transcription factor motifs according to claim 1, wherein In step S1, the multi-omics data are the scRNA-seq and scATAC-seq data of PBMC from the 10x Genomics database. By comparing the cell barcodes in the two datasets, the common cells are identified, and the number of cells in each dataset is counted to ensure the consistency of the two datasets.

3. The classification method for predicting promoter fragments in scATAC-seq based on transcription factor motifs according to claim 2, wherein In step S2, determining the top 2000 highly variable genes specifically includes the following steps: S21. Obtain the scRNA-seq dataset after quality control, standardization, and normalization; S22. Calculate the standard deviation of the expression level of each gene in all cells to measure the variability of gene expression; S23. Use the function in the Seurat software package to sort the genes according to the standard deviation of gene expression, and set parameters to select the top 2000 genes as highly variable genes.

4. The classification method for predicting promoter fragments in scATAC-seq based on transcription factor motifs according to claim 1, wherein In step S3, the classification of promoters and non-promoters specifically includes the following steps: S31. Extract the cell identifiers in the scATAC-seq data and match them with the cell identifiers in the scRNA-seq data to ensure that different omics data of the same cell are analyzed; S32. For each matched cell identifier, determine the transcription start site (TSS) of the gene in the scRNA-seq data; S33. Align each fragment in the scATAC-seq data with the TSS of the gene corresponding to the RNA expression level determined by the cell identifier; S34. If the fragment in the scATAC-seq data has a correspondence with the TSS of a certain gene, divide this fragment into the promoter region; S35. If the fragment in the scATAC-seq data has no correspondence with the TSS of any gene, divide this fragment into the non-promoter region.

5. The classification method for predicting promoter fragments in scATAC-seq based on transcription factor motifs according to claim 1, wherein Step S4 specifically includes the following steps, S41. Split the promoter and non-promoter fragments and divide them into a training set, a test set, and an independent set; S42. Extract the sequence fragments of the promoter region and non-promoter region from the scATAC-seq data; S43. Uniformly split the sequences of each promoter and non-promoter region into several 100-bp subsequences.

6. The classification method for predicting promoter fragments in scATAC-seq based on transcription factor motifs according to claim 5, wherein In step S43, for sequences with a length less than 100 bp, a zero-padding strategy is adopted for padding to ensure that all sequences have the same length and format in subsequent analysis.

7. The classification method for predicting promoter fragments in scATAC-seq based on transcription factor motifs according to claim 6, characterized in that, The number of samples of the promoter and non-promoter in the training set is 6,000 each, and the distribution numbers of the promoter and non-promoter in the test set are 2,000 each, and the distribution numbers of the promoter and non-promoter in the independent set are 2,000 each.

8. The classification method for predicting promoter fragments in scATAC-seq based on transcription factor motifs according to claim 5, wherein In step S5, motifs are used to extract features from sequences. The specific steps include: S51. Obtain the position frequency matrix PFM of motifs with a length of 12 of 600 transcription factors in the JASPAR database; S52. Convert the sequence fragment into a One-Hot encoding matrix, where each base corresponds to a unique encoding vector; S53. Use the position frequency matrix PFM as a sliding window and slide it over the One-Hot encoding matrix position by position; S54. At each position, extract the corresponding frequency value in the position frequency matrix PFM according to the position of 1 in the One-Hot encoding as the feature value at this position; S55. Through a sliding integration operation, generate a new feature matrix for each sequence fragment.

9. The classification method for predicting promoter fragments in scATAC-seq based on transcription factor motifs according to claim 1, wherein In step S6, a prediction model for promoters is constructed using a CNN model, which specifically includes the following steps: S61. According to the feature matrix extracted in step S5, determine the input dimension of the prediction model as (600, 89, 1); S62. Design the initial layer of the prediction model as a two-dimensional convolutional layer with 32 filters of size 7×7, and use batch normalization and the LeakyReLU activation function, and retain the output result; S63. Add another convolutional layer with 64 filters of size 3×3 to the prediction model, use batch normalization and the LeakyReLU activation function, and integrate the output of the initial convolutional layer into this layer through a residual connection before applying LeakyReLU; S64. Use an average pooling layer of size 2×2 to reduce the dimension of the model feature map, convert the obtained feature map into a one-dimensional vector, and add a dropout layer; S65. Add a fully connected layer with 256 units to the prediction model, and use batch normalization and the LeakyReLU activation function; S66. The output layer of the prediction model is a neuron using the sigmoid activation function to achieve binary classification, use the Adam optimizer, and use binary cross-entropy as the loss function.

10. The classification method for predicting promoter fragments in scATAC-seq based on transcription factor motifs according to claim 9, wherein In step S7, evaluate the performance of the prediction promoter model constructed by the CNN model.

Images