A method and apparatus for filtering rRNA sequences in transcriptome sequencing data

Abstract

The application discloses a method for establishing a species rRNA sequence database, and belongs to the technical field of transcriptome sequencing data processing, and comprises the following steps: obtaining first transcriptome sequencing data and reference genome data of a species; comparing the transcriptome data with the reference genome data to obtain unaligned sequences, selecting the first M sequences of the unaligned sequences, and performing comparison with a NCBI database, wherein M=8000-20000; obtaining rRNA sequences on the comparison, and establishing an rRNA sequence database. Furthermore, the application also discloses a method and device for filtering rRNA sequences in transcriptome sequencing data based on the established rRNA sequence database, which can filter rRNA sequences in a targeted manner, quickly and accurately obtain pure mRNA data, and perform subsequent analysis.

Description

Technical Field

[0001] This invention belongs to the field of transcriptome sequencing data processing technology, specifically, it relates to a method and apparatus for filtering rRNA sequences in transcriptome sequencing data. Background Technology

[0002] The transcriptome is the sum of all RNA transcribed from a species in a specific tissue or cell at a particular developmental stage or functional state. It mainly includes mRNA and non-coding RNA (ncRNA). Transcriptome research is the foundation and starting point for gene function and structure research. High-throughput sequencing allows for the comprehensive and rapid acquisition of almost all transcriptome sequence information, enabling the study of gene function and structure at a holistic level, revealing specific biological processes. It has been widely applied in fields such as plant candidate gene discovery, functional identification, and genetic improvement.

[0003] Currently, RNA sequencing (RNA-Seq) technology has become one of the important tools in transcriptomics research. It is low-cost, has few limitations, and is highly accurate. Different fragment sizes are typically selected during library preparation depending on the length of the target RNA, resulting in varying read lengths. Generally, for microRNA sequencing, the microRNA is isolated and sequenced separately. For mRNA sequencing, a 200-300 bp fragment is usually selected during library preparation, and sequencing is performed using 125PE / 150PE. Because long non-coding RNAs (lncRNAs) undergo both forward and reverse transcription, strand-specific library preparation and sequencing are often employed.

[0004] The raw data from sequencing transcriptome data consists of short sequences (also called reads) of approximately 150 bp. These short sequences cannot be directly used for data analysis. To ensure accurate and reliable analysis results, the raw data needs to be preprocessed, including removing sequencing adapters (introduced during library construction) and low-quality sequencing data (caused by errors in the sequencer itself). The resulting valid data is then used for alignment with the species' reference genome.

[0005] The above is a common procedure in transcriptome analysis. However, in prokaryotes, mRNA accounts for only 1-5% of all RNA, with the vast majority being ribosomal RNA (rRNA). Therefore, to sequence mRNA, it must first be purified. However, prokaryotic mRNA does not have the polyA structure of eukaryotic mRNA, so it cannot be directly purified using oligoT. If total RNA is used for sequencing directly, the sequencing efficiency will be very poor because a large proportion of the sequenced sequence will align with the reference genome as ribosomes, mycoplasma, or mitochondria, severely affecting alignment efficiency. Therefore, before aligning with the reference genome, ribosome, mycoplasma, and mitochondrial sequences need to be removed to obtain a pure mRNA sequence for further analysis.

[0006] Currently, the common approach in this field is to construct a ribosome sequence library for all species, which is then used for filtering in the analysis of reference transcriptomes from different species. However, existing databases have the following drawbacks:

[0007] 1. Although the database contains a large number of species, it is still not possible to guarantee that all ribosome sequences are included;

[0008] 2. Since the database selects ribosome sequences from all species, sequences that do not belong to ribosome sequences will be filtered out during the alignment process due to sequence similarity.

[0009] Therefore, there is an urgent need in this field for a rapid and accurate method to filter ribosome sequences in transcriptome sequencing data. Summary of the Invention

[0010] To solve the above-mentioned technical problems, the technical solution adopted by the present invention is as follows:

[0011] The first aspect of this invention provides a method for establishing a species rRNA sequence database, comprising the following steps:

[0012] S1, obtain the first transcriptome sequencing data and reference genome data of the species;

[0013] S2, compares the transcriptome data with the reference genome data to obtain the unaligned sequences.

[0014] S3, Select the first M uncompared sequences and compare them with the NCBI database, where M = 8000 to 20000;

[0015] S4. Obtain the aligned rRNA sequences and establish an rRNA sequence database.

[0016] In some embodiments of the present invention, after obtaining the transcriptome data of the species in step S1, the step further includes a step of preprocessing the transcriptome data: (1) filtering adapter reads; (2) filtering reads containing more than 5% of N (unknown bases); (3) filtering reads in which the number of bases with a quality value of Q ≤ 10 accounts for more than 20% of the total reads.

[0017] Furthermore, when constructing the library from the first transcriptome sequencing data, strand-specific library construction is used. Transcriptome sequencing can retain the orientation information of transcripts during transcriptome sequencing, thus it can also determine whether the transcripts originate from the positive or negative strand of the genome.

[0018] In some embodiments of the present invention, in step S2, the hisat2 comparison software is used for comparison.

[0019] In some embodiments of the present invention, step S3, comparing with the NCBI database, refers to uploading the unmatched sequence to NCBI for comparison. In other embodiments of the present invention, those skilled in the art can also download the NCBI database locally for comparison to improve comparison efficiency. Of course, those skilled in the art can also compare with other databases to enrich the rRNA sequences that can be matched. In some specific embodiments of the present invention, comparing with the NCBI database refers to comparing with the NT database within the NCBI database.

[0020] In this invention, M can also be reasonably selected based on the actual number of sequences.

[0021] Further, in step S4, "alignment" means that the alignment E value is less than a preset threshold. In some embodiments of the present invention, the preset threshold is set to 1e-5, that is, if the e value in the NCBI alignment result is less than 1e-5, it is considered an alignment.

[0022] A second aspect of the present invention provides an rRNA sequence database established based on first transcriptome sequencing data of a species using any of the methods described in the first aspect of the present invention.

[0023] Furthermore, those skilled in the art can also use multiple transcriptome sequencing data of this species to obtain a more comprehensive rRNA sequence database through the same methods described above.

[0024] A third aspect of the present invention provides a method for filtering rRNA sequences in second transcriptome sequencing data of a species, comprising the following steps:

[0025] The second transcriptome sequencing data is compared with the rRNA sequence database described in the second aspect of the present invention, and sequences that can be matched are filtered out.

[0026] The purpose of this invention in establishing an rRNA sequence database for a species is to enable rapid filtering or removal of rRNA sequences from other transcriptome data of that species. Therefore, the second transcriptome sequencing data here refers to any transcriptome data of that species from which rRNA sequences are to be removed, and of course, also includes the first transcriptome sequencing data used to establish the rRNA sequence database itself. Since only a portion of the sequences not aligned with the reference genome is used for alignment, the filtering efficiency can be greatly improved.

[0027] In one embodiment of the present invention, bowtie2 is used for alignment, and the alignment parameters are set as follows:

[0028] -I Minimum fragment length uses the default parameter 0;

[0029] -X Maximum fragment length uses the default parameter of 500;

[0030] Select the `--sensitive-local` parameter in the `--loca` mode, which is `-D 15-R 2-N 0-L 20-i S,1,0.75`. Specific parameter descriptions: `-D 15` extends a seed sequence during alignment. If no better or second-best alignment result is obtained, the alignment fails. Alignment ends after 15 consecutive failures. `-R2` determines the number of times a seed generated from a read matches too many sites on the reference sequence. If each seed matches more than 300 sites on average, a new seed is generated using a different offset. `-N 0` sets the allowed number of mismatches during seed alignment to 0. `-L 20` sets the seed length to 20. `-i` sets the distance between two adjacent seeds. This parameter is calculated using a formula and consists of three parts: (a) the calculation method, including constant (C), linear (L), square root (S), and natural logarithm (G); (b) a constant; and (c) a coefficient. S,1,0.75 is equivalent to f(x)=1+0.75*.

[0031] The -loca mode performs partial alignment of the read, omitting some bases at the ends of the read to ensure the alignment score meets the requirements. The penalty parameter –ma defaults to 2 in this mode.

[0032] A fourth aspect of the present invention provides an apparatus for filtering rRNA sequences in second transcriptome sequencing data of a species, comprising:

[0033] The data input module is used to obtain the second transcriptome sequencing data;

[0034] A database storage module is used to store the rRNA sequence database described in the second aspect of this invention;

[0035] The filtering module is connected to both the data input module and the database storage module, and is used to compare the second transcriptome sequencing data with the rRNA sequence database, filter out the sequences that can be matched, and output the results.

[0036] A fifth aspect of the present invention provides a computer device, comprising:

[0037] Memory, used to store computer programs;

[0038] A processor for executing the computer program to implement the steps of the method as described in any of the third aspects of the invention.

[0039] A sixth aspect of the present invention provides a computer-readable storage medium storing a computer program that, when executed by a processor, implements the steps of any of the methods described in the third aspect of the present invention.

[0040] Beneficial effects of the present invention

[0041] Compared with the prior art, the present invention has the following beneficial effects:

[0042] This invention first provides a method for rapidly constructing an rRNA sequence database for a species by name. Then, after preprocessing, transcriptome sequencing data is used to selectively filter rRNA sequences from this database. This allows for the rapid and accurate acquisition of clean mRNA data for subsequent analysis. Attached Figure Description

[0043] Figure 1 The results of comparing the partial reads in the mouse species of Example 2 of the present invention that were not aligned with the reference genome with the NCBI database are shown.

[0044] Figure 2 The results are shown after filtering using the rRNA sequence database constructed using this invention.

[0045] Figure 3 The results of comparing the partial reads of the Muscovy duck species in Example 3 of this invention, which were not aligned with the reference genome, with the NCBI database are shown.

[0046] Figure 4 The results are shown after filtering using the rRNA sequence database constructed in Example 3 of this invention.

[0047] Figure 5 A schematic diagram of a device for filtering rRNA sequences in transcriptome sequencing data according to an embodiment of the present invention is shown. Detailed Implementation

[0048] Unless otherwise stated, implied from the context, or as is customary in the art, all parts and percentages in this application are based on weight, and all testing and characterization methods used are concurrent with the filing date of this application. Where applicable, any patent, patent application, or disclosure relating to this application is incorporated herein by reference in its entirety, and its equivalent patent families are also incorporated herein by reference, particularly the definitions disclosed in such documents concerning the art. If any definition of a specific term disclosed in the prior art is inconsistent with any definition provided in this application, the definition provided in this application shall prevail.

[0049] The numerical ranges used in this application are approximate values ​​and therefore may include values ​​outside the range unless otherwise stated. The numerical range includes all values ​​from the lower limit to the upper limit, increasing by one unit, provided that there is an interval of at least two units between any lower and any higher value. These are merely specific examples of what is intended to be expressed, and all possible combinations of values ​​between the listed minimum and maximum values ​​are considered to be clearly stated in this application.

[0050] To make the technical problems solved by the present invention, the technical solutions and the beneficial effects of the present invention clearer, the present invention will be further described in detail below with reference to the embodiments.

[0051] Example

[0052] The following examples are used to illustrate preferred embodiments of the invention. Those skilled in the art will understand that the techniques disclosed in the examples represent techniques discovered by the inventors that can be used to implement the invention, and therefore can be considered preferred embodiments for implementing the invention. However, those skilled in the art should understand from this specification that many modifications can be made to the specific embodiments disclosed herein, still yielding the same or similar results, without departing from the spirit or scope of the invention.

[0053] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains, and all materials publicly cited herein and referenced by them are incorporated herein by reference.

[0054] Those skilled in the art will recognize, or can learn through routine experimentation, many equivalents of specific embodiments of the invention described herein. These equivalents will be included in the claims.

[0055] Unless otherwise specified, the experimental methods used in the following embodiments are conventional methods. Unless otherwise specified, the instruments and equipment used in the following embodiments are conventional laboratory instruments and equipment; unless otherwise specified, the experimental materials used in the following embodiments were purchased from conventional biochemical reagent stores.

[0056] Example 1: The reference genome of the source species was derived from transcriptome sequencing data from Ensembl, filtered by ribosome sequence.

[0057] Ensembl is a collaborative project between the Wellcome Trust at the Sanger Institute (WTSI) in the UK and the European Institute for Bioinformatics (EMBI-EBI), a division of the European Laboratory for Molecular Biology. Data obtained from the Ensembl project can be used to support comparative genomics, evolution, sequence mutation, and transcriptional regulation research.

[0058] The Ensembl database contains rRNA databases for multiple species, which can be downloaded directly for comparison. If a match is found, the rRNA can be filtered.

[0059] The inventors have established a method for downloading transcriptome data based on the species from which it originates. If the source species is a plant, download from ftp: / / ftp.ensemblgenomes.org / pub / release-version / plants / fasta / species / ncrna / *ncrna.fa.gz; if the source species is a fungus, download from ftp: / / ftp.ensemblgenomes.org / pub / release-version / fungi / fasta / species / ncrna / *ncrna.fa.gz; and if the source species is an animal, download from ftp: / / ftp.ensembl.org / pub / release-version / fasta / species / ncrna / *ncrna.fa.gz. Of course, if the URL is updated, the download address will also be updated accordingly.

[0060] If the reference genome of the source species for the transcriptome sequencing data comes from the Ensembl database, then the downloaded rRNA database can be used directly for alignment and filtering. However, for transcriptome sequencing data whose reference genome does not come from the Ensembl database, a new method needs to be developed for rRNA sequence filtering.

[0061] Example 2: Reference genome of the source species; ribosome sequence filtering of mouse transcriptome sequencing data not derived from Ensembl.

[0062] This embodiment provides a method for filtering rRNA sequences from mouse transcriptome sequencing data that are not derived from Ensembl and are used as a reference genome of the source species.

[0063] Specifically, transcriptome data of mice were obtained (using a strand-specific library construction method), with the reference genome for this species sourced from the NCBI database.

[0064] 1. Transcriptome sequencing data preprocessing

[0065] (1) Reads of the filter adapter;

[0066] (2) Filter reads containing more than 5% N (N indicates that the base information cannot be determined);

[0067] (3) Filter out low-quality reads (the number of bases with a quality value of Q≤10 accounts for more than 20% of the total read).

[0068] High-quality transcriptome sequencing data (Clean data) were obtained after preprocessing.

[0069] 2. First comparison

[0070] The clean data was aligned with the reference genome using the hisat2 alignment software. The alignment results and files containing unmapped sequences (unmapped.bam files) were obtained. The bamToFastqs software was then used to convert the bam files to fastq files (unmapped.fastq).

[0071] 3. Extraction

[0072] The converted FastQ file contains some reads with high abundance and others with low abundance, representing a small proportion of the total reads. Therefore, the inventors performed abundance statistics on the unmapped.fastq file and sorted it, resulting in the abundance-sorted unmapped.uniq.fa file. In the converted unmapped.uniq.fa file, the most abundant reads are listed first. To improve alignment efficiency, the inventors extracted the top 10,000 most abundant reads for subsequent alignments.

[0073] 4. Second comparison

[0074] The FASTA file containing 10,000 reads was uploaded to NCBI for BLAST comparison, and the results were as follows. Figure 1 As shown. After obtaining the alignment results, filter the entries that match the rRNA, select and download the FASTA files of all rRNAs, and establish an rRNA database.

[0075] 5. Filtration

[0076] The paired-end rRNA sequences of Cleandata were filtered using the rRNA database established in step 4. Specifically, version 2.2.0 of the alignment software bowtie2 was used to align the Cleandata with the rRNA database. The alignment method was selected as fr, and the alignment parameters were set as follows:

[0077] -I Minimum fragment length uses the default parameter 0;

[0078] -X Maximum fragment length uses the default parameter of 500;

[0079] Select the `--sensitive-local` parameter in the `--loca` mode, which is `-D 15-R 2-N 0-L 20-i S,1,0.75`. Specific parameter descriptions: `-D 15` extends a seed sequence during alignment. If no better or second-best alignment result is obtained, the alignment fails. Alignment ends after 15 consecutive failures. `-R2` determines the number of times a seed generated from a read matches too many sites on the reference sequence. If each seed matches more than 300 sites on average, a new seed is generated using a different offset. `-N 0` sets the allowed number of mismatches during seed alignment to 0. `-L 20` sets the seed length to 20. `-i` sets the distance between two adjacent seeds. This parameter is calculated using a formula and consists of three parts: (a) the calculation method, including constant (C), linear (L), square root (S), and natural logarithm (G); (b) a constant; and (c) a coefficient. S,1,0.75 is equivalent to f(x)=1+0.75*.

[0080] The -loca mode performs partial alignment of the read, omitting some bases at the ends of the read to ensure the alignment score meets the requirements. The penalty parameter –ma defaults to 2 in this mode.

[0081] The paired-end FASTQ data obtained after filtering the rRNA sequence through bowtie2 alignment, and the alignment results with the reference genome are as follows: Figure 2 As shown, after filtering, subsequent transcriptome analysis can be performed. The time taken for transcriptome analysis before and after rRNA removal is shown in Table 1.

[0082] Table 1. Time taken for mouse transcriptome analysis before and after rRNA removal

[0083] The time required for CleanData transcriptome alignment analysis of samples after rRNA removal 4.5 hours

[0084] Therefore, it can be seen that removing rRNA using the method of this embodiment shortened the time for mouse transcriptome alignment analysis by 10%, which is a significant effect.

[0085] In addition, because the transcriptome library is constructed using strand-specific library construction (fr-firstrand), transcriptome sequencing can retain the orientation information of transcripts during transcriptome sequencing, thus also determining whether the transcripts originate from the positive or negative strand of the genome.

[0086] Example 3: Reference genome of the source species; ribosome sequence filtering of Muscovy duck transcriptome sequencing data not derived from Ensembl.

[0087] A method for filtering rRNA sequences from Muscovy duck transcriptome sequencing data not derived from Ensembl, using the method in Example 2.

[0088] Specifically, transcriptome data of Muscovy ducks were obtained (using a strand-specific library construction method), with the reference genome for this species sourced from the NCBI database.

[0089] The clean data was aligned with the reference genome using the hisat2 alignment software. The unmapped.fastq file was converted to unmapped.fasta, and the FASTA file containing the top 10,000 reads by abundance was uploaded to NCBI for BLAST alignment. The results are as follows: Figure 3 As shown. After obtaining the alignment results, filter the entries that match the rRNA, select and download the FASTA files of all rRNAs, and establish an rRNA database.

[0090] The paired-end FASTQ data obtained after filtering the rRNA sequence through bowtie2 alignment, and the alignment results with the reference genome are as follows: Figure 4 As shown in Table 2, after filtering, subsequent transcriptome analysis can be performed. The time taken for transcriptome analysis before and after rRNA removal is shown in Table 2.

[0091] Table 2. Time taken for Muscovy duck transcriptome analysis before and after rRNA removal

[0092] The time required for CleanData transcriptome alignment analysis of samples after rRNA removal 5 hours

[0093] Therefore, it can be seen that after removing rRNA using the method of this embodiment, the time required for Muscovy duck transcriptome alignment analysis was shortened by about 17%, which is a very significant effect.

[0094] Example 4: Device for filtering rRNA sequences in transcriptome sequencing data

[0095] Based on the method of Example 2, this embodiment provides a device for filtering rRNA sequences in transcriptome sequencing data, such as... Figure 5,include:

[0096] The data input module is used to obtain the transcriptome sequencing data of this species. Here, the transcriptome sequencing data used to establish the rRNA sequence database is the same set of data.

[0097] A database storage module is used to store the rRNA sequence database established in Example 2;

[0098] The filtering module, connected to both the data input module and the database storage module, compares transcriptome sequencing data with an rRNA sequence database, filters out matching sequences, and outputs the results.

[0099] All documents mentioned in this invention are incorporated herein by reference as if each document were individually incorporated by reference. Furthermore, it should be understood that after reading the foregoing teachings of this invention, those skilled in the art can make various alterations or modifications to this invention, and these equivalent forms also fall within the scope defined by the appended claims.

Claims

1. A method for establishing a species rRNA sequence database, characterized in that, Includes the following steps: S1, Obtain the first transcriptome sequencing data and reference genome data of the species, wherein the reference genome data is not from the Ensembl database; S2, compares the transcriptome data with the reference genome data to obtain the unaligned sequences. S3, sort the uncompared sequences from highest to lowest abundance, select the top M sequences, and compare them with the NCBI database, where M = 8000~20000. The comparison with the NCBI database refers to uploading the uncompared sequences to NCBI for comparison. S4, obtain the aligned rRNA sequence and establish an rRNA sequence database. The alignment means that the alignment E value is less than a preset threshold, which is set to 1e-5.

2. The method for establishing a species rRNA sequence database according to claim 1, characterized in that, After obtaining the transcriptome data of the species, the process further includes the steps of preprocessing the transcriptome data: (1) filtering adapter reads; (2) filtering reads containing more than 5% N; and (3) filtering reads with a quality value of Q ≤ 10 accounting for more than 20% of the total reads.

3. An rRNA sequence database established using the method of claim 1 or 2 based on the first transcriptome sequencing data of a species.

4. A method for filtering rRNA sequences from second transcriptome sequencing data of a species, characterized in that, Includes the following steps: The second transcriptome sequencing data is compared with the rRNA sequence database described in claim 3, and sequences that can be matched are filtered out.

5. The method for filtering rRNA sequences from the second transcriptome sequencing data of a species according to claim 4, characterized in that, The second transcriptome sequencing data is the same as the first transcriptome sequencing data.

6. An apparatus for filtering rRNA sequences from second transcriptome sequencing data of a species, characterized in that, include: The data input module is used to obtain the second transcriptome sequencing data; A database storage module is used to store the rRNA sequence database as described in claim 3; The filtering module is connected to both the data input module and the database storage module, and is used to compare the second transcriptome sequencing data with the rRNA sequence database, filter out the sequences that can be matched, and output the results.

7. A computer device, characterized in that, include: Memory, used to store computer programs; A processor for implementing the steps of the method as described in any one of claims 4 to 5 when executing the computer program.

8. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a computer program that, when executed by a processor, implements the steps of the method as described in any one of claims 4 to 5.

Images