Method for identifying RNA (Ribonucleic Acid) from transposon element based on RNA sequencing data

Abstract

The invention relates to the technical field of bioinformatics, in particular to a method for identifying RNA (Ribonucleic Acid) from a transposon element based on RNA sequencing data. According to the invention, the exon level of the teRNA is taken as an analysis object, compared with a teRNA site, the teRNA exon not only realizes more accurate site identification, but also retains the existing form of TE in a transcription product, including whether a splicing site and a polyadenylic acid signal are provided or not, whether the TE is chimeric with a gene or other TE sites or not and the like. And compared with the level of the teRNA transcript, as a constituent part of the transcript, the teRNA exon significantly reduces the sequence assembly difficulty of the short-read-length sequencing data. According to the method, the Bayesian model, the parametric assembly algorithm and the local site de novo assembly algorithm are combined, and through mutual correction of different assembly algorithms and repeated comparison sequencing reading length before and after assembly, the assembly false positive is reduced, so that the teRNA recognition accuracy is improved.

Description

Technical Field

[0001] The present invention relates to the technical field of bioinformatics, and in particular to a method for identifying RNA derived from a transposable element based on RNA sequencing data. Background Art

[0002] Transposable elements (TEs) are a class of special DNA elements that are widely present in the genomes of eukaryotic organisms and can move freely. In the long-term evolutionary process, due to sequence variations such as mutations, deletions, homologous recombination, and the evolution of host epigenetic regulation, small RNA regulation and other defense mechanisms, most TEs have lost their normal expression and transposition functions and become silent DNA sequences in the genome, and have been considered "junk DNA" for a long time. However, with the development of high-throughput sequencing technology, more and more studies have found that some TEs can be reactivated and start transcription under external stimulation or in disease states, and then participate in the occurrence and development of various complex diseases such as cancer, neurodegenerative diseases, and autoimmune diseases. For example, the envelope protein encoded by HERVK can mediate cell fusion in melanoma; the expression of HERVW can activate host immunity and induce the expression of proinflammatory cytokines, driving the inflammatory response of multiple sclerosis. Not only that, under physiological conditions, TEs also have transcriptional capabilities in tissue cells, and their expression products even participate in multiple biological processes of growth and development. For example, HERVW and HERVFRD families can promote cell fusion and syncytiotrophoblast formation by expressing syncytin proteins, and are essential genes for placental morphogenesis. For another example, HERVH is highly expressed in embryonic stem cells and pluripotent stem cells, and can affect the differentiation process of cells by forming topologically related domain boundaries, shaping chromatin structure, and regulating the expression of surrounding genes, thereby participating in the maintenance of embryonic stemness. Therefore, systematically identifying TE-derived RNA (teRNA), analyzing its expression distribution pattern, and the factors affecting its expression will help understand the role of teRNA in the gene expression regulatory network, and lay the foundation for revealing the function of teRNA in growth and development and the occurrence and development of diseases.

[0003] Since TE is a multi-copy repeat element, TEs in the same family have highly similar sequences, and multi-site alignments are easily generated during sequence alignment, which affects the accuracy of teRNA identification. Therefore, most teRNA studies at this stage use TE family-level identification strategies for teRNA expression analysis, that is, TEs with homologous sequences are classified into one family as the basic unit for TE qualitative and quantitative analysis, thereby reducing the impact of multi-aligned reads. However, the teRNA identification method based on the TE family cannot explain whether the expression and function of TE sites in the same family have commonalities, and it is difficult to distinguish whether teRNA is derived from a specific TE site, which limits the identification of functional teRNA sites and the analysis of related regulatory mechanisms. In order to achieve site-specific identification of TE, the most direct method is to directly eliminate multi-site alignment reads or randomly assign them, but this will lead to a serious underestimation of the read counts of some TE sites. In addition, the teRNA site-specific identification method represented by Telescope proposes to repeatedly redistribute multi-aligned reads by establishing an expectation-maximizing hybrid model until the optimal site is assigned. However, existing teRNA site expression identification methods are all based on direct identification of genomic TE sites provided by the database, and TE sites are not all fully transcribed as single sites. They may produce truncated teRNAs, 3'-end read-through teRNAs, or even form co-transcriptional units with surrounding genes or TE sites to produce chimeric teRNAs. The teRNA identification method based on the transcript assembly strategy reconstructs transcripts to perform qualitative and quantitative analysis of teRNAs, which is closer to the real biological expression of teRNAs and is conducive to the analysis of the expression distribution patterns and related regulatory mechanisms of teRNAs. However, limited by short-read sequencing, the accuracy of teRNA transcript assembly is still relatively low. Summary of the invention

[0004] In order to solve the above problems, the present invention provides a method for identifying RNA derived from transposable elements based on RNA sequencing data. The identification method provided by the present invention is a teRNA identification strategy that combines a Bayesian model, a parameterized assembly algorithm, and a local site de novo assembly algorithm. By mutual correction of assembly results and quality control of secondary alignment of sequencing reads, the accuracy of teRNA identification is improved, providing a powerful tool for accurate screening and functional analysis of functional teRNA in the future.

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

[0006] The present invention provides a method for identifying RNA derived from a transposable element based on RNA sequencing data, comprising the following steps: (1) screening of teRNA sites; (2) assembly of sequencing reads derived from teRNA; (3) identification of teRNA;

[0007] The screening of teRNA sites in step (1) comprises the following steps: (1.1) obtaining sequencing data of double-ended RNA of the sample to be tested; (1.2) aligning the sequencing data to the reference genome to obtain an alignment result file; (1.3) allocating sequencing read lengths of the sequencing data, and retaining teRNA sites with sequencing read lengths > 0;

[0008] The assembly of the teRNA source sequencing reads in step (2) comprises the following steps: (2.1) extracting the alignment results corresponding to the teRNA sites in the alignment result file to obtain a teRNA alignment result file; the teRNA sites are the teRNA sites screened in step (1); (2.2) using the reference gene annotation file as a guide, using StringTie software to perform reference assembly on the teRNA alignment result file, with the parameter -G, to obtain a first assembly sequence; (2.3) using the recognition module of SERVE software to perform local site de novo assembly on the teRNA alignment result file, with the parameter --count set to 1, to obtain a second assembly sequence; (2.4) merging the first assembly sequence and the second assembly sequence to obtain a total assembly sequence;

[0009] The identification of teRNA in step (3) includes the following steps: (3.1) extracting the single exon transcript sequence in the assembled total sequence, using the merge module of BEDTools to merge the sequences, and obtaining a de-redundant single exon transcript sequence; (3.2) extracting the exon sequence of the multi-exon transcript in the assembled total sequence, and aligning and merging it with the de-redundant single exon transcript sequence to obtain a de-redundant exon sequence; the alignment and merging method includes: performing interval comparison between the exon sequence and the single exon transcript sequence, and for the exon sequence overlapping with the single exon transcript sequence, taking the interval union of the overlapping exon sequence; (3.3) extracting the teRNA exons in the de-redundant exon sequence, and using the coverage module of the BEDTools software to calculate the coverage; (3.4) eliminating teRNA exons with sequencing coverage <85% to obtain qualitatively identified teRNA.

[0010] Preferably, the identification method further comprises: (4) annotation of teRNA; (5) quantification of teRNA;

[0011] The annotation of the teRNA in step (4) comprises the following steps: (4.1) using the annotation of the transposable element database as a reference, annotating the family and site of the qualitatively identified teRNA to obtain a first annotation result; (4.2) using the annotation of the gene database as a reference, annotating the chimeric gene or adjacent gene of the qualitatively identified teRNA to obtain a second annotation result; (4.3) predicting the protein coding potential of the qualitatively identified teRNA to obtain a third annotation result; (4.4) merging the first annotation result, the second annotation result and the third annotation result to obtain an annotation set;

[0012] The quantification of the teRNA in step (5) includes the following steps: (5.1) using RSEM software to quantify the annotation set, standardizing the exon length, and obtaining an exon quantification matrix; (5.2) using Telescope software to quantify the teRNA sites to obtain a site quantification matrix; the teRNA is the teRNA site screened in step (1); (5.3) adding the expression levels of the teRNA sites from the same family in the site quantification matrix to obtain a family quantification matrix.

[0013] Preferably, the method for obtaining the sequencing data includes: obtaining from a public database or obtaining by sequencing; the public database includes one or more of GEO, GTEx, TCGA and CGGA.

[0014] Preferably, the method for allocating sequencing read lengths of the sequencing data comprises a read length allocation algorithm of Telescope software.

[0015] Preferably, the tool for aligning the sequencing data to a reference genome includes STAR software; the reference genome includes the human reference genome GRCh38 or T2T-CHM13.

[0016] Preferably, the extraction tool in step (2.1) includes SAMtools software.

[0017] Preferably, the transposable element database comprises the Dfam database.

[0018] Preferably, the gene database includes the GENCODE database.

[0019] Preferably, the tools for predicting the protein coding potential of the qualitatively identified teRNA include CPC2 software, CPAT software and PfamScan software.

[0020] Preferably, in step (5.1), the exon length is normalized using formula I.

[0021]

[0022] Among them, Ce represents the count of exons, Le represents the length of exons, Ct represents the count of transcripts, and Lt represents the length of transcripts.

[0023] Beneficial effects:

[0024] The present invention proposes to analyze the exon level of transposable element derived RNA (teRNA). Compared with teRNA sites, teRNA exons not only achieve more accurate site identification, but also retain the form of transposable elements (TE) in the transcription product, including whether splicing sites and polyadenylation signals are provided, and whether they are chimeric with genes or other TE sites. Compared with the teRNA transcript level, as a component of the transcript, teRNA exons significantly reduce the difficulty of sequence assembly of short-read sequencing data.

[0025] On the other hand, the present invention combines the Bayesian model (read length allocation in step (1.3)), the parametric assembly algorithm and the local site de novo assembly algorithm, reduces the false positive of the assembly through mutual correction of different assembly algorithms, and repeatedly compares the sequencing read lengths before and after assembly, thereby improving the accuracy of teRNA recognition.

[0026] In addition, the existing technologies mainly only identify and analyze a certain research level based on teRNA, resulting in the inability to compare the analysis results of different studies. Based on the high accuracy of teRNA qualitative identification, the present invention also provides a quantitative method for teRNA at multiple research levels (including family, site, exon), providing a convenient, unified and accurate method for future research. BRIEF DESCRIPTION OF THE DRAWINGS

[0027] 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.

[0028] Figure 1 Schematic diagram of the teRNA identification algorithm;

[0029] Figure 2 It is a technical step flow chart of the present invention;

[0030] Figure 3 Electropherogram for RT-PCR experiment to verify teRNA exon;

[0031] Figure 4 Results from RT-qPCR experiments to verify the accuracy of teRNA quantification;

[0032] Figure 5The dimensionality reduction analysis results of glioblastoma and adjacent normal samples are shown below;

[0033] Figure 6 The results of differential expression analysis of teRNA exons;

[0034] Figure 7 This is the enrichment analysis result of teRNA-related genes;

[0035] Figure 8 These are two examples of teRNAs being associated with poor prognosis in patients with glioblastoma;

[0036] Fig. 9 The results of dimensionality reduction analysis based on teRNA exon expression profiles of cell samples during neuronal differentiation;

[0037] Fig.10 The characteristic genes of the 14 co-expression modules are grouped according to their dynamic expression patterns;

[0038] Fig.11 The performance comparison results of the method of the present invention and the prior art in the simulated data set are shown;

[0039] Fig.12 The performance comparison results of the method of the present invention and the prior art in real public data sets are shown. DETAILED DESCRIPTION

[0040] The present invention provides a method for identifying RNA derived from a transposable element based on RNA sequencing data, comprising the following steps: (1) screening of teRNA sites; (2) assembling sequencing reads derived from teRNA; (3) identifying teRNA;

[0041] The screening of teRNA sites in step (1) comprises the following steps:

[0042] (1.1) obtaining sequencing data of double-ended RNA of the sample to be tested; (1.2) aligning the sequencing data to the reference genome to obtain an alignment result file; (1.3) allocating sequencing read lengths of the sequencing data, and retaining teRNA sites with sequencing read lengths > 0;

[0043] The assembly of the teRNA-derived sequencing reads in step (2) comprises the following steps:

[0044] (2.1) extracting the alignment results corresponding to the teRNA sites in the alignment result file to obtain a teRNA alignment result file; the teRNA sites are the teRNA sites screened in step (1); (2.2) using the reference gene annotation file as a guide, using StringTie software to perform reference assembly on the teRNA alignment result file, with the parameter -G, to obtain a first assembly sequence; (2.3) using the recognition module of SERVE software to perform local site de novo assembly on the teRNA alignment result file, with the parameter --count set to 1, to obtain a second assembly sequence; (2.4) merging the first assembly sequence and the second assembly sequence to obtain a total assembly sequence;

[0045] The identification of teRNA in step (3) comprises the following steps:

[0046] (3.1) extracting the single exon transcript sequence from the assembled total sequence, and using the merge module of BEDTools to merge the sequences to obtain a de-redundant single exon transcript sequence; (3.2) extracting the exon sequence of the multi-exon transcript from the assembled total sequence, and aligning and merging it with the de-redundant single exon transcript sequence to obtain a de-redundant exon sequence; the alignment and merging method comprises: performing interval comparison between the exon sequence and the single exon transcript sequence, and for the exon sequence overlapping with the single exon transcript sequence, taking the interval union of the overlapping exon sequence; (3.3) extracting the teRNA exon from the de-redundant exon sequence, and using the coverage module of the BEDTools software to calculate the coverage; (3.4) eliminating the teRNA exons with sequencing coverage <85% to obtain qualitatively identified teRNA.

[0047] The identification method provided by the present invention is a teRNA identification strategy that combines a Bayesian model, a parameterized assembly algorithm, and a local site de novo assembly algorithm. The accuracy of teRNA identification is improved by mutual correction of assembly results and quality control of secondary alignment of sequencing reads, providing a powerful tool for precise screening and functional analysis of functional teRNAs in the future.

[0048] As an embodiment, the identification method further comprises: (4) annotation of teRNA; (5) quantification of teRNA;

[0049] The annotation of the teRNA in step (4) comprises the following steps: (4.1) using the annotation of the transposable element database as a reference, annotating the family and site of the qualitatively identified teRNA to obtain a first annotation result; (4.2) using the annotation of the gene database as a reference, annotating the chimeric gene or adjacent gene of the qualitatively identified teRNA to obtain a second annotation result; (4.3) predicting the protein coding potential of the qualitatively identified teRNA to obtain a third annotation result; (4.4) merging the first annotation result, the second annotation result and the third annotation result to obtain an annotation set;

[0050] The quantification of the teRNA in step (5) includes the following steps: (5.1) using RSEM software to quantify the annotation set, standardizing the exon length, and obtaining an exon quantification matrix; (5.2) using Telescope software to quantify the teRNA sites to obtain a site quantification matrix; the teRNA is the teRNA site screened in step (1); (5.3) adding the expression levels of the teRNA sites from the same family in the site quantification matrix to obtain a family quantification matrix.

[0051] Based on the high accuracy of teRNA qualitative identification, the present invention also provides a quantitative method for teRNA at multiple research levels (including family, site, and exon), providing a convenient, unified, and accurate method for future research.

[0052] As an embodiment, the method for obtaining the sequencing data includes: obtaining from a public database or performing sequencing; the public database includes one or more of GEO, GTEx, TCGA and CGGA. As another embodiment, the format of the sequencing data is FASTQ format.

[0053] As an embodiment, the tool for aligning the sequencing data to the reference genome includes STAR software; the reference genome includes human reference genome GRCh38 or T2T-CHM13. As another embodiment, the format of the alignment result file is BAM format.

[0054] As an embodiment, the method for allocating the sequencing read length of the sequencing data includes the read length allocation algorithm of the Telescope software. As an embodiment, the tool extracted in step (2.1) includes SAMtools software. As another embodiment, the format of the teRNA comparison result file is BAM format. As an embodiment, the format of the first assembly sequence is GTF format, and the format of the second assembly sequence is GTF format. As an embodiment, the transposable element database can be a Dfam database. As an embodiment, the gene database can be a GENCODE database. As an embodiment, the tool for predicting the protein coding potential of the qualitatively identified teRNA includes CPC2 software, CPAT software and PfamScan software.

[0055] As an embodiment, in step (5.1), the exon length is normalized using formula I,

[0056]

[0057] Among them, Ce represents the count of exons, Le represents the length of exons, Ct represents the count of transcripts, and Lt represents the length of transcripts.

[0058] To further illustrate the present invention, a method for identifying RNA derived from a transposable element based on RNA sequencing data 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.

[0059] Example 1

[0060] Taking HEK293T and U251 cell lines as examples to verify the accuracy of teRNA qualitative identification, the steps are as follows:

[0061] (1) Screening steps for teRNA sites: (1.1) Perform double-end RNA sequencing on HEK293T and U251 cell lines respectively to obtain raw sequencing data; (1.2) Use STAR software to align the raw sequencing data to the human reference genome (T2T-CHM13 version) to obtain a BAM format alignment result file; (1.3) Use the assign module of Telescope software to assign sequencing reads; (1.4) Extract teRNA sites covered by sequencing reads (retain teRNA sites with sequencing reads > 0).

[0062] (2) Assembly steps of teRNA-derived sequencing reads: (2.1) According to the teRNA sites screened in step (1), the teRNA alignment result file (BAM format) was extracted using the SAMtools software; (2.2) Guided by the human reference gene annotation file (T2T-CHM13 version), StringTie was used for referenced assembly (parameter -G) to obtain the assembled sequence (GTF format); (2.3) The recognition module of the SERVE software was used for de novo assembly of local sites, with the parameter --count set to 1, to obtain the assembled sequence (GTF format); (2.4) The assembled sequences of the first two steps were merged.

[0063] (3) teRNA identification steps: (3.1) Extract the single exon transcript sequence of the assembled sequence of step (2), and use the merge module of BEDTools to merge the sequences to obtain the de-redundant single exon transcript sequence; (3.2) Extract the exon sequence of the multi-exon transcript of the assembled sequence of step (2), and compare the intervals with the de-redundant single exon transcript sequence obtained in step (3.2). For the exon sequence overlapping with the single exon transcript, take the interval union of the overlapping exons to obtain the de-redundant exon sequence; (3.3) Extract teRNA exons, and use the coverage module of BEDTools to calculate the coverage; (3.4) Eliminate teRNA exons with sequencing coverage <85%.

[0064] The teRNA specifically recognized by the present invention was selected for experimental verification, and the primer sequences are shown in Table 1.

[0065] Table 1 Amplification primers for different teRNAs

[0066]

[0067]

[0068] Results Figure 3 The selected teRNAs were all experimentally proven to be true positives, which illustrates the accuracy of the qualitative identification of teRNAs in the present invention.

[0069] Example 2

[0070] HEK293T and U251 cell lines were used as examples to verify the accuracy of teRNA quantification. The schematic diagram of the teRNA identification algorithm is shown in Figure 1 , the technical steps flow chart is shown in Figure 2 , the steps are as follows:

[0071] Based on the teRNA qualitatively identified in Example 1, the following steps were performed:

[0072] (4) The annotation of teRNA includes the following steps: (4.1) Using the TE sites annotated in the Dfam database as a reference, use BEDTools to annotate TE families and sites; (4.2) Using the GENCODE gene annotation as a reference, use BEDTools to annotate teRNA chimeric genes and teRNA adjacent genes; (4.3) Use CPC2, CPAT and PfamScan to predict the protein coding potential of teRNA.

[0073] (5) Quantification of teRNA includes the following steps: (5.1) merging the teRNA annotations obtained in step (4) with the gene reference annotations, performing quantification using RSEM software, and normalizing the exon length using formula I to obtain an exon quantification matrix;

[0074]

[0075] Wherein, Ce represents the count of exon, Le represents the length of exon, Ct represents the count of transcript, and Lt represents the length of transcript; (5.2) Use Telescope software to quantify the sites of teRNA and obtain the site quantification matrix; (5.3) Add the expression levels of teRNA sites from the same family to obtain the family quantification matrix.

[0076] teRNA was randomly selected for quantitative experimental verification, and the primer sequences are shown in Table 2.

[0077] Table 2 Amplification primers for different teRNAs

[0078]

[0079] Results Figure 4 The quantitative results of the present invention are highly consistent with the actual experimental results, proving the accuracy of the teRNA quantification of the present invention.

[0080] Example 3 Tumor Data Application

[0081] Taking the public data set of glioblastoma (ID: PRJNA613939) as an example, the present invention is applied to perform a full set of teRNA analysis processes to achieve a comprehensive analysis of tumor teRNA. The steps are as follows:

[0082] For each sample, the technical steps of Example 1 were used for analysis to obtain the teRNA annotation of each sample; then, the samples were merged using Cuffmerge software to obtain a de-redundant teRNA annotation file; and then the technical steps of Example 2 were used to quantify the teRNA of each sample.

[0083] (6) Downstream analysis includes the following steps: (6.1) Dimensionality reduction analysis was performed using the expression levels of teRNA at different research levels (cmdscale function in R language). The dimensionality reduction analysis results of glioblastoma and adjacent normal samples are shown in Figure 5 It can be seen that different teRNA research levels can clearly distinguish tumor and adjacent normal samples, and the teRNA exon level also reflects the heterogeneity of tumor samples; (6.2) teRNA differential expression analysis was performed using edgeR software. The results are shown in Figure 6 , it can be seen that the expression profiles of teRNA exons from the same family are significantly different; (6.3) clusterProfiler was used to perform teRNA-related gene enrichment analysis. The results are shown in Figure 7 , it can be seen that teRNA in glioblastoma may cause the occurrence of body immunity; (6.4) Based on the clinical information of the patients, survival analysis was performed. The analysis results of two examples of teRNA being associated with poor prognosis of patients with glioblastoma are shown in Figure 8 It can be seen that some teRNAs are significantly correlated with the poor prognosis of tumor patients.

[0084] Example 4 Application of Normal Cell Differentiation Data

[0085] Taking the neuronal differentiation public data set (ID: PRJNA596331) as an example, the data set includes various stages of neuronal differentiation (early differentiated cells, neural progenitor cells, neural progenitor cells with rosettes, and mature neuronal cells) as well as undifferentiated neural stem cells and purified neuronal cells. The present invention is applied to the full set of teRNA analysis procedures (analysis steps are the same as in Example 3), and a total of 52,838 expressed teRNA exons are identified, covering 1,075 TE families, indicating that teRNA is widely expressed during neuronal differentiation.

[0086] (7) Downstream analysis includes the following steps:

[0087] (7.1) Dimensionality reduction analysis was performed using the expression level of teRNA exons. The results of dimensionality reduction analysis based on teRNA exon expression profiles of cell samples during neuronal differentiation are shown in Fig. 9 ,The results showed that the samples were clustered according to the cell types, suggesting that teRNA may be closely related to the differentiation process of neuronal cells;

[0088] (7.2) The teRNA exon co-expression regulatory network was constructed using WGCNA software, and 14 teRNA exon co-expression modules were detected. According to the expression of the module characteristic genes at different differentiation times, Fig.10, the modules can be divided into 6 categories: neural stem cell high expression module (modules 7 and 8), early differentiated cell high expression module (modules 6, 12, 13, 14), neural progenitor cell high expression module (module 5), rosette high expression module (modules 2, 3, 10), neural progenitor cell low expression module (modules 9 and 11) and mature neuron high expression module (module 14), indicating that teRNA exons in different modules show specific dynamic expression at different differentiation times.

[0089] Comparison of simulation data of comparative example 1

[0090] By generating multiple RNA sequencing simulation data, the performance of the method provided by the present invention (Example 1) and the prior art at different research levels of teRNA was compared. The comparison results are shown in Fig.11 , the method of the prior art is as follows:

[0091] StringTie_unguided and StringTie_guided can be found in the literature [PerteaM, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg S L. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads [J]. Nat Biotechnol, 2015, 33(3):290-295.];

[0092] Cufflinks_unguided and Cufflinks_guided can be found in the literature [Trapnell C, Williams BA, Pertea G, MortazaviA, Kwan G, van Baren MJ, Salzberg SL, Wold BJ, Pachter L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation[J].NatBiotechnol,2010,28(5):511-515.];

[0093] For LIONS, please refer to the literature [BabaianA, Thompson IR, Lever J, Gagnier L, Karimi MM, Mager D L.LIONS: analysis suite for detecting and quantifying transposableelement initiated transcription from RNA-seq[J]. Bioinformatics, 2019, 35(19): 3839-3841.];

[0094] For SERVE, please refer to the literature [She J, Du M, XuZ, JinY, Zhang D, Tao C, Chen J, Wang J, YangE. The landscape of hervRNAs transcribed from human endogenous retroviruses across human body sites [J]. Genome Biol, 2022, 23(1):231.].

[0095] The results show that the method provided by the present invention has higher accuracy and sensitivity than the prior art, especially at the exon and site levels of teRNA, the present invention achieves the best balance between accuracy and sensitivity.

[0096] Comparative Example 2 Real Data Comparison

[0097] The third-generation full-length sequencing data downloaded from the GTEx (phs000424.v9.p2) and GEO (PRJNA635275) public data sets were used as the gold standard to evaluate the sensitivity of the present invention (Example 1) and the prior art in Comparative Example 1 for teRNA recognition. The comparison results are shown in Fig.12 .

[0098] It can be seen that the method provided by the present invention has a higher sensitivity in identifying teRNA in almost all samples.

[0099] In summary, the present invention provides a teRNA identification strategy that combines the Bayesian model, the parameterized assembly algorithm and the local site de novo assembly algorithm. By mutual correction of assembly results and quality control of secondary alignment of sequencing reads, the accuracy of teRNA identification is improved, providing a powerful tool for the precise screening and functional analysis of functional teRNA in the future.

[0100] 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 RNA derived from transposable elements based on RNA sequencing data, characterized in that: The following steps are involved: (1) Screening of teRNA sites; (2) Assembly of teRNA-derived sequencing reads; (3) Identification of teRNA; The screening of teRNA sites in step (1) comprises the following steps: (1.1) Obtaining sequencing data of double-ended RNA of the sample to be tested; (1.2) aligning the sequencing data to the reference genome to obtain an alignment result file; (1.3) allocating the sequencing read lengths of the sequencing data, and retaining teRNA sites with sequencing read lengths > 0; The assembly of the teRNA-derived sequencing reads in step (2) comprises the following steps: (2.1) extracting the alignment result corresponding to the teRNA site in the alignment result file to obtain a teRNA alignment result file; the teRNA site is the teRNA site screened in step (1); (2.2) Using the reference gene annotation file as a guide, the teRNA alignment result file is assembled with reference using StringTie software with the parameter -G to obtain the first assembled sequence; (2.3) Using the recognition module of SERVE software to perform local site de novo assembly on the teRNA alignment result file, the parameter --count is set to 1 to obtain the second assembly sequence; (2.4) combining the first assembled sequence and the second assembled sequence to obtain an assembled total sequence; The identification of teRNA in step (3) comprises the following steps: (3.1) extracting the single exon transcript sequence from the assembled total sequence, and merging the sequences using the merge module of BEDTools to obtain a de-redundant single exon transcript sequence; (3.2) extracting the exon sequence of the multi-exon transcript in the assembled total sequence, and aligning and merging it with the de-redundant single exon transcript sequence to obtain a de-redundant exon sequence; the alignment and merging method comprises: performing interval comparison between the exon sequence and the single exon transcript sequence, and for the exon sequence overlapping with the single exon transcript sequence, taking the interval union of the overlapping exon sequence; (3.3) extracting teRNA exons from the de-redundant exon sequence, and calculating coverage using the coverage module of BEDTools software; (3.4) Eliminate teRNA exons with sequencing coverage <85% to obtain qualitatively identified teRNAs.

2. The identification method according to claim 1, characterized in that: The identification method further comprises: (4) annotation of teRNA; (5) quantification of teRNA; The annotation of the teRNA in step (4) comprises the following steps: (4.1) using the annotation of the transposable element database as a reference, annotating the family and site of the qualitatively identified teRNA to obtain a first annotation result; (4.2) Annotating the chimeric gene or adjacent gene of the qualitatively identified teRNA with reference to the annotation of the gene database to obtain a second annotation result; (4.3) predicting the protein coding potential of the qualitatively identified teRNA to obtain a third annotation result; (4.4) merging the first annotation result, the second annotation result, and the third annotation result to obtain an annotation set; The quantification of teRNA in step (5) comprises the following steps: (5.1) quantifying the annotation set using RSEM software, standardizing the exon length, and obtaining an exon quantitative matrix; (5.2) using Telescope software to quantify the teRNA sites to obtain a site quantification matrix; the teRNA is the teRNA site screened in step (1); (5.3) Adding the expression levels of teRNA sites from the same family in the site quantitative matrix to obtain a family quantitative matrix.

3. The identification method according to claim 1, characterized in that: The method for obtaining the sequencing data includes: obtaining from a public database or performing sequencing; the public database includes one or more of GEO, GTEx, TCGA and CGGA.

4. The identification method according to claim 1, characterized in that: The method for allocating sequencing read lengths of the sequencing data includes a read length allocation algorithm of Telescope software.

5. The identification method according to claim 1, characterized in that: Tools for aligning the sequencing data to a reference genome include STAR software; the reference genome includes the human reference genome GRCh38 or T2T-CHM13.

6. The identification method according to claim 1, characterized in that: The tools extracted in step (2.1) include SAMtools software.

7. The identification method according to claim 2, characterized in that: The transposable element database includes the Dfam database.

8. The identification method according to claim 2, characterized in that: The gene database includes the GENCODE database.

9. The identification method according to claim 1, characterized in that: The tools for predicting the protein coding potential of the qualitatively identified teRNA include CPC2 software, CPAT software and PfamScan software.

10. The identification method according to claim 1, characterized in that: In step (5.1), the exon length is standardized using formula I. Among them, C e Indicates the count of exons, L e Indicates the length of the exon, C t Indicates the count of transcripts, L t Indicates the length of the transcript.

Images