A piRNA analysis method

Abstract

The present invention discloses a piRNA analysis method, comprising steps (1)-(7). The present invention analyzes the original anti-MIWI RIP-seq dataset to obtain comprehensive and ideal analysis results without the need for other processing operations on the piRNA sample through experimental means. This provides a powerful bioinformatics analysis tool for conducting piRNA tailing modification research in the reproductive and RNA fields.

Description

Technical Field

[0001] The present invention relates to the field of bioinformatics, and more particularly to a method for analyzing piRNAs (piwi-interacting RNAs), and more particularly to a bioinformatics analysis method for detecting 3'-end-tailed piRNAs (piwi-interacting RNAs) based on alignment of ribonucleic acid (RNA) sequencing reads. Background Art

[0002] Argonaute proteins, first discovered in Arabidopsis thaliana, are a class of proteins with a molecular mass of approximately 100 kDa. The Argonaute protein family is evolutionarily conserved and is divided into two subfamilies, AGO and PIWI. They bind to and mediate small non-coding RNAs (SNRs). In 2006, a new type of SNR was discovered, named PIWI-interacting RNA (piRNA) because of their specific binding to PIWI proteins. piRNAs are approximately 26 to 33 nt in length, with a U base preference at their 5' end and a 2'-O-methylation modification at their 3' end. They are specifically expressed in the animal germline, interacting with PIWI proteins to suppress transposons and regulate coding genes in animal germline cells, promoting the development and differentiation of germline cells and being essential for gametogenesis. In mice, PIWI proteins include three members: MIWI (PIWIL1), MILI (PIWIL2), and MIWI2 (PIWIL4), which are expressed in a strict temporal order during male germ cell development.

[0003] RNA tailing, the addition of non-templated nucleotides to the 3' end of RNA, is one of the most common RNA modifications. A previous study in the nematode Caenorhabditis elegans reported that the absence of 3'-terminal cleavage and 2'-O-methylation of pre-piRNAs induces 3'-terminal non-templated RNA tailing, which in turn triggers pre-piRNA degradation and severely impairs piRNA production in the nematode. However, piRNA tailing has not been studied in other organisms, such as mice. Furthermore, descriptions of tailing in nematodes are limited to experimental studies, leaving the field largely unaware of the specific tailing mechanism and base information. Therefore, a more detailed understanding of how piRNA tailing is performed in other organisms, the sequence of tailing, and how piRNA tailing influences piRNA function is warranted. To achieve this, bioinformatics analysis is required to characterize the base information of piRNAs in detail. However, traditional bioinformatics analysis can only be used to analyze untailed piRNAs and is unable to investigate the tailing of piRNAs. A systematic bioinformatics analysis approach is urgently needed to analyze piRNA tailing across diverse species and PIWI families. Summary of the Invention

[0004] Using traditional bioinformatics analysis methods, a small RNA sequencing dataset of MIWI-bound piRNAs from wild-type mouse testes was analyzed. The results showed that approximately 20% of MIWI-bound piRNAs failed to align successfully with the mouse genome; however, if one or more nucleotides were removed from their 3' end, these piRNAs would fully match the genome. This suggests that non-templated tailing also occurs in mouse testis piRNAs. Currently, there are no specific research plans for this phenomenon, nor are there relevant bioinformatics analysis processes to analyze it. To study the sequence information and specific mechanisms of piRNA tailing, traditional bioinformatics analysis needs to be updated so that the bioinformatics code can analyze the tailing of piRNAs bound to different species and different PIWI families. The present invention provides a piRNA analysis method. Compared with the traditional piRNA bioinformatics analysis process, the present invention has the following advantages: (1) The present invention establishes a new piRNA analysis method, which can comprehensively analyze untailed piRNAs and piRNAs with 3' tails, and obtain information such as the tailing abundance, form and ratio of tailed piRNAs, which is of great significance for studying the phenomenon of piRNA tailing and provides a powerful tool for piRNA analysis in the reproductive field. (2) The present invention can quickly identify the piRNA 3' tailing phenomenon (i.e., the number of tailed bases) through an iterative method, thereby achieving efficient sequence analysis. (3) piRNAs mainly have single-base tails, which are not easy to distinguish using wet experiments. The analysis method of the present invention completes the analysis completely through code, analyzes the tailing phenomenon of the entire piRNA library with high throughput, and clarifies the tailing base type, avoiding other wet experiment processing and re-sequencing analysis of the tailed piRNA samples. Therefore, the present invention can obtain comprehensive and ideal analysis results by analyzing the original anti-MIWI RIP-seq data set, without the need for other processing operations on the piRNA samples through experimental means. This provides a powerful bioinformatics analysis tool for piRNA tailing modification research in the reproductive and RNA fields.

[0005] Specifically, the first aspect of the present invention provides a piRNA analysis method, comprising the following steps:

[0006] (1) Align small RNA sequencing reads with the reference genome;

[0007] (2) The reads successfully aligned in (1) are further aligned with the pre-miRNA reference sequence; for the sequencing reads that are successfully aligned to miRNA, the normalization parameters are calculated based on the number of aligned reads of miRNA in different samples;

[0008] (3) After removing contamination from other small RNAs, reads that were not successfully aligned to miRNAs in (2) were aligned to the piRNA reference sequence;

[0009] (4) Reads that can be fully matched in (3) are defined as untailed piRNAs and normalized using the normalization parameters in step (2);

[0010] (5) Reads that were not successfully aligned in step (1) are defined as tailed piRNAs, and the 1nt base at the 3' end is removed, and then steps (1)-(3) are repeated to obtain the 1nt tailed piRNAs.

[0011] (6) Remove the 1 nt base at the 3' end of the reads that were not successfully aligned with the reference genome in (5), and then repeat steps (1)-(3) to obtain the piRNA with a 2nt tail.

[0012] (7) Iterate according to the methods in (5) and (6) until no new tailed piRNAs can be obtained; then merge all tailed piRNAs and normalize them using the normalization parameters of step (2).

[0013] The reads described in the present invention have a general definition in the field of bioinformatics, and can be defined as the read length of the downloaded sequence.

[0014] In a preferred embodiment, in step (2), a script (step2_untailed_reads.sh) is used to calculate the normalization parameter by comparing the number of reads of miRNAs in different samples.

[0015] The normalization parameters described in the present invention are obtained by calculation methods known in the art and have general definitions in the field of bioinformatics.

[0016] In a specific embodiment, the standardized parameters of the present invention are calculated by the following steps:

[0017] (i) Select a sample (e.g., wild-type mouse sample 1) as a standard sample;

[0018] (ii) calculating the total number of bases in the pre-miRNA reference sequence;

[0019] (iii) Calculate the total number of bases in the pre-miRNA reads aligned to different samples; then divide it by the total number of bases in the pre-miRNA reference sequence to obtain the sequencing depth;

[0020] (iv) The sequencing depth of different samples was divided by the standard sample to obtain the normalization parameters of each sample.

[0021] The standardization described in this invention also has a general definition in the field of bioinformatics. In a preferred embodiment, the standardization described in this invention refers to calculating and eliminating the impact of sequencing depth on the number of piRNAs in different samples, especially the impact of piRNA tailing for specific purposes (such as drug treatment, gene mutation, and environmental stress). Since different samples are treated differently, it may affect the expression of piRNAs in the organism and thus affect the sequencing depth. The normalization method of the present invention is based on the number of sequenced miRNA bases and the total number of pre-miRNA bases, and calculates and eliminates the impact of sequencing depth on the data of different samples.

[0022] In a preferred embodiment, the reference genome is the genome of the participant to be tested, such as a mouse reference genome.

[0023] In a specific embodiment, the reference genome (GRCm 38) is obtained from the Ensembl database.

[0024] In a preferred embodiment, steps (1) to (7) are performed using bioinformatics analysis and comparison software commonly used in the art, such as Bowtie, Bowtie2, BWA, STAR, and Tophat2.

[0025] In a preferred embodiment, steps (1) to (7) are implemented using Bowtie software.

[0026] In a preferred embodiment, in step (3), the other small RNA includes one or more of tRNA, rRNA or snRNA.

[0027] In a preferred embodiment, step (3) includes sequentially aligning reads that have not been successfully aligned to miRNA to tRNA, rRNA, and snRNA reference sequences to obtain reads that are not contaminated by other small RNAs.

[0028] In a preferred embodiment, step (0) is further included before step (1): performing data quality control on the sequencing reads of the small RNA.

[0029] In a preferred embodiment, the data quality control includes one or more of the following: removing 3' end sequence adapters of sequencing reads; removing PCR duplicates; removing random bases introduced before and after the sequence during library construction; setting a threshold to remove low-quality sequencing data; and filtering redundant sequence reads <18nt.

[0030] In a specific embodiment, the quality control step is as follows: remove adapters from the 3' end of reads using FASTX Toolkit (v0.0.14) subcommand fastx_clipper, followed by removing PCR duplicates (Perl script remove_duplicated_single_reads.pl), remove 4bp of random bases introduced at the beginning and end of the sequence during library construction using FASTX Toolkit subcommand fastx_trimmer, and remove low quality sequencing base data using FASTX Toolkit subcommand fastq_quality_filter with threshold (-q 20). Filter out redundant sequence reads of <18nt, then convert FASTQ data (*.fatsq) to FASTA data (*.fasta) using FASTX Toolkit subcommand fastx_collapser.

[0031] In a preferred embodiment, step (0) is further preceded by step (-1): library construction and sequencing of small RNAs to obtain sequencing reads of small RNAs.

[0032] In a preferred embodiment, step (-1) comprises: sequencing small RNAs of wild type mouse samples using Illumina HiSeq x Ten platform to screen for RNAs of 10-50nt fragment size.

[0033] In a specific embodiment, the present example is performed by GENEWIZ to sequence small RNAs of wild type mouse samples of MIWI-RIP (MIWI-RNA immunoprecipitation) using Illumina HiSeq x Ten, 2x150 paired-end sequencing to screen for RNAs of 10-50nt fragment size, with biological replicates twice.

[0034] In a preferred embodiment, the analysis method further comprises step (8): visual analysis of the un-tailed piRNAs and the tailed piRNAs.

[0035] In a preferred embodiment, the visual analysis comprises one or more of: length distribution analysis, Motif analysis, 5' end U-bases preference analysis, piRNA source analysis, and ratio analysis of tailed piRNAs to un-tailed piRNAs.

[0036] In a preferred embodiment, the method further comprises a final overall piRNA analysis, which is performed by summarizing untailed piRNAs and tailed piRNAs for: length distribution analysis, motif analysis, 5' end U base preference analysis, or piRNA source analysis.

[0037] In a specific embodiment, step (1) comprises: aligning the sequencing reads with a reference genome (Ensembl, GRCm 38) using Bowtie (v1.0.0) software (parameters: -p 20 -fv 0 -k 100 -a --best --strata).

[0038] In a specific embodiment, in step (2), the reads successfully aligned in (1) are aligned to the pre-miRNA reference sequence using Bowtie software (parameters: -p 20 -fv 0 -k 100 -a --best --strata --norc).

[0039] In a specific embodiment, the pre-miRNA reference sequence is found at ftp: / / mirbase.org / pub / mirbase / 22.1 / hairpin.fa.gz.

[0040] In a specific embodiment, in step (2), a script (step2_untailed_reads.sh) is used to calculate the normalization parameters by comparing the number of miRNA reads in different samples for subsequent piRNA normalization between different samples.

[0041] In a specific embodiment, in step (3), Bowtie software (parameters: -p 20 -fv 0 -k 100 -a --best --strata --norc) is used to align reads that fail to align to miRNA to a tRNA reference sequence to remove tRNA contamination.

[0042] In a specific embodiment, in step (3), Bowtie software (parameters: -p 20 -fv 0 -k 100 -a --best --strata --norc) is used to align reads that fail to tRNA to rRNA reference sequences to remove rRNA contamination.

[0043] In a specific embodiment, in step (3), Bowtie software (parameters: -p 20 -fv 0 -k 100 -a --best --strata --norc) is used to align reads that have not been successfully aligned to rRNA to the snRNA reference sequence to remove snRNA contamination.

[0044] In a specific embodiment, the tRNA, rRNA, and snRNA reference sequences are obtained from NCBI.

[0045] In a specific embodiment, in step (3), the reads that are successfully aligned to the snRNA are aligned to the reference sequence of the piRNA using Bowtie software (parameters: -p 20 -fv 0 -k 100 -a --best --strata --norc).

[0046] In a preferred embodiment, the reference sequence of piRNA is obtained from piRBase, v2.0, mouse, http: / / bigdata.ibp.ac.cn / piRBase / .

[0047] In a preferred embodiment, the RNA that successfully aligns with the tRNA, rRNA, and snRNA reference sequences in step (3) is a contaminating read and needs to be removed.

[0048] In a preferred embodiment, in step (8), the visualization analysis of untailed piRNA is performed by a shell script step2_untailed_reads.sh.

[0049] In a preferred embodiment, in step (8), the length distribution analysis of the untailed piRNA is performed using the R script viz_untailed.R; and / or, the motif analysis of the untailed piRNA is performed using the R script motif_code.R to draw a motif diagram; and / or, the 5'-end U base preference analysis of the untailed piRNA is performed using the R script viz_untailed.R to draw a pie chart.

[0050] In a preferred embodiment, in step (8), the piRNA source analysis of untailed piRNAs includes genomic annotation of piRNAs by converting GTF (Ensembl, GRCm 38) and RepeatMasker (ISB, v.4.0.5) annotation files to bed format, and converting the reference genome alignment results (*.genome.bam) to bed format using BEDTools (v2.25.0).

[0051] In a preferred embodiment, the piRNA source analysis includes using the BEDTools (v2.25.0) software subcommand intersect to obtain the source of the piRNA.

[0052] In a preferred embodiment, the piRNA source analysis includes visualization using the R script viz_untailed.R, for example, including visualization of intergenic regions (Intergenic), protein coding regions (Protein_coding), pseudogenes (Pseudogene), other GTF annotated regions (Others), LINE transposons (LINE), SINE transposons (SINE), LTR transposons (LTR), DNA transposons (DNA_transposon), other repetitive sequences (Other_repeat) and / or non-annotation regions (No_annotation).

[0053] In a preferred embodiment, in step (8), the visualization analysis of tailed piRNA is performed by a shell script step3_tailed_reads.sh).

[0054] In a preferred embodiment, in step (8), the code for visual analysis of tailed piRNAs refers to the code for visual analysis of untailed piRNAs.

[0055] In a preferred embodiment, in step (8), the ratio of tailed piRNA to untailed piRNA is analyzed by using the R script viz_untailed.R to calculate the ratio of each part and draw a pie chart.

[0056] In a preferred embodiment, in step (8), the tailing analysis includes using the R script viz_untailed.R to calculate the proportion of each part (1nt, 2nt, 3nt, 4nt and ≥5nt non-template tailing) and draw a pie chart.

[0057] In a preferred embodiment, step (8) includes global piRNA analysis: Shell script work.sh).

[0058] In a preferred embodiment, in step (8), the length distribution analysis of the total piRNAs is subsequently visualized using the R script viz_total.R.

[0059] In a preferred embodiment, in step (8), the motif analysis of the overall piRNA uses the R script motif_code.R to draw a motif map.

[0060] In a preferred embodiment, in step (8), the 5'-end U base preference analysis of the overall piRNA is performed using R

[0061] Script viz_total.R for visualization.

[0062] In a preferred embodiment, in step (8), the piRNA source analysis of the total piRNA uses the R script viz_untailed.R to visualize the source of the piRNA.

[0063] The second aspect of the present invention provides a universal analysis method for piRNAs of different species or different PIWI families, comprising replacing the RNA sample with a test sample and performing the analysis method as described in the first aspect of the present invention.

[0064] A third aspect of the present invention provides a computing device, comprising:

[0065] a memory for storing program instructions;

[0066] The processor is configured to call the program instructions stored in the memory and execute the analysis method as described in the first aspect of the present invention according to the obtained program instructions.

[0067] A fourth aspect of the present invention provides a computer-readable storage medium comprising computer-readable instructions. When a computer reads and executes the computer-readable instructions, the analysis method described in the first aspect of the present invention is implemented.

[0068] The fifth aspect of the present invention provides a computer program product, comprising a computer program executable by a computer device, wherein when the program is run on the computer device, the computer device executes the steps of the analysis method of the first aspect of the present invention.

[0069] On the basis of conforming to the common sense in this field, the above-mentioned preferred conditions can be arbitrarily combined to obtain the preferred embodiments of the present invention.

[0070] The reagents and raw materials used in the present invention are commercially available.

[0071] Compared with the traditional piRNA bioinformatics analysis process, the present invention has the following advantages:

[0072] (1) The present invention establishes a new piRNA analysis method, which can comprehensively analyze untailed piRNAs and piRNAs with 3'-end tailing, and obtain information such as the tailing abundance, form and ratio of tailed piRNAs. This method is of great significance for studying the phenomenon of piRNA tailing and provides a powerful tool for piRNA analysis in the field of reproduction.

[0073] (2) The present invention can rapidly identify the piRNA 3' tailing phenomenon (i.e., the number of tailing bases) through an iterative method, thereby achieving efficient sequence analysis.

[0074] (3) piRNAs primarily have single-base tails, which are difficult to distinguish using wet-label assays. The analysis method of the present invention uses code to perform analysis, enabling high-throughput analysis of tailing across the entire piRNA library and identifying the base types involved. This avoids the need for additional wet-label assay processing and subsequent resequencing of tailed piRNA samples. Therefore, the present invention can obtain comprehensive and ideal analysis results by analyzing the original anti-MIWI RIP-seq dataset without requiring additional experimental manipulation of the piRNA samples. BRIEF DESCRIPTION OF THE DRAWINGS

[0075] Figure 1 This is the technical roadmap of the present invention.

[0076] Figure 2 Length distribution analysis of untailed piRNAs.

[0077] Figure 3 Motif analysis of untailed piRNA.

[0078] Figure 4 This is an analysis of the 5' end base preference of untailed piRNAs.

[0079] Figure 5 Analysis of the genomic origin of untailed piRNAs.

[0080] Figure 6 This is the length distribution analysis of the tailed piRNA (the tailed bases of the tailed piRNA were not removed).

[0081] Figure 7 This is the Motif analysis of the tailed piRNA (the tailed bases of the tailed piRNA were not removed).

[0082] Figure 8 This is an analysis of the 5' end base preference of tailed piRNA (removing the tailing base of the tailed piRNA).

[0083] Figure 9 is the ratio of untailed piRNA to total piRNA.

[0084] Figure 10 It is the proportion of piRNAs that completely align to the piRNA reference sequence after removing one or more small RNA 3' ends.

[0085] Figure 11 This is the length distribution analysis of the total piRNA (the tailing bases of the tailed piRNA were not removed).

[0086] Figure 12 This is the Motif analysis of the overall piRNA (the tailing bases of the tailed piRNA were not removed).

[0087] Figure 13 This is an analysis of the 5' end base preference of the overall piRNA (excluding the tailing bases of the tailed piRNA).

[0088] Figure 14 This is the source analysis of the total piRNA (excluding the tailing bases of the tailed piRNA). DETAILED DESCRIPTION

[0089] Overview

[0090] Overall, the present invention provides a bioinformatics analysis method for studying piRNA 3' end tailing modification, which provides a powerful bioinformatics analysis tool for conducting piRNA tailing modification research in the reproductive and RNA fields.

[0091] Specifically, this application proposes the following technical solutions for piRNA analysis (see Figure 1 ):

[0092] (1) RNA samples were sequenced using the Illumina HiSeq×Ten platform, with two biological replicates.

[0093] (2) Data quality control (QC): FASTX Toolkit was used to remove sequence adapters from the 3' end of the reads, remove PCR duplicates, and remove 4 bp of random bases introduced before and after the sequence during library construction. A threshold (phred quality of 20) was then set to remove low-quality sequencing data. Redundant sequence reads <18 nt were filtered, and FASTQ data (*.fatsq) were converted to FASTA data (*.fasta).

[0094] (3) Alignment with the reference genome: Sequencing reads were aligned with the mouse (Mus_musculus) reference genome using Bowtie software. Reads that failed to be successfully aligned to the genome were used for further analysis of 3' end-tailed piRNAs.

[0095] (4) Data normalization: The reads that successfully aligned to the reference genome were aligned to the pre-miRNA reference sequence to obtain sequencing reads that successfully aligned to the miRNA, and normalization parameters were calculated for subsequent piRNA normalization between different samples.

[0096] (5) Removal of tRNA contamination: reads that failed to successfully map to miRNA were aligned to the tRNA reference sequence to obtain reads that were not contaminated by tRNA.

[0097] (6) Removal of rRNA contamination: reads that failed to successfully align to tRNA were aligned to the rRNA reference sequence to obtain reads that were not contaminated by rRNA.

[0098] (7) Removal of snRNA contamination: reads that failed to successfully map to rRNA were aligned to the snRNA reference sequence to obtain reads that were not contaminated by snRNA.

[0099] (8) Obtaining untailed piRNAs: Reads that fail to successfully align to snRNAs will be used to align to the piRNA reference sequence. If these reads completely match the reference sequence (no base mismatch), they will be defined as untailed piRNAs and normalized using the normalization factor obtained in step (4) for further piRNA visualization analysis.

[0100] Visual analysis of untailed piRNAs:

[0101] Length distribution analysis;

[0102] Motif analysis;

[0103] Analysis of 5'-end U base preference.

[0104] PiRNA source analysis: For piRNA genomic annotation, we converted GTF and RepeatMasker annotation files to bed format and converted the reference genome alignment results (*.genome.bam) to bed format using BEDTools. We then used the BEDTools subcommand intersect to obtain the piRNA source.

[0105] (9) Obtaining tailed piRNAs: Remove the 1 nt base at the 3' end of the reads that failed to successfully map to the genome, and then repeat steps (3-8) with the reads that removed the 1 nt base to obtain 1-nt-tailed piRNAs. Then, remove the 1 nt base at the 3' end of the reads that still failed to map to the reference genome, and then repeat steps (3-8) to obtain 2-nt-tailed piRNAs. Repeat the above method until no new tailed piRNAs can be obtained. Then, merge the tailed piRNA reads into a FASTA file (*.fasta) and normalize them using the normalization factor obtained in step (4) for further visualization analysis of tailed piRNAs.

[0106] Visual analysis of tailed piRNAs:

[0107] Length distribution analysis;

[0108] Motif analysis;

[0109] Analysis of 5'-end U base preference;

[0110] Analysis of the proportion of tailed piRNAs and untailed piRNAs.

[0111] (10) Overall piRNA analysis: Summarize untailed piRNAs and tailed piRNAs for visualization analysis of overall piRNAs:

[0112] Length distribution analysis;

[0113] Motif analysis;

[0114] Analysis of 5'-end U base preference;

[0115] piRNA source analysis: same as step (8).

[0116] The present invention is further illustrated by way of examples below, but the present invention is not limited to the scope of the examples. Experimental methods in the following examples where specific conditions are not specified were performed according to conventional methods and conditions, or selected according to the product specifications.

[0117] Example 1

[0118] The present application is described below with reference to a specific example (small RNA sequencing combined with MIWI):

[0119] Using traditional bioinformatics analysis methods to analyze the MIWI-bound small RNA sequencing dataset in wild-type mouse testis (PRJNA1103657), we found that approximately 20% of MIWI-bound piRNAs failed to successfully align with the mouse genome; however, if one or more nucleotides (1-10nt) were removed at their 3' end, these piRNAs would completely match the genome.

[0120] Therefore, the present invention proposes a new piRNA analysis method, which can analyze piRNA more comprehensively and specifically. The technical route is as follows (see Figure 1 ):

[0121] (1) All small RNA sequencing platforms are applicable. In this case, GENEWIZ Corporation prepared the small RNA library using Illumina HiSeq×Ten, 2×150 paired-end sequencing, and performed small RNA sequencing on wild-type mouse samples from MIWI-RIP (MIWI-RNA immunoprecipitation) to screen for RNA fragments of 10-50 nt in size. The biological results were repeated twice.

[0122] (2) Data quality control (QC, step1_qc_fastx.sh): The FASTX Toolkit (v0.0.14) subcommand fastx_clipper was used to remove adapters from the 3' end of the reads. PCR duplicates were then removed (Perl script remove_duplicated_single_reads.pl). The FASTX Toolkit subcommand fastx_trimmer was used to remove 4 bp of random bases introduced before and after the sequence during library construction. Low-quality sequencing base data were removed using the FASTX Toolkit subcommand fastq_quality_filter with a threshold of -q 20. Redundant sequence reads <18 nt were then filtered, and the FASTX Toolkit subcommand fastx_collapser was used to convert the FASTQ data (*.fatsq) to FASTA data (*.fasta).

[0123] (3) Alignment to the reference genome: Sequencing reads were aligned to the reference genome (Ensembl, GRCm 38) using Bowtie (v1.0.0) software (parameters: -p 20 -fv 0 -k 100 -a --best --strata). Reads that failed to be successfully aligned to the genome were used for further analysis of tailed piRNAs.

[0124] (4) Data normalization: Reads successfully aligned to the reference genome were aligned to the pre-miRNA reference sequence (ftp: / / mirbase.org / pub / mirbase / 22.1 / hairpin.fa.gz) using Bowtie software (parameters: -p 20 -fv 0 -k 100 -a--best--strata--norc) to obtain sequencing reads that successfully aligned to miRNAs. The script (step2_untailed_reads.sh) was then used to calculate normalization parameters based on the number of aligned reads for miRNAs in different samples for subsequent piRNA normalization between different samples.

[0125] (5) Removal of tRNA contamination: Bowtie software (parameters: -p 20 -fv 0 -k 100 -a --best --strata --norc) was used to align the reads that failed to successfully map to miRNA with the tRNA reference sequence (mouse tRNA obtained from NCBI) to obtain reads that were not contaminated by tRNA.

[0126] (6) Removal of rRNA contamination: Reads that failed to successfully align to tRNA were aligned to the rRNA reference sequence (NCBI) using Bowtie software (parameters: -p 20 -fv 0 -k 100 -a --best --strata --norc) to obtain reads that were not contaminated by rRNA.

[0127] (7) Removal of snRNA contamination: Bowtie software (parameters: -p 20 -fv 0 -k 100 -a --best --strata --norc) was used to align the reads that failed to successfully map to rRNA with the snRNA reference sequence (NCBI) to obtain reads that were not contaminated by snRNA.

[0128] (8) Obtaining untailed piRNAs: Reads that failed to successfully align to snRNAs will be aligned to the piRNA reference sequence (piRBase, v2.0, mouse, http: / / bigdata.ibp.ac.cn / piRBase / ) using Bowtie software (parameters: -p20 -fv 0 -k 100 -a --best --strata --norc). If these reads completely match the reference sequence (no base mismatch), they will be defined as untailed piRNAs and normalized using the normalization parameters obtained in step (4) for further piRNA visualization analysis. Visualization analysis of untailed piRNAs (Shell script step2_untailed_reads.sh):

[0129] Length distribution analysis (see Figure 2 ): Combined with sequencing reads, the overall length distribution of untailed piRNAs present in MIWI-IP small RNAs was analyzed and subsequently visualized using the R script viz_untailed.R.

[0130] Motif analysis (see Figure 3 ): Combined with the sequence information of untailed piRNAs present in MIWI-IP small RNAs, the R script motif_code.R was used to draw a Motif diagram to display the sequence base distribution of piRNAs;

[0131] 5' end U preference analysis (see Figure 4 ): Based on the first base sequence of the untailed piRNA present in the MIWI-IP small RNA, the R script viz_untailed.R was used to draw a pie chart showing the U preference of the first base of the piRNA;

[0132] · piRNA source analysis (participate in Figure 5 ): For piRNA genomic annotation, GTF (Ensembl, GRCm38) and RepeatMasker (ISB, v.4.0.5) annotation files were converted to bed format, and the reference genome alignment results (*.genome.bam) were converted to bed format using BEDTools (v2.25.0). The BEDTools (v2.25.0) software subcommand intersect was used to obtain the source of piRNAs. The R script viz_untailed.R was then used to visualize the genomic sources of untailed piRNAs present in MIWI-IP small RNAs, including intergenic regions (Intergenic), protein coding regions (Protein_coding), pseudogenes (Pseudogene), other GTF annotated regions (Others), LINE transposons (LINE), SINE transposons (SINE), LTR transposons (LTR), DNA transposons (DNA_transposon), other repetitive sequences (Other_repeat), and non-annotated regions (No_annotation).

[0133] (9) Obtain tailed piRNAs: Remove the 1nt base at the 3' end of the reads that failed to successfully map to the genome, and then repeat steps (3-8) with the reads that removed the 1nt base to obtain 1nt-tailed piRNAs; then remove the 1nt base at the 3' end of the reads that still failed to map to the reference genome, and then repeat steps (3-8) to obtain 2nt-tailed piRNAs. Iterate in the above manner until no new tailed piRNAs can be obtained, then merge the tailed piRNA reads into a FASTA file (*.fasta), and normalize them using the normalization factor obtained in step (4) for further visualization analysis of tailed piRNAs. Visualization analysis of tailed piRNAs (Shell script step3_tailed_reads.sh):

[0134] Length distribution analysis (see Figure 6 ) : Combined with sequencing reads, the overall length distribution of the tailed piRNAs (without removing the tailing bases) present in the MIWI-IP small RNA was analyzed to verify whether there is a length preference in the tailing;

[0135] Motif analysis (see Figure 7 ): Combined with the sequence information of the tailed piRNA (without removing the tailed base) present in the small RNA of MIWI-IP, the base distribution of the entire piRNA sequence is displayed;

[0136] • 5' end U bias analysis (see Figure 8 ): The first base sequence of the sequence of the tailed piRNA (remove the tailing base) existing in the small RNA combined with MIWI-IP shows the U bias of the first base of piRNA;

[0137] • Tailing analysis (see Figure 9 ): According to the reads of the different degrees of non-template tailing piRNA existing in the small RNA combined with MIWI-IP, the proportion of each part (1 nt, 2 nt, 3 nt, 4 nt and ≥5 nt non-template tailing) is calculated using the R script viz_untailed.R, and a pie chart is drawn to show that the tailing of the piRNA of MIWI-IP is biased towards single base tailing;

[0138] • Tailing analysis (see Figure 10 ): According to the reads of the different degrees of non-template tailing piRNA existing in the small RNA combined with MIWI-IP, the proportion of each part (1 nt, 2 nt, 3 nt, 4 nt and ≥5 nt non-template tailing) is calculated using the R script viz_untailed.R, and a pie chart is drawn to show that the tailing of the piRNA of MIWI-IP is biased towards single base tailing;

[0139] (10) Overall piRNA analysis: The sequencing reads of untailed piRNA and tailed piRNA are summarized for the visualization analysis of overall piRNA. Visualization analysis of overall piRNA (Shell script work.sh):

[0140] • Length distribution analysis (see Figure 11 ): Combined with the sequencing reads, the overall length distribution of all piRNAs (tailed piRNA and untailed piRNA, without removing the tailing base of tailed piRNA) existing in the small RNA combined with MIWI-IP is analyzed, and the subsequent visualization is performed using the R script viz_total.R;

[0141] • Motif analysis (see Figure 12 ): Combined with the sequence information of all piRNAs (tailed piRNA and untailed piRNA, without removing the tailing base of tailed piRNA) existing in the small RNA combined with MIWI-IP, the Motif graph is drawn using the R script motif_code.R to show the sequence base distribution of piRNA;

[0142] • 5' end U bias analysis (see Figure 13): The first base information of all piRNAs (tailing piRNAs and un-tailing piRNAs, removing the tailing base of tailing piRNAs) present in the small RNAs bound by MIWI-IP is visualized using R script viz_total.R, showing the U-bias of the first base of piRNAs.

[0143] • piRNA source analysis (see Figure 14 ): Similar to step (8), the post-alignment sequencing reads of all piRNAs (tailing piRNAs and un-tailing piRNAs, removing the tailing base of tailing piRNAs) present in the small RNAs bound by MIWI-IP is visualized using R script viz_untailed.R, showing the source of piRNAs.

Claims

1. A method for analyzing piRNA, characterized in that: The following steps are involved: (0) performing data quality control on the sequencing reads of small RNA of wild-type mouse samples, wherein the data quality control includes removing the 3' end sequence adapter of the sequencing reads; (1) Align small RNA sequencing reads with the reference genome; (2) The reads successfully aligned in (1) are further aligned with the pre-miRNA reference sequence; for the sequencing reads that are successfully aligned to miRNA, the normalization parameters are calculated based on the number of aligned reads of miRNA in different samples; (3) After removing contamination from other small RNAs, the reads that were not successfully aligned to miRNAs in (2) were aligned to the piRNA reference sequence. The other small RNAs included tRNA, rRNA, and snRNA. (4) Reads that can be fully matched in (3) are defined as untailed piRNAs and normalized using the normalization parameters in step (2); (5) Reads that were not successfully aligned in step (1) are defined as tailed piRNAs, and the 1nt base at the 3' end is removed, and then steps (1)-(3) are repeated to obtain the 1nt tailed piRNAs. (6) Remove the 1 nt base at the 3' end of the reads that were not successfully aligned with the reference genome in (5), and then repeat steps (1)-(3) to obtain the piRNA with a 2nt tail. (7) Iterate according to the methods in (5) and (6) until no new tailed piRNAs can be obtained; then merge all tailed piRNAs and normalize them using the normalization parameters of step (2); Wherein, the reference genome is the mouse reference genome.

2. The analysis method according to claim 1, wherein Steps (1) to (7) are implemented using Bowtie software.

3. The method according to claim 2, wherein Step (3) includes sequentially aligning the reads that failed to align to miRNA to the tRNA, rRNA, and snRNA reference sequences to obtain reads that are not contaminated by other small RNAs.

4. The analysis method according to any one of claims 1 to 3, wherein In step (0), the data quality control also includes: removing PCR duplicates; removing random bases introduced before and after the sequence during library construction; setting a threshold to remove low-quality sequencing data; and filtering redundant sequence reads <18nt.

5. The analysis method according to any one of claims 1 to 3, wherein Before step (0), step (-1) is also included: constructing a library and sequencing the small RNA to obtain sequencing reads of the small RNA.

6. The analysis method according to claim 5, wherein Step (-1) includes: performing small RNA sequencing using the Illumina HiSeq×Ten platform to screen RNA with a fragment size of 10-50 nt.

7. The analysis method according to any one of claims 1 to 3, wherein The analysis method further comprises step (8): visually analyzing the untailed piRNA and the tailed piRNA.

8. The analysis method according to claim 7, wherein The visualization analysis includes one or more of the following: length distribution analysis, motif analysis, 5' end U base preference analysis, piRNA source analysis, and ratio analysis of tailed piRNA and untailed piRNA.

9. A computing device, characterized in that include: a memory for storing program instructions; A processor is configured to call the program instructions stored in the memory, and execute the analysis method according to any one of claims 1 to 8 according to the obtained program instructions.

10. A computer-readable storage medium, characterized in that The method comprises computer-readable instructions, which, when read and executed by a computer, enable the analysis method according to any one of claims 1 to 8 to be implemented.

11. A computer program product, characterized in that The invention comprises a computer program executable by a computer device, and when the program is run on the computer device, the computer device is caused to execute the steps of the analysis method according to any one of claims 1 to 8.

Images