A method for predicting non-coding region-derived neoantigens by integrating multi-omics data

Abstract

The application discloses a method for predicting new antigens from non-coding regions by integrating multi-omics data. The method uses whole genome data and transcriptome data of tumor and normal samples, and combines small peptide databases from non-coding regions to predict new antigens from non-coding regions. The method provides reliable clues for functional research for experimenters, and further improves the research and development efficiency of new antigen vaccines, and has important significance for improving the clinical immunotherapy effect.

Description

Technical Field

[0001] This invention relates to the field of tumor immunotherapy, and more specifically to a method for predicting neoantigens for immunotherapy originating from non-coding regions by integrating whole-genome data of DNA and transcriptome data of RNA. Background Technology

[0002] The human immune system plays a crucial role in the healing of tumors. In recent years, tumor immunotherapy technologies developed around the immune system have repeatedly achieved breakthroughs in efficacy. Mechanistically, immunotherapy technologies modify T cells to activate the immune system or inhibit its inhibitory pathways, thereby enhancing the immune system's ability to recognize and kill tumor cells.

[0003] In recent years, although researchers from various countries have been actively engaged in the research of tumor vaccines, they have focused their attention on the coding regions of tumor cell genomes in order to find tumor-specific antigens, but have not yet made any breakthrough progress.

[0004] Researchers at the Institute for Immunology and Cancer (IRIC) at the University of Montreal published an article in Science Translational Medicine entitled “Noncoding regions are the main source of targetable tumor-specific antigens.” In this study, researchers used a proteogenomic approach to identify that approximately 90% of tumor-specific antigens in two mouse tumor lines and seven primary human samples originated from noncoding DNA.

[0005] Therefore, effective screening of tumor-specific neoantigens originating from non-coding regions will be of great significance for improving the efficiency of antigen vaccine development and the efficacy of immunotherapy. Summary of the Invention

[0006] To address the aforementioned issues, this invention comprehensively considers the possibility that peptides generated by mutations in non-coding regions of tumor tissue can become neoantigens, and develops a bioinformatics method for obtaining tumor-specific neoantigens.

[0007] This invention provides a method for predicting neoantigens originating from non-coding regions by integrating multi-omics data, comprising the following steps:

[0008] (1) Collection, analysis, and result organization of small peptide data, and construction of a small peptide database;

[0009] (2) Comparison of WGS and RNA-seq data between tumor and normal samples;

[0010] (3) Detection, annotation, and eigenvalue calculation of tumor somatic mutations;

[0011] (4) Obtaining the modified sequences of non-coding small peptides from tumors;

[0012] (5) Screening of variant small peptides expressed in tumors;

[0013] (6) Human leukocyte antigen typing in specific cases;

[0014] (7) Assessment of the affinity of tumor-associated non-coding small peptides for HLA, and selection of potential neotumor antigens based on the results.

[0015] In a preferred embodiment of the invention, the method further includes deleting tumor neoantigens derived from non-coding regions from small peptides, which includes one or more of the following steps:

[0016] A. Use cutadapt, STAR, and Ribo-TISH to preprocess, align, and analyze ribosome profiling data to obtain translated open reading frames, and filter small peptides from non-coding sources based on sequence length, genomic location of non-coding regions, and gene annotation;

[0017] B. Use fastqc and trimmomatic to perform quality control on WGS and RNA-seq data of tumor-normal samples to remove sequencing adapters and low-quality sequences; use bwa to process WGS data and STAR to process RNA-seq data to align to the reference genome.

[0018] C. Based on the small peptide annotation information obtained in step A, use Varscan to analyze the WGS and RNA-seq data of tumor-normal samples aligned to the reference genome, and detect tumor somatic mutations and allele mutation frequencies (VAF); use VEP to annotate tumor somatic mutations; use Kallisto, StringTie, etc. to calculate gene expression values ​​(TPM) in tumor-normal samples respectively.

[0019] D. Based on the genomic location of the non-coding small peptide obtained in step A, filter the tumor somatic mutations obtained in step C to obtain tumor somatic mutations that overlap with the non-coding small peptide; according to the reading frame of the small peptide, filter non-synonymous mutations according to the codon table, and obtain the amino acid sequence of the mutated small peptide.

[0020] E. Based on the expression data obtained in step C, filter the non-coding variant peptides obtained in step D to obtain variant peptides expressed in tumors (TPM>0);

[0021] F. For the sequences aligned to the reference genome in step B, use Optitype, seq2HLA, etc. to perform HLA typing;

[0022] G. Using ACME, NetCTL, NetMHCpan 4.0, Neopesee, etc., predict the affinity of the variant peptides in the tumor obtained in step E for each HLA type, and screen variant peptides with high affinity for HLA as potential non-coding region tumor neoantigens.

[0023] In a preferred embodiment of the present invention, the method further includes the step of finding small peptides encoded by non-coding regions, which comprises the following steps:

[0024] a. Collect ribosome profiling sequencing data from the GEO and ENA databases, call cutadapt to remove sequencing adapters, and retain sequences with a remaining length of 25-35 base pairs;

[0025] b calls STAR in EndToEnd mode to align the sequence to the human hg38 reference genome, allowing a maximum of 2 mismatches;

[0026] c calls the Ribo-TISH quality module to analyze the quality of ribosome profiling data and peptidyl site offset;

[0027] d calls Ribo-TISH's predict module to predict the open reading frame (ORF) of the translation. For different start codons of the same ORF, only the best test result (--framebest) is selected.

[0028] ORFs with a length of less than 100 amino acid residues were retained as small peptides. Small peptides from untranslated regions (UTRs) were obtained based on genome annotation. Genomic locations of noncoding RNAs (ncRNAs) were obtained from the NONCODE database to filter and obtain small peptides from noncoding regions. Detailed Implementation

[0029] The present invention will be further described in detail below through specific embodiments. Through these descriptions, the features and advantages of the present invention will become clearer and more apparent. The specific embodiments described herein are for illustration and explanation only and are not intended to limit the present invention.

[0030] Example 1: Integrating multi-omics data to predict neoantigens originating from non-coding regions

[0031] 1. Collection, analysis, and result organization of small peptide data; construction of a small peptide database.

[0032] Step 1.1. From the GEO database ( https: / / www.ncbi.nlm.nih.gov / gds / ) and ENA database ( https: / / www.ebi.ac.uk / ena / browser / home Collect ribosome profiling sequencing data, call cutadapt to remove sequencing adapters, and retain the remaining sequence length of 25-35 base pairs.

[0033] Step 1.2. Call STAR in EndToEnd mode to align the sequence to the human hg38 reference genome, allowing up to 2 mismatches.

[0034] Step 1.3. Call the Ribo-TISH quality module to analyze the quality of ribosome profiling data and peptidyl site offset (P-site offset).

[0035] Step 1.4. Call the predict module of Ribo-TISH to predict the open reading frame (ORF) of the translation. For different start codons of the same ORF, only the best test result (--framebest) is selected.

[0036] Step 1.5. Retain ORFs with a length of less than 100 amino acid residues as small peptides. Obtain small peptides from untranslated regions (UTRs) based on genome annotation. Obtain the genomic location of noncoding RNA (ncRNA) from the NONCODE database (http: / / www.noncode.org / ) to filter and obtain small peptides from noncoding regions.

[0037] Step 1.6. Conduct literature mining to obtain small peptides recorded in published studies, and construct a small peptide database by combining the results of ribosome profiling analysis.

[0038] Table 1. Partial data from the small peptide database

[0039]

[0040]

[0041]

[0042]

[0043]

[0044]

[0045] 2. Alignment of WGS and RNA-seq data between tumor and normal samples

[0046] Step 2.1. Use FastQC and TrimMomatic to perform quality control on the WGS and RNA-seq data of the samples to remove sequencing adapters and low-quality base sequences.

[0047] Step 2.2. Use bwa to align the WGS data to the human hg38 reference genome; use STAR to align the RNA-seq data to the human hg38 reference genome.

[0048] 3. Detection, annotation, and eigenvalue calculation of tumor somatic mutations

[0049] Step 3.1. Add the small peptide annotation information obtained in Step 1 to the human reference genome annotation file as the annotation file required for subsequent analysis.

[0050] Step 3.2. Call Varscan, using WGS and RNA-seq data of tumor-normal samples aligned to the reference genome as input, to detect tumor somatic mutations and allele mutation frequencies (VAF).

[0051] Step 3.3. Call VEP to annotate the detected tumor somatic mutations.

[0052] Step 3.4. Call Kallisto or StringTie to calculate the gene expression value TPM in the tumor-normal samples respectively.

[0053] Table 2. Some somatic mutations detected in colon adenocarcinoma samples from TCGA.

[0054]

[0055]

[0056]

[0057]

[0058]

[0059]

[0060]

[0061] Table 3. Gene expression levels in TCGA colon cancer tumor-normal samples

[0062]

[0063]

[0064]

[0065]

[0066] 4. Obtaining the post-mutated sequence of non-coding small peptides from tumors

[0067] Step 4.1. Based on the genomic location of the non-coding small peptide obtained in Step 1, filter the tumor somatic mutations obtained in Step 3.2 to obtain the tumor somatic mutations that overlap with the non-coding small peptide.

[0068] Step 4.2. For the mutations obtained in Step 4.1, filter according to the reading frame of the small peptide using the codon table to obtain non-synonymous mutations, and obtain the amino acid sequence of the mutated small peptide, for example...

[0069] normal:MKKLPSPYITNLCCAPTRTCFKFPWQFTTPILYQARYQHQDDLW

[0070] mutant:MKKLPIPYITNLCCAPTRTCFKFPWQFWTPILYQARYQHQDDLW

[0071] 5. Screening of variant small peptides expressed in tumors

[0072] Step 5.1. Based on the expression data obtained in Step 3.4, filter the non-coding variant peptides obtained in Step 4 to obtain variant peptides expressed in tumors (TPM>0).

[0073] 6. Human leukocyte antigen (HLA) typing of the case

[0074] Step 6.1. Based on the IMGT / HLA database (https: / / www.ebi.ac.uk / ipd / imgt / hla / ), perform HLA typing on the sequencing data aligned to the reference genome in Step 2 (including sequences not matched to the reference genome) using Optitype, seq2HLA, etc.

[0075] Table 4. Individual HLA typing results

[0076] HLA-A 24:02:01:03 24:02:43 HLA-B 46:01:10 51:01:01:02 HLA-C 01:03 15:02:01:01 HLA-E 01:03:02:02 01:03:01:01 HLA-F 01:01:02:02 01:01:02:02 HLA-G 01:04:01 01:04:01 MICA 010:02 009:01 MICB 002:01:01 005:02:04 HLA-DMA 01:01:01:01 01:01:01:01 HLA-DMB 01:01:01:01 01:01:01:01 HLA-DOA 01:01:03 01:01:01 HLA-DOB 01:01:01:01 01:01:01:01 HLA-DPA1 01:03:01:02 02:02:02 HLA-DPB1 04:02:01:02 04:02:01:02 HLA-DQA1 03:01:01 05:05:01:02 HLA-DQB1 03:01:01:01 03:03:02:04 HLA-DRA 01:01:01:03 01:01:01:03 HLA-DRB1 09:01:02 11:01:01 HLA-DRB5 . . TAP1 02:01:01 02:01:02 TAP2 01:01:03:03 01:04

[0077] 7. Assessment of HLA affinity for tumor-associated non-coding small peptides and selection of potential neoantigens.

[0078] Step 7.1. Based on known research (Maria, Bonsack, Stephanie, et al. Performance Evaluation of MHC Class-I Binding Prediction Tools Based on an Experimentally Validated MHC-Peptide Binding Data Set. [J]. Cancer Immunology Research, 2019), different affinity prediction algorithms such as ACME, NetCTL, NetMHCpan 4.0, and Neopesee are selected for different HLA types and peptide lengths.

[0079] Step 7.2. Select different affinity thresholds for different HLA types and HLA typing algorithms.

[0080] Step 7.3. For the tumor-associated non-coding small peptides obtained in Step 5, predict their affinity for each HLA type before and after mutation based on the algorithm and parameters determined in Steps 7.1 and 7.2.

[0081] Table 5. Partial results of HLA affinity prediction between variant small peptides in tumors (taking the variant small peptides in step 4 as an example) and HLA.

[0082]

[0083]

[0084]

[0085]

[0086]

[0087]

[0088]

[0089] Explanation of each column in the table:

[0090] Peptide: The peptide fragment sequence of a potential ligand.

[0091] Core: The core sequence consisting of 9-11 amino acids that are in direct contact with the MHC.

[0092] Of: The start position of the core sequence in the peptide fragment (if >0, this method can predict N-terminal overhangs).

[0093] Gp: The location to delete (if any).

[0094] Gl: The length to be deleted (if any).

[0095] Icore: The core sequence of interactions.

[0096] Score_EL: Raw predicted score.

[0097] %Rank_EL: The rank of the predicted binding fraction compared to a random set of natural peptides. This metric is not affected by the inherent bias of some molecules towards higher or lower average predicted affinities.

[0098] BindLevel: (SB: Strong affinity, WB: Weak affinity). If %Rank is below the specified threshold for strong affinity (default is 0.5%), the peptide will be identified as having strong affinity. If %Rank is above the threshold for strong affinity but below the threshold for weak affinity (default is 2%), the peptide will be identified as having weak affinity.

[0099] Step 7.4. Based on the affinity prediction results of Step 7.3, filter out the tumor-associated non-coding small peptides that are most likely to be recognized and bound by HLA, which are potential non-coding tumor neoantigens.

[0100] Taking the mutant peptide in step 4 as an example, the core fragment TRFKFPWQF of this peptide has a strong affinity for HLA-B*27:05. If the HLA typing of the case includes HLA-B*27:05, then this peptide can serve as a potential neotumor antigen for the case.

[0101] While many embodiments of the invention have been disclosed above, including those currently preferred, many other embodiments and variations are possible within the scope of the disclosure and in the appended claims. Therefore, the details of the preferred embodiments and examples provided should not be construed as limiting. It should be understood that the terminology used herein is merely descriptive and not limiting, and various variations and many equivalents may be obtained without departing from the spirit or scope of the claimed invention.

Claims

1. A method for predicting neoantigens originating from non-coding regions by integrating multi-omics data, comprising the following steps: (1) Collection, analysis, and result organization of small peptide data, and construction of a small peptide database; (2) Comparison of WGS and RNA-seq data between tumor and normal samples; (3) Detection, annotation, and eigenvalue calculation of tumor somatic mutations; (4) Obtaining the modified sequences of non-coding small peptides from tumors; (5) Screening of variant small peptides expressed in tumors; (6) Human leukocyte antigen typing for cases; (7) Assessment of HLA affinity for tumor-associated non-coding variant peptides. Step (1) also includes screening for tumor neoantigens derived from non-coding regions in small peptides, which includes the following steps: A. Use cutadapt, STAR, and Ribo-TISH to preprocess, align, and analyze ribosome profiling data to obtain translated open reading frames, and filter small peptides from non-coding sources based on sequence length, genomic location of non-coding regions, and gene annotation.

2. The method according to claim 1, characterized in that, It also includes the following steps: B. Use fastqc and trimmomatic to perform quality control on WGS and RNA-seq data of tumor-normal samples to remove sequencing adapters and low-quality sequences; use bwa to process WGS data and STAR to process RNA-seq data to align to the reference genome.

3. The method according to claim 1 or 2, characterized in that, It also includes the following steps: C. Based on the small peptide annotation information obtained in step A, use Varscan to analyze the WGS and RNA-seq data of tumor-normal samples aligned to the reference genome, and detect tumor somatic mutations and allele mutation frequencies (VAF); use VEP to annotate tumor somatic mutations; use Kallisto and StringTie to calculate gene expression values ​​(TPM) in tumor-normal samples respectively.

4. The method according to claim 3, characterized in that, It also includes the following steps: D. Based on the genomic location of the non-coding small peptide obtained in step A, filter the tumor somatic mutations obtained in step C to obtain tumor somatic mutations that overlap with the non-coding small peptide; according to the reading frame of the small peptide, filter non-synonymous mutations according to the codon table, and obtain the amino acid sequence of the mutated small peptide.

5. The method according to claim 4, characterized in that, It also includes the following steps: E. Based on the expression data obtained in step C, filter the non-coding variant peptides obtained in step D to obtain variant peptides expressed in tumors (TPM>0).

6. The method according to claim 5, characterized in that, It also includes the following steps: F. For the sequences aligned to the reference genome in step B, use Optitype and seq2HLA to perform HLA typing.

7. The method according to claim 6, characterized in that, It also includes the following steps: G. Using ACME, NetCTL, NetMHCpan 4.0, and Neopesee, predict the affinity of the variant peptides in the tumor obtained in step E for each HLA type, and screen variant peptides with high affinity for HLA as potential non-coding region tumor neoantigens.

8. The method according to claim 7, characterized in that, It also includes the search for small peptides encoded by non-coding regions, which involves the following steps: a. Collect ribosomal imprint sequencing data from the GEO and ENA databases, call cutadapt to remove sequencing adapters, and retain sequences with a remaining length of 25-35 base pairs; b calls STAR in EndToEnd mode to align the sequence to the human hg38 reference genome, allowing a maximum of 2 mismatches; c calls the Ribo-TISH quality module to analyze the quality of ribosomal blot data and the offset of peptidyl sites; d calls Ribo-TISH's predict module to predict the open reading frame of the translation, and selects only the best test result for different start codons of the same ORF.

9. The method according to claim 8, characterized in that, ORFs with a length of less than 100 amino acid residues were retained as small peptides. Small peptides from untranslated regions of genes were obtained based on genome annotation, and the genomic locations of non-coding RNAs were obtained from the NONCODE database to filter and obtain small peptides from non-coding regions.

Images