A method for identifying transposable element-derived RNAs based on RNA sequencing data

Abstract

The present invention relates to the technical field of bioinformatics, and particularly to a method for identifying transposable element-derived RNAs based on RNA sequencing data. The present invention proposes to take the exon level of teRNAs as the analysis object. Compared with the teRNA locus, the teRNA exon not only achieves more accurate locus identification, but also retains the form of the TE in the transcript, including whether it provides splicing sites and polyadenylation signals, and whether it is chimeric with genes or other TE loci. And compared with the teRNA transcript level, as a component of the transcript, the teRNA exon significantly reduces the sequence assembly difficulty of short-read sequencing data. The present invention combines Bayesian models, reference-based assembly algorithms, and local de novo assembly algorithms, and through the mutual correction of different assembly algorithms and repeated alignment of sequencing reads before and after assembly, reduces the false positives of assembly, thereby improving the accuracy of teRNA identification.

Description

Technical Field

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

[0002] Transposable elements (TEs) are a class of special DNA elements that can move freely and are ubiquitously present in eukaryotic genomes. During the long-term evolution process, due to sequence variations such as mutations, deletions, and homologous recombination of their own sequences, as well as the evolution of defense mechanisms such as host epigenetic regulation and small RNA regulation, most TEs have lost their normal expression and transposition functions and become silent DNA sequences in the genome. Therefore, they were 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 initiate transcription under external stimuli 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 drive the inflammatory response in multiple sclerosis by activating the host immune system and inducing the expression of pro-inflammatory cytokines. Moreover, under physiological conditions, TEs also have the ability to transcribe in tissue cells, and their expression products even participate in multiple biological processes of growth and development. For example, the HERVW and HERVFRD families can promote cell fusion and syncytiotrophoblast formation by expressing syncytin proteins and are essential genes for placental morphogenesis. Another example is that HERVH is highly expressed in embryonic stem cells and pluripotent stem cells. It can shape chromatin structure by forming topological associated domain boundaries, regulate the expression of surrounding genes, and affect the cell differentiation process, thus participating in the maintenance of embryonic stemness. Therefore, systematically identifying TE-derived RNAs (teRNAs), analyzing their expression distribution patterns, and expression influencing factors will help to understand the role of teRNAs in the gene expression regulatory network and lay a foundation for revealing the functions of teRNAs in growth and development and disease occurrence and development.

[0003] Since TE is a multi-copy repetitive element, TEs within the same family have highly similar sequences, which are prone to multi-site alignment during sequence alignment, thus affecting the accuracy of teRNA identification. Therefore, at present, most teRNA studies adopt the identification strategy at the TE family level for teRNA expression analysis, that is, TEs with homologous sequences are grouped into a 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 TE families cannot explain whether the expression and function of TE sites within the same family are common, nor can it easily distinguish whether teRNA originates from specific TE sites, restricting the identification of functional teRNA sites and the analysis of related regulatory mechanisms. To achieve site-specific identification of TE, the most direct method is to directly remove multi-site alignment reads or randomly allocate them, but this will lead to a serious underestimation of the read count at some TE sites. In addition, the teRNA site-specific identification method represented by Telescope proposes to repeatedly reallocate multi-aligned reads by establishing an expectation-maximization mixture model until the best site is allocated. However, existing teRNA site expression identification methods directly identify based on genomic TE sites provided by the database, and not all TE sites are transcribed in full length with their own single sites as units. They may produce truncated teRNA, 3'-end read-through teRNA, or even form co-transcription units with surrounding genes or TE sites to generate chimeric teRNA. The teRNA identification method based on the transcript assembly strategy performs qualitative and quantitative analysis of teRNA by reconstructing transcripts, which is closer to the true biological expression of teRNA and is conducive to the analysis of the expression distribution pattern of teRNA and related regulatory mechanisms. However, limited by short-read sequencing, the accuracy of teRNA transcript assembly is still relatively low at present. Summary of the Invention

[0004] To solve the above problems, the present invention provides a method for identifying transposable element-derived RNA based on RNA sequencing data. The identification method provided by the present invention is a teRNA recognition strategy that combines a Bayesian model, a reference-based assembly algorithm, and a local de novo assembly algorithm. By mutually correcting the assembly results and performing quality control on the secondary alignment of sequencing reads, the accuracy of teRNA recognition is improved, providing a powerful tool for the precise screening and functional analysis of future functional teRNAs.

[0005] To achieve the above object, the present invention provides the following technical solutions:

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

[0007] The screening of the teRNA locus described in step (1) includes the following steps: (1.1) Obtain the sequencing data of the double-stranded RNA of the sample to be tested; (1.2) Align the sequencing data to the reference genome to obtain an alignment result file; (1.3) Allocate the sequencing read lengths of the sequencing data, and retain the teRNA loci with the number of sequencing read lengths > 0.

[0008] The assembly of the teRNA-derived sequencing read lengths described in step (2) includes the following steps: (2.1) Extract the alignment results corresponding to the teRNA loci in the alignment result file to obtain a teRNA alignment result file; the teRNA loci are the teRNA loci screened in step (1); (2.2) Guided by the reference gene annotation file, use the StringTie software to perform parametric assembly on the teRNA alignment result file with the parameter -G to obtain the first assembled sequence; (2.3) Use the recognition module of the SERVE software to perform local locus de novo assembly on the teRNA alignment result file with the parameter --count set to 1 to obtain the second assembled sequence; (2.4) Merge the first assembled sequence and the second assembled sequence to obtain the total assembled sequence.

[0009] The identification of the teRNA described in step (3) includes the following steps: (3.1) Extract the single-exon transcript sequences in the total assembled sequence, and use the merge module of BEDTools to merge the sequences to obtain the non-redundant single-exon transcript sequences; (3.2) Extract the exon sequences of the multi-exon transcripts in the total assembled sequence, and perform alignment and merging with the non-redundant single-exon transcript sequences to obtain the non-redundant exon sequences; the method of alignment and merging includes: performing interval comparison on the exon sequences and the single-exon transcript sequences, and for the exon sequences overlapping with the single-exon transcript sequences, taking the union of the intervals of the overlapping exon sequences; (3.3) Extract the teRNA exons in the non-redundant exon sequences, and use the coverage module of the BEDTools software to calculate the coverage; (3.4) Eliminate the teRNA exons with a sequencing coverage < 85% to obtain the qualitatively identified teRNA.

[0010] Preferably, the identification method further includes: (4) Annotation of teRNA; (5) Quantification of teRNA.

[0011] The annotation of the teRNA described in step (4) includes the following steps: (4.1) Using the transposable element database annotation as a reference, perform family and locus annotation on the qualitatively identified teRNA to obtain a first annotation result; (4.2) Using the gene database annotation as a reference, perform chimeric gene or adjacent gene annotation on the qualitatively identified teRNA to obtain a second annotation result; (4.3) Perform protein-coding potential prediction on the qualitatively identified teRNA to obtain a third annotation result; (4.4) Combine the first annotation result, the second annotation result, and the third annotation result to obtain an annotation set;

[0012] The quantification of the teRNA described in step (5) includes the following steps: (5.1) Use the RSEM software to quantify the annotation set and standardize the exon lengths to obtain an exon quantification matrix; (5.2) Use the Telescope software to perform locus quantification of the teRNA to obtain a locus quantification matrix; the teRNA is the teRNA locus screened in step (1); (5.3) Add up the expression levels of the teRNA loci from the same family in the locus quantification matrix to obtain a family quantification matrix.

[0013] Preferably, 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.

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

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

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

[0017] Preferably, the transposable element database includes 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 the CPC2 software, the CPAT software, and the PfamScan software.

[0020] Preferably, in step (5.1), the exon lengths are standardized 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 take the exon level of transposable element-derived RNA (teRNA) as the analysis object. Compared with the teRNA locus, the teRNA exon not only achieves more accurate locus identification, but also retains the form of the transposable element (TE) in the transcript, including whether it provides splicing sites and polyadenylation signals, and whether it is chimeric with genes or other TE loci. And compared with the teRNA transcript level, as a component of the transcript, the teRNA exon significantly reduces the sequence assembly difficulty 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 reference-based assembly algorithm and the local site de novo assembly algorithm. By mutual correction of different assembly algorithms and repeatedly aligning the sequencing reads before and after assembly, the false positives of the assembly are reduced, thereby improving the accuracy of teRNA recognition.

[0026] In addition, the existing technologies mainly only conduct identification and analysis on 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 recognition, the present invention also provides quantitative methods for multiple research levels of teRNA (including family, locus, exon), providing a convenient, unified and accurate method for future research. Description of the drawings

[0027] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required to be used in the embodiments.

[0028] Figure 1 It is a schematic diagram of the teRNA recognition algorithm;

[0029] Figure 2 It is a flowchart of the technical steps of the present invention;

[0030] Figure 3 It is an electrophoresis diagram for verifying teRNA exons by RT-PCR experiment;

[0031] Figure 4 It is the result for verifying the accuracy of teRNA quantification by RT-qPCR experiment;

[0032] Figure 5Results of dimensionality reduction analysis for glioblastoma multiforme and para-cancer samples;

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

[0034] Figure 7 Results of enrichment analysis of teRNA-related genes;

[0035] Figure 8 Two examples of the association between teRNA and poor prognosis in glioblastoma multiforme patients;

[0036] Figure 9 Results of dimensionality reduction analysis based on the teRNA exon expression profile of cell samples during neuronal differentiation;

[0037] Figure 10 Results of grouping the characteristic genes of 14 co-expression modules according to the dynamic expression pattern;

[0038] Figure 11 Results of performance comparison between the method of the present invention and the prior art in the simulated dataset;

[0039] Figure 12 Results of performance comparison between the method of the present invention and the prior art in the real public dataset. Detailed implementation manners

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

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

[0042] (1.1) Obtaining the sequencing data of the paired-end 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) Assigning the sequencing reads of the sequencing data and retaining the teRNA sites with the number of sequencing reads > 0.

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

[0044] (2.1) Extract the alignment results corresponding to the teRNA sites from the alignment result file to obtain a teRNA alignment result file; the teRNA sites are the teRNA sites screened in step (1); (2.2) Guided by the reference gene annotation file, use the StringTie software to perform reference-based assembly on the teRNA alignment result file with the parameter -G to obtain a first assembled sequence; (2.3) Use the recognition module of the 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 assembled sequence; (2.4) Merge the first assembled sequence and the second assembled sequence to obtain an overall assembled sequence;

[0045] The recognition of teRNA in step (3) includes the following steps:

[0046] (3.1) Extract the single-exon transcript sequences from the overall assembled sequence, and use the merge module of BEDTools to merge the sequences to obtain non-redundant single-exon transcript sequences; (3.2) Extract the exon sequences of the multi-exon transcripts from the overall assembled sequence, and perform alignment and merging with the non-redundant single-exon transcript sequences to obtain non-redundant exon sequences; the method of alignment and merging includes: performing interval comparison between the exon sequences and the single-exon transcript sequences, and for the exon sequences overlapping with the single-exon transcript sequences, taking the union of the intervals of the overlapping exon sequences; (3.3) Extract the teRNA exons from the non-redundant exon sequences, and use the coverage module of the BEDTools software to calculate the coverage; (3.4) Eliminate the teRNA exons with a sequencing coverage <85% to obtain qualitatively recognized teRNA.

[0047] The identification method provided by the present invention is a teRNA recognition strategy that combines a Bayesian model, a reference-based assembly algorithm, and a local site de novo assembly algorithm. By mutually correcting the assembly results and performing quality control on the secondary alignment of sequencing reads, the accuracy of teRNA recognition is improved, providing a powerful tool for the precise screening and functional analysis of future functional teRNA.

[0048] As an implementation, the identification method further includes: (4) Annotation of teRNA; (5) Quantification of teRNA.

[0049] The annotation of the teRNA described in step (4) includes the following steps: (4.1) Using the transposable element database annotation as a reference, perform family and locus annotation on the qualitatively identified teRNA to obtain a first annotation result; (4.2) Using the gene database annotation as a reference, perform chimeric gene or neighboring gene annotation on the qualitatively identified teRNA to obtain a second annotation result; (4.3) Predict the protein-coding potential of the qualitatively identified teRNA to obtain a third annotation result; (4.4) Combine the first annotation result, the second annotation result, and the third annotation result to obtain an annotation set;

[0050] The quantification of the teRNA described in step (5) includes the following steps: (5.1) Use the RSEM software to quantify the annotation set and standardize the exon lengths to obtain an exon quantification matrix; (5.2) Use the Telescope software to perform locus quantification of the teRNA to obtain a locus quantification matrix; the teRNA is the teRNA locus screened in step (1); (5.3) Add up the expression levels of the teRNA loci from the same family in the locus quantification matrix to obtain a family quantification matrix.

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

[0052] As an implementation method, 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 implementation method, the format of the sequencing data is the FASTQ format.

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

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

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

[0056]

[0057] where 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 transposable element-derived RNA based on RNA sequencing data provided by the present invention will be described in detail below with reference to the accompanying drawings and embodiments, but they should not be construed as limiting the protection scope of the present invention.

[0059] Example 1

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

[0061] (1) Screening steps for teRNA sites: (1.1) Perform paired-end RNA sequencing on the HEK293T and U251 cell lines respectively to obtain the original sequencing data; (1.2) Use the STAR software to align the original 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 the Telescope software to allocate the sequencing read lengths; (1.4) Extract the teRNA sites covered by the sequencing read lengths (retain the teRNA sites with the sequencing read count > 0).

[0062] (2) Assembly steps of teRNA-derived sequencing reads: (2.1) According to the teRNA sites screened in step (1), use the SAMtools software to extract the alignment result file (BAM format) of teRNA; (2.2) Guided by the human reference gene annotation file (T2T-CHM13 version), use StringTie for reference-based assembly (parameter -G) to obtain the assembled sequence (GTF format); (2.3) Use the recognition module of the SERVE software for local site de novo assembly, set the parameter --count to 1, and obtain the assembled sequence (GTF format); (2.4) Merge the assembled sequences of the previous two steps.

[0063] (3) Recognition steps of teRNA: (3.1) Extract the single-exon transcript sequences of the assembled sequences in step (2), and use the merge module of BEDTools to merge the sequences to obtain the non-redundant single-exon transcript sequences; (3.2) Extract the exon sequences of the multi-exon transcripts of the assembled sequences in step (2), and through interval comparison with the non-redundant single-exon transcript sequences obtained in step (3.2), for the exon sequences overlapping with the single-exon transcripts, take the union of the overlapping exon intervals to obtain the non-redundant exon sequences; (3.3) Extract teRNA exons, and use the coverage module of BEDTools to calculate the coverage; (3.4) Remove teRNA exons with a sequencing coverage < 85%.

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

[0065] Table 1 Amplification primers for different teRNAs

[0066]

[0067]

[0068] The results are shown in Figure 3 , and all the selected teRNAs are experimentally proven to be true positives, indicating the accuracy of the qualitative recognition of teRNA by the present invention.

[0069] Example 2

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

[0071] Based on the teRNA qualitatively recognized in Example 1, perform the following steps:

[0072] (4) Annotation of teRNA includes the following steps: (4.1) Using BEDTools for TE family and locus annotation with the TE loci annotated in the Dfam database as a reference; (4.2) Using BEDTools to annotate teRNA chimeric genes and teRNA neighboring genes with the GENCODE gene annotation as a reference; (4.3) Using CPC2, CPAT, and PfamScan to predict the protein-coding potential of teRNA.

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

[0074]

[0075] where Ce represents the count of the exon, Le represents the length of the exon, Ct represents the count of the transcript, and Lt represents the length of the transcript; (5.2) Using the Telescope software to perform locus quantification of teRNA to obtain a locus quantification matrix; (5.3) Summing up the expression levels of teRNA loci from the same family to obtain a family quantification matrix.

[0076] Randomly select teRNA for quantitative experimental verification, and the primer sequences are shown in Table 2.

[0077] Table 2 Amplification primers for different teRNAs

[0078]

[0079] The results are shown in Figure 4 , and the quantification results of the present invention show a high degree of consistency with the real experimental results, proving the accuracy of the teRNA quantification of the present invention.

[0080] Example 3 Application to tumor data

[0081] Taking the glioblastoma public dataset (ID: PRJNA613939) as an example, applying the full analysis process of teRNA of the present invention to achieve a comprehensive analysis of tumor teRNA, the steps are as follows:

[0082] For each sample, perform analysis using the technical steps of Example 1 to obtain the teRNA annotation of each sample; then, use the Cuffmerge software to merge the samples to obtain a non-redundant teRNA annotation file; and then perform teRNA quantification on each sample using the technical steps of Example 2.

[0083] (6) Downstream analysis includes the following steps: (6.1) Perform dimensionality reduction analysis using the expression levels of teRNA at different research levels (the cmdscale function in R language). The results of dimensionality reduction analysis for glioblastoma and adjacent cancer samples are shown in Figure 5 , indicating that teRNA at different research levels can clearly distinguish tumor and adjacent cancer samples. Among them, the teRNA exon level also reflects the heterogeneity of tumor samples; (6.2) Perform differential expression analysis of teRNA using the edgeR software. The results are shown in Figure 6 , indicating that there are obvious differences in the expression profiles of teRNA exons from the same family; (6.3) Perform enrichment analysis of teRNA-related genes using clusterProfiler. The results are shown in Figure 7 , indicating that teRNA in glioblastoma may trigger the occurrence of the body's immunity; (6.4) According to the clinical information of the patients, perform survival analysis. The analysis results of two examples related to the poor prognosis of teRNA and glioblastoma patients are shown in Figure 8 , indicating that some teRNA are significantly associated with the poor prognosis of tumor patients.

[0084] Example 4 Application of normal cell differentiation data

[0085] Taking the public dataset of neuron differentiation (ID: PRJNA596331) as an example, this dataset includes various stages of neuron differentiation (early differentiated cells, neural progenitor cells, neural progenitor cells with rosettes, and mature neuron cells), as well as undifferentiated neural stem cells and purified neuron cells. Apply the full set of analysis processes of teRNA of the present invention (the analysis steps are the same as in Example 3). A total of 52,838 expressed teRNA exons were identified, covering 1,075 TE families, indicating the extensive expression of teRNA during neuron differentiation.

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

[0087] (7.1) Perform dimensionality reduction analysis using the expression levels of teRNA exons. The results of dimensionality reduction analysis based on the teRNA exon expression profile for cell samples during neuron differentiation are shown in Figure 9 , and the results show that the samples are clustered according to cell type, suggesting that teRNA may be closely related to the process of neuron cell differentiation;

[0088] (7.2) Use the WGCNA software to construct a co-expression regulatory network of teRNA exons, and detect 14 co-expression modules of teRNA exons. According to the expression of module characteristic genes at different differentiation times, such as Figure 10, the modules can be divided into 6 categories: the modules highly expressed in neural stem cells (Modules 7 and 8), the modules highly expressed in early differentiated cells (Modules 6, 12, 13, 14), the modules highly expressed in neural progenitor cells (Module 5), the modules highly expressed in rosettes (Modules 2, 3, 10), the modules lowly expressed in neural progenitor cells (Modules 9 and 11), and the modules highly expressed in mature neurons (Module 14), indicating that the teRNA exons in different modules show specific dynamic expression at different differentiation times.

[0089] Comparison of simulation data in Comparative Example 1

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

[0091] StringTie_unguided and StringTie_guided refer to the literature [Pertea M, Pertea G M, Antonescu C M, Chang T C, Mendell J T, 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 refer to the literature [Trapnell C, Williams B A, Pertea G, Mortazavi A, Kwan G, van Baren M J, Salzberg S L, Wold B J, Pachter L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation[J]. Nat Biotechnol, 2010, 28(5): 511 - 515.];

[0093] LIONS See reference

Babaian A, Thompson I R, Lever J, Gagnier L, Karimi M M, Mager D L. LIONS: analysis suite for detecting and quantifying transposable element initiated transcription from RNA-seq[J]. Bioinformatics, 2019, 35(19): 3839-3841.

[0094] SERVE See reference

She J, Du M, Xu Z, Jin Y, Zhang D, Tao C, Chen J, Wang J, Yang E. The landscape of hervRNAs transcribed from human endogenous retroviruses across human body sites[J]. Genome Biol, 2022, 23(1): 231.

[0095] As can be seen from the results, the method provided by the present invention is more accurate and sensitive than the prior art. Especially at the exon and site levels of teRNA, the present invention achieves the best balance between precision and sensitivity.

[0096] Comparison of real data in Comparative Example 2

[0097] Using the third-generation full-length sequencing data downloaded from the GTEx (phs000424.v9.p2) and GEO (PRJNA635275) public datasets as the gold standard, the sensitivity of the present invention (Example 1) and the prior art in Comparative Example 1 for teRNA recognition was evaluated. The comparison results are shown in Figure 12 .

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

[0099] In summary, the present invention provides a teRNA recognition strategy that combines a Bayesian model, a reference-based assembly algorithm, and a local site de novo assembly algorithm. By mutually correcting the assembly results and performing quality control on the secondary alignment of sequencing reads, the accuracy of teRNA recognition is improved, providing a powerful tool for the precise screening and functional analysis of future functional teRNAs.

[0100] Although the above embodiments have described the present invention in detail, they are only a part of the embodiments of the present invention, rather than all embodiments. People can also obtain other embodiments based on this embodiment without creative efforts, and these embodiments all fall within the protection scope of the present invention.

Claims

1. A method for identifying transposable element-derived RNAs based on RNA sequencing data, characterized in that, Including the following steps: (1) Screening of teRNA sites; (2) Assembly of teRNA-derived sequencing reads; (3) Identification of teRNA; The screening of the teRNA sites described in step (1) includes the following steps: (1.1) Obtaining the sequencing data of the double-stranded 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 reads of the sequencing data and retaining the teRNA sites with the number of sequencing reads > 0; The assembly of the teRNA-derived sequencing reads described in step (2) includes 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 StringTie software to perform parameterized assembly on the teRNA alignment result file with reference to the reference gene annotation file, with the parameter -G, to obtain a first assembled sequence; (2.3) Using the identification module of the 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 assembled sequence; (2.4) Merging the first assembled sequence and the second assembled sequence to obtain an assembled total sequence; The identification of teRNA described in step (3) includes the following steps: (3.1) Extracting the single-exon transcript sequences in the assembled total sequence and using the merge module of BEDTools to merge the sequences to obtain non-redundant single-exon transcript sequences; (3.2) Extracting the exon sequences of the multi-exon transcripts in the assembled total sequence and performing alignment and merging with the non-redundant single-exon transcript sequences to obtain non-redundant exon sequences; the method of alignment and merging includes: performing interval comparison on the exon sequences and the single-exon transcript sequences, and for the exon sequences overlapping with the single-exon transcript sequences, taking the union of the overlapping exon sequence intervals; (3.3) Extracting the teRNA exons in the non-redundant exon sequences and using the coverage module of the BEDTools software to calculate the coverage; (3.4) Removing the teRNA exons with a sequencing coverage < 85% to obtain qualitatively identified teRNA.

2. The identification method according to claim 1, wherein The identification method further includes: (4) Annotation of teRNA; (5) Quantification of teRNA; The annotation of teRNA described in step (4) includes the following steps: (4.1) Taking the transposable element database annotation as a reference, performing family and site annotation on the qualitatively identified teRNA to obtain a first annotation result; (4.2) Taking the gene database annotation as a reference, performing chimeric gene or adjacent gene annotation on the qualitatively identified teRNA to obtain a second annotation result; (4.3) Performing protein-coding potential prediction on the qualitatively identified teRNA to obtain a third annotation result; (4.4) Combine the first annotation result, the second annotation result, and the third annotation result to obtain an annotation set; The quantification of the teRNA in step (5) includes the following steps: (5.1) Use the RSEM software to quantify the annotation set, standardize the exon lengths, and obtain an exon quantification matrix; (5.2) Use the Telescope software to perform site quantification of the teRNA to obtain a site quantification matrix; the teRNA is the teRNA site screened in step (1); (5.3) Add up the expression levels of the teRNA sites from the same family in the site quantification matrix to obtain a family quantification 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, wherein The method for allocating the sequencing reads of the sequencing data includes the read length allocation algorithm of the Telescope software.

5. The identification method according to claim 1, wherein The tool for aligning the sequencing data to the reference genome includes the 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 tool for extraction in step (2.1) includes the SAMtools software.

7. The identification method according to claim 2, wherein 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 the CPC2 software, the CPAT software, and the PfamScan software.

10. The identification method according to claim 1, characterized in that, In step (5.1), the exon lengths are standardized using Equation Ⅰ, Among them, C e represents the count of exons, and L e represents the length of exons. C t represents the count of transcripts, and L t represents the length of transcripts.

Images