A screening method for atypical new antigens

Abstract

The present invention discloses a method for screening atypical neoantigens, comprising: an RNA-seq screening step: an HLA typing prediction step: a gene transcript expression quantification step: a transcript variable splicing mutation detection step: an RNA editing mutant peptide generation step; a reference mutant peptide data construction step; a neoantigen prediction step: combining the results obtained in the HLA typing prediction step to obtain the binding affinity of candidate neoantigens to MHC class I and the binding affinity to MHC class II; obtaining atypical neoantigens by sorting candidate neoantigens; the screening method, through RNA-seq data of tumor patients, combined with bioinformatics tools, utilizes the characteristics that effective neoantigens have similar sequences, and shares neoantigens between patients with the same disease.

Description

Technical Field

[0001] The present invention relates to a method for screening atypical neoantigens, and belongs to the technical field of biomedicine. Background Art

[0002] In recent years, one of the research focuses in the field of tumor immunotherapy is how to effectively identify and utilize tumor neoantigens, that is, antigens that are only expressed in tumor cells but not in normal tissues, as targets for stimulating specific immune responses. At present, most of the neoantigens screened by the neoantigen recognition process are mainly produced by non-synonymous single nucleotide variations (SNVs) and small insertions and deletions (InDels) in somatic DNA, which are called typical neoantigens. However, more and more studies have shown that in addition to neoantigens produced by DNA mutations, mutations in the transcriptome and proteome are also important sources of tumor neoantigens, and these neoantigens can also efficiently induce T cells to produce strong immune responses. These non-DNA mutations produce neoantigens called atypical neoantigens.

[0003] In this context, researchers have used new technologies and methods to discover a variety of non-DNA variations that can produce new antigens. These non-DNA variations include but are not limited to: 1. Transcript alternative splicing variation: Different splicing methods of genes can produce mRNAs with different sequences, which are then translated into different proteins, which provides an additional pathway for the generation of new antigens; 2. Post-translational modification: After translation, proteins undergo chemical modifications (such as phosphorylation, glycosylation, etc.), which may change their immunogenicity and generate new epitopes; 3. Transposable element movement: Such as the insertion of retroviral elements, which can lead to changes in genome structure and produce unique antigens; 4. Non-coding region variation: Although it does not directly change the coding sequence, it may affect gene expression regulation and indirectly contribute to the emergence of new antigens; these atypical variations may cause the original non-coding sequence to be translated into protein, or affect the translation and modification of protein, thereby generating new immunogenic epitopes. Compared with classical neoantigens, unconventional neoantigens are shared more frequently among patients, which provides hope for the development of universal immunotherapy drugs suitable for a wider population. For example, abnormal splicing factor mutations, such as splicing factor 3b subunit 1 (SF3B1), are common in many cancers and can cause changes in splicing patterns, thereby generating new antigens. In addition, viruses integrate into the host genome and generate shared antigens through atypical pathways.

[0004] However, as a new direction with great potential in tumor immunotherapy, the screening methods of atypical neoantigens still have several unresolved problems. For example, the understanding of the mechanism is not comprehensive: the biological mechanism of the formation of atypical neoantigens, such as alternative splicing, non-coding region variation, post-translational modification and other processes. Although we have some understanding, the specific mechanism and its differences in different tumor types still need to be studied in depth. This ambiguity at the mechanism level limits the design and implementation of targeted treatment strategies. Identification of non-DNA variations: Related studies have shown that atypical variations such as selective splicing variation, transposable element movement, non-coding region variation, and post-translational modification can produce effective neoantigens. Due to insufficient research on the mechanism, there is a serious lack of relevant bioinformatics tools that can accurately and efficiently identify these atypical variations. Optimization of atypical neoantigens: Atypical neoantigens, such as antigens produced by abnormal splicing, are often more immunogenic than antigens formed by non-synonymous mutations. However, the current immunogenicity assessment methods for atypical neoantigens are not perfect, making the accurate selection of atypical neoantigens challenging. Summary of the invention

[0005] In order to overcome the shortcomings of the prior art, the purpose of the present invention is to provide a method for screening atypical neoantigens, which uses RNA-seq data of tumor patients, combined with bioinformatics tools, and utilizes the characteristics that effective neoantigens have similar sequences. There are shared neoantigens among patients with the same disease.

[0006] The purpose of the present invention can be achieved by adopting the following technical scheme: A method for screening atypical neoantigens, comprising:

[0007] RNA-seq screening steps: RNA-seq and proteomic data of tumor tissues of tumor patients were obtained, and RNA-seq data were screened: reads containing sequencing adapters were removed, reads with a base quality score of <15 and a base ratio of >40% and reads containing N bases ≥5 were removed, the total GC base content of the reads was: 45-55%, and the proportion of the base quality value of the reads >30 exceeded 85%;

[0008] HLA typing prediction step: predicting HLA-I class genotyping and HLA-II class genotyping;

[0009] Steps for quantifying gene transcript expression: TPM quantitative formula (Transcripts Per Million) = Quantify transcript expression;

[0010] Transcript variable splicing mutation detection step: determine the exon-to-intron splicing event and exon-exon junctions event;

[0011] RNA editing mutant peptide generation steps: After filtering non-synonymous events through annotation information, RNA editing sites are obtained; then candidate peptides are selected according to the RNA editing sites, which can contain up to ten amino acids before and after the editing sites. Based on the candidate peptides, RNA editing mutant peptides are obtained through sequence similarity;

[0012] Reference mutant peptide data construction step: construct a reference mutant peptide database using the data obtained in the transcript variable splicing mutation detection step, the data obtained in the RNA editing mutant peptide generation step, the human normal protein sequence and the sequence in the cRAP data;

[0013] Neoantigen prediction step: Based on the reference mutant peptide database, the maximum length of the peptide is 30 amino acids and the minimum length of the peptide is 7 amino acids to screen candidate neoantigens; combined with the results obtained in the HLA typing prediction step, the binding affinity of the candidate neoantigen to MHC class I and the binding affinity to MHC class II are obtained; by sorting the candidate neoantigens, atypical neoantigens are obtained.

[0014] Furthermore, in the RNA-seq screening step, the amount of RNA-seq and proteomic data of tumor tissues of tumor patients is >10G.

[0015] Furthermore, in the HLA typing prediction step, the HLA class I genotyping was predicted using the optimal parameter 0.009 of homology judgment by Optitype as a condition.

[0016] Furthermore, in the HLA typing prediction step, class II genotyping is predicted using a prediction method based on fragment assembly.

[0017] Furthermore, in the transcript variable splicing mutation detection step, the conditions for determining the Exitron splicing event are: the minimum number of reads observed for transcription is 3, the alignment score of reads to the reference genome is ≥50, and the minimum threshold of PSO is 0.05.

[0018] Furthermore, in the transcript variable splicing mutation detection step, the conditions for determining exon-exon junctions events are: the number of reads supporting the junction (junction-supportive read counts) ≥ 10, the sequencing depth at the junction (read depth at the junction) ≥ 20 and the read frequency supporting the junction (junction-supportive read frequency) ≥ 0.01.

[0019] Furthermore, in the neoantigen prediction step, NetMHCpan was used to predict the binding affinity of the candidate neoantigen to MHC class I, and then NetMHCIIpan was used to predict the binding affinity of the candidate neoantigen to MHC class II.

[0020] Furthermore, in the neoantigen prediction step, the candidate neoantigens are ranked by scores.

[0021] Score = Rank_bf*0.5+Rank_pr*0.3+Rank_bs*0.2; where bf is the binding affinity, bs is the binding stability, and pr is the percentile rank score.

[0022] Furthermore, the screening method also includes a neoantigen optimization step: the atypical neoantigens obtained in the neoantigen prediction step are compared for sequence similarity with the experimentally verified neoantigens in the dbPepNeo v2.0 database to obtain atypical neoantigens with a similarity ≥ 90%.

[0023] Furthermore, the screening method also includes a new antigen optimization step: after the new antigen prediction step, the candidate new antigens binding to MHC class I molecules are screened with conditional binding strength ≤34nM and tumor gene expression abundance ≥33TPM.

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

[0025] 1. The present invention uses RNA-seq data of tumor patients, combined with bioinformatics tools, and takes advantage of the fact that effective new antigens have similar sequences. There are shared new antigens among patients with the same disease.

[0026] 2. The present invention is conducive to the accurate selection of atypical neoantigens by predicting the binding affinity and stability of mutant peptides to MHC molecules. DETAILED DESCRIPTION

[0027] Below, in conjunction with specific implementation modes, the present invention is further described:

[0028] 1. RNA-seq screening steps: RNA-seq data of tumor tissues are obtained by high-throughput sequencing, and the data volume is required to be greater than 10G; proteomic data of tumor tissues are obtained by LC-MS / MS mass spectrometry; then the RNA-seq data are screened: reads containing sequencing adapters are removed, reads with a base quality score of <15 and a ratio of >40% and reads containing N bases ≥5 are removed, the total GC base content of the reads is: 45-55%, and the proportion of reads with a base quality value of >30 exceeds 85%. The proportion of reads with a base quality value of >30 must exceed 85%; Fastp (v0.23.2) software can be used to achieve the above detection.

[0029] 2. HLA typing prediction steps: For the prediction of HLA class I genotyping, Optitype is a tool that can use RNA-seq data for accurate prediction. When using it, using the optimal parameter of 0.009 for homology judgment, an optimal prediction solution can be obtained. For the prediction of HLA class II genotyping, HLAminer (v1.4) is used to infer HLA-II class genotyping. When using it, choosing a prediction method based on fragment assembly can obtain more accurate prediction results. The important function of class I molecules of the histocompatibility complex (MHC) is to participate in the process of presenting antigens to CD8+ cell lysing T cells. MHC class I molecules have a peptide binding groove with two closed ends, which limits the size of their ligands to about 8-11 amino acids. The antigen binding site is mainly targeted at the antigen peptide backbone with relatively conservative amino acid sequences, which is conducive to the binding of variable amino acid side chains in the free state on the antigen peptide to TCR. MHC class II molecules mainly present processed exogenous antigen fragments to CD4+T cells in the initial stage of the immune response. In contrast to MHC class I, class II molecules have an open binding groove that enables them to bind longer peptides, typically 12-20 amino acid residues in length. The complex open binding groove makes it difficult to assess the binding affinity of antigens to MHC class II. Predicting HLA typing is essential for the subsequent identification of new antigens. Different MHC haplotypes have different binding affinities to peptides, so accurate genotyping of HLA alleles is essential before predicting peptide-MHC binding affinity.

[0030] 3. Quantification steps of gene transcript expression: TPM quantitative formula (Transcripts Per Million) = Quantify transcript expression;

[0031] There are three common methods for quantifying gene expression:

[0032] ①RPKM, which is the number of reads mapped to the gene divided by the number of all reads mapped to the genome (in millions) and the length of RNA (in KB), is a formula for measuring the relative expression of genes and is applicable to single-end sequencing;

[0033] ②FPKM, the number of fragments mapped to the gene divided by the number of all reads mapped to the genome (in millions) and the length of the RNA (in KB), applicable to single-end and double-end sequencing;

[0034] ③TPM, first correct the number of reads of each gene with the length of the gene, and then take the quotient of the corrected number of reads of this gene and the corrected number of reads of all the corrected reads of this sample, which is suitable for single-end and double-end sequencing.

[0035] The results of RPKM and FPKM are affected by gene length, while TPM takes gene length and sequencing depth into consideration to make its results more robust.

[0036] Therefore, we chose to use TPM as the quantification method for transcripts. Kallisto is a transcript expression quantification tool that does not require alignment and uses a pseudo-alignment algorithm. First, each transcript sequence is divided into kmer, and the kmer sequences of all transcripts are used to construct the de Bgujin Graph, referred to as T-DBG. In the T-DBG graph, each node is a kmer, and each path represents a transcript. Due to the redundancy of transcript sequences, each kmer actually corresponds to multiple paths, that is, multiple transcripts; then the sequenced reads are also divided into kmers and mapped to T-DBG. Compared with other quantitative algorithms, pseudo-alignment quantifies transcripts through the kmer features of reads with higher accuracy and faster calculation speed. When using kallisto, the calculation speed can be accelerated by increasing the kmer length, for example, the read length is 150 and the kmer is 31. This setting can have more advantages than other similar tools.

[0037] 4. Transcript variable splicing mutation detection steps: Determine Exitron splicing events and exon-exon junctions events; related studies have shown that variable splicing is not only the core mechanism of protein expression, but also found that the incidence in tumor samples is about 1.2 times that of matched healthy samples. Variable splicing abnormalities are highly correlated with the generation of new antigens. There are many variable splicing events, but the current detection methods and means are not yet mature. Therefore, here we only detect two relatively accurate variable splicing mutations. The advantage of this is to reduce the occurrence of false positive new antigens.

[0038] Exitron splicing is an atypical alternative splicing mechanism that creates a hidden intron, or exitron, inside a protein-coding exon. This term comes from the combination of "exon" and "intron", meaning an intron-like splicing event that occurs inside an exon. Incorrect Exitron splicing may lead to abnormal protein function, is associated with the occurrence of certain cancers, and is an important source of new antigens. ScanExitron (v1.3.1) identifies exitrons based on splicing junctions and gene annotations, and filters out low-confidence candidates through a splicing percentage (percent spliced ​​out, PSO) indicator, which measures the percentage of transcripts that are spliced ​​by a specific exitron. Therefore, when determining Exitron splicing events, the minimum number of reads observed in the transcription is 3, and the alignment score of the reads to the reference genome is not less than 50. The minimum threshold of PSO is 0.05. These parameters can effectively reduce false positive detection results.

[0039] Exon-exon junctions event: This mutation event refers to the mutation that occurs at the junction of exons during gene splicing. Such mutations may affect the normal splicing pattern, thereby causing changes or loss of protein function. These mutation events are common in tumor diseases because they can cause significant changes in protein structure and function, sometimes even completely destroying protein activity, and are also another important source of new antigens. These mutations can be divided into several categories: 1) splice site mutations, which occur at the 5' or 3' splice site; 2) exon skipping, when the mutation affects the recognition of the splicing mechanism for the exon, it may cause the entire exon to be skipped. In this case, two exons that were not directly connected will be directly connected together to form a new exon-exon junction; 3) intron retention: On the contrary, if the mutation prevents the normal splicing of the intron, it may cause the intron to be incorrectly retained in the mature mRNA, which also changes the original exon-exon junctions. SplAdder (v3.0.5) is a tool for detecting and quantifying alternative splicing events in RNA-Seq data. The tool uses given annotations and standardized RNA-Seq read alignments, converts the annotations into variable splicing map representations, extracts additional information from the read data to expand the variable splicing map, extracts alternative variable splicing events from the map based on the alignment data, and quantifies the events. The quantified events can then be used for differential analysis. When using the SplAdder software to detect this event, the variation determination conditions that need to be focused on are: (1) junction-supportive readcounts ≥ 10, (2) read depth at the junction ≥ 20, and (3) junction-supportive read frequency ≥ 0.01 (1%). These conditions can effectively remove false positive data in the test results.

[0040] 5. RNA editing mutant peptide generation steps: After screening non-synonymous events through annotation information, RNA editing sites are obtained; then candidate peptides are selected according to the RNA editing sites, which can contain up to ten amino acids before and after the editing sites. Based on the candidate peptides, RNA editing mutant peptides are obtained through sequence similarity;

[0041] RNA editing refers to the process of chemically modifying RNA sequences by specific enzymes after RNA molecules are transcribed, thereby changing the information encoded by the original gene. RNA editing can occur on a variety of RNA molecules such as mRNA, tRNA, and rRNA. The most common types include adenosine to inosine (A-to-I) editing and cytidine to uridine (C-to-U) editing. Peptides derived from edited RNA transcripts - edited peptides - may be presented on human leukocyte antigens (HLA) and become a source of tumor antigens. The DARNED database contains RNA editing site information for multiple species such as human. In addition to the most common A->I type of RNA editing, it also includes some C->U RNA editing sites. The RADAR version 2 database contains 2,576,459 records, which were downloaded and annotated using ANNOVAR based on the RefSeq annotations hg19_refGeneMrna.fa and hg19_refGene.txt. After screening non-synonymous events through annotation information, 1,369 RNA editing sites were obtained. The R package sapFinder (v1.18.0) was further used to infer the corresponding amino acid sequence. Due to the existence of different protein isomers, these editing sites corresponded to 2516 records, which in turn corresponded to 1387 unique and different candidate peptides (21-mer peptides), which can contain up to ten amino acids before and after the editing site. Based on these candidate peptides, RNA editing mutant peptides were obtained by sequence similarity. The method of using data to form reference RNA editing mutant peptides is more accurate and has lower false positives than the detection method of conventional tools. Due to the lack of reliable model training data, the current detection tools for RNA editing variations are not accurate, and the use of such tools will produce a large number of false positive neoantigens.

[0042] 6. Reference mutant peptide data construction step: Construct a reference mutant peptide database using the data obtained in the transcript variable splicing mutation detection step, the data obtained in the RNA editing mutant peptide generation step, the human normal protein sequence and the sequence in the cRAP data;

[0043] By collecting existing protein and peptide sequence data, a reference amino acid sequence library was constructed for the analysis of subsequent proteome data. An R script was written to process the amino acid sequence, with the goal of making the reference amino acid sequence into a fasta format file. By integrating the results in steps 4 and 5, the human normal protein sequence, and the sequences in the cRAP data, a customized Uniprot+cRAP+mutant peptide database was constructed using the written R script. In addition, cRAP is a database containing protein sequences found as contaminants in proteomics experiments.

[0044] 7. Neoantigen prediction step: Based on the reference mutant peptide database, the maximum length of the peptide is 30 amino acids and the minimum length of the peptide is 7 amino acids, and the candidate neoantigens are screened; combined with the results obtained in the HLA typing prediction step, the binding affinity of the candidate neoantigen to MHC class I and the binding affinity to MHC class II are obtained; by sorting the candidate neoantigens, atypical neoantigens are obtained;

[0045] The data processing process includes two processes: protein identification and quantification. First, protein identification is performed through MS / MS data, and then the corresponding peptide signal intensity or specific quantitative tags are used to quantify the expression level of the protein. However, in practice, protein identification and quantification are usually performed jointly, and there are many methods and tools for analysis. We chose the MS-GF+ (v2020.07.02) tool, which is a protein identification and quantification tool that is both sensitive and versatile. Compared with most other search tools, it can identify more peptides and is compatible with different types of mass spectrometry data, mass spectrometry instruments with different configurations, and different experimental strategies. During the analysis, combined with the customized database obtained in step 6, the maximum length of the peptide is 30 amino acids, and the minimum length of the peptide is 7 amino acids, which is also the minimum effective length of the new antigen. Personalized candidate new antigens are obtained through this step. By using a customized database, protein mass spectrometry data can not only screen candidate new antigens produced by RNA mutations, but also screen out new antigens produced by abnormal protein translation modification.

[0046] The prerequisite for whether a candidate neoantigen is immunogenic is that it can bind to MHC molecules. Therefore, the prediction results of the binding of neoantigens to MHC molecules are very important and can be used as a basis for screening candidate neoantigens. MHC molecules have different molecular structures and different binding modes, and different prediction methods are required. Related studies have shown that NetMHCpan (v4.1) based on artificial neural network algorithms has a good prediction effect for MHC class I molecules, and NetMHCIIpan (v4.0) based on deep learning algorithms has a good prediction effect for MHC class II molecules. Therefore, the obtained candidate neoantigens are combined with the HLA allele typing results in step 2, and NetMHCpan is used to predict the binding affinity of the candidate neoantigen to MHC class I, and then NetMHCIIpan is used to predict the binding affinity of the candidate neoantigen to MHC class II. There are three indicators in the prediction results to measure binding affinity: 1) binding affinity (bf) (bf < 500uM), 2) binding stability (bs) (> 1h), 3) pr (% rank) percentile rank score (HLA-I: pr ≤ 0.5; HLA-II: pr ≤ 2). pr (% rank) contains information about peptide MHC binding events and previous steps in the biological antigen presentation pathway. The three indicators represent a condition for binding, and their importance varies. Based on the results of related studies, the weight coefficients based on the three indicators were determined, and the following formula was obtained by calculation:

[0047] Score=Rank_bf*0.5+Rank_pr*0.3+Rank_bs*0.2

[0048] The score represents the binding score between the candidate neoantigen and the MHC molecule determined by combining the three indicators. The score indicator is more balanced and more robust than using a single indicator. Finally, the candidate neoantigens are ranked according to the comprehensive binding score.

[0049] 8. Neoantigen Optimization Steps:

[0050] The atypical neoantigens obtained in the neoantigen prediction step were compared with the experimentally verified neoantigens in the dbPepNeo v2.0 database for sequence similarity, and atypical neoantigens with similarity ≥ 90% were obtained. For candidate neoantigens that bind to MHC class I molecules, they were screened based on the conditional binding strength ≤ 34nM and tumor gene expression abundance ≥ 33TPM.

[0051] Using the sequence alignment tool BLASTp, the atypical neoantigens in step 7 were compared with the experimentally verified neoantigens in the dbPepNeo v2.0 database for sequence similarity to obtain atypical neoantigens with a sequence similarity of ≥90%. We used a total of 801 high-confidence neoantigens that were directly verified by experiments to induce immune responses.

[0052] Strict threshold filtering of new antigens bound to MHC class I molecules: For new antigens bound to MHC class I molecules, a more stringent binding threshold is used, with binding strength ≤34nM and tumor gene expression abundance ≥33TPM. The purpose of setting this condition is to reduce false positive new antigens bound to MHC class I molecules, because most of the current prediction tools are biased towards MHC I molecules, which will introduce more errors. Compared with the method of selecting new antigens based solely on the binding score to MHC molecules, our optimization strategy can obtain fewer false positive candidate new antigens.

[0053] Embodiment 1:

[0054] Obtain tumor tissue from a lung cancer patient, extract RNA from the tumor tissue to construct a sequencing library, perform RNA-seq sequencing, obtain 10G of RNA data, process the tumor tissue, and obtain high-quality LC-MS / MS protein mass spectrometry data. Data storage location:

[0055] RNAseq data:

[0056] https: / / www.ncbi.nlm.nih.gov / geo / query / acc.cgi?acc=GSM1104129

[0057] Proteomic data:

[0058] https: / / db.systemsbiology.net / sbeams / cgi / PeptideAtlas / PASS_View?identifier=PASS00215

[0059] The data analysis is as follows:

[0060] (1) Quality control: Fastp (v0.23.2) software was used to perform quality control on the patient RNA-seq sequencing data, including removing the sequencing adapters of each read, removing sequences with a base quality score of <15 and a base ratio of >40% and reads with N bases ≥5, the total GC base content of the reads was 45-55%, and the proportion of reads with a base quality value greater than 30 was more than 85%. If any of the above indicators was not met, the analysis was terminated.

[0061] (2) Prediction of patient HLA genotype: Based on the RNA data after quality control in (1), HLA-I alleles were predicted using Optitype (v1.3.5); HLA-II alleles were predicted using HLAminer (v1.4). Optitype (v1.3.5) and HLA-HD (v1.0) used default parameters.

[0062] (3) Use kallisto (v0.50.0) software to perform quantitative analysis on the data processed in (1), using TPM quantification with parameters: t8-b100 and the rest as default values.

[0063] (4) Detect transcript variable splicing mutations based on the results of (1) and (3): 1. ScanExitron (v1.3.1) was used to detect Exitron splicing events. The parameter configuration was: --ao3 --pso0.05 -m50. The other parameters were default. The screening conditions for mutations were: TPM>1, PSI (efficiency of exon splicing into the transcript population)>0.8. 2. SplAdder (v3.0.5) was used to detect exon-exon junctions events. The parameter configuration was: -T exon_skip, -c2. The other parameters were default. The screening conditions for mutations were: (1) junction-supportive read counts ≥ 10, (2) read depth at the junction ≥ 20, and (3) junction-supportive read frequency ≥ 0.01 (1%).

[0064] (5) Based on the hg38 reference genome, the RADAR version 2 database containing 2,576,459 records was annotated using ANNOVAR to screen RNA editing mutation sites of non-synonymous events. The R package sapFinder (v1.18.0) was used to infer the amino acid sequence corresponding to the mutation site and obtain the reference mutant peptide.

[0065] (6) A customized Uniprot+cRAP+mutant peptide amino acid sequence database was constructed using a self-written R script to process the results of 4 and 5.

[0066] (7) Based on the 6 databases, MSGF+ (v2020.07.02) software was used to process the LC-MS / MS mass spectrometry data to obtain candidate neoantigens. Parameter settings: (6) Generated database, PrecursorMassTolerance=20ppm; IsotopeError: 0, 1; NumMods=3; FragmentationMethodID=0; InstrumentID=1; EnzymeID=1; NTT=2; TDA=1; NumThreads=All; MinPepLength=7; MaxPepLength=30; MinCharge=2; MaxCharge=5; NumMatchesPerSpec=1; MaxMissedCleavages=1; StaticMod=229.1629, *, fix, N-term, TMT11plex; StaticMod=229.1629, K, fix,any, TMT11plex; StaticMod=C2H3N1O1, C, fix, any, Carbamidomethyl; DynamicMod=O1, MP, opt, any, Oxidation. Java environment configuration: java -Xmx3500M -jar MSGFPlus.jar -s$MZML -d $FASTA -conf $CONFIG -o $MZID.

[0067] (8) Based on the results in (2) and (7), the NetMHCpan (v4.1) tool was used to predict the binding affinity of candidate neoantigens to MHC class I molecules, and the NetMHCIIpan (v4.0) tool was used to predict the binding affinity of candidate neoantigens to MHC class II. Parameter configuration: NetMHCpan: -BA-l 8,9,10,11, the effective length of the neoantigen is 8-11 amino acids; NetMHCIIpan: -BA-length 13,14,15,16,17, the effective length of the neoantigen is 12-20 amino acids.

[0068] Each candidate neoantigen has three binding indicators: 1) binding affinity (bf) (bf < 500uM), 2) binding stability (bs) (> 1h), 3) pr (%rank) percentile rank score (HLA-I: pr ≤ 0.5; HLA-II: pr ≤ 2). According to the formula: score = Rank_bf*0.5+Rank_pr*0.3+Rank_bs*0.2, the score of each candidate neoantigen is calculated by a self-written data processing script, and then ranked.

[0069] After screening according to the above three indicators, there are 77 candidate new antigens, as shown in Table 1:

[0070] Table 1 Score ranking.

[0071]

[0072] (10) Neoantigen optimization: 1) Using the sequence alignment tool BLASTp, the 9 candidate neoantigens were compared with the high-confidence neoantigens in the dbPepNeov2.0 database for sequence similarity to obtain ≥90% neoantigens; parameter settings: –num_threads 11 –max_target_seqs 6, and the rest were default values. The results after screening are shown in Table 2. databaseID is the ID of the high-confidence peptide in the database, and Pos-percentage is the comparison score.

[0073] Table 2 BLASTp screening results

[0074]

[0075] 2) Strict threshold filtering of new antigens bound to MHC class I molecules: For new antigens bound to MHC class I molecules, a more stringent binding threshold is used, with binding strength ≤34nM and tumor gene expression abundance ≥33TPM. The results after screening are shown in Table 3; 3) The results after combining the first two steps of filtering are shown in Table 4.

[0076] Table 3 More stringent HLA class I neoantigens

[0077]

[0078] Table 4 Neoantigens that meet two conditions

[0079]

[0080] Embodiment 2:

[0081] By analyzing the RNA-seq sequencing data and protein mass spectrometry data of the patient's tumor tissue, we screen for new antigens produced by RNA mutations and abnormal protein post-translational modifications. The conventional new antigen screening program is to screen for new antigens produced by RNA variable splicing by analyzing the patient's blood NDA sequencing data, tumor tissue DNA, and RNA-seq sequencing data. The conventional method and this program are very different in terms of the data used, the means of variation detection, the optimization of new antigens, and the consumption of resources and time, as shown in Table 5:

[0082] Table 5 Differences between Example 1 and conventional method

[0083]

[0084] Embodiment 3:

[0085] The same RNA-seq sequencing data as in Example 1 was used, and the conventional neoantigen screening scheme in Example 2 was adopted. After analysis, 53 candidate neoantigens generated by RNA mutations were obtained. This candidate neoantigen was compared with the top 200 candidate neoantigens ranked by scores in Example 1, and it was found that 10 neoantigens were the same as the neoantigens in Example 1, as shown in Table 6, of which 3 (MYO1B, SYNJ2, CNOT1) met the screening criteria of Example 1, and the rest did not meet the criteria mainly due to poor binding stability.

[0086] Table 6 Same neoantigens

[0087]

[0088] This scheme was further optimized and BLASTp was used to compare the 53 candidate neoantigens with the 801 high-confidence neoantigens in dbPepNeo. The results are shown in Table 7:

[0089] Table 7 Comparison of neoantigens

[0090]

[0091] The above results show that compared with conventional schemes, this technical scheme has obvious advantages in detecting atypical variations: RNA variations and abnormal protein post-translational modifications, and also has obvious superiority in screening atypical new antigens.

[0092] For those skilled in the art, various other corresponding changes and modifications can be made according to the technical solutions and concepts described above, and all of these changes and modifications should fall within the protection scope of the claims of the present invention.

Claims

1. A method for screening atypical neoantigens, characterized in that: include: RNA-seq screening steps: RNA-seq and proteomic data of tumor tissues of tumor patients were obtained, and RNA-seq data were screened: reads containing sequencing adapters were removed, reads with a base quality score of <15 and a base ratio of >40% and reads containing N bases ≥5 were removed, the total GC base content of the reads was: 45-55%, and the proportion of the base quality value of the reads >30 exceeded 85%; HLA typing prediction step: predicting HLA-I class genotyping and HLA-II class genotyping; Steps for quantifying gene transcript expression: TPM quantitative formula = Quantify transcript expression; Transcript variable splicing mutation detection step: determine the splicing event of exon conversion to intron and the exon-exon connection event; RNA editing mutant peptide generation steps: After filtering non-synonymous events through annotation information, RNA editing sites are obtained; then candidate peptides are selected according to the RNA editing sites, which can contain up to ten amino acids before and after the editing sites. Based on the candidate peptides, RNA editing mutant peptides are obtained through sequence similarity; Reference mutant peptide data construction step: construct a reference mutant peptide database using the data obtained in the transcript variable splicing mutation detection step, the data obtained in the RNA editing mutant peptide generation step, the human normal protein sequence and the sequence in the cRAP data; Neoantigen prediction step: Based on the reference mutant peptide database, the maximum length of the peptide is 30 amino acids and the minimum length is 7 amino acids, and the candidate neoantigens are screened; combined with the results obtained in the HLA typing prediction step, the binding affinity of the candidate neoantigen to MHC class I and the binding affinity to MHC class II are obtained; the candidate neoantigens are ranked by scores to obtain atypical neoantigens; Score = rank_bf*0.5+rank_pr*0.3+rank_bs*0.2; where bf is binding affinity, bs is binding stability, and pr is percentile rank score; Neoantigen optimization step: The atypical neoantigens obtained in the neoantigen prediction step were compared with the experimentally verified neoantigens in the dbPepNeo v2.0 database for sequence similarity to obtain atypical neoantigens with a similarity ≥ 90%; for candidate neoantigens binding to MHC class I molecules, screening was performed based on the conditional binding strength ≤ 34nM and tumor gene expression abundance ≥ 33TPM.

2. The method for screening atypical neoantigens according to claim 1, characterized in that: In the RNA-seq screening step, the amount of RNA-seq and proteomic data of tumor tissues of tumor patients is >10G.

3. The method for screening atypical neoantigens according to claim 1, characterized in that: In the HLA typing prediction step, the HLA class I genotyping was predicted using the optimal parameter 0.009 of homology judgment by Optitype as a condition.

4. The method for screening atypical neoantigens according to claim 1, characterized in that: In the HLA typing prediction step, class II genotyping is predicted using a fragment assembly-based prediction method.

5. The method for screening atypical neoantigens according to claim 1, characterized in that: In the transcript variable splicing mutation detection step, the conditions for determining the splicing event of exon to intron conversion are: the minimum number of reads observed in transcription is 3, the alignment score of reads to the reference genome is ≥50, and the minimum threshold of PSO is 0.

05.

6. The method for screening atypical neoantigens according to claim 1, characterized in that: In the transcript alternative splicing mutation detection step, the conditions for determining exon-exon junction events are: the number of reads supporting the junction is ≥ 10, the sequencing depth of the junction is ≥ 20, and the frequency of reads supporting the junction is ≥ 0.

01.

7. The method for screening atypical neoantigens according to claim 1, characterized in that: In the neoantigen prediction step, NetMHCpan was used to predict the binding affinity of candidate neoantigens to MHC class I, and then NetMHCIIpan was used to predict the binding affinity of candidate neoantigens to MHC class II.

Images