Method for identifying elncRNAs regulation target gene in single cell based on single cell ATAC-seq

Abstract

The invention relates to the technical field of bioinformatics, in particular to a method for identifying an elncRNAs regulation target gene in a single cell based on single cell ATAC-seq. According to the method provided by the invention, a regulatory network of an elncRNAs-target gene on a single cell whole genome is identified by utilizing chromatin co-accessibility based on scATAC-seq data. Meanwhile, the method is applied to immune cell differentiation, and the elncRNAs-target gene regulation dynamic state of immune cell differentiation is accurately recognized. In addition, the biological functions of the elncRNAs in different cell types are further analyzed. The method provided by the invention can help related researchers in the field of biomedicine to more elaborately and accurately predict the regulation effect and function of the elncRNAs under the single cell resolution.

Description

Technical Field

[0001] The present invention relates to the technical field of bioinformatics, and in particular to a method for identifying elncRNAs regulated target genes in single cells based on single-cell ATAC-seq. Background Art

[0002] The proposal of the human genome project and the development of high-throughput sequencing technology have prompted researchers to discover a new type of non-coding RNA, long non-coding RNAs (lncRNAs), which are a class of RNA transcript molecules with a length of more than 200 bases (nucleotides, nt) and usually do not encode proteins, which are widely present in both animals and plants. Enhancer lncRNAs (elncRNAs) are a class of lncRNAs transcribed from enhancer regions. They can promote chromatin looping between enhancers and promoters to enhance gene expression, and the chromatin openness of elncRNAs has an important influence on their transcriptional regulation. More and more studies have shown that elncRNAs play a key role in various cell biological processes, including cell differentiation, immune response, and disease progression. However, early studies on lncRNAs mainly focused on their participation in the competitive endogenous RNA (ceRNA) regulatory network, while studies on the regulatory mechanism of elncRNAs are still less. In addition, previous studies mainly used bulk data to perform co-expression analysis of lncRNAs and genes to identify lncRNAs-target gene regulatory relationships. Although co-expression analysis provides valuable insights into potential functional associations, it often ignores the importance of the spatial distance between lncRNAs and target genes. Summary of the invention

[0003] In order to solve the above problems, the present invention provides a method for identifying elncRNAs regulated target genes in single cells based on single cell ATAC-seq. The present invention identifies the elncRNAs-target gene regulatory network at the single cell level based on single cell chromatin accessibility high-throughput sequencing (scATAC-seq) of transposase, reveals the regulatory role and function of elncRNAs in single cells, and provides a new reference method and application tool for the study of gene regulation mechanism in cells from the perspective of non-coding RNA regulation.

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

[0005] The present invention provides a method for identifying elncRNAs regulated target genes in single cells based on single cell ATAC-seq, comprising the following steps:

[0006] Step 1: Collect single-cell elncRNAs-related data, including: lncRNAs, protein-coding gene annotation information, enhancer annotation information and scATAC-seq data;

[0007] Step 2: Map the data related to single-cell elncRNAs in step 1 through genomic position to identify the elncRNAs with open chromatin in single cells, including:

[0008] 1) Extract the transcripts of lncRNAs and the annotation information of protein-coding genes;

[0009] 2) The ±2 kb window of the transcription start site of the protein-coding gene was used as the promoter region of the protein-coding gene, and the lncRNAs overlapping with the promoter region of the protein-coding gene were removed to obtain deduplicated lncRNAs;

[0010] 3) extracting active enhancers from the Roadmap database; taking the deduplicated lncRNAs whose transcription start sites are located in the active enhancer regions as primary enhancer lncRNAs;

[0011] 4) remove primary enhancer lncRNAs overlapping with ENCODE hg19 blacklisted genomic regions to obtain secondary enhancer lncRNAs;

[0012] 5) Comparing the secondary enhancer lncRNAs obtained in step 4) with the scATAC-seq data to obtain the elncRNAs with open chromatin in the single cell; the transcription start site of the elncRNAs with open chromatin in the single cell falls within the open peak region of the single cell;

[0013] Step 3: For the single-cell chromatin-accessible elncRNAs, protein-coding genes, and scATAC-seq data obtained in step 2, use Cicero to analyze the single-cell chromatin co-accessibility and identify the single-cell elncRNAs-target gene interactions, including:

[0014] a) Single cells were divided into multiple cell types based on cell markers, and the chromatin co-accessibility elncRNAs-target gene interactions of each cell type were calculated separately; for each human cell type, the scATAC-seq data was first reduced in dimension to obtain the binary UMAP coordinates of each single cell, and then the binary single cell peak matrix M was normalized using the word frequency-inverse document frequency weighting method to obtain the TF-IDF matrix. The calculation formula is as follows:

[0015] TF=t[t(M) / colSums(M)] Formula I;

[0016] IDF=log[1+ncol(M) / rowSums(M)] Formula II;

[0017] TF-IDF=TF×IDF Formula III;

[0018] Where colSums(M) and rowSums(M) are the sums of each column and row of the matrix M, respectively, ncol(M) is the number of columns of the matrix M, t represents the T function, and t(M) is a function related to M (the number of columns);

[0019] b) using the TF-IDF matrix in step a) to perform matrix singular value decomposition for dimensionality reduction, retaining 1 to 50 dimensional principal components, and obtaining a binary single-cell peak matrix;

[0020] c) Generate a cell dataset using the binarized single-cell peak matrix in step b), and then perform detectedGenes() and estimateSizeFactors(); Generate an integrated cell dataset cicero_cds object using the KNN method and the UMAP dimension reduction coordinates of each single cell in each cell type as the object input into the Cicero software to obtain an integrated chromatin accessibility matrix;

[0021] d) The peak-to-peak interactions were calculated using the integrated chromatin accessibility matrix, the Cicero co-accessibility threshold was set to 0.1, and the interactions between peaks within a center distance of 250 kb were retained;

[0022] e) using the protein-coding genes obtained in step 1 and the elncRNAs obtained in step 2 to annotate the interactions between peaks, and obtain single-cell elncRNAs-target gene interaction data; the promoter region of the protein-coding gene overlaps with the peak; the first column of the interaction data frame represents elncRNAs, the second column represents target genes, and the third column represents the Cicero co-accessibility size of the elncRNAs-target gene interaction;

[0023] Step 4: Perform functional enrichment analysis on the single-cell elncRNAs-target gene interaction data obtained in step 3, including: using clusterProfiler4.0 software in R language to perform Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis and gene ontology functional enrichment analysis. Pathway entries with a p-value less than 0.01 are considered to be significantly enriched;

[0024] Step 5: Integrate the elncRNAs chromatin accessibility matrix obtained in steps 1 to 2 with the protein-coding gene activity score matrix calculated from scATAC-seq data, and use the gene activity score to prune the regulatory dynamics of related elncRNAs in the cell differentiation trajectory, including:

[0025] A) ENCODE hg19 blacklisted genomic regions were filtered out and protein-coding genes were annotated using the addGeneScoreMatrix() function in ArchR software to create a gene activity score matrix in scATAC-seq data;

[0026] B) The addTrajectory() function in ArchR software was used to calculate the cell differentiation trajectory. The Euclidean distance between the reduced-dimensional coordinates of each cell and the average coordinates of the next cell type cluster was used to sort the single cells in the 2D subspace.

[0027] C) Use the getTrajectory() function in ArchR software to visualize the activity changes of elncRNAs and protein-coding genes along the cell differentiation trajectory;

[0028] D) Based on the genomic position overlap method, the correlateTrajectories() function in ArchR software was used to perform correlation analysis between the elncRNAs chromatin accessibility matrix and the protein-coding gene activity score matrix within 250 kb, thereby trimming the elncRNAs chromatin accessibility data and retaining the elncRNAs chromatin accessibility data with Pearson correlation > 0.35 and FDR < 0.001 to identify elncRNAs that are more likely to be involved in target gene regulation in cell differentiation trajectories.

[0029] Preferably, the sources of the lncRNAs and protein-coding gene annotation information include: GENCODE database; the sources of the enhancer annotation information include: Roadmap database; the sources of the scATAC-seq data include: scEnhancer and GEO databases.

[0030] Preferably, the active enhancers include Active enhancer 1, Genic enhancer 1, Activeenhancer2 and Genic enhancer2.

[0031] Preferably, the binary single-cell peak matrix is: input 2 to 50 dimensional principal components and perform two-dimensional visualization projection in UMAP.

[0032] Preferably, step 3 e) further comprises: for each cell type, constructing an elncRNAs-target gene regulatory network using elncRNAs-target gene interactions in Cytoscape software.

[0033] Preferably, the normalization parameter of the getTrajectory() function in step 5C) is set to log2Norm="TRUE".

[0034] Preferably, step D) in step 5 comprises:

[0035] For the dynamic regulatory association between a pair of elncRNAs and protein-coding gene candidates, the Pearson correlation coefficient on the pseudo-time vector of cell differentiation is calculated as follows:

[0036]

[0037] Where n is the length of the pseudo-time vector for cell differentiation, X i is the expression value of protein-coding gene X in time period i, X u is the mean expression value of gene X during the pseudo-time of cell differentiation, Y i is the expression value of elncRNAsY in time period i, Y u is the mean expression value of elncRNAsY in the pseudo-time of cell differentiation;

[0038] The correlation p-value was calculated based on the t-test using Equation IV, and the Benjamini-Hochberg correction was further performed using the p.adjust() function in the R language.

[0039] Preferably, the method parameter of the Benjamini-Hochberg correction is set to method="fdr".

[0040] Preferably, the variance quantiles of the activity matrix of elncRNAs and protein-coding genes are calculated to identify the top 75% of the most variable (ie, more variable) elncRNAs and protein-coding genes.

[0041] Beneficial effects:

[0042] The method provided by the present invention uses chromatin co-accessibility based on scATAC-seq data to identify the regulatory network of elncRNAs-target genes on the whole genome of a single cell. At the same time, the method is applied to immune cell differentiation, accurately identifying the elncRNAs-target gene regulatory dynamics of immune cell differentiation. In addition, the biological functions of these elncRNAs in different cell types are further analyzed. The method provided by the present invention can help researchers in the biomedical field to more finely and accurately predict the regulatory effects and functions of elncRNAs at single-cell resolution. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required to be used in the embodiments are briefly introduced below.

[0044] Figure 1 Schematic diagram of the transcriptional regulation mechanism of elncRNAs and the identification of elncRNAs in scATAC-seq;

[0045] Figure 2 The results of constructing the regulatory network of elncRNAs-target genes in different cell types;

[0046] Figure 3 The figure shows the elncRNAs regulatory trajectory of the immune cell lineage; Baso is basophil; NK is mature natural killer cell; ImmNK is immature natural killer cell; N CD8 is initial CD8 T cell; CD8 CM is intermediate memory CD8 T cell; CD8 EM is effector memory CD8 T cell; N CD4 is initial CD4 T cell; M CD4 is memory CD4 T cell; the lower right corner shows the regulatory relationship between elncRNAs and target genes;

[0047] Figure 4 is the pseudo-time projection of the B cell lineage differentiation trajectory; the smooth line is the result of fitting along the pseudo-time;

[0048] Figure 5 Single-cell elncRNAs-target gene regulatory dynamics during human B cell lineage differentiation. DETAILED DESCRIPTION

[0049] The present invention provides a method for identifying elncRNAs regulated target genes in single cells based on single cell ATAC-seq, comprising the following steps:

[0050] Step 1: Collect single-cell elncRNAs-related data, including: lncRNAs, protein-coding gene annotation information, enhancer annotation information and scATAC-seq data;

[0051] Step 2: Map the data related to single-cell elncRNAs in step 1 through genomic position to identify the elncRNAs with open chromatin in single cells, including:

[0052] 1) Extract the transcripts of lncRNAs and the annotation information of protein-coding genes;

[0053] 2) The ±2 kb window of the transcription start site of the protein-coding gene was used as the promoter region of the protein-coding gene, and the lncRNAs overlapping with the promoter region of the protein-coding gene were removed to obtain deduplicated lncRNAs;

[0054] 3) extracting active enhancers from the Roadmap database; taking the deduplicated lncRNAs whose transcription start sites are located in the active enhancer regions as primary enhancer lncRNAs;

[0055] 4) remove primary enhancer lncRNAs overlapping with ENCODE hg19 blacklisted genomic regions to obtain secondary enhancer lncRNAs;

[0056] 5) Comparing the secondary enhancer lncRNAs obtained in step 4) with the scATAC-seq data to obtain the elncRNAs with open chromatin in single cells; the transcription start sites of the elncRNAs with open chromatin in single cells fall within the open peak region of the single cell, so as to identify the elncRNAs with open chromatin in single cells that can be transcribed and play a role;

[0057] Step 3: For the single-cell chromatin-accessible elncRNAs, protein-coding genes, and scATAC-seq data obtained in step 2, use Cicero to analyze the single-cell chromatin co-accessibility and identify the single-cell elncRNAs-target gene interactions, including:

[0058] a) Single cells were divided into multiple cell types based on cell markers, and the chromatin co-accessibility elncRNAs-target gene interactions of each cell type were calculated separately; for each human cell type, the scATAC-seq data was first reduced in dimension to obtain the binary UMAP coordinates of each single cell, and then the binary single cell peak matrix M was normalized using the word frequency-inverse document frequency weighting method to obtain the TF-IDF matrix. The calculation formula is as follows:

[0059] TF=t[t(M) / colSums(M)] Formula I;

[0060] IDF=log[1+ncol(M) / rowSums(M)] Formula II;

[0061] TF-IDF=TF×IDF Formula III;

[0062] Where colSums(M) and rowSums(M) are the sums of each column and row of the matrix M, respectively, ncol(M) is the number of columns of the matrix M, t represents the T function, and t(M) is a function related to M (the number of columns);

[0063] b) using the TF-IDF matrix in step a) to perform matrix singular value decomposition for dimensionality reduction, retaining 1 to 50 dimensional principal components, and obtaining a binary single-cell peak matrix;

[0064] c) Generate a cell dataset using the binarized single-cell peak matrix in step b), and then perform detectedGenes() and estimateSizeFactors(); Generate an integrated cell dataset cicero_cds object using the KNN method and the UMAP dimension reduction coordinates of each single cell in each cell type as the object input into the Cicero software to obtain an integrated chromatin accessibility matrix;

[0065] d) The peak-to-peak interactions were calculated using the integrated chromatin accessibility matrix, the Cicero co-accessibility threshold was set to 0.1, and the interactions between peaks within a center distance of 250 kb were retained;

[0066] e) using the protein-coding genes obtained in step 1 and the elncRNAs obtained in step 2 to annotate the interactions between peaks, and obtain single-cell elncRNAs-target gene interaction data; the promoter region of the protein-coding gene overlaps with the peak; the first column of the interaction data frame represents elncRNAs, the second column represents target genes, and the third column represents the Cicero co-accessibility size of the elncRNAs-target gene interaction;

[0067] Step 4: Perform functional enrichment analysis on the single-cell elncRNAs-target gene interaction data obtained in step 3, including: using clusterProfiler4.0 software in R language to perform Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis and gene ontology functional enrichment analysis. Pathway entries with a p-value less than 0.01 are considered to be significantly enriched;

[0068] Step 5: Integrate the elncRNAs chromatin accessibility matrix obtained in steps 1 to 2 with the protein-coding gene activity score matrix calculated from scATAC-seq data, and use the gene activity score to prune the regulatory dynamics of related elncRNAs in the cell differentiation trajectory, including:

[0069] A) ENCODE hg19 blacklisted genomic regions were filtered out and protein-coding genes were annotated using the addGeneScoreMatrix() function in ArchR software to create a gene activity score matrix in scATAC-seq data;

[0070] B) The addTrajectory() function in ArchR software was used to calculate the cell differentiation trajectory. The Euclidean distance between the reduced-dimensional coordinates of each cell and the average coordinates of the next cell type cluster was used to sort the single cells in the 2D subspace.

[0071] C) Use the getTrajectory() function in ArchR software to visualize the activity changes of elncRNAs and protein-coding genes along the cell differentiation trajectory;

[0072] D) Based on the genomic position overlap method, the correlateTrajectories() function in ArchR software was used to perform correlation analysis between the elncRNAs chromatin accessibility matrix and the protein-coding gene activity score matrix within 250 kb, thereby trimming the elncRNAs chromatin accessibility data and retaining the elncRNAs chromatin accessibility data with Pearson correlation > 0.35 and FDR < 0.001 to identify elncRNAs that are more likely to be involved in target gene regulation in cell differentiation trajectories.

[0073] As a specific embodiment of the present invention, the sources of the lncRNAs and protein-coding gene annotation information include: GENCODE database; the sources of the enhancer annotation information include: Roadmap database; the sources of the scATAC-seq data include: scEnhancer and GEO databases.

[0074] As a specific embodiment of the present invention, the active enhancers include Active enhancer 1, Genic enhancer 1, Active enhancer 2 and Genic enhancer 2.

[0075] As a specific embodiment of the present invention, the binary single-cell peak matrix is: input 2 to 50 dimensional principal components to perform two-dimensional visualization projection in UMAP.

[0076] As a specific embodiment of the present invention, step e) in step 3 further includes: for each cell type, constructing an elncRNAs-target gene regulatory network using elncRNAs-target gene interactions in Cytoscape software.

[0077] As a specific implementation of the present invention, the normalization parameter of the getTrajectory() function in step C) of step 5 is set to log2Norm="TRUE".

[0078] As a specific embodiment of the present invention, step D) in step 5 includes:

[0079] For the dynamic regulatory association between a pair of elncRNAs and protein-coding gene candidates, the Pearson correlation coefficient on the pseudo-time vector of cell differentiation is calculated as follows:

[0080]

[0081] Where n is the length of the pseudo-time vector for cell differentiation, X i is the expression value of protein-coding gene X in time period i, X u is the mean expression value of gene X during the pseudo-time of cell differentiation, Y i is the expression value of elncRNAsY in time period i, Y u is the mean expression value of elncRNAsY in the pseudo-time of cell differentiation;

[0082] The correlation p-value was calculated based on the t-test using Equation IV, and the Benjamini-Hochberg correction was further performed using the p.adjust() function in the R language.

[0083] As a specific implementation of the present invention, the method parameter of the Benjamini-Hochberg correction is set to method = "fdr".

[0084] As a specific embodiment of the present invention, the variance quantiles of the activity matrix of elncRNAs and protein-coding genes are calculated to identify the top 75% of the most variable elncRNAs and protein-coding genes.

[0085] The present invention identifies the elncRNAs-target gene regulatory network at the single-cell level based on scATAC-seq, reveals the regulatory effects and functions of elncRNAs in single cells, and provides new reference methods and application tools for the study of gene regulatory mechanisms in cells from the perspective of non-coding RNA regulation. This is the main content of the technology of the present invention.

[0086] The present invention is mainly divided into three parts:

[0087] Part I Data Collection and Preprocessing:

[0088] Obtain scATAC-seq data of human cells from scEnhancer and GEO data resources. Obtain annotation information of active enhancers (marked by H3K27ac) from the Roadmap database, and obtain lncRNAs annotation information from GENCODE. Obtain elncRNAs with chromatin openness in single cells through genomic position mapping. Integrate biological pathway data such as KEGG to predict the biological functions regulated by elncRNAs, and use them as the source of biological pathway data for the present invention.

[0089] Part II Identification of regulatory relationships between single-cell elncRNAs and target genes:

[0090] ElncRNAs are lncRNAs transcribed from enhancer regions, which can enhance gene expression by promoting chromatin looping between enhancers and promoters. In the present invention, based on scATAC-seq data, the chromatin co-accessibility of elncRNAs and protein-coding genes is used to identify protein-coding genes regulated by elncRNAs. Then, based on the "guilt by association" method, the present invention performs functional enrichment analysis on the target genes regulated by elncRNAs, and then infers the function of elncRNAs in cells.

[0091] Part III Network Platform Construction:

[0092] The present invention has developed a user-friendly web server website, scElncTar, for analyzing the interaction of elncRNAs and elncRNAs-target genes at single-cell resolution. This database uses MySQL (v.5.7) for data management. Among them, the web program is written in Java language (v.1.8.0) server page, and the web service is built with ApacheTomcat web server (v.6.0). In order to develop and visualize data, the present invention uses many Javascript plug-ins, including Datatable (v.1.10.10), jQuery (v.1.10.2) and ECharts (v.5.0). In addition, the present invention uses R language (v.3.6.3) for statistical summary analysis. Most importantly, the present invention has also developed several database analysis tools for users to explore the regulatory mechanism of elncRNAs at the single-cell level. Finally, the scElncTar database can be run on Windows systems, Linux systems and Mac systems using Chrome, Firefox, Safari and Microsoft Edge browsers. In addition, the scElncTar database data resources can be obtained at http: / / bio-bigdata.hrbmu.edu.cn / scElncTar.

[0093] The method provided by the present invention can identify the cell biological functions regulated by elncRNAs in various cell types, and at the same time develop a database platform to provide flexible query and convenient result acquisition, which is of great directional guidance for studying the functions of elncRNAs in cells, analyzing the complex regulatory mechanisms in cells, and revealing the elncRNAs-target gene regulatory dynamics in the cell differentiation process, while providing relevant candidate data for laboratory verification.

[0094] In order to further illustrate the present invention, a method for identifying elncRNAs regulated target genes in single cells based on single-cell ATAC-seq provided by the present invention is described in detail below in conjunction with the accompanying drawings and examples, but they should not be construed as limiting the scope of protection of the present invention.

[0095] Example 1

[0096] A method for identifying enhancer lncRNAs regulated target genes in single cells, comprising the following steps:

[0097] Step 1: Collect data related to single-cell elncRNAs

[0098] 1) lncRNAs and protein-coding gene annotation information source: GENCODE database;

[0099] 2) Enhancer annotation information source: Roadmap database;

[0100] 3) scATAC-seq data source: scEnhancer and GEO database;

[0101] Step 2: elncRNAs are lncRNAs transcribed from enhancer regions (see Figure 1 ). Therefore, the present invention identifies the elncRNAs with open chromatin in single cells by mapping the enhancers, lncRNAs, and scATAC-seq data collected in step 1 to the genome position. The specific method is as follows:

[0102] 1) Extract the annotation information of lncRNAs transcripts and protein-coding genes from the GENCODE database annotation file.

[0103] 2) The transcription start site (TSS) ±2kb window of the protein-coding gene was used as the promoter region of the protein-coding gene, and the lncRNAs overlapping with the promoter region of the protein-coding gene were removed to obtain deduplicated lncRNAs; active enhancers were extracted from the Roadmap database, including Active enhancer1, Genic enhancer 1, Active enhancer 2 and Genic enhancer2.

[0104] 3) Deduplicated lncRNAs whose transcription start sites are located in the active enhancer regions are used as primary enhancer lncRNAs.

[0105] 4) Primary enhancer lncRNAs overlapping with ENCODE hg19 blacklisted genomic regions were removed to identify elncRNAs transcribed from active enhancer regions and obtain secondary enhancer lncRNAs.

[0106] 5) Compare the secondary enhancer lncRNAs obtained in step 4) with the scATAC-seq data to obtain the elncRNAs with open chromatin in single cells; the transcription start sites of the elncRNAs with open chromatin in single cells fall within the "open" peak region of the single cell, so as to identify the elncRNAs in which the chromatin is opened and can be transcribed and function.

[0107] Step 3: elncRNAs transcribed from enhancers work together with enhancers to enhance the expression of nearby protein-coding genes by promoting the looping of enhancers and promoters (see Figure 1). Therefore, the present invention uses Cicero (a software that uses scATAC-seq data to predict cis-regulatory effects on the genome by analyzing chromatin co-accessibility) to identify single-cell elncRNAs-target gene interactions for the elncRNAs, protein-coding genes and scATAC-seq data of single-cell chromatin accessibility obtained in step 2. The Cicero software uses the k-nearest neighbor (KNN) method to integrate the chromatin accessibility data between each group of cells to calculate the co-accessibility interactions, specifically:

[0108] 1) The present invention divides single cells into 532 cell types based on cell markers, and calculates the chromatin co-accessibility elncRNAs-target gene interactions of each cell type. For each human cell type, the present invention first reduces the dimension of the scATAC-seq data to obtain the UMAP coordinates of each single cell in a two-dimensional subspace (i.e., binary), and then uses the term frequency-inverse document frequency (TF-IDF) weighting method to standardize the binary single cell peak matrix M to obtain the TF-IDF matrix, and the calculation formula is as follows:

[0109] TF=t[t(M) / colSums(M)] Formula I;

[0110] IDF=log[1+ncol(M) / rowSums(M)] Formula II;

[0111] TF-IDF=TF×IDF Formula III;

[0112] Among them, colSums(M) and rowSums(M) are the sums of each column and each row of the matrix M, ncol(M) is the number of columns of the matrix M, t represents the T function, and t(M) is a function related to M (the number of columns). This TF-IDF method is a two-step data normalization process that not only corrects for differences in sequencing depth between all cells, but also gives higher weights to rarer peaks among all chromatin accessibility peaks.

[0113] 2) Use the TF-IDF matrix in step 1) to perform matrix singular value decomposition (SVD) for dimensionality reduction, retaining 1 to 50 dimensional principal components. This method of combining TF-IDF and SVD is called latent semantic indexing (LSI) method. The first LSI principal component usually captures technical differences (sequencing depth differences) rather than biological differences, and the sequencing fragment (reads) depth is highly correlated with the first dimension. In this case, the present invention will perform subsequent analysis after removing the first dimension. Finally, the present invention inputs 2 to 50 dimensional principal components for two-dimensional visualization projection in UMAP.

[0114] 3) Generate a cell dataset "cellDataSet" using the binary single-cell peak matrix, and then perform "detectedGenes" (calculate the number of cells in the "cellDataSet" object that can detect expression above the minimum threshold for each feature) and "estimateSizeFactors" (calculate the Size Factor for each cell to normalize the expression level of each cell). Since single-cell ATAC-seq data is very sparse, the present invention uses the KNN (k = 50) method and the UMAP dimensionality reduction coordinates of each single cell in each cell type to generate an integrated cell dataset "cicero_cds" object as an object input into the Cicero software.

[0115] 4) The peak-to-peak interactions were calculated using the integrated chromatin accessibility matrix, and the Cicero co-accessibility threshold was set to 0.1. In addition, if the center distance between the peaks was within 250 kb, the present invention retained the interaction between the peaks to identify the co-accessible interactions within a certain range that the chromatin could form loops.

[0116] 5) Use protein-coding genes (whose promoter regions overlap with the peaks retained in step 4) and elncRNAs (whose transcription start sites fall within the open peak regions) to annotate the interactions between peaks. Thus, the present invention generates an interaction data frame, in which the first column represents elncRNAs, the second column represents target genes, and the third column represents the Cicero co-accessibility size of elncRNAs-target gene interactions. For each cell type, the present invention uses elncRNAs-target gene interactions in Cytoscape software to construct elncRNAs-target gene regulatory networks ( Figure 2). The present invention identified 444,925 pairs of elncRNAs-target gene interaction networks from 532 cell types.

[0117] Step 4: elncRNAs can regulate the corresponding target genes and thus participate in the biological functions of cells. The present invention performs functional enrichment analysis on the single-cell elncRNAs-target gene interaction data obtained in step 3 and elncRNAs-related target genes to infer the functions of elncRNAs in 532 different cell types. The present invention uses clusterProfiler4.0 software in R language to perform Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis and Gene Ontology (GO) functional enrichment analysis. For the results of KEGG and GO functional enrichment analysis, pathway entries with p values ​​less than 0.01 are considered to have significant enrichment.

[0118] Step 5: According to the results of steps 1 to 4 above, the regulatory network of single-cell elncRNAs-target genes is a static regulatory mechanism for a specific cell type. This regulatory relationship may change dynamically during cell differentiation. At the same time, elncRNAs may participate in the dynamic regulation of immune cell differentiation by regulating the expression of protein-coding genes ( Figure 3 ). Therefore, the present invention performs a dynamic analysis of the single-cell elncRNAs-target gene regulatory trajectory. The present invention integrates the elncRNAs chromatin accessibility matrix obtained in steps 1 to 2 with the protein-coding gene activity score matrix calculated from scATAC-seq data, and uses the gene activity score to trim the regulatory dynamics of the relevant elncRNAs in the cell differentiation trajectory. This method enables the present invention to identify elncRNAs that are more likely to participate in the regulation of target genes in the cell differentiation trajectory, as follows:

[0119] 1) The present invention filters out ENCODE hg19 blacklisted genomic regions and uses the addGeneScoreMatrix() function in ArchR software to annotate protein-coding genes and create a gene activity score matrix in scATAC-seq data.

[0120] 2) The present invention uses the addTrajectory() function in ArchR software to calculate the cell differentiation trajectory. In order to create a single-cell pseudo-temporal differentiation trajectory, the present invention calculates the Euclidean distance between the reduced-dimensional coordinates of each cell and the average coordinates of the next cell type cluster, and sorts the single cells in the 2D subspace (i.e., the UMAP reduced-dimensional space) ( Figure 4 ).

[0121] 3) The present invention uses the getTrajectory() function in ArchR software (with the normalization parameter set to log2Norm = "TRUE") to visualize the activity changes of elncRNAs and protein-coding genes along the cell differentiation trajectory. This method sorts single cells along the identified pseudo-temporal vectors and creates pseudo-temporal heat maps for elncRNAs and protein-coding genes using the standardized elncRNAs chromatin accessibility matrix and protein-coding gene activity score matrix, respectively.

[0122] 4) Based on the genomic position overlap method, the present invention uses the correlateTrajectories() function in the ArchR software to perform correlation analysis between the elncRNAs chromatin accessibility matrix and the protein-coding gene activity score matrix within 250 kb, thereby trimming the elncRNAs chromatin accessibility data.

[0123] For the dynamic regulatory association between a pair of elncRNAs and protein-coding gene candidates, the Pearson correlation coefficient on the pseudo-time vector of cell differentiation is calculated as follows:

[0124]

[0125] Where n is the length of the pseudo-time vector for cell differentiation, X i is the expression value of protein-coding gene X in time period i, X u is the mean expression value of gene X during the pseudo-time of cell differentiation, Y i is the expression value of elncRNAsY in time period i, Y u is the mean expression value of elncRNAsY during the pseudo-time of cell differentiation.

[0126] The present invention uses this method to calculate the correlation p value based on the t-test, and further uses the p.adjust() function in the R language to perform Benjamini-Hochberg correction, and the method parameter is set to method = "fdr". At the same time, the present invention also calculates the variance quantiles of the elncRNAs and protein-coding gene activity matrix to identify the top 75% of the most variable elncRNAs and protein-coding genes (i.e., identifying 75% of the elncRNAs and protein-coding genes with greater changes).

[0127] 5) The present invention trimmed the regulatory trajectories of elncRNAs and genes according to their correlation (Pearson correlation>0.35 and FDR<0.001), and identified elncRNAs that are more likely to participate in the regulation of target genes in the cell differentiation trajectory ( Figure 5 ).

[0128] Although the above embodiment describes the present invention in detail, it is only a part of the embodiments of the present invention, not all of the embodiments. People can also obtain other embodiments based on this embodiment without creativity, and these embodiments all fall within the protection scope of the present invention.

Claims

1. A method for identifying elncRNAs regulated target genes in single cells based on single-cell ATAC-seq, characterized in that: The following steps are involved: Step 1: Collect single-cell elncRNAs-related data, including: lncRNAs, protein-coding gene annotation information, enhancer annotation information and scATAC-seq data; Step 2: Map the data related to single-cell elncRNAs in step 1 through genomic position to identify the elncRNAs with open chromatin in single cells, including: 1) Extract the transcripts of lncRNAs and the annotation information of protein-coding genes; 2) The ±2 kb window of the transcription start site of the protein-coding gene was used as the promoter region of the protein-coding gene, and the lncRNAs overlapping with the promoter region of the protein-coding gene were removed to obtain deduplicated lncRNAs; 3) extracting active enhancers from the Roadmap database; taking the deduplicated lncRNAs whose transcription start sites are located in the active enhancer regions as primary enhancer lncRNAs; 4) remove primary enhancer lncRNAs overlapping with ENCODE hg19 blacklisted genomic regions to obtain secondary enhancer lncRNAs; 5) comparing the secondary enhancer lncRNAs obtained in step 4) with the scATAC-seq data to obtain the elncRNAs with open chromatin in the single cell; the transcription start site of the elncRNAs with open chromatin in the single cell falls within the open peak region of the single cell; Step 3: For the single-cell chromatin-accessible elncRNAs, protein-coding genes, and scATAC-seq data obtained in step 2, use Cicero to analyze the single-cell chromatin co-accessibility and identify the single-cell elncRNAs-target gene interactions, including: a) Single cells were divided into multiple cell types based on cell markers, and the chromatin co-accessibility elncRNAs-target gene interactions of each cell type were calculated separately; for each human cell type, the scATAC-seq data was first reduced in dimension to obtain the binary UMAP coordinates of each single cell, and then the binary single cell peak matrix M was normalized using the word frequency-inverse document frequency weighting method to obtain the TF-IDF matrix. The calculation formula is as follows: TF=t[t(M) / colSums(M)] Formula I; IDF=log[1+ncol(M) / rowSums(M)] Formula II; TF-IDF=TF×IDF Formula III; Where colSums(M) and rowSums(M) are the sums of each column and row of the matrix M, respectively, ncol(M) is the number of columns of the matrix M, t represents the T function, and t(M) is a function related to M (the number of columns); b) using the TF-IDF matrix in step a) to perform matrix singular value decomposition for dimensionality reduction, retaining 1 to 50 dimensional principal components, and obtaining a binary single-cell peak matrix; c) Generate a cell dataset using the binarized single-cell peak matrix in step b), and then perform detectedGenes() and estimateSizeFactors(); Generate an integrated cell dataset cicero_cds object using the KNN method and the UMAP dimension reduction coordinates of each single cell in each cell type as the object input into the Cicero software to obtain an integrated chromatin accessibility matrix; d) The peak-to-peak interactions were calculated using the integrated chromatin accessibility matrix, the Cicero co-accessibility threshold was set to 0.1, and the interactions between peaks within a center distance of 250 kb were retained; e) using the protein-coding genes obtained in step 1 and the elncRNAs obtained in step 2 to annotate the interactions between peaks, and obtain single-cell elncRNAs-target gene interaction data; the promoter region of the protein-coding gene overlaps with the peak; the first column of the interaction data frame represents elncRNAs, the second column represents target genes, and the third column represents the Cicero co-accessibility size of the elncRNAs-target gene interaction; Step 4: Perform functional enrichment analysis on the single-cell elncRNAs-target gene interaction data obtained in step 3, including: using clusterProfiler4.0 software in R language to perform Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis and gene ontology functional enrichment analysis. Pathway entries with a p-value less than 0.01 are considered to be significantly enriched; Step 5: Integrate the elncRNAs chromatin accessibility matrix obtained in steps 1 to 2 with the protein-coding gene activity score matrix calculated from scATAC-seq data, and use the gene activity score to prune the regulatory dynamics of related elncRNAs in the cell differentiation trajectory, including: A) ENCODE hg19 blacklisted genomic regions were filtered out and protein-coding genes were annotated using the addGeneScoreMatrix() function in ArchR software to create a gene activity score matrix in scATAC-seq data; B) The addTrajectory() function in ArchR software was used to calculate the cell differentiation trajectory. The Euclidean distance between the reduced-dimensional coordinates of each cell and the average coordinates of the next cell type cluster was used to sort the single cells in the 2D subspace. C) Use the getTrajectory() function in ArchR software to visualize the activity changes of elncRNAs and protein-coding genes along the cell differentiation trajectory; D) Based on the genomic position overlap method, the correlateTrajectories() function in ArchR software was used to perform correlation analysis between the elncRNAs chromatin accessibility matrix and the protein-coding gene activity score matrix within 250 kb, thereby trimming the elncRNAs chromatin accessibility data and retaining the elncRNAs chromatin accessibility data with Pearson correlation > 0.35 and FDR < 0.001 to identify elncRNAs that are more likely to be involved in target gene regulation in cell differentiation trajectories.

2. The method according to claim 1, characterized in that: The sources of the lncRNAs and protein-coding gene annotation information include: GENCODE database; the sources of the enhancer annotation information include: Roadmap database; the sources of the scATAC-seq data include: scEnhancer and GEO databases.

3. The method according to claim 1, characterized in that The active enhancers include Active enhancer 1, Genic enhancer 1, Active enhancer2 and Genic enhancer2.

4. The method according to claim 1, characterized in that: The binary single-cell peak matrix is ​​as follows: 2 to 50 dimensional principal components are input and two-dimensional visualization projection is performed in UMAP.

5. The method according to claim 1, characterized in that Step 3 e) also includes: for each cell type, constructing an elncRNAs-target gene regulatory network using elncRNAs-target gene interactions in Cytoscape software.

6. The method according to claim 1, characterized in that The normalization parameter of the getTrajectory() function in step 5 (C) is set to log2Norm = "TRUE".

7. The method according to claim 1, characterized in that Step 5 D) includes: For the dynamic regulatory association between a pair of elncRNAs and protein-coding gene candidates, the Pearson correlation coefficient on the pseudo-time vector of cell differentiation is calculated as follows: Where n is the length of the pseudo-time vector for cell differentiation, X i is the expression value of protein-coding gene X in time period i, X u is the mean expression value of gene X during the pseudo-time of cell differentiation, Y i is the expression value of elncRNAsY in time period i, Y u is the mean expression value of elncRNAsY in the pseudo-time of cell differentiation; The correlation p-value was calculated based on the t-test using Equation IV, and the Benjamini-Hochberg correction was further performed using the p.adjust() function in the R language.

8. The method according to claim 7, characterized in that The method parameters of the Benjamini-Hochberg correction are set to method="fdr".

9. The method according to claim 7, characterized in that: The variance quantiles of the activity matrices of elncRNAs and protein-coding genes were calculated to identify the top 75% of the most variable elncRNAs and protein-coding genes.

Images