Methods of predicting tumor neoantigens of endogenous retroviral origin

Abstract

This invention provides a method for mining tumor antigens derived from endogenous retroviruses. The method involves bulk RNA-seq sequencing of tumor tissue and adjacent normal tissue samples using a splicing strategy to obtain DNA sequences of human endogenous retrovirus (HERV)-derived transcripts (hervRNA). Subsequently, a specific screening strategy is employed to obtain candidate open reading frames (ORF) sequences on tumor-specific hervRNAs. Finally, homology analysis, prediction of binding energy to human leukocyte antigen class I (HLA-I) molecules, and multi-site analysis are used to predict and mine neoantigens derived from HERV in tumors.

Description

Technical Field

[0001] This application relates to the field of tumor immunotherapy technology, and specifically to a method for predicting tumor neoantigens originating from endogenous retroviruses. Background Technology

[0002] Tumor neoantigens are polypeptides expressed exclusively in tumor cells, presented to the cell surface by major histocompatibility complex (MHC), and specifically binding to the T cell receptor (TCR) on the surface of T cells, thereby activating T cells to generate an immune response. Given that tumor neoantigens can activate cytotoxic T cells to kill tumor cells, therapies developed based on tumor neoantigens, such as tumor vaccines and adoptive T cell transplantation, have achieved promising clinical results. In conclusion, the discovery of tumor neoantigens is a fundamental and crucial step in tumor immunotherapy.

[0003] However, directly discovering tumor neoantigens through experimental methods is currently a labor-intensive, time-consuming, and costly process. Computational prediction of tumor neoantigens can narrow down the range of candidate antigen peptides, thereby reducing the workload of experimental verification and accelerating the discovery of tumor neoantigens. However, current patents related to computational prediction mainly focus on tumor neoantigens caused by single nucleotide polymorphisms (see Chinese patent application CN116580771A, which discloses a "method and apparatus for predicting tumor neoantigens"). There are no patent reports on the prediction of tumor neoantigens from other sources, such as human endogenous retroviruses (HERV).

[0004] Human endogenous retroviruses (HERVs) refer to DNA sequences in the human genome derived from exogenous retroviruses. During human evolution, exogenous retroviruses infected the germ cells of human ancestors and incorporated their genetic sequences into the human genome through reverse transcription, thus becoming part of the human genome. In normal tissues, HERVs are in a "silent" state and are not expressed. However, under certain conditions, such as in tumor cells, due to changes in epigenetic status, HERVs begin to express and encode certain proteins. Since these proteins are not expressed in normal tissues, they have the potential to act as tumor neoantigens, activating T cells to generate an immune response and thus exerting anti-cancer functions. Compared with tumor neoantigens derived from single nucleotide polymorphisms (SNPs), HERV-derived tumor neoantigens have greater universality and can exert therapeutic effects in a wider range of patient populations. However, there are currently no patents related to the discovery of HERV-derived tumor neoantigens. Only some research articles have reported methods based on computational prediction to discover tumor neoantigens in tumors. However, these prediction methods have problems such as not taking into account alternative splicing during transcription and incomplete reference HERV sequences, resulting in the need to improve the efficiency of related antigen discovery and prediction. Summary of the Invention

[0005] Therefore, the purpose of this invention is to provide a system and method for predicting tumor neoantigens of endogenous retrovirus origin.

[0006] This invention processes bulk RNA-seq sequencing data using a splicing strategy to obtain DNA sequences of human endogenous retrovirus (HERV)-derived transcripts (hervRNA); then, a specific screening strategy is employed to obtain candidate open reading frames (ORF) sequences on tumor-specific hervRNAs; finally, homology analysis, prediction of binding energy to human leukocyte antigen class I (HLA-I) molecules, and multi-site analysis are used to predict and mine neoantigens of HERV-derived tumors in tumors.

[0007] The objective of this invention is achieved through the following technical solution:

[0008] <First Aspect>

[0009] A method for mining tumor antigens derived from endogenous retroviruses, comprising the following steps:

[0010] S0, bulk RNA-seq sequencing of tumor tissue and adjacent normal tissue samples;

[0011] S1. Perform quality control on the bulk RNA-seq sequencing data from step S0 to obtain filtered bulk RNA-seq data;

[0012] S2. Based on the splicing strategy, process the bulk RNA-seq data filtered in S1 to obtain the DNA sequence of hervRNA, its hervRNA expression matrix in each sample, and the gene expression matrix of each sample, and predict the human leukocyte antigen (HLA) type of the sample.

[0013] S3. Based on the hervRNA expression matrix of each sample in step S2, tumor-specific hervRNAs are screened out according to the specificity of hervRNA expression in tumor tissues.

[0014] S4. Based on the gene expression matrix of each sample in step S2, calculate the exhaustion score and cytotoxicity score of each sample. After excluding samples with high exhaustion scores, perform correlation analysis on the expression level of the tumor-specific hervRNA and the cytotoxicity score based on the remaining samples to obtain candidate tumor-specific hervRNAs related to cytotoxicity.

[0015] S5. Based on the candidate tumor-specific hervRNA and the corresponding hervRNA DNA sequence described in S4, determine the possible ORF (open reading frame) sequence, and obtain the candidate ORF sequence of the candidate tumor-specific hervRNA source with high confidence through homology analysis.

[0016] S6. Based on the human leukocyte antigen (HLA) type of the sample described in S2 and the candidate ORF sequence described in step S5, use NetMHCpan to predict the binding energy of 8-12 amino acids in the ORF with HLA class I molecules, and screen out candidate antigen peptides with binding energy less than 500 nM.

[0017] S7. Perform multi-site analysis on the candidate antigen peptides to obtain the predicted HERV-derived tumor antigen peptides.

[0018] In step S1, quality control refers to using the software FastQC and MultiQC to check the FastQ files of the bulk RNA-seq sequencing data of each sample in step S1 for low-quality sequencing fragments (such as adapter sequences, contamination from low-quality sequencing sequences, and samples with abnormal GC content); and using the software FastP to remove low-quality fragments, as well as removing sample data of poor quality. The samples referred to are valid samples that meet sequencing standards, including tumor samples and normal tissue samples.

[0019] In S2, the splicing strategy involves using bioinformatics tools to assemble sequencing data to obtain the complete DNA sequence of hervRNA.

[0020] S2 includes the following steps:

[0021] S21. Using the 2-pass mode of the STAR software, align the reads in the fastq file after quality control in step S1 to the human reference genome GRCh38.p13 sequence to obtain a BAM file sorted according to the alignment position to the reference genome GRCh38.p13. Then, use the sambamba software to extract the sequencing fragments in the BAM file that are aligned to the HERV region annotated in the HERVd database (human endogenous retrovirus database) as HERV regions.

[0022] S22. Using the Trinity software, a genome-guided de novo transcription splicing mode was used to splice the DNA sequences of hervRNAs from each sample based on the fragments extracted in step S21. The align_and_estimate_abundance.pl module of the Trinity software was used for quantification. Then, the spliced ​​hervRNA sequences were compared with the reference genome GRCh38.p13 sequence using the GMAP software. HervRNAs with a count value greater than 5 and a similarity greater than 96% were retained.

[0023] S23. Using the software TACO, hervRNAs from each sample are integrated based on their presence and expression in each sample. For the integrated hervRNAs, those longer than 200 bp are retained. Then, the software GMAP is used to re-align the HERVs of the filtered hervRNAs to the human GRCh38.p13 reference genome (these hervRNAs also undergo screening in S22). Only HERVs that are individually aligned to the HERV region are retained as candidate HERVs. (The hervRNAs obtained in S22 whose DNA sequences encoding these hervRNAs are individually aligned to the HERVd database and annotated as HERV regions are retained.) This is to obtain more accurate hervRNA sequence information.

[0024] S24. Based on the hervRNA DNA sequence obtained in S23, the occurrence of the DNA sequence in different samples was analyzed. The hervRNA was then fused using the software TACO to obtain the hervRNA gene annotation file (GTF file).

[0025] S25. Merge the hervRNA gene annotation file (GTF file) obtained in step S24 with the human reference genome GTF annotation file, and then use the software RSEM for quantification to obtain the hervRNA expression matrix and gene expression matrix of each sample.

[0026] S26. Using the fastq file after quality control in step S1 as the input file, predict the HLA-I molecular type of each sample using seq2hla.

[0027] include:

[0028] HLA-A: A gene locus in the human leukocyte antigen class I molecule, located within the major histocompatibility complex (MHC) gene group.

[0029] HLA-B: A gene locus in the human leukocyte antigen class I molecule, located within the major histocompatibility complex (MHC) gene group.

[0030] HLA-C: The gene is located on human chromosome 6.

[0031] Step S3 includes the following steps:

[0032] S31. Based on the hervRNA TPM expression (Transcripts Per Million) matrix of each cohort in step S2, retain hervRNAs with TPM greater than 1 in at least 10% of tumor samples.

[0033] S32. Convert TPM values ​​less than 1 in the hervRNA TPM expression matrix to 1, and then take the logarithm of the entire matrix to the base 2 to obtain the logarithmic hervRNA TPM expression matrix.

[0034] S33, Based on Formula Specificity is calculated, where specificity refers to the TAU value τ of a hervRNA. HervRNAs with a mean expression level in tumor samples greater than that in normal tissue samples and a TAU value greater than 0.8 are considered tumor-specific hervRNAs.

[0035] xi represents the mean expression level of the hervRNA in tissue i; the mean expression level is derived from the logarithmic hervRNA TPM (Transcripts Per Million) matrix.

[0036] n represents the number of organization types;

[0037] Representative: Values ​​calculated based on the expression status of hervRNA in all tissues (expression status of hervRNA in tumor tissues and normal tissues).

[0038] The specific calculation process is as follows: First, the TPM values ​​less than 1 in the TPM matrix of the hervRNA are converted to 1. Then, the logarithm of the converted TPM matrix is ​​taken with the base 2 to obtain the log-converted TPM matrix. Next, the mean value of the log-converted TPM of the hervRNA in tumor tissue and adjacent normal tissue is calculated. Finally, the TAU value of the hervRNA is calculated. Only hervRNAs with TPM greater than a set threshold and TAU value greater than a set threshold in tumor tissue are retained as tumor-specific hervRNAs.

[0039] In S3, screening for tumor-specific HERVs involves comparing the expression levels of HERVs in tumor and normal tissues and selecting those HERVs that are specifically highly expressed in tumor tissues.

[0040] In S4, the exhaustion score is calculated as follows: First, the geometric mean of TPM for the CD8A and CD8B genes in each tumor sample is calculated to characterize the CD8+ T cell level of the sample. Then, the geometric mean of TPM for the five genes PDCD1, CTLA4, LAG3, HAVCR2, and KLRG1 in each sample is calculated. Finally, the geometric mean of TPM for the CD8A and CD8B genes is used as the independent variable, and the geometric mean of TPM for the five genes PDCD1, CTLA4, LAG3, HAVCR2, and KLRG1 is used as the dependent variable. The residuals of each tumor sample are calculated through linear regression, and the residuals are used as the exhaustion score of the sample to characterize the exhaustion degree of the tumor sample. It is worth noting that the exhaustion score is calculated by performing linear regression independently for each cohort.

[0041] The cytotoxicity score is calculated by calculating the geometric mean of the TPM of the two genes PRF1 (NCBI: 5551) and GZMB (NCBI: 3002) in each tumor sample in the hervRNA TPM expression matrix, which is used as the cytotoxicity score of the sample.

[0042] Tumor samples with high exhaustion scores are those with exhaustion scores greater than 0; the association analysis refers to the calculation of Spearman correlation coefficients. Specifically, after removing tumor samples with exhaustion scores greater than 0, each cohort independently calculates the Spearman correlation coefficient between the cytotoxicity score of each sample in its cohort and hervRNA TPM, and retains hervRNAs with Spearman correlation coefficients greater than 0.

[0043] In S5: Determining possible ORF sequences refers to first determining the DNA sequence of hervRNA based on the gene annotation information of hervRNA using the software gffread, and then using the software ORFfinder to determine all possible ORF sequences on hervRNA.

[0044] The homology analysis refers to comparing the ORF sequence with HERV-related protein sequences in the gEVE database and NCBI database using BLASTp, and retaining the ORF sequences with a similarity greater than a set threshold.

[0045] The homology analysis also includes BLASTp comparison of the ORF sequence with the ORF of hervRNA expressed in normal tissue, and annotation of the ORF of tumor-specific hervRNA with 100% similarity and all alignments.

[0046] It is worth noting that if there is data in multiple independent queues, the ORF sequence that appears in at least two queues is preserved based on sequence consistency.

[0047] In S6, the screening of candidate antigenic peptides also includes considering their accessibility in tumor tissues.

[0048] In S7, multi-site analysis includes analyzing the differences in binding ability of candidate antigen peptides across different HLA alleles in order to select tumor antigen peptides with broad coverage.

[0049] Furthermore, in S7, the multi-site analysis refers to first using the software BLAT to align the candidate antigen peptide to all loci on the human reference genome (GRCh38.p13) that may encode the antigen peptide, then generating a saf file and using the software featureCount to calculate the expression values ​​of these locus regions, and finally screening based on the expression in normal tissues.

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

[0051] (1) Through a specific analysis strategy (through splicing strategy and specific analysis, relevant information of hervRNA in tumor tissue was obtained, then through homology analysis and multi-cohort analysis, the ORF sequence of highly reliable hervRNA was obtained, then through binding energy prediction, candidate peptides were obtained, and finally through multi-site analysis, tumor neoantigens derived from HERV in tumor tissue were obtained), the prediction of tumor neoantigens derived from HERV was achieved, which is different from the current mainstream mining of tumor neoantigens derived from gene mutations.

[0052] (2) By adopting a splicing strategy-based approach, the prediction of HERV-derived antigens, alternative splicing, and specific loci of hervRNA can be obtained;

[0053] (3) By adopting the TAU value algorithm, association analysis, multi-cohort joint analysis, homology analysis and multi-site analysis, the problem of a large number of HERV-derived antigens obtained by prediction was solved, and the reliability of the predicted antigens was improved. Attached Figure Description

[0054] Other features, objects, and advantages of the present invention will become more apparent from the following detailed description of non-limiting embodiments with reference to the accompanying drawings:

[0055] Figure 1 A graph showing the number of hervRNAs screened from five independent cohorts;

[0056] Figure 2 A distribution map of cytotoxicity scores in different tissues across different cohorts;

[0057] Figure 3 Distribution of the number of hervRNAs after filtering and screening for each cohort;

[0058] Figure 4 A graph showing the number of ORFs determined for the five queues;

[0059] Figure 5 An ORF relationship diagram between different queues;

[0060] Figure 6 This is a histogram showing the number of alignment sites for the peptide.

[0061] Figure 7 This is a distribution map showing the number of candidate antigenic peptides expressed and not expressed in normal tissues due to their multi-site coding.

[0062] Figure 8 The distribution map of the 10 HLA molecules with the most strongly binding antigenic peptides among different HLA molecules;

[0063] Figure 9 This is a flowchart of the method for predicting tumor neoantigens of endogenous retrovirus origin as described in Example 1. Detailed Implementation

[0064] The present invention will be described in detail below with reference to embodiments. These embodiments will help those skilled in the art to further understand the present invention, but do not limit the invention in any way. It should be noted that those skilled in the art can make several adjustments and improvements without departing from the concept of the present invention. These all fall within the protection scope of the present invention.

[0065] Example 1

[0066] This embodiment provides a method for predicting tumor neoantigens of endogenous retrovirus origin. Bulk RNA-seq data were analyzed from 924 samples of tumor tissue and adjacent normal tissue from 728 patients with clear cell renal cell carcinoma (ccRCC). The samples involved are shown in Table 1.

[0067] Table 1

[0068] queue name Tumor sample number Number of normal tissue samples Queue 1 55 11 Queue 2 98 61 Queue 3 31 35 Queue 4 11 8 Queue 5 542 72

[0069] Note: In Table 1, the samples for each cohort came from different hospitals.

[0070] Wherein: the sample can be any valid sample that meets the sequencing standards, the tumor sample is a tumor sample that has been removed from the human body, and the control sample is a normal sample that has been removed from the human body.

[0071] Specifically, it includes the following steps (see flowchart). Figure 9 ):

[0072] S0. Bulk RNA-seq sequencing was performed on tumor tissue and adjacent normal tissue samples from tumor patients in Table 1.

[0073] S1. Perform quality control on the bulk RNA-seq sequencing data from step S0; obtain the filtered bulk RNA-seq data.

[0074] In step S1, quality control refers to using the software FastQC and MultiQC to check whether the FastQ files of the bulk RNA-seq sequencing data of each sample in step S1 contain low-quality sequencing fragments (such as adapter sequences, contamination from low-quality sequences, and samples with abnormal GC content); using the software FastP to remove low-quality fragments and the sample data after removing low-quality fragments. The samples mentioned are valid samples that meet sequencing standards, including tumor samples and normal tissue samples.

[0075] S2. Obtain the hervRNA DNA sequence, expression matrix, gene expression matrix, and sample HLA type:

[0076] The bulk RNA-seq data filtered by S1 was processed using a splicing strategy to obtain the DNA sequence of hervRNA, its hervRNA expression matrix in each sample, and the gene expression matrix of each sample, and to predict the human leukocyte antigen (HLA) type of the sample.

[0077] Specifically, the splicing strategy employs a bioinformatics data analysis pipeline called SERVE ("The landscape of hervRNAs transcribed from human endogenous retroviruses across human body sites") to predict the gene annotation files of hervRNAs in tumor samples and obtain hervRNA expression matrices and gene expression matrices. This specifically includes using the "2-pass mode" of the STAR software to align the bulk RNA-seq sequencing data fastq files to the human reference genome GRCh38.p13 (NCBI). The RefSeqassembly was located on the GCF_000001405.39 sequence. Sequencing reads annotated as HERV regions in the HERVd database (human endogenous retrovirus database) were then extracted using the Sambamba software. Genome-guided de novo assembly was performed using the Trinity software to assemble the DNA sequences of hervRNAs for each sample. Subsequently, the TACO software was used to fuse the DNA sequences based on their occurrence in different samples to obtain the hervRNA gene annotation file. Finally, RSEM software was used to quantify expression and obtain the hervRNA expression matrix and gene expression matrix.

[0078] The hervRNA expression matrix consists of hervRNA in each row and a sample in each column.

[0079] The gene expression matrix consists of a row of genes and a column representing a sample.

[0080] The two matrices contain TPM values ​​obtained quantitatively by RSEM.

[0081] Specifically, the steps include the following:

[0082] S21. Alignment of sequencing data with the reference genome and extraction of HERV region reads (using STAR and Sambamba software)

[0083] Using the 2-pass mode of the STAR software, the reads in the fastq file after quality control in step S1 were aligned to the human reference genome GRCh38.p13 sequence to obtain a BAM file sorted according to the alignment position to the reference genome GRCh38.p13. Then, the sequencing fragments (hereinafter referred to as reads) aligned to the HERV region annotated in the HERVd database (human endogenous retrovirus database) in the BAM file were extracted using the sambamba software.

[0084] S22. HervRNA transcript sequence assembly and preliminary screening (using Trinity and GMAP software)

[0085] The Trinity software was used for genome-guided de novo transcription splicing. Based on the fragments extracted in step S21, the DNA sequences of hervRNAs for each sample were spliced. The alignment_and_estimate_abundance.pl module of the Trinity software was used for quantification. Then, the spliced ​​hervRNA sequences were compared with the reference genome GRCh38.p13 sequence using the GMAP software. HervRNAs with a count value greater than 5 and a similarity greater than 96% were retained.

[0086] S23. Using the software TACO, integrate the hervRNAs from each sample based on their presence and expression. For the integrated hervRNAs, those longer than 200 bp are retained. Then, using the software GMAP, the HERV sequence of the filtered hervRNAs is re-aligned to the human GRCh38.p13 reference genome, and the screening process in S22 is repeated. The DNA sequences encoding the hervRNAs obtained in S22 are individually aligned to the HERVd database and annotated as HERV regions. Those that meet the criteria are retained.

[0087] S24. Obtain the hervRNA gene annotation file (using TACO software).

[0088] Based on the hervRNA DNA sequence obtained from S23, the occurrence of the hervRNA DNA sequence in different samples was analyzed. The hervRNA gene annotation file (GTF file) was obtained by performing fusion and other operations using the software TACO.

[0089] S25. Quantitative expression analysis: Obtaining the expression matrix (using RSEM software)

[0090] The hervRNA gene annotation file (GTF file) obtained in step S24 was merged with the human reference genome GTF annotation file, and then the hervRNA expression matrix and gene expression matrix of each sample were obtained by using RSEM software.

[0091] S26. Predict the HLA type of the sample (using seq2hla software)

[0092] Using the FASTQ file after quality control in step S1 as the input file, the HLA-I molecule type of each sample is predicted using seq2hla. In this example, the input file for seq2hla is the FASTQ file of the quality-controlled bulk RNA-seq data, and the other parameters are the default parameters. The distribution of HLA-I molecule types in tumor samples from five independent cohorts is shown in Table 2.

[0093] Table 2 HLA-I molecular information of tumor samples

[0094]

[0095]

[0096] In Table 2:

[0097] HLA-A refers to a gene locus in the human leukocyte antigen class I molecule, located within the major histocompatibility complex (MHC) gene group.

[0098] It is a gene locus in the HLA-B human leukocyte antigen class I molecule, located within the major histocompatibility complex (MHC) gene group.

[0099] The HLA-C gene is located on human chromosome 6.

[0100] S3. Screen for tumor-specific hervRNAs.

[0101] Based on the hervRNA expression matrix of each sample in step S2, tumor-specific hervRNAs are screened out according to the specificity of hervRNA expression in tumor tissues.

[0102] Specifically, it includes:

[0103] S31. Based on the hervRNA TPM expression matrix of each cohort in step S2, retain hervRNAs with TPM greater than 1 in at least 10% of tumor samples.

[0104] 10% means that in S25, the expression matrix of hervRNA will be obtained, the TPM value of a certain hervRNA in all tumor samples will be obtained, and then, in order to reduce the influence of random error, the hervRNA will be screened based on the proportion of the number of samples with a TPM greater than 1 relative to the total number of tumor samples in the cohort, and hervRNAs with a proportion of at least 10% will be retained.

[0105] For example, if a cohort of tumor samples is found to contain a hervRNA, and the cohort contains data from 30 tumor samples, and then it is found that 3 tumor samples in the cohort have a TPM value greater than 1 for the hervRNA, and 3 is 10% relative to 30, then the hervRNA is retained.

[0106] TPM (transcripts per million) is a commonly used standardization method for representing gene expression levels.

[0107] S32. Convert TPM values ​​less than 1 in the hervRNA TPM expression matrix to 1, and then take the logarithm of the entire matrix to the base 2 to obtain the logarithmic hervRNA TPM expression matrix.

[0108] S33, Based on Formula Specificity is calculated, where specificity refers to the TAU value τ of a hervRNA. HervRNAs with a mean expression level in tumor samples greater than that in normal tissue samples and a TAU value greater than 0.8 are considered tumor-specific hervRNAs.

[0109] xi represents the mean expression level of the hervRNA in sample tissue i, which is derived from the logarithmic hervRNA TPM (Transcripts Per Million) matrix.

[0110] n represents the number of tissue types. In this example, there are only tumor tissue and adjacent normal tissue, so the value of n is 2.

[0111] Representative: Values ​​calculated based on the expression status of hervRNA in all tissues (expression status of hervRNA in tumor tissues and normal tissues).

[0112] The purpose of calculating TAU values ​​is to screen for tumor-specific hervRNAs.

[0113] There are two specific screening criteria:

[0114] First, the mean expression level of the hervRNA in the tumor sample must be greater than that in the normal tissue sample. This initial screening, based on the relative expression levels, means that only hervRNAs that are relatively highly expressed in tumor tissues are likely to meet the subsequent criteria.

[0115] Secondly, in addition to this, the calculated TAU value must be greater than 0.8. Only hervRNAs that meet both of these conditions can be considered tumor-specific hervRNAs. The number of hervRNAs screened from the five independent cohorts is shown in [link to table]. Figure 1 ;

[0116] S4. Obtain candidate tumor-specific hervRNAs

[0117] HervRNAs are screened based on their relationship with cytotoxicity. The process of cytotoxic T cells killing target cells involves the specific binding of the peptide-HLA complex (pHLA) on the target cell surface to specific T cell surface TCR receptors, thereby activating cytotoxic T cells to produce effector molecules such as perforin and granzymes, which then kill the target cells. Therefore, if a hervRNA encodes a tumor neoantigen and activates cytotoxic T cells to generate an immune response, the expression level of this hervRNA will be positively correlated with the proportion of activated cytotoxic T cells and the effector molecule levels of cytotoxic T cells. Based on this principle, hervRNAs can be screened by analyzing the correlation between hervRNAs and the proportion of effector T cells or the effector molecule levels of cytotoxic T cells.

[0118] However, in tumor cells, due to factors such as the tumor microenvironment, T cells experience "T cell exhaustion." Specifically, even when T cells encounter the pHLA complex, which specifically binds to their TCRs, they fail to activate and function. Therefore, samples with high levels of exhaustion need to be removed before association analysis to avoid interference with the results. Association analysis specifically refers to:

[0119] S41. Calculate the cytotoxicity score.

[0120] The cytotoxicity score for each sample was calculated based on the gene expression matrix of the samples. Specifically, the cytotoxicity score was calculated by taking the geometric mean of the TPM values ​​of the PRF1 (NCBI: 5551) and GZMB (NCBI: 3002) genes in each tumor sample within the hervRNA TPM expression matrix. This geometric mean was used as the cytotoxicity score for each sample. In this example, the distribution of cytotoxicity scores in different tissues across different cohorts is shown in [reference needed]. Figure 2 ;

[0121] S42. Calculate the sample exhaustion score and remove highly exhausted samples.

[0122] First, the geometric mean of TPM for the CD8A (Entrez ID: 925) and CD8B (Entrez ID: 926) genes in each tumor sample was calculated to characterize the CD8+ T cell level of the sample. Then, the geometric mean of TPM for the five genes PDCD1 (Entrez ID: 5133), CTLA4 (Entrez ID: 1493), LAG3 (Entrez ID: 3902), HAVCR2 (Entrez ID: 84868), and KLRG1 (Entrez ID: 10219) in each sample was calculated. Finally, using the geometric mean of TPM for CD8A and CD8B genes as independent variables and the geometric mean of TPM for PDCD1, CTLA4, LAG3, HAVCR2, and KLRG1 genes as dependent variables, the residuals for each tumor sample were calculated using linear regression. These residuals were used as the exhaustion score for that sample, characterizing the degree of exhaustion in the tumor sample. It is worth noting that the exhaustion score was calculated by performing linear regression independently for each cohort.

[0123] Samples with exhaustion scores greater than 0 are considered to have a high degree of exhaustion and are removed. The number of samples removed in each queue is as follows: 26 samples are removed in the first queue, 35 samples are removed in the second queue, 6 samples are removed in the third queue, 225 samples are removed in the fourth queue, and 10 samples are removed in the fifth queue.

[0124] S43. Association analysis: Calculate correlation coefficients and screen hervRNAs.

[0125] The association analysis refers to calculating the Spearman correlation coefficient. Specifically, after removing tumor samples with exhaustion scores greater than 0, each cohort independently calculates the Spearman correlation coefficient between the cytotoxicity score and the hervRNA TPM for each sample within that cohort, retaining hervRNAs with Spearman correlation coefficients greater than 0. The distribution of the remaining hervRNA numbers in each cohort is shown in the figure. Figure 3 .

[0126] In calculating exhaustion and cytotoxicity scores: It should be noted that the genes used in the above calculations can be derived from single-cell transcriptome data and published literature. Single-cell transcriptome data refers to genes significantly highly expressed in cytotoxic T cells and exhausted T cells, identified through differential gene analysis after cell annotation. If gene function can be combined with screening, the selected genes will be more ideal. Furthermore, published literature refers to genes significantly highly expressed in cytotoxic T cells and exhausted T cells in the corresponding tissues, or genes that have been shown to be closely related to processes such as cell killing or exhaustion, based on currently published relevant research results.

[0127] S5. Obtain the ORF sequence of the candidate hervRNA with high confidence.

[0128] Based on the candidate tumor-specific hervRNA and its corresponding DNA sequence in S4, its possible open reading frame (ORF) sequence is determined, and a high-confidence candidate ORF sequence from the candidate tumor-specific hervRNA is obtained through homology analysis; specifically, the following steps are included:

[0129] S51. Determine possible ORF sequences

[0130] Determining possible ORF sequences involves first using the software gffread to determine the DNA sequence of the hervRNA based on its gene annotation information, and then using the software ORFfinder to determine all possible ORF sequences on the hervRNA.

[0131] Since the fastq data comes from chain-nonspecific library construction, the assembled hervRNA direction can be either forward or reverse. Therefore, the "strand" parameter of ORFfinder is set to both; the "s" parameter of ORFfinder is set to 0 (start codon is only ATG), and "g" is set to 1 (standard codon is used). The number of ORFs determined for each cohort is shown in [link to documentation]. Figure 4 ;

[0132] S52. Ensure homology through BLASTp alignment.

[0133] The homology analysis refers to comparing the ORF sequence with HERV-related protein sequences in the gEVE database and NCBI database using BLASTp, and retaining the ORF sequences with a similarity greater than a set threshold.

[0134] To further ensure the homology between the screened ORF and HERV, the ORF sequences screened in step S51 were compared with the "MetORF" sequence of Homo sapiens in the gEVE database and all protein sequences with taxonomy of 206037 in the NCBI database using BLASTp. ORFs with a similarity greater than 85% were retained to ensure the homology between the ORF and HERV.

[0135] The gEVE database contains ORF sequences on HERVs in the human genome based on computational predictions. "MetORF" refers to an ORF sequence with the first amino acid of Met (methionine) or ATG as the start codon. All protein sequences with a taxonomy of 206037 in the NCBI database represent protein sequences of all human endogenous retroviruses, and 206037 is the taxonomy number of human endogenous retroviruses.

[0136] S53. Screening for commonly existing ORF sequences

[0137] Previously, each queue was operated independently. Due to the existence of systematic and random errors, the analysis and mining of data from a particular queue may uncover ORF sequences unique to that queue. However, these unique ORF sequences may be caused by systematic and random errors, and these ORF sequences themselves are not universal, or may not even exist in reality.

[0138] Considering the systematic and random errors present in the cohort analyses, and to avoid interference from cohort-specific ORFs in subsequent antigen peptide mining, ORF sequences appearing in at least two independent cohorts were selected based on sequence consistency. A total of 873 such ORF sequences were selected from the five cohorts. The ORF relationships between different cohorts are shown in [link to ORF analysis]. Figure 5 .

[0139] S54. Alignment and annotation of ORF sequences of tumor-specific hervRNAs with ORF sequences of hervRNAs expressed in normal tissues.

[0140] The homology analysis also includes BLASTp comparison of the ORF sequence with the ORF of hervRNA expressed in normal tissue, and annotation of the ORF of tumor-specific hervRNA with 100% similarity and all alignments.

[0141] Referring to the hervRNA information in normal tissues in the article "The landscape of hervRNAs stranded from human endogenous retroviruses across human body sites" published by Professor Yang Enze's research group at Peking University in 2022, BLASTp was used to compare the ORF sequences of tumor-specific hervRNAs with the ORF sequences of hervRNAs expressed in normal tissues, and the ORF sequences of tumor-specific hervRNAs with 100% similarity and all matched were labeled.

[0142] It is worth noting that if there is data in multiple independent queues, the ORF sequence that appears in at least two queues is preserved based on sequence consistency.

[0143] S6. Obtain candidate antigen peptides

[0144] Based on the ORF sequences obtained in step S53 that appear in at least two independent cohorts and the types of HLA-I molecules in each sample in S2, the binding energy of peptides with a length of 8-12 amino acids in the ORF to HLA-I molecules is predicted using the software NetMHCpan 4.1. Based on the BA (Bind Affinity) in the prediction results, candidate antigen peptides are screened, and peptides with BA ≤ 500 nM are selected as candidate antigen peptides. In this example, the number of candidate antigen peptides is 45616.

[0145] S7. Perform multi-site analysis on the candidate antigen peptides to obtain the predicted tumor antigen peptides derived from hervRNA.

[0146] Multisite analysis includes analyzing the differences in binding ability of candidate antigen peptides across different HLA alleles in order to select tumor antigen peptides with broad coverage.

[0147] The multi-site analysis refers to first using the software BLAT to align candidate antigen peptides to all loci on the human reference genome (GRCh38.p13) that may encode the antigen peptide, then generating a saf file and using the software featureCount to calculate the expression values ​​of these locus regions, and finally screening based on the expression in normal tissues.

[0148] The predicted antigenic peptide was obtained. In this example, the predicted candidate antigenic peptide had 8-12 amino acids. Since 8-12 amino acids correspond to 24-36 nucleotides according to codons, and considering that HERV is a type of transposon capable of transposition, leading to multiple copies of some HERV sequences in the genome, and given that the number of nucleotides in the DNA sequence encoding the antigenic peptide is much smaller than the approximately 3 billion base pairs in the human genome, even though previous steps using the TAU value algorithm and other methods have ensured the tumor specificity of the HERV RNA, it cannot be guaranteed that the copy number of the DNA sequence encoding the candidate antigenic peptide in the human genome is 1, or that the antigenic peptide can be encoded by DNA sequences at other loci, and that the DNA at that locus can be expressed and encode proteins in normal tissues, thus affecting the specificity of the antigenic peptide. In short, to ensure that the predicted candidate antigenic peptide will not lose specificity due to expression by DNA at other loci in the genome, further determination of its specificity is needed, specifically including:

[0149] S71. Use BLAT alignment to determine peptide coding.

[0150] BLAT was used to align previously predicted candidate antigen peptides to the GRCh38.p13 sequence of the human reference genome, thereby determining whether the peptide could be encoded by DNA sequences at multiple loci in the genome, and the SAF file of the peptide based on the BLAT results.

[0151] In this example, the BLAT parameters are set as follows: t = dnax, q = prot, minScore = 0, tileSize = 5, stepSize = 1, minIdentity = 5, repMatch = 10000000. The reference sequence is GRCh38.p13, and the sequence to be aligned is the peptide sequence obtained in the previous analysis. After obtaining the BLAT psl results, the alignment information that is completely aligned to the genome is extracted (i.e., the parameters "matches" and "Q size" are the same). Then, the peptide ID, chromosome number, coordinates of the alignment position, and whether it is aligned to the positive or negative strand are extracted. The above information is then converted into an SAF file, and the number of alignment sites for each peptide is obtained. The histogram of the number of peptide alignment sites is shown below. Figure 6 ,Depend on Figure 6 It is known that most polypeptides can only be encoded by sequences at one or a few loci;

[0152] S72. Use FeatureCount to quantitatively determine the expression status in normal tissues.

[0153] FeatureCount is used to quantify candidate antigenic peptides based on SAF and BAM files to determine whether they are expressed in normal tissues. In this example, if an antigenic peptide has a CPM (Cost Per Million Mapping Reads) greater than 1 in more than 50% of normal tissues, it is removed. In this example, the FeatureCount parameter is set to the default. The input files are the BAM file of normal tissue and the SAF file generated in the above steps for quantification. The BAM file is obtained by aligning the quality-controlled fastq file of normal tissue to the human reference genome GRCh38.p13 sequence using STAR's "2-pass mode," and then sorting the BAM file according to the alignment position to the reference genome. After obtaining the normal tissue expression matrix of the encoded peptide DNA sequence, this matrix is ​​converted into a CPM matrix (Counts Per Million), a standardized method used to measure the expression level of raw gene or transcript reads across different samples. Peptides with a CPM greater than 1 in more than 50% of normal tissues are removed. The number of candidate antigenic peptides expressed in normal tissues due to multi-site encoding is 1172, and the number of candidate antigenic peptides not expressed in normal tissues is 44444. The relevant distribution is shown in [reference needed]. Figure 7 ;

[0154] S73, Final Prediction Results

[0155] In this example, a total of 44,444 HERV-derived tumor neoantigen peptides were predicted in ccRCC. These 44,444 neoantigens ranged in length from 8 to 12 amino acids and strongly bound to 194 HLA class I molecules (see Table 2) (BA < 500 nM). The distribution of the 10 HLA molecules with the most strongly bound antigenic peptides among different HLA molecules is shown in Table 2. Figure 8Of these, 4773 antigenic peptides strongly bind to HLA-A*02:01. Among these 4773 candidate antigenic peptides, four are neoantigens derived from HERV in ccRCC reported in the literature. Related literature experimentally found that these four neoantigens derived from HERV in ccRCC can be presented by HLA-A*02:01 and activate T cells to generate an immune response. These four antigenic peptides are antigen 1 (MVGPWPRPV), antigen 2 (NSWQEMVPV), antigen 3 (SLNITSCYV), and antigen 4 (LLLQIMRSFV). All four neoantigens passed multi-site screening. However, two neoantigens, antigen 3 (SLNITSCYV) and antigen 4 (LLLQIMRSFV), were marked because their ORF sequences are identical to the ORF sequences on HERV RNA expressed in normal tissues. Table 3 shows the relevant information of neoantigens of HERV origin captured from ccRCCs reported in the literature. The first column is the predicted antigen peptide ID, the second column is the peptide sequence, the third column is the HLA-I molecules that are predicted to strongly bind to the peptide based on binding energy, the fourth column is the number of strongly binding HLA-I molecules, the fifth column is the ORF id of the peptide in the prediction results, the sixth column is whether the ORF of the peptide is completely consistent with the ORF sequence of hervRNA in normal tissue, and the seventh column is whether the peptide will be encoded and expressed in normal tissue due to multi-site coding.

[0156] Table 3

[0157]

[0158]

[0159] The specific embodiments of the present invention have been described above. It should be understood that the present invention is not limited to the specific embodiments described above, and those skilled in the art can make various modifications or variations within the scope of the claims, which do not affect the essence of the present invention.

Claims

1. A method for mining tumor antigens derived from endogenous retroviruses, characterized in that, Includes the following steps: S0. Bulk RNA-seq sequencing was performed on tumor tissue and adjacent normal tissue samples. S1. Perform quality control on bulk RNA-seq sequencing data; S2. Process the bulk RNA-seq data filtered in S1 based on the splicing strategy to obtain the DNA sequence of hervRNA, its hervRNA expression matrix in each sample, and the gene expression matrix of each sample, and predict the human leukocyte antigen (HLA) type of the sample. S3. Based on the hervRNA expression matrix of each sample in step S2, tumor-specific hervRNAs are screened out according to the specificity of hervRNA expression in tumor tissues. S4. Based on the gene expression matrix of each sample in step S2, calculate the exhaustion score and cytotoxicity score of each sample. After excluding samples with high exhaustion scores, perform correlation analysis on the expression level of the tumor-specific hervRNA and the cytotoxicity score based on the remaining samples to obtain candidate tumor-specific hervRNAs related to cytotoxicity. S5. Based on the candidate tumor-specific hervRNA and the corresponding hervRNA DNA sequence described in S4, determine the possible ORF sequence, and obtain the candidate ORF sequence of the candidate tumor-specific hervRNA source with high confidence through homology analysis. S6. Based on the human leukocyte antigen HLA type of the sample described in S2 and the candidate ORF sequence described in step S5, predict the binding energy of 8-12 amino acids in the ORF with HLA class I molecules, and screen out candidate antigen peptides with binding energy less than 500 nM. S7. Perform multi-site analysis on the candidate antigen peptides to obtain the predicted HERV-derived tumor antigen peptides.

2. The method for mining tumor antigens derived from endogenous retroviruses according to claim 1, characterized in that, In S1, quality control refers to using the software FastQC and MultiQC to check whether the FastQ files of the bulk RNA-seq sequencing data of each sample in step S1 contain low-quality sequencing fragments; using the software FastP to remove low-quality fragments, and removing sample data containing low-quality fragments.

3. The method for mining tumor antigens derived from endogenous retroviruses according to claim 1, characterized in that, S2 includes the following steps: S21. The reads in the fastq file after quality control in step S1 are compared with the human reference genome GRCh38.p13 sequence to obtain the BAM file sorted according to the position of the reference genome GRCh38.p13. The sequencing fragments that are aligned to the HERVd database and annotated as HERV regions are extracted from the BAM file using the software sambamba. S22. Based on the fragments extracted in step S21, the DNA sequences of hervRNAs of each sample are spliced ​​together, and the align_and_estimate_abundance.pl module of the software Trinity is used for quantification. Then, the spliced ​​hervRNA sequences are compared with the reference genome GRCh38.p13 sequence using the software GMAP, and hervRNAs with a count value greater than 5 and a similarity greater than 96% are retained. S23. Using the software TACO, integrate the hervRNAs from each sample based on their presence and expression. For the integrated hervRNAs, those longer than 200 bp are retained. Then, using the software GMAP, the HERV sequence of the filtered hervRNAs is re-aligned to the human GRCh38.p13 reference genome, and the screening process in S22 is repeated. The DNA sequences encoding the hervRNAs obtained in S22 are individually aligned to the HERVd database and annotated as HERV regions. Those that meet the criteria are retained. S24. Based on the hervRNA DNA sequence obtained in S23, we analyzed its appearance in different samples and performed fusion operations using the software TACO to obtain the hervRNA gene annotation file GTF file. S25. Merge the hervRNA gene annotation file GTF file obtained in step S24 with the human reference genome GTF annotation file, and then use the software RSEM for quantification to obtain the hervRNA expression matrix and gene expression matrix of each sample. S26. Predict the HLA type of the sample.

4. The method for mining tumor antigens derived from endogenous retroviruses according to claim 1, characterized in that, S3 includes the following steps: S31. Based on the hervRNA TPM expression matrix of each cohort in step S2, retain hervRNAs with TPM greater than 1 in at least 10% of tumor samples. S32. Convert TPM values ​​less than 1 in the hervRNA TPM expression matrix to 1, and then take the logarithm of the entire matrix to the base 2 to obtain the logarithmic hervRNA TPM expression matrix. S33, Based on Formula Specificity is calculated, where specificity refers to the TAU value τ of a hervRNA. HervRNAs with a mean expression level in tumor samples greater than that in normal tissue samples and a TAU value greater than 0.8 are considered tumor-specific hervRNAs. xi represents the mean expression level of the hervRNA in sample tissue i, which is derived from the logarithmic hervRNA TPM (Transcripts Per Million) matrix. n represents the number of tissue types. In this example, there are only tumor tissue and adjacent normal tissue, so the value of n is 2. Representative: Values ​​calculated based on the expression of hervRNA in tumor and normal tissues.

5. The method for mining tumor antigens derived from endogenous retroviruses according to claim 1, characterized in that, In S4, the cytotoxicity score is calculated as follows: The geometric mean of the TPM values ​​of the PRF1 and GZMB genes in each tumor sample in the hervRNA TPM expression matrix was calculated and used as the cytotoxicity score of the sample.

6. The method for mining tumor antigens derived from endogenous retroviruses according to claim 1, characterized in that, In S4, the exhaustion score is calculated as follows: First, the geometric mean of TPM for the CD8A and CD8B genes in each tumor sample was calculated to characterize the CD8+ T cell level in the sample. Then, the geometric mean of TPM for the five genes PDCD1, CTLA4, LAG3, HAVCR2 and KLRG1 in each sample was calculated. Finally, the geometric mean of TPM for the CD8A and CD8B genes was used as the independent variable. The geometric mean of TPM of five genes, PDCD1, CTLA4, LAG3, HAVCR2 and KLRG1, was used as the dependent variable. The residuals of each tumor sample were calculated by linear regression. The residuals were used as the exhaustion score of the sample to characterize the exhaustion degree of the tumor sample. The exhaustion score was calculated by performing linear regression on each cohort independently.

7. The method for mining tumor antigens derived from endogenous retroviruses according to claim 1, characterized in that, In S5, homology analysis refers to comparing ORF sequences with HERV-related protein sequences in the gEVE and NCBI databases using BLASTp, and retaining ORF sequences with similarity greater than a set threshold.

8. The method for mining tumor antigens derived from endogenous retroviruses according to claim 7, characterized in that, In S5, the homology analysis also includes comparing the ORF sequence with the ORF of hervRNA expressed in normal tissue using BLASTp, and annotating the ORF of tumor-specific hervRNA with 100% similarity and all alignments.

9. The method for mining tumor antigens derived from endogenous retroviruses according to claim 7, characterized in that, Multisite analysis includes analyzing the differences in binding ability of candidate antigen peptides across different HLA alleles in order to select tumor antigen peptides with broad coverage.

10. A tumor antigen derived from endogenous retroviruses obtained by a method for mining tumor antigens derived from endogenous retroviruses as described in any one of claims 1-9.

Images