Neonatal protein detection method and system based on full-length and short-segment transcriptome high-throughput sequencing technology

Abstract

The invention develops a sample newborn protein detection and analysis method, the method provides a one-key detection and analysis process, full-length and short-fragment RNA-seq data can be comprehensively analyzed, all variations which may induce potential new proteins, including large-scale variations and small-scale variations, are covered, and new transcripts are identified based on the various variations. Finally, all types of new transcripts are integrated and converted into a fasta-format protein sequence for subsequent protein analysis, such as protein structure prediction, neoantigen immunoepitope analysis and the like.

Description

Technical Field

[0001] The present invention relates to the field of bioinformatics and nascent protein identification and detection, and in particular, to a nascent protein detection method and system based on full-length and short-fragment transcriptome high-throughput sequencing technology. Background Art

[0002] Many fundamental problems in the field of human diseases, including cancer, are closely related to new proteins generated by base substitutions, deletions, insertions, frameshifts, intron retention, alternative splicing, and translation of new unannotated open reading frames (nuORFs) at the gene and transcription levels. Since new proteins do not exist in normal cells, they will not be immune-tolerant and are therefore easily recognized as foreign antigens. New proteins can be digested into short peptides, which are regarded as T cell antigens by the host system, leading to the stimulation of cellular immune responses and the elimination of cells expressing new proteins. Based on this principle, pharmaceutical companies or research institutions have developed new antigen-based immunotherapy as a promising cancer treatment method. In addition, new proteins are also related to other complex diseases, such as Alzheimer's disease, type II diabetes, liver disease, kidney disease, obesity, and cardiovascular disease. Therefore, the effective identification of new proteins is crucial to understanding human diseases and developing corresponding treatments.

[0003] The traditional method for identifying and characterizing proteins is to combine tandem mass spectrometry (MS) with database searching after mass spectrometry. MS relies on a reference database containing regular sequences in the canonical proteome, so it is almost impossible to identify proteins outside the canonical proteome, especially those arising from patient-specific individual variations. The MS peptidomics approach has other technical drawbacks: the detection rate of Cys-rich peptides is significantly reduced, and the detection sensitivity of low molecular weight (<500Da) or high molecular weight (>3kDa) peptides is low. Obviously, the bias caused by technical problems limits the application of MS in identifying new proteins. In the past decade, next-generation sequencing (NGS) has made significant progress, achieving sample-specific whole genome sequencing or total RNA sequencing, and the two methods of MS and NGS can be combined for new protein analysis. However, the following shortcomings still exist: (1) Due to the short read length of next-generation sequencing, there are inherent technical defects in detecting large-scale variants; (2) Standard protein identification protocols are usually redundant and require detection at the protein level by MS and at the RNA level by NGS, which has high requirements for input samples; (3) Most of the existing bioinformatics analysis workflows based on nucleotide sequencing only analyze general routine transcripts or focus on single or a few variant categories, and therefore cannot cover all sample-specific variants.

[0004] In summary, there is an urgent need in the art for a method and system for detecting new proteins that can cover multiple types of variants and is efficient and accurate. Summary of the invention

[0005] The purpose of the present invention is to provide a new protein identification method and system based on full-length and short-fragment transcriptome sequencing data, specifically, it involves using RNA-seq data from NGS and HiFi to finally generate new protein sequences through quality control, variant detection, annotation and hierarchical filtering processes.

[0006] The first aspect of the present invention provides a method for detecting a new protein, the method comprising the steps of:

[0007] (A) providing raw data and constructing a library, wherein the data includes second-generation (NGS) transcriptome sequencing data and third-generation transcriptome sequencing data, wherein the third-generation transcriptome sequencing data is HiFi data;

[0008] (B) Preprocessing and quality control of raw data, including:

[0009] The preprocessing comprises the steps of: polishing the third-generation transcriptome sequencing data; aligning the second-generation transcriptome sequencing data to the reference genome to obtain a first BAM file; trimming the third-generation transcriptome sequencing data to obtain a full-length non-chimeric sequence (FLNC), aligning the FLNC to the reference genome to obtain a second BAM file; the two BAM files are used for small-scale variation analysis; clustering the FLNC to obtain a third BAM file, and the BAM file is used for large-scale variation analysis;

[0010] The quality control includes: calculating the number of zero-mode waveguide holes (ZMWs); calculating the number of reads aligned to a unique position in the reference genome; calculating the number of reads aligned to multiple positions in the reference genome; calculating the number of reads that cannot be aligned to the reference genome; using the BAM file to calculate the number and proportion of ribosome-associated RNA bases, intron bases, intergenic region bases and UTR region bases; calculating the sequencing depth coverage of the 5' end and the 3' end; calculating the proportion of reads aligned to the positive and negative chains of the reference genome;

[0011] (C) Detection of new variant isoforms (NVI);

[0012] (D) performing hierarchical filtering and annotation on the detection results to obtain filtered NVI, wherein the filtered NVI includes NVI in nucleotide sequence format and NVI in transcript format;

[0013] (E) Generation of new proteins using filtered NVI.

[0014] In another preferred example, the second-generation transcriptome sequencing data includes Illumina sequencing data.

[0015] In another preferred embodiment, the original data includes one or more selected from the following group: human brain reference RNA samples spiked with 2% SIRV-Set4, and tumor allogeneic transplant RNA samples.

[0016] In another preferred embodiment, the method for obtaining the tumor xenograft RNA sample comprises: subcutaneously injecting MC38 cell line into C57BL / 6 mice.

[0017] In another preferred embodiment, the third-generation transcriptome sequencing data is polished by circular consensus sequencing (CCS).

[0018] In another preferred embodiment, the FLNC is clustered by a method selected from the group consisting of hierarchical alignment, iterative clustering, or a combination thereof.

[0019] In another preferred embodiment, the third-generation transcriptome sequencing data is pruned by a method selected from the following group:

[0020] (a) trimming the barcodes of the third-generation transcriptome sequencing data using Lima software;

[0021] (b) trimming the ployA tail and adapter sequence of the third-generation transcriptome sequencing data using Iso-Seq3 software;

[0022] or a combination thereof.

[0023] In another preferred example, whether the library has residual ribosome-associated RNA is determined by calculating the number and ratio of ribosome-associated RNA bases.

[0024] In another preferred embodiment, the number and proportion of the ribosome-associated RNA bases is 0%.

[0025] In another preferred embodiment, whether the library is degraded is determined by calculating the sequencing depth coverage of the 5' and 3' ends.

[0026] In another preferred embodiment, the sequencing depth coverage of the 5' end is ≥0.7.

[0027] In another preferred embodiment, the sequencing depth coverage of the 3' end is ≥0.7, preferably ≥0.9.

[0028] In another preferred example, the properties of the library are determined by calculating the ratio of reads aligned to the positive and negative strands of the reference genome.

[0029] In another preferred embodiment, the library properties include non-strand-specific library and strand-specific library.

[0030] In another preferred example, the ratio of reads of the non-strand-specific library mapped to the positive and negative chains of the reference genome remains balanced, and the "maintaining balance" means that the number of reads N1 mapped to the positive chain of the reference genome, compared with the number of reads N2 mapped to the negative chain of the reference genome, satisfies N1 / N2=100±5%.

[0031] In another preferred embodiment, the ratio of reads of the chain-specific library mapped to the positive and negative chains of the reference genome is significantly unbalanced, and the "significant imbalance" means that the number of reads N1 mapped to the positive chain of the reference genome, compared with the number of reads N2 mapped to the negative chain of the reference genome, satisfies N1 / N2≥1900%, or N2 / N1≥1900%.

[0032] In another preferred embodiment, in step (C), the step is included:

[0033] (C1) using collapse_isoforms_by_sam.py of cDNA_cupcake to remove redundancy in the HiFi data and detect NVI; using fusion_finder.py of cDNA_cupcake to detect gene fusion; using PBsv software to detect NVI; removing the results obtained from multiple software to obtain the detection results of HiFi data;

[0034] (C2) Using Vardict software to detect NVI in the NGS data, and obtaining a VCF file of the detection results of the NGS data.

[0035] In another preferred embodiment, in step (C1), redundancy is removed by folding identical transcript isoforms into unique transcript isoforms.

[0036] In another preferred embodiment, in step (C1), the number of full-length reads of the unique transcript isoform is quantified to obtain a quantified value for use in step (D).

[0037] In another preferred embodiment, in step (D), the layered filtration comprises:

[0038] (D1) filtering isoforms / fusions;

[0039] (D2) Filtering of SNV / INDELs;

[0040] Wherein, in step (D1), the steps are included:

[0041] (d1.1) remove NVI with abundance below the threshold;

[0042] (d1.2) removal of NVI lacking a 5' exon;

[0043] (d1.3) noise removal NVI;

[0044] (d1.4) annotating the remaining NVIs using a database, wherein the database is the chimerDB 4.0 database;

[0045] In step (D2), comprising:

[0046] (d2.1) Removing NVIs whose sequencing depth is lower than a threshold. The threshold of the sequencing depth is calculated as follows:

[0047]

[0048] Among them, RPK refers to reads per kilobase; TPM refers to transcripts per million; mean_read_length refers to the average length of reads;

[0049] (d2.2) Remove NVI with noise mark in the FILTER column of the VCF file;

[0050] (d2.3) Removal of NVIs from SNPs or germline mutations;

[0051] (d2.4) The remaining NVIs are annotated using databases, including the COSMIC database.

[0052] In another preferred embodiment, the abundance threshold is 2.

[0053] In another preferred example, the TPM=35.

[0054] In another preferred example, the noise indicators include P8, Q10, and bias.

[0055] In another preferred example, NVIs from SNPs or germline mutations are removed based on dbSNP annotation information and allele frequency (AF).

[0056] In another preferred embodiment, in step (E), the step of:

[0057] (E1) converting the NVI in the nucleotide sequence format into the NVI in the transcript format;

[0058] (E2) merging the NVIs of the two transcript formats and removing redundancy, and storing the NVIs without redundancy in a fasta file;

[0059] (E3) converting the redundantly removed NVI into a protein sequence using GeneMarkS-T software;

[0060] (E4) Using BLAST software, the protein sequence is compared with the UniProt protein database to determine whether the protein sequence is a new protein.

[0061] The second aspect of the present invention provides a software or system for identifying a nascent protein, the software or system comprising:

[0062] An input unit, wherein the input unit is configured to input sequence data, wherein the input sequence data is NGS transcriptome sequencing data and third-generation transcriptome sequencing data;

[0063] A detection module, wherein the detection module is configured to perform a predetermined detection on the input sequence data to obtain an analysis result, and the detection module comprises:

[0064] (Z1) a preprocessing and quality control module, wherein the preprocessing and quality control module is configured to perform step (B) of the method according to the first aspect of the present invention, thereby obtaining a result of preprocessing and quality control;

[0065] (Z2) a detection module, wherein the detection module is configured to perform step (C) of the method according to the first aspect of the present invention, thereby obtaining a detection result;

[0066] (Z4) a filtering and annotation module, wherein the preprocessing and quality control module is configured to execute step (D) of the method according to the first aspect of the present invention, thereby obtaining filtering and annotation results;

[0067] An output unit, wherein the output unit is configured to execute step (E) in the method described in the first aspect of the present invention, thereby outputting the detection result of the detection module.

[0068] The third aspect of the present invention provides a protein generated by the method of the first aspect of the present invention, or the software or system of the second aspect of the present invention, wherein the protein is encoded by a nucleotide sequence as shown in SEQ ID NO: 1-4.

[0069] The fourth aspect of the present invention provides a polynucleotide comprising a nucleotide sequence as shown in SEQ ID NO: 1-4.

[0070] The fifth aspect of the present invention provides a composition comprising the protein described in the fourth aspect of the present invention.

[0071] In another preferred embodiment, the composition further comprises a pharmaceutically acceptable carrier, diluent or excipient.

[0072] The sixth aspect of the present invention provides a use of the protein described in the third aspect of the present invention, the polynucleotide described in the fourth aspect of the present invention, or the composition described in the fifth aspect of the present invention for preparing an agent or drug for treating tumors.

[0073] In the seventh aspect of the present invention, there is provided a method according to the first aspect of the present invention or a use of the software or system according to the second aspect of the present invention for protein analysis, wherein the analysis includes protein structure prediction and new antigen immune epitope analysis.

[0074] It should be understood that within the scope of the present invention, the above-mentioned technical features of the present invention and the technical features specifically described below (such as embodiments) can be combined with each other to form a new or preferred technical solution. Due to space limitations, they will not be described one by one here. BRIEF DESCRIPTION OF THE DRAWINGS

[0075] Figure 1 The sample novel protein detection and analysis process of the present invention is shown. By comprehensively analyzing NGS and HiFi sequencing data, all new variants and isomers are detected in an integrated manner, including single nucleotide variants (SNVs), insertion-deletion mutations (INDELs), gene fusions (FUSION), intron retention, and ISOFORM of alternative splicing. The intermediate steps of the analysis include multiple strict filtering, and finally obtain results with low background noise. Finally, the converted amino acid sequence is compared with the protein database to select new or classic proteins.

[0076] Figure 2 Categories of novel variant isoform (NVI) detection are shown. For small-scale variants, including (A) point mutations and (C) gene fusions, analysis was performed using short-read sequencing reads. For large-scale variants, including (B) RNA alternative splicing, (C) gene fusions, and (D) NuORFs, analysis was performed using long-read sequencing data. RNA fusions were analyzed using both long-read and short-read data.

[0077] Figure 3 (A) illustrates the folding process of long read sequences. Because the introns of sequences a and b, and sequences c and d are exactly the same, they are folded into unique potential transcript candidates, while sequences a and c cannot be further folded due to the differences in introns. (B) shows that many sequences cannot be aligned at the 5' end, but can be aligned at the 3' end. Such sequences may be biologically 5' degradations of the same original transcript, so they are folded into a unique transcript isomer.

[0078] Figure 4 The hierarchical filtering of the preliminary detection results of all new variant isomers (NVIs) is shown. The inverted triangles represent the logic of filtering. Each layer represents a step of filtering operation.

[0079] Figure 5 The results of the MC38 analysis are shown. (A) The analysis results show the nucleotide and corresponding protein sequences of the variants studied under the Deliverables section. (B) The results of the PB1596.5 variant displayed in IGV using the PacBio BAM file. (C) The final new protein sequence obtained from the integrated analysis. (D) Comparison of the protein sequence of the PB1596.5 variant with Q1KYM0 using BLAST. (E) The electrophoresis results of PB1596.5, from left to right, are the Marker group, MC38 cells, normal tissue control group, and MC38 allogeneic transplants. DETAILED DESCRIPTION

[0080] After extensive and in-depth research, the inventors have creatively developed a method for detecting and analyzing new proteins in samples. The method provides a one-click detection and analysis process that can comprehensively analyze full-length and short-fragment RNA-seq data, covering all possible mutations that may induce potential new proteins, such as large-scale mutations including large gene fusions, transcript exon skipping, intron retention, RNA alternative splicing, and abnormal expression of non-coding regions, as well as small-scale mutations including insertions, deletions, duplications, substitutions, and small-scale gene fusions, and identifies new transcripts based on the above-mentioned multiple mutations. Finally, all types of new transcripts are integrated and converted into protein sequences in fasta format for subsequent protein analysis, such as protein structure prediction, new antigen immune epitope analysis, etc.

[0081] It should be understood that the specific methods and experimental conditions of the present invention are described below in various degrees of detail to provide a substantial understanding of the present invention. The following provides definitions of certain terms used in this specification. Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by those of ordinary skill in the art to which the present invention belongs.

[0082] the term

[0083] Where a numerical range is provided, it is understood that every intermediate integer of the value 20, every tenth of every intermediate integer of the value, between the upper and lower limits of the range, and any other intermediate values ​​in the specified range are included in the invention, unless the context clearly indicates otherwise. The upper and lower limits of these smaller ranges may be independently included in the smaller ranges and are also encompassed in the invention, subject to any explicitly excluded limitations in the specified range. For example, "1 to 50" includes "2 to 25", "5 to 20", "25 to 50", "1 to 10", etc.

[0084] As used herein, the term "comprising" or "including (comprising)" may be open, semi-closed and closed. In other words, the term also includes "consisting essentially of" or "consisting of".

[0085] As used herein, the terms "nascent protein", "new protein" and "novel protein" are used interchangeably and refer to non-naturally occurring proteins. In the present invention, the nascent protein can be obtained from long-length data and short-read data by the detection method provided by the present invention.

[0086] As used herein, the terms "comparison" and "alignment" are used interchangeably and refer to arranging two or more sequences together to compare similarities. The present invention compares transcriptome data with a reference genome to obtain a BAM file.

[0087] As used herein, the term "barcode" refers to a unique identifier for a single sample in a mixed sample, which is used to distinguish different samples. Since the sequencing capacity of the test instrument is much greater than the number of test samples, multiple samples are often measured in the same lane to avoid waste. Specific identifiers are added to different samples to distinguish different sample data during subsequent data analysis.

[0088] As used herein, the term "gene" refers to a nucleotide sequence in DNA or RNA that can encode and synthesize gene products including RNA and proteins. A gene is generally regarded as a basic genetic unit, i.e., a functional DNA or RNA.

[0089] As used herein, the term "transcript" refers to one or more mature RNA molecules formed by transcription of a gene. The RNA formed by gene transcription can form different transcripts through different processing processes. The different transcripts are "transcript isomers".

[0090] As used herein, the terms "new variant isomer" and "NVI" are used interchangeably and refer to the results obtained based on genetic variation data using the detection method or detection system of the present invention.

[0091] As used herein, the term "single nucleotide variation (SNV)" refers to a variation occurring in a single nucleotide, and its variation frequency is not limited.

[0092] As used herein, the term "insertion-deletion mutations (INDELs)" refers to the insertion or deletion of small fragments in a genomic sequence, with a length of between 1-50 bp.

[0093] As used herein, the terms "short fragment" and "short read length" are used interchangeably and refer to sequences obtained using second-generation sequencing technology. Second-generation sequencing technology has low sequencing costs, high sequencing speeds, and high accuracy. The read length of second-generation sequencing is usually no more than 500bp. Due to the short read length, short read length data usually contains small-scale variations. In the present invention, short fragment data is used for small-scale variation analysis.

[0094] As used herein, the terms "full length" and "long read length" are used interchangeably and refer to sequences obtained using third-generation sequencing technology. The read length of third-generation sequencing can usually reach several thousand or even base pairs, so it can contain information on large-scale variations. In the present invention, the data from third-generation sequencing is mainly used for large-scale variation analysis, and can also be used for small-scale variation analysis.

[0095] Sample nascent protein detection and analysis method

[0096] The sample nascent protein detection and analysis method (hereinafter referred to as the detection method) uses Snakemake as its workflow management system to achieve controlled and scalable process execution in a high-performance computing environment. Snakemake is an efficient and general workflow management tool that provides many advantages for scientific data analysis, including excellent memory management, strong portability, modular design and repeatable results. The overall process of the detection method of the present invention is shown in Figure 1 .

[0097] As a comprehensive bioinformatics pipeline, the assay is compatible with a variety of data types, including NGS and HiFi RNA sequencing data. It provides a series of analysis steps, including quality control, variant detection, annotation and stratified filtering, and new protein sequence generation. The assay provides a wide range of variant detection categories, such as insertions, deletions, frameshifts, gene fusions, intron retention, and alternative splicing, to enable comprehensive and detailed analysis of transcriptome data. These variant categories help identify new transcripts or proteins, which are referred to herein as novel variant isoforms (NVIs).

[0098] To use this detection method flow, raw data from a subreads.bam file generated by HiFi sequencing is required, as well as a configuration file containing optional parameters and necessary reference and annotation file paths. Raw data generated by the NGS platform can also be included if available. Users can use the default parameters or adjust the parameters according to specific needs before the initial test run. After local deployment, users only need to specify sample-specific raw data paths for subsequent analysis. The detection method flow contains more than twenty steps and can be divided into four different parts (such as Figure 1 ):

[0099] 1) Quality control and preprocessing of raw data,

[0100] 2) Detection of new variant isomers (NVI),

[0101] 3) Hierarchical filtering and annotation of detection results,

[0102] 4) Generate protein sequences from the filtered NVI in fasta format and identify novel proteins not present in the canonical proteome.

[0103] Circularization consensus sequencing

[0104] As used herein, the terms "circularization consensus sequencing", "circularization sequencing", and "CCS" can be used interchangeably, and high-precision single-molecule consistency HiFi reads can be obtained through circularization consensus sequencing. The CCS process includes cutting the "polymerase reads" that appear alternately between the linker and the insert sequence to remove the linker, thereby obtaining multiple insert sequences called "sub-reads". These sub-read sequences are generated by circularization sequencing of the insert sequence. Sub-reads from the same source insert sequence are largely identical, but sequencing errors introduce random errors in these sub-reads from the same insert sequence. To solve this problem, the sub-reads are aligned to each other and polished. In general, as long as three or more circularization sequencing is performed, that is, three or more sub-read sequences are obtained, CCS polishing can produce HiFi reads with an accuracy of up to 99.9%. The CCS process is the most time-consuming step in the present invention.

[0105] Genetic variation

[0106] Genetic variation in the present invention refers to all variations in genes that may generate new proteins, including small-scale variations and large-scale variations. The small-scale variations include short insertions, deletions, duplications, substitutions and small-scale gene fusions; the large-scale variations include large-scale gene fusions, transcript exon skipping, intron retention, RNA alternative splicing and abnormal expression of non-coding regions, as well as long insertions, deletions and duplications.

[0107] The main advantages of the present invention include:

[0108] (1) The detection method of the present invention integrates long-read and short-read sequencing technologies and combines them with innovative bioinformatics analysis processes to comprehensively detect genetic variations of various scales, including point mutations, insertions and deletions, gene fusions, intron retention, etc., thereby identifying more potential new proteins and solving the limitations of detecting large-scale and small-scale variations and analyzing complex samples.

[0109] (2) The detection method of the present invention adopts strict quality control and filtering steps to reduce the false positive rate and improve the reliability and accuracy of the detection results.

[0110] (3) The detection method of the present invention accelerates the analysis process and improves detection efficiency by optimizing algorithms and parallel computing.

[0111] The present invention will be further described below in conjunction with specific examples. It should be understood that these examples are intended to illustrate the present invention only and are not intended to limit the scope of the present invention. The experimental methods in the following examples where specific conditions are not specified are usually performed under conventional conditions, such as those described in Sambrook et al., Molecular Cloning: A Laboratory Manual (New York: Cold Spring Harbor Laboratory Press, 1989), or under conditions recommended by the manufacturer. Unless otherwise indicated, percentages and parts are weight percentages and weight parts.

[0112] Example 1: Quality control and preprocessing of raw data

[0113] This embodiment involves quality control and preprocessing of raw HiFi sequencing data.

[0114] This embodiment uses the open source software provided by PacBio to preprocess the data. Specifically, the circular consensus sequencing (CCS) software is used to obtain high-precision single-molecule consensus HiFi reads. This process involves cutting the "polymerase reads" that appear alternately between the linker and the inserted sequence to remove the linker, thereby obtaining multiple inserted sequences called "sub-reads". These sub-read sequences are generated by circular sequencing of the inserted sequence. Sub-reads from the same source insert sequence are largely identical, but sequencing errors introduce random errors in these sub-reads from the same insert sequence. To solve this problem, the sub-reads are aligned to each other and polished. In general, as long as three or more circular sequencings are performed, that is, three or more sub-read sequences are obtained, CCS polishing can produce HiFi reads with an accuracy of up to 99.9%. Among them, the CCS step is the most time-consuming step in the entire process.

[0115] After generating HiFi reads, barcodes were trimmed using Lima software to distinguish multiple samples constructed in the same library pool. This step can be omitted if the library contains only one sample. Subsequently, the polyA tails and adapter sequences at both ends of the HiFi reads were trimmed using the refine module of Iso-Seq3 software to obtain full-length non-chimeric (FLNC) sequences, while also removing residual tandem repeats. The obtained FLNC sequences were aligned to the reference genome to generate a BAM file. At the same time, the NGS data were aligned to the same reference genome to generate another BAM file. These two BAM files will be further used for small-scale variation analysis. In addition, the FLNC reads were clustered using Iso-Seq3 software, using hierarchical n*log(n) alignment and iterative cluster merging to identify other large-scale variations.

[0116] Data preprocessing involves multiple processing steps and analysis modules, which may introduce anomalies and affect subsequent analysis. Therefore, multidimensional quality control is implemented at each step to ensure downstream analysis and the final results about new proteins. For example, the number of holes used for single-molecule sequencing, called Zero-Mode Waveguide (ZMWs), represents the amount of data. In addition, the polished HiFi reads and the FLNC reads remaining after removing the polyA tail must reach the required accuracy. In addition, indicators such as "uniq_mapped_reads (reads aligned to a unique position in the reference genome)", "NonUniq_mapped_reads (reads aligned to multiple positions)" and "Unmapped_reads (unmatched reads)" are also used to evaluate the quality of sequencing data. These quality control indicators can be used for analysis at each step and provide feedback on analysis problems. For example, if the number of ZMWs is too low, it indicates that the input library for sequencing is insufficient, prompting researchers to increase sample input to increase data volume. Preferably, ZMWs ≥ 1M (million); if the proportion of Unmapped_reads is too high, it may indicate contamination from other species.

[0117] The number and ratio of RIBOSOMAL_BASES, INTRONIC_BASES, INTERGENIC_BASES, and UTR_BASES were calculated using BAM files to assess library quality. A high ratio of INTERGENIC_BASES may indicate DNA contamination, while a high ratio of RIBOSOMAL_BASES indicates problems with ribosomal RNA removal, which can lead to reduced detection sensitivity. The geneBody_coverage.py module in the RSeQC software was used to analyze the sequencing depth coverage of the 5' and 3' ends of the gene to obtain MEDIAN_5PRIME_BIAS and MEDIAN_3PRIME_BIAS. A low MEDIAN_5PRIME_BIAS value indicates a significant reduction in the coverage depth of the 5' end of the transcript, which may indicate severe library degradation. Combined with the RIN value for assessing RNA quality, the above analysis helps to fully determine whether resequencing is necessary.

[0118] In addition, the ratio of reads that aligned to the "+" and "-" strands on the reference genome was determined using the infer_experiment.py module in the RSeQC software. A balanced ratio indicates a non-strand-specific library construction, while a significantly unbalanced ratio indicates a strand-specific library construction. Strand-specific libraries can be used to analyze gene expression and isoform detection of two genes that overlap.

[0119] Example 2: Detection of new variant isomer (NVI)

[0120] This embodiment involves using HiFi data as well as NGS data to perform large-scale variation and small-scale variation analysis to detect new variant isoforms.

[0121] After the raw data was preprocessed as in Example 1, new variant isomers (NVI) were detected. All genetic variations that may generate new proteins were divided into two categories: small-scale variations (such as short insertions, deletions, duplications, substitutions, and small-scale gene fusions) and large-scale variations (such as large-scale gene fusions, transcript exon skipping, intron retention, RNA alternative splicing, abnormal expression of non-coding regions, and long insertions, deletions, and duplications). The types are shown in Figure 2 .

[0122] In order to simultaneously analyze variants of different scales, this embodiment integrates various variant analysis software into the detection method to meet the analysis requirements of various variants and improve the overall analysis performance. This embodiment optimizes the analysis method by utilizing the characteristics of the data. Due to the advantages of long read length sequences, HiFi reads are mainly used to analyze large-scale variants; NGS reads are mainly used to analyze small-scale variants due to their higher accuracy at the single base level and easier acquisition of high-depth sequencing data.

[0123] For HiFi reads, collapse_isoforms_by_sam.py module of cDNA_cupcake was used to collapse identical transcript isoforms into unique transcript isoforms. The number of full-length reads corresponding to different unique isoforms was quantified to be used as parameters in subsequent filtering steps. Gene fusions were detected using fusion_finder.py of cDNA_cupcake, as they can be an important source of new or chimeric proteins. HiFi reads can also be used to analyze small-scale variants, for example, various structural variants including deletions of about 10 bp were analyzed using PBsv, which may also appear in NGS variant analysis results. Redundancy in variant results obtained from multiple software was removed.

[0124] For NGS reads, Vardict software was used to analyze small-scale variants and generate VCF files that store all possible variants. Variants that do not generate new proteins, such as single base substitutions in non-coding regions or synonymous, nonsense, or silent mutations, were filtered out.

[0125] Large-scale variation based on NGS data analysis was compared and integrated with the results of HiFi data analysis to further improve the reliability of the results.

[0126] Example 3: Layered Filtration

[0127] 3.1 Filtering of isomers / fusions

[0128] Since fusions are also considered a specific type of isomer, the filtering process and steps are basically the same as for isomers. Filtering isomers for detection usually involves 4 steps (e.g. Figure 4 ):

[0129] (S1) Remove isoforms below a preset threshold (default is 2, high-quality isoforms should be supported by at least 2 FLNC reads) based on abundance. Typically, the abundance of each detected isoform can be quantified directly from the PacBio HiFi sequencing data. However, if additional NGS sequencing data is available from another aliquot of the same sample, the abundance of the isoform can be quantified more accurately.

[0130] (S2) Removal of some isomers (see Figure 3 ), these isoforms have the same exons at the 3' end as the longest isoform, but lack some 5' exons, indicating that these 5' truncated transcripts may be caused by 5' end degradation of RNA during library construction. The more severely the library is degraded, the more transcripts are filtered out due to 5' truncation.

[0131] (S3) Remove "artificial" isomers introduced by PCR chimerism, restriction enzyme ligation or sequencing errors. These artificial isomers are real signals detected from sequencing data, but they are noise to the sample and must be removed. Since it is difficult to identify the random interference that leads to artificial isomers with a single feature, this embodiment uses the SQANTI3 tool constructed based on a random forest machine learning method to solve this problem. SQANTI3 extracts up to 47 features, including coverage, exon coordinates, number of supported reads, etc., and effectively implements filtering based on these features.

[0132] (S4) Annotation with database. This step can identify meaningful isoforms / fusions associated with clinical diseases and potential therapeutic targets. The level of detail and accuracy of gene annotation results are highly dependent on the quality and quantity of the annotation database. For human fusion gene annotation, the chimerDB 4.0 database is used, which summarizes nearly 100,000 fusion candidate genes of patients in the Cancer Genome Atlas (TCGA) project, including tens of thousands of fusion genes with experimental evidence support.

[0133] 3.2 SNV / INDEL Filtering

[0134] SNV / INDEL filtering involves four steps (see Figure 4 ):

[0135] (X1) Remove variants whose sequencing depth is below the confidence threshold. The depth threshold is calculated using the following formula:

[0136]

[0137] Among them, RPK refers to reads per kilobase, which is calculated by dividing the read count by the kilobase length of each gene; TPM refers to transcripts per million, which is calculated by dividing the PRK value by the "per million" scaling factor; mean_read_length refers to the average read length.

[0138] Since peptides translated from variant transcripts with a TPM higher than 35 are considered immunogenic, the TPM is set to 35 and the threshold for variant depth can be calculated using the formula.

[0139] (X2) Remove background noise, mainly using the FILTER column in the VCF file to remove variants with P8, Q10, bias and other flags.

[0140] (X3) Remove variants from SNPs or germline mutations by comprehensive analysis of dbSNP annotation information and allele frequency (AF).

[0141] (X4) Use databases such as COSMIC to annotate and identify mutations associated with clinical diseases.

[0142] Example 4: New protein generation

[0143] This embodiment involves integrating the results of different analysis software, converting the file format and removing redundancy to obtain the final new protein generation results.

[0144] All kinds of variant results from different analysis software were filtered in layers to obtain their respective final results. To facilitate subsequent analysis, the variant nucleotide sequences were first converted to transcript format. Transcript isoform and fusion analysis results are already in nucleotide sequence format and do not require further conversion. However, small-scale variants in VCF format of NGS data obtained by Vardict and structural variants in HiFi reads obtained by PBsv need to be converted to transcript format. After integrating the result variant files and removing redundant variants, the integrated variant transcripts are stored in fasta files.

[0145] After transcript analysis, the expressed protein sequences were checked using GeneMarkS-T software to convert transcript sequences to corresponding protein sequences. Due to the presence of multiple open reading frames, longer transcripts will produce multiple protein sequences. The variants obtained through this analysis usually constitute new proteins, including tumor-specific proteins. The above proteins were compared with the protein database downloaded from UniProt using BLAST software to determine whether the proteins produced by the variants have been studied and published, so as to facilitate manual judgment and selection of the final results in a specific research context. Published protein sequences are not eliminated arbitrarily, but possible homologous protein IDs and consistency values ​​are annotated in the fasta file. In addition, some published proteins may be specifically expressed in tumor tissues and are therefore also considered relevant to the research purpose.

[0146] Example 5: Verification of model performance

[0147] This embodiment relates to verifying the working principle of the detection method of the present invention, including using two kinds of data or samples to verify it.

[0148] 5.1 Dataset

[0149] The first sample was a human brain reference RNA sample spiked with 2% SIRV-Set 4. SIRV-Set 4 is a biological reference standard containing a set of artificial isoforms. PacBio sequencing data of this RNA sample was obtained from NCBI in the format of FLNC.read.bam, indicating that it had been preprocessed. SIRV-Set 4 contains 69 artificial transcript variants that simulate the splicing characteristics of 7 human model loci, fully reflecting the diversity of alternative splicing, alternative transcription start and stop sites, overlapping genes, and antisense transcripts. In addition, SIRV-Set 4 also includes 15 artificial transcript variants in five length categories (4kb, 6kb, 8kb, 10kb, and 12kb, respectively) to simulate the length complexity of human transcripts.

[0150] The second sample was obtained by subcutaneously injecting the MC38 cell line into C57BL / 6 mice to form tumor allogeneic grafts. After killing the mice, the tumor allogeneic grafts and paired normal small intestinal tissues were collected and divided into small pieces of about 300 mg. These tissue pieces were immediately washed with 1xPBS, placed in separate centrifuge tubes, and transported to GeneDenovo Co., Ltd. at low temperature for RNA extraction, PacBio HiFi long read sequencing, and Illumina short read sequencing. After obtaining the sequencing data, the detection method of the present invention was used for analysis.

[0151] Since there are significant differences between tumor allogeneic transplants and normal tissues in transcriptome and proteome, the purpose of validation is to use this detection method to identify mutations and abnormal transcript expression unique to tumor allogeneic transplants and obtain corresponding novel proteins.

[0152] 5.2 Verification Results

[0153] The detection method of the present invention uses PacBio HiFi long-read data to analyze human brain reference RNA samples, and uses two types of read data to analyze MC38 allogeneic transplant samples.

[0154] Since the human brain reference RNA samples have been pre-processed, including a CCS polishing step, the validation process starts with the alignment of the pre-polished HiFi reads using MINIMAP2. Since the validation dataset lacks short-read data, this validation only focuses on analyzing large-scale variants, especially the different isoforms incorporated into SIRV-Set 4. The MC38 allogeneic transplant samples contain both PacBio HiFi long reads and Illumina short reads, so both large-scale and small-scale mutation types need to be analyzed.

[0155] Since both data sets were raw data, the analysis of MC38 allogeneic transplantation started with data preprocessing and quality control steps, and all analysis modules for nascent protein detection were executed. The nascent protein detection was run on a high-performance computer cluster equipped with Intel(R) Xeon(R) Platinum processors (8269CY CPU, 3.10GHz) and Ubuntu 20.04.4LTS Linux operating system. After about 1 hour, the detection method of the present invention completed the analysis of the human brain reference RNA sample data; while the analysis of the MC38 allogeneic transplantation data took nearly 5 hours. By statistically analyzing the running time of each module, it was found that the CCS step was the most time-consuming. The human brain reference RNA samples that did not require the CCS step saved a lot of time.

[0156] For the human brain reference RNA sample, the quality control indicators were first evaluated, and the total number of HiFi reads obtained was 1.9 million, which almost occupied half of the SMRT Cell chip data capacity. The PCT_RIBOSOMAL_BASES indicator was 0%, indicating that the presence of tRNA was negligible. The MEDIAN_5PRIME_BIAS and MEDIAN_3PRIME_BIAS values ​​were 0.73 and 0.99, respectively, indicating that the RNA library was relatively intact and no obvious degradation occurred. The "++--" and "+-,-+" values ​​were 0.859 and 0.0042, respectively, indicating that the library construction was strand-specific. The above quality control indicators together indicate that the quantity and quality of the data meet the required standards.

[0157] Since SIRV-Set 4 spiked with human brain reference RNA samples is a biological reference standard with accurate isoform information, this example confirms the detection results of all true isoforms, and the detection rate is calculated to determine whether the new protein detection method of the present invention will mistakenly miss true positives. For the 12 long transcripts of SIRV-set 4, the detection method of the present invention achieved a detection rate of 100% (12 / 12) (see Tables 1 and 2). However, for the detection of 69 short SIRV transcripts, the detection method of the present invention only achieved a detection rate of 89.9% (62 / 69), missing 7 isoforms. Afterwards, the undetected isoforms were explored, and it was found that the undetected isoforms were due to their short introns or exons, resulting in high mapping scores in the global alignment, and were subsequently misclassified as soft clipping (Soft clip) or deletions. Statistics of the human RefSeq transcriptome show that the vast majority of exons and introns are more than 20bp in length. The exons and introns of these three misclassified isoforms are only 9 bp, 31 bp, and 20 bp, respectively, representing extreme cases. The other four isoforms missing in the original data indicate that they were lost during library construction or nucleotide sequencing.

[0158] Table 1 The detection methods correctly detected the short SIRV isomers

[0159]

[0160]

[0161] Table 2 The detection methods correctly detected the isomers of long SIRV

[0162]

[0163] For the MC38 tumor allogeneic graft, this example first performed a quality control assessment on it. The results showed that the quantity and quality of the data were sufficient for subsequent analysis. A large number of high-confidence small-scale variants were detected on short-read data using Vardict and PBsv; large-scale variants, such as intron retention, exon skipping, and abnormal expression of non-coding regions, were mainly detected by analyzing long-read HiFi data using cDNA_cupcake, including gene fusions identified by fusion_finder. These large-scale variants were integrated with small-scale variants after SQANTI and hierarchical filtering and stored in separate fasta files (such as Figure 5 A). Regarding large-scale genetic variation, no clear evidence of gene fusions or structural alterations was found, which is consistent with previous observations of the MC38 cell line. However, this example confirms the existence of five transcript variants of a specific transcript isoform, and these variants are strongly supported by evidence. Four of these variants, named PB.349.2, PB.458.1, PB.1517.1 and PB.1528.1, are located in the Cep131, F2RL1, Gm5345 and KLHL26 loci of the mouse genome, respectively. All four variants express an additional segment in the intronic region as a new exon. The remaining transcript variant PB.1596.5 is located in the intergenic region of the mouse genome and is abnormally expressed as a new transcript (such as Figure 5 B). All of these variants were initially identified using HiFi long reads and subsequently validated by NGS short reads at their respective splice sites. None of these variants were detected in the corresponding paired normal tissues, indicating that they are tumor specific.

[0164] In addition, this example uses GMST to generate corresponding protein sequences for these variants (such as Figure 5C), and homology analysis was performed by BLAST with the UniProt database tr_mouse_canon_isoform.fasta. The results showed that the identity values ​​of the highest homologous proteins of the first four variants (including PB.349.2, PB.458.1, PB.1517.1, and PB.1528.1) were less than 60%, indicating their novelty. However, PB.1596.5 showed a different pattern, and the alignment results showed that its identity value with the most homologous protein Q1KYM0 was 99% (667 / 669) (e.g. Figure 5 D). Q1KYM0 is an endogenous retroviral (ERV) protein that has been reported to be expressed in Neuro-2a tumor cells and can lead to immune escape. Knockdown of the gene expressing Q1KYM0 with siRNA delayed tumor growth and increased mouse survival. These findings strengthen the evidence supporting the tumor specificity of PB.1596.5 and highlight the importance of combining large-scale and small-scale variant detection and subsequent analysis of novel proteins.

[0165] Finally, to further confirm the reliability of the identified variant results, primers for PB.1596.5 were designed and RNA was extracted for PCR verification. The results showed that there was significant overexpression in MC38 cells and MC38 allogeneic grafts, while there was no expression in control normal tissues (such as Figure 5 E) All results indicate that the nascent protein detection method can analyze the full range of variants and shows high performance in identifying long isoforms.

[0166] discuss

[0167] In the new era of cancer immunotherapy, the concept of "neoantigens" has become widely popular. Neoantigens can be defined as immunogenic proteins / epitopes derived from abnormal genetic changes in tumor cells. Since normal cells do not have such abnormal genetic changes, neoantigens can be regarded as "absolute" specific biomarkers for tumor cells. Therefore, in recent years, biotechnology companies and research institutions have been actively developing immunotherapies targeting neoantigens. Today, how to effectively identify neoantigens from clinical samples has become an important issue that needs to be addressed in translational research.

[0168] The "classical" approach to identifying neoantigens is based on searching for nonsynonymous mutations in canonical open reading frames (ORFs). This "classical" approach is easy to understand, but its limitations are also obvious: only tumors with a high mutation load have enough nonsynonymous mutations within ORFs to induce new proteins; in addition, immunogenic neoantigens are only a small subset of new proteins. Since many cancer types have low tumor mutation loads (<5 mutations / Mb), such as prostate cancer, glioblastoma, acute myeloid leukemia, etc. For these cancer types, it is difficult to find enough neoantigens in canonical ORFs for immunotherapy.

[0169] However, due to abnormal genetic changes, most types of cancer cells involve alternative RNA splicing that is not present in normal cells. RNA-seq analysis of The Cancer Genome Atlas data has provided compelling evidence for the prevalence of tumor-specific alternative splicing events, which generate neoantigens that are expected to be more immunogenic than missense mutations. In addition, some studies have shown that alternative RNA splicing provides another source of new proteins that are translated from RNA sequences that are not translated in normal cells. Therefore, identifying new proteins and new antigens based on large-scale mutations (such as alternative RNA splicing) will be an effective new approach to break through the limitations of "classical" methods, thereby greatly expanding the application of immunotherapy targeting neoantigens.

[0170] As a powerful tool for studying the peptidome or proteome, mass spectrometry (MS) can theoretically identify novel proteins derived from alternative RNA splicing. In fact, the performance of MS depends heavily on the characteristics of the reference database. As discussed in the introduction, the canonical human proteome is an unfavorable reference database for identifying novel proteins that are completely outside the canonical human proteome. To resolve the above logical paradox, the ubiquitous use of HiFi and NGS technologies provides a scientifically sound approach to establish sample-specific reference databases containing hypothetical protein sequences derived from RNA translation.

[0171] When analyzing RNA sequencing data, comprehensive detection of all possible variants that may lead to protein changes is a key requirement. In this regard, multiple mechanisms including fusions, single nucleotide polymorphisms (SNPs), and insertions / deletions (INDELs) can lead to the generation of new RNA transcripts and the subsequent induction of new proteins. Given the complexity of the task, a powerful bioinformatics tool is needed to process RNA sequencing data and cover all types of RNA transcript variations.

[0172] All documents mentioned in the present invention are cited as references in this application, just as each document is cited as reference individually. In addition, it should be understood that after reading the above teachings of the present invention, those skilled in the art can make various changes or modifications to the present invention, and these equivalent forms also fall within the scope defined by the claims attached to this application.

Claims

1. A method for detecting a new protein, characterized in that: The method comprises the steps of: (A) providing raw data and constructing a library, wherein the data includes second-generation (NGS) transcriptome sequencing data and third-generation transcriptome sequencing data, wherein the third-generation transcriptome sequencing data is HiFi data; (B) Preprocessing and quality control of raw data, including: The preprocessing comprises the steps of: polishing the third-generation transcriptome sequencing data; aligning the second-generation transcriptome sequencing data to the reference genome to obtain a first BAM file; trimming the third-generation transcriptome sequencing data to obtain a full-length non-chimeric sequence (FLNC), aligning the FLNC to the reference genome to obtain a second BAM file; the two BAM files are used for small-scale variation analysis; clustering the FLNC to obtain a third BAM file, and the BAM file is used for large-scale variation analysis; The quality control includes: calculating the number of zero-mode waveguide holes (ZMWs); calculating the number of reads aligned to a unique position in the reference genome; calculating the number of reads aligned to multiple positions in the reference genome; calculating the number of reads that cannot be aligned to the reference genome; using the BAM file to calculate the number and proportion of ribosome-associated RNA bases, intron bases, intergenic region bases and UTR region bases; calculating the sequencing depth coverage of the 5' end and the 3' end; calculating the proportion of reads aligned to the positive and negative chains of the reference genome; (C) Detection of new variant isoforms (NVI); (D) performing hierarchical filtering and annotation on the detection results to obtain filtered NVI, wherein the filtered NVI includes NVI in nucleotide sequence format and NVI in transcript format; (E) Generation of new proteins using filtered NVI.

2. The method according to claim 1, characterized in that In step (C), the steps include: (C1) Use collapse_isoforms_by_sam.py of cDNA_cupcake to remove redundancy in the HiFi data and detect NVI; use fusion_finder.py of cDNA_cupcake to detect gene fusion; Use PBsv software to detect NVI; remove the results obtained from multiple software to obtain the detection results of HiFi data; (C2) Using Vardict software to detect NVI in the NGS data, and obtaining a VCF file of the detection results of the NGS data.

3. The method according to claim 1, characterized in that In step (D), the layered filtering includes: (D1) filtering isoforms / fusions; (D2) Filtering of SNV / INDELs; Wherein, in step (D1), the steps are included: (d1.1) remove NVI with abundance below the threshold; (d1.2) removal of NVI lacking a 5' exon; (d1.3) noise removal NVI; (d1.4) annotating the remaining NVIs using a database, wherein the database is the chimerDB 4.0 database; In step (D2), comprising: (d2.1) Removing NVIs whose sequencing depth is lower than a threshold, the threshold of the sequencing depth is calculated as follows: Among them, RPK refers to reads per kilobase; TPM refers to transcripts per million; mean_read_length refers to the average length of reads; (d2.2) Remove NVI with noise mark in the FILTER column of the VCF file; (d2.3) Removal of NVIs from SNPs or germline mutations; (d2.4) The remaining NVIs are annotated using databases, including the COSMIC database.

4. The method according to claim 1, characterized in that In step (E), the steps include: (E1) converting the NVI in the nucleotide sequence format into the NVI in the transcript format; (E2) merging the NVIs of the two transcript formats and removing redundancy, and storing the NVIs without redundancy in a fasta file; (E3) converting the redundantly removed NVI into a protein sequence using GeneMarkS-T software; (E4) Using BLAST software, the protein sequence is compared with the UniProt protein database to determine whether the protein sequence is a new protein.

5. A software or system for identifying nascent proteins, characterized in that: The software or system includes: An input unit, wherein the input unit is configured to input sequence data, wherein the input sequence data is NGS transcriptome sequencing data and third-generation transcriptome sequencing data; A detection module, wherein the detection module is configured to perform a predetermined detection on the input sequence data to obtain an analysis result, and the detection module comprises: (Z1) a preprocessing and quality control module, wherein the preprocessing and quality control module is configured to execute step (B) of the method according to claim 1, thereby obtaining a preprocessing and quality control result; (Z2) a detection module, wherein the detection module is configured to execute the method of claim 2 to obtain a detection result; (Z4) a filtering and annotation module, wherein the preprocessing and quality control module is configured to execute the method of claim 3 to obtain filtering and annotation results; An output unit, wherein the output unit is configured to execute the method of claim 4, thereby outputting the detection result of the detection module.

6. A protein generated by the method of any one of claims 1 to 4, or the software or system of claim 5, characterized in that: The protein is encoded by the nucleotide sequence shown in SEQ ID NO: 1-4.

7. A polynucleotide, characterized in that The polynucleotide comprises the nucleotide sequence shown in SEQ ID NO: 1-4.

8. A composition, characterized in that The composition comprises the protein of claim 6 or the polynucleotide of claim 7.

9. Use of the protein according to claim 6, the polynucleotide according to claim 7 or the composition according to claim 8, characterized in that: Used for preparing reagents or drugs for treating tumors.

10. A method according to any one of claims 1 to 4, or a use of the software or system according to claim 5, characterized in that: Used for protein analysis, including protein structure prediction and new antigen immune epitope analysis.

Images