Pinus massoniana miRNA identification method based on sRNA sequencing data

Abstract

The invention discloses a masson pine miRNA recognition method based on sRNA sequencing data, and belongs to the technical field of plant bioinformatics. The method mainly comprises the following steps: acquiring original offline data of masson pine sRNA sequencing; performing quality control processing on reads in the original offline data of the sRNA sequencing; converting a file of the high-quality sRNA sequencing data from an FASTQ format into an FASTA format, and combining the repetitive sequences to obtain a combined sRNA sequence; finding out all non-miRNA sequences in the Rfam database, comparing the merged sRNA sequences with all the non-miRNA sequences in the Rfam database by using bowtie according to the standard that one base is mismatched, and reserving the sequences which cannot be compared to form sRNA data in an FASTA format; and continuously performing reference genome sequence comparison, secondary structure extraction, structure filtration, abundance filtration and the like on the sRNA data to finally obtain the miRNA sequence.

Description

Technical Field

[0001] The present invention belongs to the technical field of plant bioinformatics, and in particular relates to a method for identifying Pinus massoniana miRNA based on sRNA sequencing data. Background Art

[0002] MicroRNAs (miRNAs), a class of non-coding RNA molecules approximately 20-24 nt in length, play a central role in regulating biological processes such as plant growth and development, metabolic pathways, and stress responses. Precise identification of miRNA forms and expression patterns is a prerequisite for in-depth research into miRNA function. Early identification relied primarily on techniques such as clone sequencing, gene chips, and real-time fluorescence quantitative PCR. However, these methods generally suffer from low throughput and limited sensitivity, making them inadequate for large-scale, genome-wide miRNA studies.

[0003] In recent years, with the widespread adoption of high-throughput sequencing technology, miRNA screening and identification based on sRNA sequencing has become a mainstream method for studying miRNA function. sRNA sequencing primarily utilizes the unique structure of the 3' and 5' ends of sRNA molecules (an intact phosphate group at the 5' end and a hydroxyl group at the 3' end). First, total RNA from the sample or tissue is extracted and RNA fragments of 18 to 30 nt are isolated. T4 ligase is then used to attach adapters to both ends of the sRNA molecule, followed by reverse transcription to synthesize cDNA. After PCR amplification, the target DNA fragments are separated using PAGE gel electrophoresis. The gel is then excised and recovered to obtain a cDNA library. Finally, sequencing is performed using an Illumina high-throughput instrument. After analysis of the sequencing data, sRNA sequencing data for the sample or tissue is obtained.

[0004] There are two main analytical paradigms for miRNA identification based on sRNA sequencing data:

[0005] 1) "Homology-first" process

[0006] First, sequencing reads are rigorously aligned to reference databases such as miRBase (typically requiring zero mismatches) to identify known miRNAs. Subsequently, known miRNA reads are removed, and the remaining sequences are screened for novel miRNAs using features such as secondary structure prediction and free energy calculation. Representative tools include miRDeep-P2.

[0007] 2) De novo-first process

[0008] First, all reads are globally profiled to predict the characteristics, and sequences meeting miRNA criteria (such as hairpin structure, Dicer cleavage site, etc.) are marked as candidate miRNAs; then known / new miRNAs are distinguished by homology alignment. The mainstream tools include miRDeep2, ShortStack, etc.

[0009] The core contradictions in current plant miRNA research mainly include three aspects. The first aspect is technical compatibility barriers. Different software prediction results have poor consistency, which hinders data integration and cross-study comparison. The second aspect is insufficient process reliability. The existing two analysis paradigms both have methodological defects, which affects the identification accuracy. The third aspect is standard iteration lag. New annotation standards (such as plant-specific hairpin structure, Dicer processing characteristics, etc.) have not been fully integrated into the analysis process.

[0010] Especially noteworthy is that in 2018, Meyers et al. updated the '08 standard and proposed the '18 standard suitable for the big data era. Therefore, the miRNA data set based on the previous loose standard needs to be systematically re-evaluated. And developing a new analysis platform that can integrate the latest biological knowledge, overcome the defects of the existing process, and realize automatic standard verification has become an urgent technical need to break through the bottleneck of plant miRNA research. SUMMARY

[0011] In view of the problems existing in the prior art, the technical problem to be solved by the present application is to provide a miRNA identification method based on sRNA sequencing data.

[0012] To solve the problems in the prior art, the technical scheme of the present application is as follows:

[0013] A Pinus massoniana miRNA identification method based on sRNA sequencing data, comprising the following steps:

[0014] S1, obtaining the original sRNA sequencing data of Pinus massoniana;

[0015] S2, performing quality control processing on the reads in the original sRNA sequencing data to obtain high-quality sRNA sequencing data;

[0016] S3, converting the file of high-quality sRNA sequencing data from FASTQ format to FASTA format, and merging repeated sequences to obtain merged sRNA sequences;

[0017] S4, finding all non-miRNA sequences in the Rfam database, using bowtie to align the merged sRNA sequences with all non-miRNA sequences in the Rfam database with a 1-base mismatch allowed standard, and retaining the sequences that cannot be aligned to form FASTA format sRNA data;

[0018] S5, pre-filtering: pre-filtering the sRNA data in FASTA format, directly retaining sequences of 20-22 nt, but only retaining sequences of 23-24 nt length with RPM≥20 and abundance≥10, to obtain pre-filtered sRNA data;

[0019] S6, using bowtie to align the pre-filtered sRNA data with reference sequences to obtain a sam file, merging sequences with positions close to each other; removing sRNA sequences with alignment positions on the reference sequences exceeding 20; the reference sequences are genome sequences of plants or reference sequences of second-generation transcriptome assembly;

[0020] S7, using bowtie to extract potential miRNA precursors at the alignment positions on the reference sequences and predict local stable secondary structures thereof, with the length of the precursors being limited to 300 nt;

[0021] S8, performing structure filtering on the miRNA of the precursors limited to 300 nt;

[0022] S9, performing abundance filtering on the miRNA sequences obtained in step S8 to obtain miRNA data of Pinus massoniana.

[0023] In some embodiments, the quality control processing in step S1 includes: removing 3p adapters of reads, retaining only reads of 18-30 nt, removing reads not containing 3p adapters, removing reads containing unknown bases, and removing reads with too low quality values.

[0024] In some embodiments, the sam file in step S6 records information of alignment of each sRNA sequence with reference sequences, including positions, directions, and positive and negative strands of sequence alignment to the reference sequences.

[0025] In some embodiments, step S7 specifically includes: extending 150 nt sequences before and after each position obtained by bowtie alignment on the reference sequences, and extracting them as potential miRNA precursors using bedtools software, predicting local stable secondary structures of each potential miRNA precursor using RNALfold software, and limiting the length of the precursors to 300 nt.

[0026] In some embodiments, the structure filtering criteria in step S8 are: 1) retaining the miRNA sequence which can be predicted according to the miRNA and precursor sequence; 2) the total number of miRNA / miRNA* duplex mismatches cannot exceed 5, and the maximum continuous mismatch in the duplex region cannot exceed 3; 3) the length of miRNA is between 20-24 nt, and the length of miRNA* is between 18-24 nt; 4) the length difference between miRNA and miRNA* is less than or equal to 2 nt.

[0027] In some embodiments, the structure filtering criteria in step S8 for miRNA with length of 23-24 nt further include that the miRNA* sequence must be present in the sequencing data.

[0028] In some embodiments, the abundance filtering criteria in step S9 are: 1) the minimum distance between miRNA and miRNA* is 5 nt, and the maximum distance is 300 nt; 2) the maximum length of miRNA precursor is 300 nt; 3) the abundance is calculated by considering the isomer, and the isomer is selected according to the criterion of only one base difference with the miRNA or miRNA*; 4) the proportion of precursor chain reads in the total number of double-stranded reads is greater than 0.9; 5) the proportion of the total number of miRNA and miRNA* and single nucleotide isomer reads in the total number of precursor reads is greater than 0.7; 6) the duplex with higher expression abundance is temporarily defined as miRNA, and the duplex with lower expression abundance is temporarily defined as miRNA*, and the abundance high-low criterion is self-defined.

[0029] In some embodiments, the method further comprises S10: comparing the miRNA data of Pinus massoniana with the miRNA database to distinguish known miRNA and new miRNA.

[0030] The miRNA sequence of Pinus massoniana obtained by the method.

[0031] The miRNA sequence of Pinus massoniana obtained by the method.

[0032] Compared with the prior art, the method has the following beneficial effects:

[0033] The application develops the above miRNA identification method based on sRNA sequencing data. The software implemented in particular is written in Python, Shell and R, and some commonly used bioinformatics software is built in. After the sRNA sequencing data is preprocessed, the secondary structure meeting the miRNA production conditions is screened according to the characteristics of plant miRNA production, and the miRNA information is obtained. In order to enable the software to process more different types of data, 24 parameters are additionally set, and the person skilled in the art can flexibly adjust the threshold according to the characteristics of the research species. The miRNA identification method and software developed by the application have more than 80% consistency with the identification results of other software such as sRNAminer. In addition to being able to better meet the miRNA annotation standard, only a single command is needed to conveniently, accurately and quickly identify miRNA, and the running results are detailed, including many information of miRNA sequence and expression amount, and only the corresponding data needs to be extracted for subsequent analysis. BRIEF DESCRIPTION OF DRAWINGS

[0034] Figure 1 It is a flow chart of miRNA identification method and program;

[0035] Figure 2 It is a schematic diagram of the original data format;

[0036] Figure 3 It is a pinus massoniana miRNA identification and filtering standard diagram;

[0037] Figure 4 It is a diagram of miRNA precursor secondary structure, wherein the left diagram is a secondary structure meeting the standard, and the right diagram is a secondary structure not meeting the standard;

[0038] Figure 5 It is a diagram of maximum continuous mismatch of miRNA, wherein the left diagram has a mismatch of 2, meeting the standard, and the right diagram has a mismatch of 4, not meeting the standard;

[0039] Figure 6 It is a schematic diagram of the content of the structure_filtering.pri2mature.input file, note: this diagram only shows part of the file content, and the length of miRNA and miRNA* is adjusted;

[0040] Figure 7 It is a diagram of comparison between miRNA production site and siRNA production site, wherein A is a diagram of miRNA production site, and B is a diagram of siRNA production site;

[0041] Figure 8The 2-7 miRNA sequences and isomers (isomiRs); Note: Red is the miRNA sequence, green is the miRNA* sequence, black is their single nucleotide isomers, the number in the brackets of the sequence is the expression level, the gray sequence is not the single nucleotide isomer of miRNA and miRNA*, which is not considered when calculating the abundance;

[0042] Figure 9 The complete help document content diagram of python miRNApipeline.py;

[0043] Figure 10 The command execution diagram for the identification of Pinus massoniana miRNAs;

[0044] Figure 11 The number distribution diagram of known and new miRNAs of Pinus massoniana, wherein each point represents one miRNA. DETAILED DESCRIPTION

[0045] In order to make the purpose, technical scheme and advantages of the present application clearer, the present application will be further described below in combination with specific examples. If no detailed description is given in the following examples, the technical means used are all conventional means familiar to those skilled in the art. Or according to the reagent kit and product instruction. The materials, reagents and the like used in the following examples, if no special description, can be obtained from commercial channels.

[0046] Example 1

[0047] Test data

[0048] The program test data is 19 sRNA libraries of diploid Fragaria vesca (Table 1). The above sequencing raw data FTP address was obtained using SRAExplorer (https: / / sra-explorer.info / #) and downloaded using axel (https: / / github.com / axel-download-accelerator / axel). The reference genome is F. vesca V4.0.

[0049] Table 1 Test data

[0050] Run number Accession number Library type Developmental stage SRR1586419 GSE61798 sRNA-Seq 4DPA ovary wall SRR1586420 GSE61798 sRNA-Seq 10DPA ovary wall SRR1586421 GSE61798 sRNA-Seq 4DPA seeds SRR1586422 GSE61798 sRNA-Seq 10DPA seeds SRR1586423 GSE61798 sRNA-Seq 10DPA receptacles SRR1586424 GSE61798 sRNA-Seq Unopen flowers SRR1586425 GSE61798 sRNA-Seq Unexpanded leaves SRR1586426 GSE61798 sRNA-Seq 4W seedling SRR771535 GSE44930 sRNA-Seq Open flowers of 1D SRR771536 GSE44930 sRNA-Seq Open flowers of 1D SRR771537 GSE44930 sRNA-Seq Open flowers of 1D SRR7250917 SRP149613 sRNA-Seq Anthers of stage_6_S1 SRR7250918 SRP149613 sRNA-Seq Anthers of stage_7_S2 SRR7250919 SRP149613 sRNA-Seq Anthers of stage_8_S3 SRR7250920 SRP149613 sRNA-Seq Anthers of stage_9_S4 SRR2003391 PRJNA282518 sRNA-Seq Stolons of 8WP SRR2003393 PRJNA282518 sRNA-Seq Open flowers of 8WP SRR2003394 PRJNA282518 sRNA-Seq Fruits of 8WP SRR2003396 PRJNA282518 sRNA-Seq Leaves of 8WP

[0051] Note: Days post anthesis (DPA); Week-old plants (WP); Months plant (MP); Day-opening (D); Week old plants (W).

[0052] System software and hardware and environment requirements:

[0053] The operating system was Ubuntu 18.04, with 96GB of memory and 10TB of hard disk storage. The Python environment was v3.7.9, and software such as dnapi.py, fastx_collapser (https: / / github.com / agordon / fastx_toolkit), cutadapt, RNALfold, bowtie, samtools, bedtools, starseeker, multiqc, fastqc (http: / / www.bioinformatics.babraham.ac.uk / projects / fastqc / ), axel, and TBtools were installed. Furthermore, the Python modules click, os, subprocess, and sys were used.

[0054] miRNA identification methods and program design

[0055] In response to the problems existing in the existing miRNA identification process, this application independently constructs plant miRNA identification methods and procedures based on the latest miRNA annotation standards ( Figure 1 The main process is to obtain sRNA sequencing data, preprocess, and filter the structure and abundance to complete plant miRNA screening and identification, and obtain plant miRNA data. The specific principles and methods include the following:

[0056] 1. Preprocessing module

[0057] (1) Data quality control

[0058] The raw data of sRNA sequencing is usually stored in FASTQ format, which consists of four lines per sequence. Figure 2 As shown, the first line is the sequence ID and description information of the sequence. A00881 in the first sequence ID represents the sequencer number; 291 is the number of the sequencing reaction; HJVY7DRXX is the flowcell number; the next four numbers are the sequence position information, 1 represents the first lane in the flowcell, 2101 represents the 2101th tile in the first lane, and 28818:1016 represents the x:y coordinate information of the sequence in the tile. The second line is the sequence details. The third line must start with "+". The fourth line is the quality value, each character corresponds to the quality of the base at the corresponding position in the second line. The character can be converted into a base quality score according to certain rules, and the base quality score can reflect the error rate of the base. The correspondence between the quality value and the error rate can be obtained by the formula,

[0059] Q = -10 log 10 P

[0060] Where Q represents the quality value, i.e. the ASCII value corresponding to the fourth line character, and P represents the base error rate.

[0061] Second-generation sequencing requires adding adapters to both ends of the sequencing fragments to enable sequencing, resulting in the presence of adapter sequences in the original sequencing data, which can affect subsequent analysis. In addition, some reads in the original data may have poor quality values, which can also adversely affect subsequent analysis. Therefore, data quality control can remove the adapters (3p adapters) of reads, remove too long and too short reads (retain 18-30 nt reads), remove reads without adapters (3p adapters), remove reads containing N (N represents unknown bases), and remove reads with too low quality values.

[0062] (2) Merging of identical reads

[0063] The file obtained in the previous step is a FASTQ file, which needs to be converted to a FASTA file for subsequent analysis. Unlike the FASTQ format, the FASTA format has only two lines per sequence, and the first line starts with “>”, not “@”. The first line, similar to FASTQ, is a description of the sequence. After converting the FASTQ file to FASTA, the next step can be performed. At this time, the file contains several million sequences, and a significant portion of these sequences are repetitive. If the repetitive sequences are not merged, the original file will be large, and the temporary files generated during analysis will also be very large, consuming a lot of server resources and significantly increasing the running time. After merging identical reads, the data storage space is only 1 / 50 of the original, and the processing time is greatly reduced. Therefore, in this step, the fastx_collapser software is used to merge identical sequences, and the resulting file size is approximately 20M. It is worth noting that if the file size is significantly higher than this threshold, it may indicate a problem with the data.

[0064] (3) Removal of other RNA fragments

[0065] It has been found that rRNA and tRNA can also derive sRNA, which will interfere with the analysis and need to be filtered out. First, Rfam_withoutMIR.fa is obtained by merging the non-miRNA sequences in the Rfam database (https: / / www.jianshu.com / p / cb75b0bf134e); the sequences after merging are aligned with Rfam_withoutMIR.fa using bowtie with one base mismatch standard, and the sequences that cannot be aligned are retained. Removing these sequences does not mean that this data is meaningless, but only means that in the analysis of miRNA, this part of the sequence may affect the results. In recent years, more and more researchers have focused their attention and research on rRNA and tRNA-derived sRNA, and have developed corresponding databases. It is believed that with the continuous deepening of the understanding of sRNA, they will be the next research hotspot.

[0066] 2. Structure filtering

[0067] Pre-filtering

[0068] After filtering out other RNA fragments in the previous step, a FASTA file is obtained again. For 20-22 nt sequences, directly start the subsequent analysis, and for longer sequences (23-24 nt), further filtering by expression is required, and sequences with RPM (Reads Per Million) ≥ 20 are selected for subsequent analysis. The most abundant sRNA in plants is siRNA, which has a peak at 24 nt, and longer sRNA is likely not miRNA. In order to reduce false positives, these sequences need to be additionally filtered ( Figure 3 ). sRNA sequencing will produce a large number of low-abundance sequences, which may only have 1 repeat in an sRNA library, and such sequences account for a considerable part of the entire library. Some newly evolved younger miRNA may have low abundance, but this miRNA may not have function at this stage, and most importantly, low-abundance sequences often have low reliability, so sequences with abundance less than 10 need to be filtered out. The above steps will obtain a file named Structure_filtering.input.fa. This file can be used as a structure filtering input file. All pre-filtering operations can be completed by the Pre_filtering function.

[0069] Structure filtering

[0070] The pre-filtered result file is aligned with the reference sequence using bowtie with zero mismatch standard to get a sam file, which contains detailed information of sequence alignment with the reference sequence, including sequence alignment position, direction, positive and negative strand, etc. Different positions on the genome can produce the same miRNA, so the same sequence can be aligned to several positions of the reference sequence. If the alignment times are too many, it indicates that these data come from repetitive sequences. Filter out the sequences with more than 20 alignment positions and write the position information into structure_filtering.extrac_seq.input. The merged_position function merges sequences with close positions (overlap). Then the extrac_seq function extends 150 nt before and after each position after merging to get potential precursor position information. Finally, the bedtools software getfasta command is used to extract the sequence to get potential miRNA precursor Figure 3

[0071] After obtaining the miRNA precursor, Structure_filtering.input.fa is aligned with the potential precursor using bowtie, and the sam file part information is stored in a dictionary named sam_information, in the form of:

[0072] {ref_name__ref_strand:[reads1_name,reads1_strand,reads1_position,reads2_name,reads2_strand,reads2_position,……]},

[0073] Structure_filtering.input.fa is stored in a dictionary named reads_dict, in the form of:

[0074] {reads_name:reads_sequence},

[0075] The potential precursor information is stored in a dictionary named flank150, in the form of:

[0076] {ref_name:ref_sequence},

[0077] The RNAfold function calls the RNALfold software to predict the local stable secondary structure of the precursor, and limits the length of the precursor to 300 nt. The RNALfold result is arranged, and the result is stored in the RNALfold_information dictionary, in the form of:

[0078] ​{secondary_structure: [MFE, start_position, sequence]},

[0079] The above information is passed to the structure_filtering function for structure filtering.

[0080] The filtering criteria for miRNAs of 20-22 nt in length are: 1) the miRNA* sequence must be able to be predicted from the miRNA and precursor sequence (2) the total number of mismatches in the miRNA / miRNA* duplex must not exceed 5, and the maximum number of consecutive mismatches in the duplex region must not exceed 3 (3) the miRNA length must be between 20-24 nt, and the miRNA* length must be between 18-24 nt (4) the length difference between the miRNA and miRNA* must be controlled to within 2 nt (5) for miRNAs of 23-24 nt in length, an additional criterion is added, which requires that the miRNA* sequence must be present in the sequencing data. Figure 4 Figure 5

[0081] The structure_filtering function outputs a structure_filtering.pri2mature.input file. The structure of this file is complex, with each sequence consisting of 6 lines (1) the pri-miRNA sequence (2) the miRNA sequence (3) the pre-miRNA sequence (4) the miRNA* sequence (5) the pre-miRNA secondary structure information (6) a split line, which has no actual meaning. For each pre-miRNA, there can be multiple local stable secondary structures, and the secondary structure with the minimum free energy, which is the most stable, and which has the highest read expression abundance, is selected as the pre-miRNA. The above steps can all be accomplished by the pri2mature function. Figure 6 The pri2mature function generates an Abundance_filtering.input file, which has 7 columns per line (1) the miRNA ID number (2) the miRNA sequence (3) the miRNA* ID number (4) the miRNA* sequence (5) the pre-miRNA sequence (6) the precursor secondary structure information (7) the precursor minimum free energy. The entire structure filtering module is packaged in the Structure_filtering function.

[0082] 3. Abundance filtering module

[0083]

[0084] ​​​Abundance_filtering.input can be used as an abundance filtering input file. Plant sRNA is divided into miRNA and siRNA. The production location of miRNA is relatively fixed ( Figure 7 A), while the location of siRNA production is not fixed ( Figure 7 (B) miRNAs can be screened based on this characteristic.

[0085] The specific criteria for abundance filtering are: 1) the minimum distance between miRNA and miRNA* is 5 nt and the maximum distance is 300 nt; 2) the maximum length of miRNA precursor is 300 nt; 3) isomers are also considered in abundance calculation, and the criterion for isomer screening is that they differ from miRNA or miRNA* by only one base ( Figure 8 ); 4) The ratio of precursor strand reads to total double-strand reads must be greater than 0.9; 5) The ratio of miRNA and miRNA* reads and their single-nucleotide isomers to total precursor reads must be greater than 0.7; 6) The dimer with the highest expression abundance is tentatively designated as a miRNA, while the dimer with the lowest expression abundance is tentatively designated as a miRNA*. The abundance criteria are customizable. The above process is encapsulated in the ABundance_filtering function.

[0086] After the above steps, relatively reliable miRNAs have been obtained. If there are many sRNA libraries, you can screen at least two libraries for miRNAs supported by them for downstream analysis. If there are fewer libraries, you can skip this step directly.

[0087] 4. Documentation and program operation

[0088] (1) Documentation

[0089] After rigorous screening, relatively reliable miRNA sequences are obtained. All information, including miRNA ID and sequence, miRNA* ID and sequence, pre-miRNA sequence, and expression level, is stored in a file named after the library. Manual correction is performed based on the four-mismatch principle to distinguish between known miRNAs (using the miRBase and sRNAanno databases as references), and the resulting sequences are uniformly named. Known miRNAs are named pma-miR319-mature, where pma represents the species abbreviation, miR319 represents the miRNA family, and mature indicates mature form. New miRNAs are named pma-miRN001-mature, with miRN used to denote novel miRNAs.

[0090] The above steps complete the first step of the study of miRNA, i.e. accurate identification of miRNA. Subsequent analysis of miRNA can include renaming, calculation of expression, differential analysis, identification of miRNA target genes, etc.

[0091] (2) Program running

[0092] The program is very simple to use, and complete help documentation can be obtained by typing python miRNApipeline.py--help on the command line, with details as shown in Figure 9

[0093] The --srna parameter specifies sRNA sequencing data, the --rfam parameter specifies the Rfam database file, and the --genome parameter specifies the reference sequence. The above three parameters are mandatory parameters and have no default values. To facilitate data processing, 24 optional parameters are also set to allow more personalized processing of sRNA data. The parameter information is shown in Table 2.

[0094] Table 2. Software parameter list

[0095]

[0096]

[0097] In view of the many problems existing in current miRNA data processing, the above software for miRNA screening and identification is developed. The software is written using Python, Shell and R, and some commonly used bioinformatics software is built in. After preprocessing of sRNA sequencing data, secondary structures that meet the conditions for miRNA production are screened according to the characteristics of plant miRNA production, and miRNA information is obtained. In order to enable the software to process more different types of data, 24 additional parameters are set. Researchers can flexibly adjust the thresholds according to the characteristics of the research species. The miRNA identification method and software developed in the present application have more than 80% consistency with the identification results of other software such as sRNAminer. In addition to being able to better meet the miRNA annotation standards, miRNA can be identified conveniently, accurately and quickly by a single command line. The running results are detailed and contain a lot of information such as miRNA sequences and expression, and only the corresponding data needs to be extracted for subsequent analysis.

[0098] Example 2: Identification of Pinus massoniana miRNA

[0099] 1. Construction and sequencing of Pinus massoniana sRNA library

[0100] ​Experimental materials were derived from different ramets of elite Masson pine clones. Three 10-year-old grafted ramets served as three biological replicates. Samples included megasporophyll (female cones that were not fully opened, PF), microsporophyll (male cones that were not fully opened, SF), young cones (12 months old, BF), bark (phloem and immature xylem, TB), young stems (3 months old, TS), and current-year needles (NL). Field samples were snap-frozen in liquid nitrogen and transported on dry ice.

[0101] Total RNA was extracted from 18 samples of 6 tissues of Pinus massoniana using the RNAprep Pure Polysaccharide and Polyphenol Plant Total RNA Extraction Kit (Tiangen, Beijing, China). The method and steps are detailed in the kit instructions. RNA integrity and concentration were determined according to the method of Chen et al.; RNA was used. Multiplex Small RNA Library Prep Set for Seventeen sRNA libraries were constructed using a DNA sequencing kit (NEB, CA, USA) (of which library construction failed for one TB sample). The specific principles and procedures were referred to the method of Wang et al.; the prepared sRNA libraries were sequenced at 1 × 50 bp using the Illumina HiSeq 2500 platform (Novogene, Beijing, China).

[0102] 2. Use miRNApipeline.py prepared in Example 1 to identify Pinus massoniana miRNA

[0103] Since the reference genome of Pinus massoniana has not yet been published, the reference sequence assembled using the second-generation transcriptome may result in fewer miRNA identification results. Therefore, this embodiment appropriately relaxes the Pinus massoniana miRNA identification standards, and modifies the new standards as follows: 1) The maximum mismatch of miRNA / miRNA* dimer is adjusted to 6, and the maximum length of the protrusion structure continues to be 3nt; 2) The maximum length of the miRNA precursor is adjusted from 300nt to a maximum distance of 300nt between miRNA and miRNA*, and the minimum distance between miRNA and miRNA* remains unchanged at 5; 3) When performing abundance filtering, the ratio of the total abundance of miRNA and miRNA* with its 8 single nucleotide isomers to the total precursor reads is adjusted from 0.75 to 0.7; the remaining standards refer to the "18 new standards". Specific identification commands are as follows Figure 10 shown.

[0104] After manually sorting the analysis results, 447 miRNAs were obtained, including 110 known miRNAs and 337 novel miRNAs. These miRNAs all met the latest annotation standards well, with known miRNAs and novel miRNAs accounting for 24.61% and 75.39% respectively ( Figure 11 , Table 3).

[0105] Table 3 Number of known miRNAs and new miRNAs in 17 libraries

[0106] Sample Known miRNAs New miRNAs BF1 90 242 BF2 90 247 BF3 99 263 NL1 79 213 NL2 81 246 NL3 83 242 PF1 97 249 PF2 99 262 PF3 98 264 SF1 91 220 SF2 94 243 SF3 93 233 TB1 83 242 TB2 82 232 TS1 94 249 TS2 91 252 TS3 87 211 Total 110 337

[0107] The identification results show that the number of new miRNAs of Pinus massoniana is much larger than that of known miRNAs. There are 78 known miRNAs expressed in all tissues, accounting for 70.91%, and there are 1 known miRNA only in BF, PF and SF, respectively, and no known miRNA exists only in NL, TB and TS. There are 194 new miRNAs in all six tissues, accounting for 57.57%, and the number of new miRNAs only existing in BF, NL, SF, TB and TS is 2, 11, 5, 4 and 1, respectively. The proportion of known miRNAs expressed in all tissues is much larger than that of new miRNAs, which is likely because these new miRNAs are newly evolved miRNAs with obvious tissue specificity. In addition, the length distribution of miRNAs is between 20-22 nt, and the number of miRNAs with a length of 21 nt is the largest, which is consistent with the research results of other species. The relative expression of miRNA and miRNA* is not fixed, if the expression of miRNA* is greater than 30% of the sum of the expression of miRNA and miRNA*, then this miRNA* will be included in the subsequent analysis. After analysis, 69 miRNA* meet the conditions.

[0108] In summary, the program of the miRNA identification method based on sRNA sequencing data provided in the present application has more than 80% consistency with other software such as sRNAminer. In addition to meeting the miRNA annotation standard well, it only needs a single command to identify miRNA conveniently, accurately and quickly, and the running results are detailed, including many information of miRNA sequence and expression, and only the corresponding data needs to be extracted for subsequent analysis.

Claims

1. A method for identifying Pinus massoniana miRNA based on sRNA sequencing data, characterized in that: The steps include: S1. Obtain the raw data of sRNA sequencing of Pinus massoniana; S2. Perform quality control on the reads in the original data of sRNA sequencing to obtain high-quality sRNA sequencing data; S3. Convert the high-quality sRNA sequencing data file from FASTQ format to FASTA format, merge the duplicate sequences, and obtain the merged sRNA sequence; S4. Find all non-miRNA sequences in the Rfam database, use bowtie to align the merged sRNA sequences with all non-miRNA sequences in the Rfam database, allowing a maximum of one base mismatch, and retain the sequences that cannot be aligned to form sRNA data in FASTA format; S5. Pre-filtering: Pre-filter the FASTA format sRNA data, directly retaining sequences of 20-22 nt, but only retaining sequences of 23-24 nt with RPM ≥ 20 and abundance ≥ 10, to obtain pre-filtered sRNA data; S6. Use bowtie to align the pre-filtered sRNA data with the reference sequence to obtain a sam file, merge overlapping sequences, and remove sRNA sequences that align more than 20 positions on the reference sequence; the reference sequence is the genome sequence of Pinus massoniana or the reference sequence of the second-generation transcriptome assembly; S7. Use bowtie to extract potential miRNA precursors and predict their local stable secondary structures at the positions compared to the reference sequence, limiting the precursor length to 300 nt. S8, structural filtering of miRNAs limited to 300 nt precursors; S9. Perform abundance filtering on the miRNA sequence obtained in step S8 to obtain miRNA data of Pinus massoniana.

2. The method according to claim 1, characterized in that The quality control process described in step S1 includes: removing 3p adapters from reads, retaining only reads of 18 to 30 nt, removing reads without 3p adapters, removing reads containing unknown bases, and removing reads with too low quality values.

3. The method according to claim 1, characterized in that The sam file described in step S6 records the information of the alignment of each sRNA sequence with the plant genome sequence, and the information includes the position, direction, positive and negative strands of the sequence alignment to the reference sequence.

4. The method according to claim 1, wherein Step S7 specifically includes extending the sequence before and after 150 nt of each position obtained by bowtie alignment on the reference sequence, extracting them as potential miRNA precursors using bedtools software, and predicting the local stable secondary structure of each potential miRNA precursor using RNALfold software, and limiting its length to 300 nt.

5. The method according to claim 1, wherein The criteria for structure filtering in step S8 are: 1) retaining miRNA sequences that can predict miRNA* sequences based on miRNA and precursor sequences; 2) the total number of miRNA / miRNA* dimer mismatches cannot exceed 5, and the maximum continuous mismatches in the dimer region cannot exceed 3; 3) the miRNA length is between 20 and 24 nt, and the miRNA* length is between 18 and 24 nt; 4) the difference in length between miRNA and miRNA* is less than or equal to 2 nt.

6. The method according to claim 1, characterized in that The structure filtering criteria in step S8 for miRNAs with a length of 23 to 24 nt also include: the miRNA* sequence must appear in the sequencing data.

7. The method according to claim 1, characterized in that The abundance filtering criteria described in step S9 are as follows: 1) the minimum distance between miRNA and miRNA* is 5 nt, and the maximum distance is 300 nt; 2) the maximum length of miRNA precursor is 300 nt; 3) isoforms are included in the abundance calculation, and the criterion for isoform selection is that they differ from miRNA or miRNA* by only one base; 4) the ratio of the number of precursor chain reads to the total number of double-stranded reads is greater than 0.9; 5) the ratio of the total number of miRNA and miRNA* and their single nucleotide isoform reads to the total number of precursor reads is greater than 0.7; 6) the dimer with higher expression abundance is tentatively designated as miRNA, and the dimer with lower expression abundance is tentatively designated as miRNA*. The abundance standards are self-determined.

8. The method according to claim 1, characterized in that The method further comprises step S10: comparing the Pinus massoniana miRNA data with the miRNA database to distinguish known miRNAs from novel miRNAs.

9. The miRNA sequence of Pinus massoniana obtained by any one of the methods described in claims 1-8.

10. Use of the Pinus massoniana miRNA sequence according to claim 9 in the study of Pinus massoniana gene regulatory network.

Images