Methods for eRNA identification, regulatory target prediction and functional annotation based on high-throughput transcriptome sequencing data

Abstract

The application discloses a method for eRNA transcription identification, regulation target prediction and function annotation based on high-throughput transcriptome sequencing data, and is characterized in that the method comprises the following steps: identifying part of non-coding RNA in which a transcription start site is located in an enhancer region as eRNA; obtaining eRNA-related protein coding genes which simultaneously exist in an eRNA-protein coding gene co-expression network and an eRNA-centered regulation network, constructing an eRNA-protein coding gene relationship network, and extracting protein coding genes which are directly connected or closely connected to the eRNA, so as to predict potential regulation targets of the eRNA; and finally, performing function enrichment analysis on the potential regulation targets of the eRNA, so as to obtain the results of eRNA function annotation, and the method has the advantages of wider application range, and can be applied to all eRNAs and more accurately obtain the action forms between the eRNA and the protein coding gene.

Description

TECHNICAL FIELD

[0001] The present application relates to the fields of molecular biology and bioinformatics, and in particular to a method for eRNA identification, regulatory target prediction and functional annotation based on high-throughput transcriptome sequencing data. BACKGROUND

[0002] Enhancer is a cis-regulatory element in the genome for regulating the expression of target genes. Recent studies have shown that the activated enhancer region can transcribe some non-coding RNAs, which are called enhancer RNA (eRNA). The activity of enhancer varies in different tissues and cells, which leads to the transcription of eRNA showing high tissue and cell specificity. Therefore, accurately capturing the transcription of eRNA in cells is an important prerequisite for studying its function.

[0003] In vivo, eRNA can participate in the regulation of promoter-enhancer looping, histone modification, transcription elongation and pause, etc. by interacting with other biological molecules (such as adhesion molecules, histone acetyltransferase, transcription factors, etc.). The regulatory role of eRNA can affect the occurrence and development of various diseases, such as tumors, rheumatoid arthritis, cardiovascular diseases, etc. Therefore, exploring the function of eRNA is of great significance for studying the mechanism of disease occurrence and development. However, the functions of most eRNAs are still unclear. Due to the large scale of eRNA transcription in the genome, time-consuming and laborious experimental verification cannot meet the needs of its functional research, which makes the computer technology-based functional prediction method particularly important.

[0004] Currently, there are several databases (such as HeRA, TCeA, Animal-eRNAdb and eRic) for recording the transcription and potential targets of eRNA. However, the transcription and regulatory relationship of eRNA has tissue and cell specificity, and the database is difficult to meet the personalized research needs. In addition, the existing non-coding RNA functional prediction platform is not suitable for eRNA (enhancer RNA). For example, ncFANsv2.0 requires known gene names or IDs as input, but eRNA does not have fixed names or IDs. AnnoLnc2 can predict the function of non-coding RNA through co-expression network, but it does not consider tissue and cell specificity, nor does it provide the characteristics of eRNA (such as histone modification, chromatin structure, etc.). Currently, there is no public method for eRNA transcription identification, regulatory target prediction and functional annotation based on high-throughput transcriptome sequencing data. SUMMARY

[0005] The technical problem to be solved by the present application is to provide a method for eRNA transcription identification, regulatory target prediction and functional annotation based on high-throughput transcriptome sequencing data, which has a wider application range and can be applied to all eRNAs and more accurately obtain the action form between eRNA and protein coding genes.

[0006] The technical solution adopted by the present application to solve the above technical problem is: a method for eRNA identification, regulatory target prediction and functional annotation based on high-throughput transcriptome sequencing data, comprising the following steps:

[0007] (1) eRNA identification

[0008] A. Collecting transcriptome RNA sequencing (RNA-seq) data or nascent RNA sequencing (such as GRO-seq) data of cell or tissue samples for de novo assembly to obtain assembled transcripts;

[0009] B. Filtering the transcripts obtained in step 1(A), using CPC2 software to predict the protein coding ability of the transcripts, and screening to obtain non-coding RNAs;

[0010] C. Defining the region with at least one of the following characteristics in the genomic database as an enhancer candidate region: H3K27ac histone modification characteristics, H3K4me1 histone modification characteristics, chromatin open characteristics and RNA polymerase II binding characteristics, defining the region 3000bp upstream to 3000bp downstream of the center of the enhancer candidate region as an enhancer region, or defining one or more annotated regions in the enhancer database as an enhancer region;

[0011] D. Identifying the part of the non-coding RNA obtained in step (1) B whose transcription start site is located in the enhancer region obtained in step (1) C as eRNA, and obtaining the chromosomal location information of the eRNA;

[0012] (2) Regulatory target prediction of eRNA

[0013] A. Expression quantification of eRNA and protein coding genes

[0014] Collecting tissue or cell samples for RNA-seq sequencing, or obtaining RNA-seq data from transcriptome sequencing databases (such as TCGA, GTEx), calculating the expression of eRNA and the expression of protein-coding genes annotated in the reference genome in each sample, if the obtained sample is a single type of sample, then the annotated protein-coding genes are used for subsequent analysis, if the obtained sample is a control sample containing normal tissue and disease tissue, then the protein-coding genes with differential expression are screened from the annotated protein-coding genes for subsequent analysis, and the expression of eRNA and protein-coding genes is quantified by using featureCount software or a calculation formula;

[0015] B. Constructing eRNA-protein-coding gene co-expression network

[0016] Pearson correlation analysis is performed on the eRNA expression of all samples obtained in step (2) A and the protein-coding gene expression used for subsequent analysis to construct an eRNA expression correlation matrix of eRNA and protein-coding genes with expression correlation; Pearson correlation analysis is performed between the protein-coding gene expression obtained in step (2) A for subsequent analysis to construct a protein expression correlation matrix between protein-coding genes with expression correlation; the eRNA expression correlation matrix and the protein expression correlation matrix are combined to obtain an eRNA-protein-coding gene co-expression network;

[0017] C. Constructing eRNA-centered regulatory network

[0018] Transcription factors bound to the eRNA chromosome position, eRNA-related RNA binding proteins, and eRNA-mediated loop-related target genes are obtained to construct an eRNA-centered regulatory network;

[0019] D. Integrating the eRNA-protein-coding gene co-expression network constructed in step (2) B and the eRNA-centered regulatory network constructed in step (2) C, extracting the connection relationship of eRNA and protein-coding genes that appear in both networks to construct a more reliable eRNA-protein-coding gene relationship network;

[0020] E. The eRNA-protein-coding gene relationship network constructed in step (2) D is extracted by hub method to extract protein-coding genes directly connected to eRNA; or based on clustering algorithm software (such as MCODE, SPICi) for module method (module) analysis, extracting gene modules closely connected to eRNA, and the protein-coding genes in the gene module are used as potential regulatory targets of eRNA;

[0021] (3) Functional annotation of eRNA

[0022] The potential regulatory targets of the eRNA obtained in step (2) are subjected to functional enrichment analysis by using software based on hypergeometric distribution algorithm (such as clusterProfiler, https: / / github.com / zhangyw0713 / FunctionEnrichment, etc.), and the results of eRNA functional annotation are predicted.

[0023] Further, step (1) A is specifically as follows: the fastq file of the transcriptome RNA sequencing of the collected cell or tissue sample is aligned to the reference genome to obtain the aligned BAM file, and then Stringtie software is used for de novo assembly; or the fastq file of the collected nascent RNA sequencing data is aligned to the reference genome to obtain the aligned BAM file, and then Homer software is used for de novo assembly to obtain the assembled transcripts.

[0024] Further, the filtering in step (1) B is specifically as follows:

[0025] a. If the transcript overlaps with the chromosomal location of the annotated protein-coding gene in the reference genome, it is discarded;

[0026] b. If the transcript overlaps with the blacklist region in the ENCODE database and the simple repeat region in the UCSC database, it is discarded;

[0027] c. If the transcript overlaps with the promoter region of the annotated gene in the reference genome, it is discarded;

[0028] Further, the genomic database in step (1) C includes Cistrome, ENCODE database, and the enhancer database includes SEdb v2.0, EnhancerAtlas v2.0, FANTOM5 and SCREEN / ENCODE database.

[0029] Further, the expression calculation formula in step (2) A is as follows:

[0030] Wherein FPKM represents the expression amount of eRNA or protein-coding gene in each sequencing sample, ∑(Cov) represents the read coverage in the chromosomal location interval of eRNA or protein-coding gene, R represents the read length of the base sequence obtained by single sequencing of the sequencer, L is the length of eRNA or protein-coding gene, T is the library size, i.e. the sum of all reads in each sequencing sample, and the sample size N>10.

[0031] Further, step (2)B is specifically as follows:

[0032] a. Perform Pearson correlation analysis on the eRNA expression of all samples obtained in step (2)A and the protein-coding gene expression for subsequent analysis to obtain correlation coefficient r e If the expression of a certain eRNA and a certain protein-coding gene has p<0.05 and the correlation coefficient |r e |>0.4, then the protein-coding gene and the eRNA have expression correlation, and an eRNA expression correlation matrix of eRNA and protein-coding genes with expression correlation is constructed;

[0033] b. Perform Pearson correlation analysis on the protein-coding gene expression of all samples obtained in step (1)A for subsequent analysis between each other to obtain correlation coefficient r p If the expression of a certain two protein-coding genes has p<0.05 and the correlation coefficient |r p |>0.7, then the two protein-coding genes have expression correlation, and a protein expression correlation matrix between protein-coding genes with expression correlation is constructed;

[0034] c. Merge the eRNA expression correlation matrix and the protein expression correlation matrix to obtain an eRNA-protein-coding gene co-expression network.

[0035] Further, step (2)C is specifically as follows:

[0036] a. Obtain transcription factors having binding sites in the chromosomal position region of eRNA from transcription factor binding map data (such as ChIP-seq data in Cistrome and ENCODE databases) to construct the regulatory relationship between eRNA and transcription factors;

[0037] b. Obtain RNA binding proteins having mutual binding relationship with eRNA from RNA binding protein binding map data (such as CLIP-seq data in POSTAR3 database) to construct the regulatory relationship between eRNA and RNA binding proteins;

[0038] c. Obtain genes corresponding to promoters having loop formation with the enhancer region where eRNA is located from chromatin interaction data (such as HiChIP data in HiChIPdb database), and the gene is a loop-related target gene regulated by eRNA through mediation of loop formation, and the regulatory relationship between eRNA and the loop-related target gene is constructed;

[0039] d. Integrate the regulatory relationship of eRNA with transcription factors, RNA binding proteins, and enhancer-promoter loop-related target genes to construct an eRNA-centered regulatory network.

[0040] Further, the functional enrichment analysis in step (3) includes gene ontology (GO), KEGG pathway, and MSigDB Hallmark feature gene analysis.

[0041] Compared with the prior art, the present application has the advantages that:

[0042] (1) The present application is the only one-stop solution for eRNA identification, regulatory target prediction, and functional annotation at present; the present application defines enhancer regions using a method combining multiple enhancer markers. Compared with the traditional method based on a single enhancer marker, this method has stronger flexibility.

[0043] (2) A large number of eRNAs in organisms have not been identified, and the past RNA functional annotation platform (such as ncFANs-NET) can only serve the RNAs that have been annotated and have official IDs or names, and cannot serve most eRNAs. The present application performs regulatory target prediction and functional annotation based on the chromosomal position information of eRNAs, and can be applied to all eRNAs, and has stronger biological practicability.

[0044] (3) The present application uses a method combining co-expression networks and eRNA-centered regulatory networks to construct an eRNA-protein coding gene relationship network. Compared with the method based on co-expression networks alone, the present application can more accurately obtain the action forms between eRNAs and protein coding genes.

[0045] (4) In the construction of eRNA-protein coding regulatory relationship, tissue or cell-specific data can be used to obtain eRNA regulatory targets and functions in a specific tissue or cell. Compared with tools such as AnnoLnc2 that do not consider biological specificity, the present application can provide more accurate clues for biological research of eRNAs.

[0046] (5) The expression quantification formula used in the present application is 400 times faster than the traditional featureCounts expression quantification method, and the expression quantification results have high similarity Figure 2 ) to the traditional method. BRIEF DESCRIPTION OF DRAWINGS

[0047] Figure 1 The left graph shows the process of identifying eRNAs from transcriptome sequencing data, and the right graph shows the flowchart of eRNA regulatory target prediction and functional annotation;

[0048] Figure 2 The chromosomal position of eRNA STRG.10047.2;

[0049] Figure 3 To compare the expression quantification formula used in this invention with the traditional featureCounts method, (A) the correlation between gene expression levels in the GTEx database obtained by the two methods, and (B) the correlation between gene expression levels in the TCGA database obtained by the two methods. In (A) and (B), the x-axis represents the tissue or cancer type, and the y-axis represents the Pearson correlation coefficient. (C) A comparison of the processing speed of the two methods on the same expression profile (GTEX-ZYFC-2626-SM-5NQ6S). The task was performed by a Dell Precision T7920 workstation.

[0050] Figure 4 Comparison of eRNA STRG.10047.2 expression levels in colon adenocarcinoma and normal colon samples;

[0051] Figure 5 This represents the co-expression network between eRNA STRG.10047.2 and protein-coding genes;

[0052] Figure 6 It is a regulatory network centered on eRNA STRG.10047.2;

[0053] Figure 7 This is a high-reliability eRNA-protein-coding gene-related network formed by integrating co-expression networks and regulatory networks;

[0054] Figure 8 To use the modular approach from Figure 7 Tightly connected subnetworks extracted from the network;

[0055] Figure 9 The results of functional enrichment analysis of protein-coding genes closely linked to eRNA STRG.10047.2 are as follows: (A) Gene ontology enrichment analysis, (B) KEGG pathway enrichment analysis, and (C) MSigDB characteristic gene functional enrichment analysis. Detailed Implementation

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

[0057] The following section will use RNA-seq data (sample ID: TCGA-5M-AAT6-01A) of human colon adenocarcinoma from the TCGA database as an example to illustrate the method used for eRNA identification, regulatory target prediction, and functional annotation. Specific Implementation Example 1

[0059] The steps for eRNA identification are as follows: Figure 1 As shown on the left, it specifically includes:

[0060] 1. Download the fastq files of RNA-seq data of colon adenocarcinoma sample TCGA-5M-AAT6-01A from TCGA database, align them to human hg38 reference genome using HISAT2 software, obtain the aligned BAM files, and use Stringtie software to perform de novo assembly on the BAM files to obtain the assembled transcripts;

[0061] 2. Filter the transcripts obtained in step 1, and use CPC2 software to predict the protein coding ability of the transcripts to screen and obtain non-coding RNAs; wherein the filtering is as follows:

[0062] A. If the transcript overlaps with the chromosomal location of the human protein-coding genes annotated in the GENCODE database, it is discarded;

[0063] B. If the transcript overlaps with the blacklist region in the ENCODE database and the simple repeat region in the UCSC database, it is discarded;

[0064] C. If the transcript overlaps with the promoter region (defined as the 2kb region upstream and downstream of the transcription start site) of the annotated human genes in the GENCODE database, it is discarded;

[0065] 3. Obtain the information of H3K27ac and H3K4me1 histone modification regions of human colon tissue from Cistrome database, and obtain the regions in the genome that simultaneously appear H3K27ac and H3K4me1 histone modification. The region from 3000bp upstream to 3000bp downstream of the center position of these regions is defined as the enhancer region. In this example, a total of 5654 enhancer regions are obtained;

[0066] 4. Identify the part of the non-coding RNAs in step 2 whose transcription start site is located in the enhancer region obtained in step 3 as eRNA. In this example, a total of 124 non-coding RNAs with transcription start site in the enhancer region are obtained, i.e. the eRNAs identified from the colon adenocarcinoma RNA-seq data, and the chromosomal position information of the eRNAs is obtained, which is shown in Table 1.

[0067] Table 1: List of eRNAs identified from colon adenocarcinoma sequencing samples.

[0068]

[0069]

[0070]

[0071]

[0072] Among the eRNAs identified in Table 1 above, STRG.10047.2 has a large range of overlap with the known colon cancer related eRNA CCAT1 in the chromosome position, as shown in Figure 2 . This suggests that the identified eRNA STRG.10047.2 may be CCAT1. Existing literature reports that CCAT1 participates in the regulation of the occurrence and development of colon cancer by affecting the cell cycle. The following examples will continue to predict the regulatory targets and potential functions of STRG.10047.2, and observe whether they are consistent with the functions of CCAT1. Specific embodiment two

[0074] The steps of eRNA target prediction and function annotation are as shown in Figure 1 , and specifically include:

[0075] 1. Obtain the RNA-seq data of 454 colon adenocarcinoma samples in the TCGA database and the RNA-seq data of 374 normal colon tissue (Colon Sigmoid) samples in the GTEx database, calculate the expression amount of eRNA STRG.10047.2 and the expression amount of the protein-coding genes already annotated in the GENCODE database in each sample, and obtain the protein-coding genes with differential expression by comparing the expression amounts of the cancer tissues and the normal tissues, wherein the expression amount calculation formula is as follows:

[0076] wherein FPKM represents the expression amount of the eRNA or the protein-coding gene in each sequencing sample, ∑(Cov) represents the read coverage in the chromosome position interval of the eRNA or the protein-coding gene, R represents the read length of the base sequence obtained by single sequencing of the sequencer, L is the length of the eRNA or the protein-coding gene, T is the library size, i.e., the sum of all reads in each sequencing sample, and the sample amount N > 10.

[0077] Figure 3 For comparison of the eRNA quantification method of the present application with the traditional featureCounts, (A) the correlation between the gene expression amounts obtained by the two methods of the GTEx database, (B) the correlation between the gene expression amounts obtained by the two methods of the TCGA database, and (C) the processing speed comparison of the same expression profile (GTEX-ZYFC-2626-SM-5NQ6S) by the two methods. The task is performed by a Dell Precision T7920 workstation. The correlation between the gene expression amounts obtained by the two methods of the GTEx database is shown in Figure 3 (A), Figure 3 (B) and Figure 3(C)It can be known that the eRNA expression quantification method is 400 times faster than the traditional featureCounts expression quantification method, and the expression quantification results are highly similar.

[0078] By Figure 4 It can be known that the expression amount of eRNA STRG.10047.2 in cancer tissues is significantly higher than that in normal tissues, which indicates that STRG.10047.2 plays an important role in cancer. At the same time, 4234 protein-coding genes annotated in the GENCODE database and differentially expressed in colon adenocarcinoma are obtained by differential analysis.

[0079] 2. Constructing eRNA-protein-coding gene co-expression network in colon adenocarcinoma:

[0080] According to the expression amount of eRNA STRG.10047.2 and the differentially expressed and annotated protein-coding genes in the TCGA colon adenocarcinoma sample in step 2, Pearson correlation analysis is performed on the expression amount of eRNA STRG.10047.2 and the differentially expressed and annotated protein-coding genes; the expression amount of the differentially expressed and annotated protein-coding genes is pairwise correlated; according to the threshold of p<0.05 and the correlation coefficient |r e |>0.4, 233 protein-coding genes with expression correlation with STRG.10047.2 are screened, and an eRNA expression correlation matrix is constructed; according to the threshold of p<0.05 and the correlation coefficient |r e |>0.7, the expression correlation matrix between the 233 protein-coding genes is obtained, and a protein-coding gene expression correlation matrix is constructed; the expression correlation matrix of STRG.10047.2 and the protein-coding genes and the expression correlation matrix between the protein-coding genes are combined to obtain a STRG.10047.2-protein-coding gene co-expression network. As Figure 5 shown, the network consists of 234 nodes (1 eRNA and 233 protein-coding genes) and 3229 edges.

[0081] 3. Constructing STRG.10047.2-centered regulatory network:

[0082] A. Obtain transcription factors with binding sites at the chromosomal position of eRNA STRG.10047.2 from the Cistrome database, and construct the regulatory relationship between eRNA and transcription factors;

[0083] B. Obtain RNA binding proteins with mutual binding relationship with STRG.10047.2 from CLIP-seq data in the POSTAR3 database, and construct the regulatory relationship between eRNA and RNA binding proteins;

[0084] C. Obtain the gene corresponding to the promoter existing in the loop from the HiChIP data in the HiChIPdb database in the enhancer region where STRG.10047.2 is located, which is the loop-related target gene regulated by the eRNA mediated loop, and construct the regulatory relationship between the eRNA and the loop-related target gene;

[0085] D. Integrate the regulatory relationship of eRNA STRG.10047. and transcription factors, RNA binding proteins, and sub-loop-related target genes, and construct the regulatory network centered on STRG.10047.2. As shown in Figure 6 , the network consists of 147 nodes (1 eRNA and 146 regulatory targets) and 3272 edges.

[0086] 4. Integrate the eRNA-protein coding gene co-expression network obtained in step 2 and the regulatory network centered on STRG.10047.2 obtained in step 3, extract the eRNA STRG.10047.2 related protein coding genes (Table 2) that appear in both networks, and construct the integrated STRG.10047.2 related network Figure 7 );

[0087] Table 2 eRNA STRG.10047.2 related protein coding genes in the integrated network.

[0088]

[0089] 5. Use SPICi software for module method analysis, extract the gene module closely connected to STRG.10047.2 from the integrated STRG.10047.2 related network obtained in step 4, and the protein coding genes in the gene module are the potential regulatory targets of STRG.10047.2. The closely connected module needs to meet the requirements of node density greater than 0.5 and support greater than 0.5.

[0090] As shown in Figure 8 , it can be seen from Figure 8 that multiple CCAT1 regulatory targets appear in the regulatory network of STRG.10047.2 (such as SOX4 and EZH2), further proving that STRG.10047.2 and CCAT1 may be the same eRNA. Specific embodiment three

[0092] Functional annotation

[0093] The protein-coding genes obtained in Specific Embodiment Two that have regulatory relationship with eRNA STRG.10047.2 were subjected to Gene ontology (GO) enrichment analysis and KEGG pathway analysis using clusterProfiler, respectively, and it was found that these protein-coding genes were mainly involved in the regulation of cell cycle. The protein-coding genes obtained in Specific Embodiment Two were subjected to MSigDB Hallmark enrichment analysis using https: / / github.com / zhangyw0713 / FunctionEnrichment.

[0094] Figure 9 The results of functional enrichment analysis of protein-coding genes closely connected with eRNA STRG.10047.2. (A) Gene ontology enrichment analysis. (B) KEGG pathway enrichment analysis. (C) MSigDB Hallmark enrichment analysis. The results were obtained by using clusterProfiler and https: / / github.com / zhangyw0713 / FunctionEnrichment, respectively. Figure 9 It can be seen that these protein-coding genes have significant overlap with MYC, E2F-related genes and G2M cell cycle checkpoint-related genes, which more powerfully proves that these genes are involved in the regulation of cell cycle. Since STRG.10047.2 exists in the same closely connected subnetwork with these protein-coding genes, it can be considered that STRG.10047.2 also has similar functions, i.e. the potential function of STRG.10047.2 is predicted to be affecting the occurrence and development of colon adenocarcinoma by regulating cell cycle.

[0095] In summary, Specific Embodiment One found that STRG.10047.2 overlaps with eRNA CCAT1 in chromosomal location, Specific Embodiment Two found that STRG.10047.2 and CCAT1 have many consistent regulatory targets, and Specific Embodiment Three found that STRG.10047.2 and CCAT1 have consistent functions of regulating cell cycle. Therefore, it can be judged that they are the same eRNA, and at the same time, the reliability of the method for predicting eRNA identification, regulatory targets and functions is also proved.

[0096] The above description is not a limitation of the present application, and the present application is not limited to the above examples. Changes, modifications, additions or substitutions made by those of ordinary skill in the art within the spirit and scope of the present application should also be within the protection scope of the present application.

Claims

1. A method for eRNA identification, regulatory target prediction and functional annotation based on high-throughput transcriptome sequencing data, characterized in that The method comprises the following steps: (1) eRNA identification A. Collecting transcriptome RNA sequencing data or nascent RNA sequencing data of a cell or tissue sample for de novo assembly to obtain assembled transcripts; B. Filtering the transcripts obtained in step (1) A, and using CPC2 software to predict the protein coding ability of the transcripts to screen non-coding RNAs; C. Defining a region having at least one of the following characteristics in a genome database as a candidate enhancer region: H3K27ac histone modification characteristics, H3K4me1 histone modification characteristics, chromatin opening characteristics and RNA polymerase II binding characteristics, defining a region 3000 bp upstream to 3000 bp downstream of the center of the candidate enhancer region as an enhancer region, or defining one or more regions annotated in an enhancer database as an enhancer region; D. Identifying the part of the non-coding RNAs obtained in step (1) B, in which the transcription start site is located in the enhancer region obtained in step (1) C, as eRNA, and obtaining the chromosomal location information of the eRNA; (2) Prediction of the regulatory target of eRNA A. Expression quantification of eRNA and protein coding genes Collecting tissue or cell samples for RNA-seq sequencing, or obtaining RNA-seq data from a transcriptome sequencing database, calculating the expression of eRNA and the expression of protein coding genes annotated in the reference genome in each sample, if the obtained sample is a single type of sample, using the annotated protein coding genes for subsequent analysis, if the obtained sample is a control sample containing normal tissue and disease tissue, screening the protein coding genes with differential expression from the annotated protein coding genes for subsequent analysis, and quantifying the expression of eRNA and protein coding genes using featureCount software or a calculation formula; B. Construction of eRNA-protein coding gene co-expression network Performing Pearson correlation analysis on the eRNA expression of all samples obtained in step (2) A and the expression of protein coding genes used for subsequent analysis to construct an eRNA expression correlation matrix of eRNA and protein coding genes having expression correlation with the eRNA; performing Pearson correlation analysis between the expression of protein coding genes used for subsequent analysis obtained in step (2) A to construct a protein expression correlation matrix between protein coding genes having expression correlation; and combining the eRNA expression correlation matrix and the protein expression correlation matrix to obtain an eRNA-protein coding gene co-expression network; C. Construction of eRNA-centered regulatory network Obtaining transcription factors combined at the chromosomal location of eRNA, RNA binding proteins related to eRNA, and eRNA-mediated loop-related target genes to construct an eRNA-centered regulatory network. D. Integrating the eRNA-protein coding gene co-expression network constructed in step (2)B and the eRNA-centered regulatory network constructed in step (2)C, the connection relationship of eRNA and protein coding gene appearing in both networks is extracted to construct a more reliable eRNA-protein coding gene relationship network; E. The eRNA-protein coding gene relationship network constructed in step (2)D is extracted by the hub method to extract protein coding genes directly connected to eRNA as potential regulatory targets of eRNA; or the gene module method is used for module analysis based on clustering algorithm software to extract protein coding genes in the gene module as potential regulatory targets of eRNA; (3) Functional annotation of eRNA The potential regulatory targets of eRNA obtained in step (2) are subjected to functional enrichment analysis by software based on hypergeometric distribution algorithm to predict the results of eRNA functional annotation.

2. The method for eRNA identification, regulatory target prediction and functional annotation based on high-throughput transcriptome sequencing data according to claim 1, characterized in that Step (1)A is as follows: The fastq file of the collected transcriptome RNA sequencing of the cell or tissue sample is aligned to the reference genome to obtain the aligned BAM file, and then the Stringtie software is used for de novo assembly; or the fastq file of the collected nascent RNA sequencing data is aligned to the reference genome to obtain the aligned BAM file, and then the Homer software is used for de novo assembly to obtain the assembled transcripts.

3. The method for eRNA identification, regulatory target prediction and functional annotation based on high-throughput transcriptome sequencing data according to claim 1, characterized in that The filtering in step (1)B is as follows: a. If the transcript overlaps with the protein coding gene annotated in the GENCODE database in terms of chromosomal location, it is discarded; b. If the transcript overlaps with the blacklist region in the ENCODE database and the simple repeat region in the UCSC database, it is discarded; c. If the transcript overlaps with the promoter region of the GENCODE database annotated gene, it is discarded.

4. The method for eRNA identification, regulatory target prediction and functional annotation based on high-throughput transcriptome sequencing data according to claim 1, characterized in that The genomic database in step (1)C includes the Cistrome database, and the enhancer database includes SEdb v2.0, EnhancerAtlas v2.0, FANTOM5 and SCREEN / ENCODE databases.

5. The method for eRNA identification, regulatory target prediction and functional annotation based on high-throughput transcriptome sequencing data according to claim 1, characterized in that The calculation formula in step (2)A is as follows: where FPKM represents the expression amount of eRNA or protein coding gene in each sequencing sample, represents the read coverage in the chromosome position interval of eRNA or protein coding gene, R represents the read length of the base sequence obtained by single sequencing of the sequencer, L is the length of eRNA or protein coding gene, T is the library size, that is, the sum of all reads in each sequencing sample, and the sample amount N>10.

6. The method for eRNA identification, regulatory target prediction and functional annotation based on high-throughput transcriptome sequencing data according to claim 1, characterized in that Step (2)B is as follows: a. Perform Pearson correlation analysis on the eRNA expression of all samples obtained in step (2) A and the protein-coding gene expression for subsequent analysis to obtain correlation coefficient r e If the correlation coefficient |r e |>0.4, then the protein-coding gene and the eRNA have expression correlation, and an eRNA expression correlation matrix of eRNA and protein-coding genes with expression correlation is constructed. b. Perform Pearson correlation analysis between the expression levels of protein-coding genes for subsequent analysis in all samples obtained in step (1) A to obtain correlation coefficient r p If the p value between the expression levels of two protein-coding genes is less than 0.05 and the absolute value of the correlation coefficient |r p | is greater than 0.7, then the two protein-coding genes have expression correlation, and a protein expression correlation matrix between protein-coding genes with expression correlation is constructed. c. Combine the eRNA expression correlation matrix and the protein expression correlation matrix to obtain an eRNA-protein coding gene co-expression network.

7. The method for eRNA identification, regulatory target prediction and functional annotation based on high-throughput transcriptome sequencing data according to claim 1, characterized in that Step (2)C is as follows: a. From the transcription factor binding map data, obtain transcription factors having binding sites with eRNA in terms of chromosomal location to construct the regulatory relationship between eRNA and transcription factors; b. From the RNA binding protein binding map data, obtain RNA binding proteins having mutual binding relationship with eRNA to construct the regulatory relationship between eRNA and RNA binding proteins; c. obtaining the gene corresponding to the promoter existing in the enhancer region where the eRNA is located from the chromatin interaction data, which is the enhancer-promoter loop related target gene regulated by the eRNA mediated loop, and constructing the regulatory relationship between the eRNA and the loop related target gene; d. integrating the regulatory relationship between the eRNA and the transcription factor, the RNA binding protein and the loop related target gene, and constructing the regulatory network with the eRNA as the center.

8. The method for eRNA identification, regulatory target prediction and functional annotation based on high-throughput transcriptome sequencing data according to claim 1, characterized in that: The functional enrichment analysis in step (3) includes gene ontology, KEGG pathway and MSigDB Hallmark feature gene analysis.

Images