Method for optimizing molecular marker algorithms based on machine learning and immune escape mechanisms

Abstract

The application discloses a kind of based on machine learning and immune escape mechanism optimization molecular marker algorithm method, the method includes: obtaining the NGS sequencing data of tumor tissue and normal tissue, the NGS sequencing data of tumor tissue and matched normal tissue are quality controlled, alignment, obtain original alignment file, the original alignment file is pretreated to obtain final alignment file;Tumor somatic single nucleotide variation and small fragment insertion and deletion detection are carried out to alignment file by software, obtain somatic variation information;Construct multivariate linear regression model, calculate the prediction ability of neoantigen of SNV site and INDEL site;The TMB calculation method is optimized in combination with the two aspects of factors of different neoantigen prediction ability and tumor patient immune escape mechanism.The application more truly reflects the real immune state of patient, and also has better prompting effect on patient immunotherapy response and prognosis.

Description

Technical Field

[0001] This invention relates to the field of molecular biomarker algorithm optimization, and in particular to a method for optimizing molecular biomarker algorithms based on machine learning and immune escape mechanisms. Background Technology

[0002] The study of tumor mutation burden (TMB) as a prognostic biomarker for immunotherapy in non-small cell lung cancer (NSCLC) originated from a 2015 Science study, which found that NSCLC patients with higher median TMB had longer progression-free survival (PFS). Subsequently, several large studies, including Check Mate-026 and Check Mate-227, confirmed that TMB can serve as an independent biomarker for immunotherapy. In 2020, pembrolizumab, a drug targeting programmed cell death protein-1 (PD-1), one of the immune checkpoints, was approved by the FDA for the treatment of patients with solid tumors with high mutation burden, and TMB has been widely used as a biomarker for tumor immunotherapy.

[0003] In theory, a higher mutation burden can increase the likelihood of neoantigen generation, thereby increasing the likelihood of immune recognition and killing of tumor cells, thus affecting the efficacy of immunotherapy. Therefore, current research on TMB-related scientific research and clinical testing products focuses on optimizing sequencing strategies and mutation detection algorithms to accurately detect patients' mutation status. The number of high-quality mutations detected is statistically analyzed, and combined with the coverage of the testing product, to assess the tumor mutation burden of each patient, ultimately providing appropriate immunotherapy plans for cancer patients.

[0004] The primary pathway by which mutated tumor neoantigens exert their immune effects is through the delivery of neoantigen peptides to CD8+ T cells by HLA class I molecules, thereby triggering a downstream immune response. However, in advanced stages of cancers such as lung cancer, a large number of cancer patients experience loss of heterozygosity in their HLA genes. The missing HLA type loses its ability to bind to the original mutated peptides, meaning not all predicted neoantigen peptides are presented to CD8+ T cells. Therefore, traditional methods that only consider the number of neoantigens generated by mutations without taking into account the patient's HLA loss-of-heterozygosity status are insufficient. Based on current technological advancements, we can correct biomarkers based on immune escape mechanisms. Summary of the Invention

[0005] The purpose of this invention is to provide a method for optimizing molecular biomarker algorithms based on machine learning and immune escape mechanisms, which more realistically reflects the patient's true immune status and provides better indications for the patient's response to immunotherapy and prognosis.

[0006] To achieve the above objectives, the present invention provides the following solution:

[0007] A method for optimizing molecular biomarker algorithms based on machine learning and immune escape mechanisms includes the following steps:

[0008] S1. Obtain NGS sequencing data of tumor tissue and normal tissue, including whole genome sequencing data or whole exon capture sequencing data based on next-generation sequencing.

[0009] S2. Perform quality control and alignment on the NGS sequencing data of tumor tissue and paired normal tissue to obtain the original alignment file. Preprocess the original alignment file to obtain the final alignment file.

[0010] S3. Perform tumor somatic cell single nucleotide variant and small fragment insertion / deletion detection on the final alignment file, and filter the original variant list obtained from the detection to obtain somatic cell variant information.

[0011] S4. Perform HLA typing on NGS sequencing data of normal tissues from tumor patients to obtain the patient's HLA four-position typing results;

[0012] S5. Based on the somatic cell variation information filtered in S3, the NetMHCpan deep learning algorithm is used to predict the tumor neoantigen peptides of patients and screen the tumor neoantigen peptides.

[0013] S6. Based on the final alignment file obtained in S2, assess the HLA heterozygosity loss of tumor samples;

[0014] S7. Statistically count the number of SNV and INDEL mutations in the TCGA database and the number of predicted neoantigen peptides to construct a multiple linear regression model and calculate the neoantigen prediction capabilities of SNV and INDEL sites.

[0015] S8. The TMB calculation method was optimized by combining two factors: the different predictive abilities of neoantigens and the immune escape mechanism of tumor patients.

[0016] Preferably, in S2, the NGS sequencing data of tumor tissue and paired normal tissue are quality controlled and compared to obtain the original alignment file. The original alignment file is then preprocessed to obtain the final alignment file, specifically including:

[0017] The NGS sequencing data of tumor tissue and paired normal tissue were subjected to quality control, and low-quality sequences, low-complexity sequences and adapter sequences in the NGS sequencing data were removed using FastP software.

[0018] Then, the filtered clean reads were aligned to the human reference genome hg19 using BWA-MEN software to obtain the original alignment file.

[0019] The original alignment file was preprocessed as follows: (1) Sorting was performed using samtools software; (2) Repeated sequences were marked using Picard software; (3) Local re-alignment and base quality correction were performed using GATK software to obtain the final alignment file.

[0020] Preferably, in S3, the final alignment file is subjected to tumor somatic single nucleotide variant and small fragment insertion / deletion detection, and the original variant list obtained is filtered to obtain somatic variant information, specifically including:

[0021] Based on the final alignment file, MuTect2 software was used to detect single nucleotide variants and small fragment insertions / deletions in tumor somatic cells;

[0022] The original variant list obtained from the detection is filtered by the following steps to obtain somatic variant information. The filtering steps are as follows: (1) Retain non-silent variants, including missense mutations, nonsense mutations, frameshift mutations, and splice site variants; (2) Retain variant sites with high confidence: the sequence alignment quality of the site is greater than 30, the base sequencing quality is greater than 20, the mutation positive and negative strand readings are unbiased, the variant site depth is greater than 50X, and the mutant allele frequency is greater than 0.03; (3) Remove mutation sites with high mutation frequency and retain sites with a frequency of less than 0.01 in the TCGA database.

[0023] Preferably, in S5, based on the somatic cell variation information filtered in S3, the NetMHCpan deep learning algorithm is used to predict the patient's tumor neoantigen peptides, and the screening of tumor neoantigen peptides specifically includes:

[0024] Based on the filtered somatic mutation information in S3, the NetMHCpan deep learning algorithm is used to predict the tumor neoantigen peptides of patients. The binding affinity inhibition concentrations of site mutants and paired wild-types with the predicted peptides are calculated respectively. Based on the affinity inhibition concentrations, the differential affinity fold for each site is calculated, and then the predicted tumor neoantigen peptides are screened.

[0025] Preferably, in S7, the number of SNV and INDEL mutations in the TCGA database and the number of predicted neoantigen peptides are statistically analyzed to construct a multiple linear regression model. The ability to predict neoantigens at SNV and INDEL sites specifically includes:

[0026] (1) Obtain TCGA mutation data, classify mutations according to TCGA rules, filter mutation data, and retain mutations listed as PASS in the FILTER column; then use vcf2maf software to convert the maf file in the mutation data into a standard vcf file, and sort the mutation data in the vcf file according to chromosome position.

[0027] (2) Obtain the HLA type classification data of patients from the TCGA database;

[0028] (3) After matching the mutation information and HLA typing data of the sample, the tumor neoantigen peptides of the sample are predicted and filtered;

[0029] (4) Count the number of SNV and INDEL mutations in patients in the TCGA database and the number of predicted tumor neoantigen peptides; construct a multiple linear regression model to calculate the ability of SNV sites and INDEL sites to predict neoantigens.

[0030] Count Neo =β+α1*Count SNV +α2*Count INDEL

[0031] Where, Count Neo Count is the number of neoantigens predicted for patient samples included in the model. SNV To include patient samples in the model for predicting the number of SNV sites for neoantigens, Count INDEL The number of INDEL sites for neoantigens was used to predict the number of neoantigens in the patient samples included in the model. α1 is the weight of the ability of SNV to generate neoantigens, α2 is the weight of the ability of INDEL to generate neoantigens, and β is the intercept value of the multiple linear regression model.

[0032] Preferably, in S8, the optimization of the TMB calculation method by combining two factors—different neoantigen prediction capabilities and the immune escape mechanism of tumor patients—specifically includes:

[0033] The original TMB calculation formula is as follows:

[0034]

[0035] Where, Count snv+indel For the mutation sites ultimately detected in the patient, including SNV and INDEL, Countbed The size of the genome covered by the detection kit;

[0036] Considering both the varying predictive abilities of neoantigens and the immune escape mechanisms in cancer patients, the following optimization method for TMB calculation is proposed:

[0037] (1) Combining the list of neoantigen predicted peptides obtained in S5 and the HLA heterozygosity loss status obtained in S6, the proportion of the number of neoantigens bound by the missing HLA types in the total number of neoantigens bound by all alleles in the patient was calculated, and the ratio was used as the neoantigen failure status caused by HLA heterozygosity loss in the patient.

[0038] (2) Based on the neoantigen prediction ability of SNV and INDEL obtained by linear regression, the TMB value after mutation type correction was obtained by calculating SNV and INDEL of tumor patients with different production weights.

[0039] (3) Combining the above two steps, the bias in neoantigen assessment caused by immune escape and mutation type is corrected, and the calculation formula is as follows:

[0040]

[0041] Where, r hla The correction coefficient for TMB based on the immune escape mechanism, TMB adj Neo represents the corrected tumor mutation burden value. g* This represents the number of high-quality neoantigen peptides from the screened patients that have bound to the g* type HLA alleles. Here, hhlaloh is the list of HLA alleles identified as lacking heterozygosity in step S6, and gi is the allele in the hhlaloh list; total is the list of all HLA alleles from the patient obtained in step S4, and gj is the allele in the total list. Count snv The number of single nucleotide mutations ultimately detected in the patient, Count indel Count represents the number of insertion / deletion mutations ultimately detected in the patient. bed The size of the genome covered by the detection kit.

[0042] Preferably, step S8 further includes optimizing the TMB calculation method based on the different neoantigen prediction capabilities of SNV and INDEL:

[0043] Based on the neoantigen predictive ability of SNV and INDEL obtained from linear regression, the TMB value after mutation type correction was calculated for SNV and INDEL in tumor patients with different production weights. The calculation formula is as follows:

[0044]

[0045] Where, Count snv The number of single nucleotide mutations ultimately detected in the patient, Count indel Count represents the number of insertion / deletion mutations ultimately detected in the patient. bed The size of the genome covered by the detection kit.

[0046] According to specific embodiments provided by the present invention, the present invention discloses the following technical effects:

[0047] This invention combines machine learning and other methods to mine tumor neoantigen data from public databases, constructs a multiple linear regression model, calculates the predictive power weights of neoantigens for two different mutation types, SNV and INDEL, and corrects the TMB calculation process based on these weights. This invention also uses cutting-edge algorithms to assess the patient's immune escape status, reflecting the patient's true neoantigen status. By optimizing the algorithm from two perspectives, this invention more realistically reflects the patient's true immune status compared to existing algorithms, and provides better indications for the patient's immunotherapy response and prognosis. Attached Figure Description

[0048] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0049] Figure 1 A flowchart of a method for optimizing molecular marker algorithms based on machine learning and immune escape mechanisms provided by this invention;

[0050] Figure 2 Forest graph showing the total lifetime of the three algorithms;

[0051] Figure 3a Original TMB of 88 non-small cell lung cancer patients raw Overall survival KM curve;

[0052] Figure 3b For 88 patients with non-small cell lung cancer, TMB was calculated based on combined mutation type and immune escape correction. adj Overall survival KM curve;

[0053] Figure 3c TMB after mutation type correction in 88 non-small cell lung cancer patients madj Overall survival KM curve;

[0054] Figure 4 The graph shows the no-progression-survival forest for the three algorithms.

[0055] Figure 5a Original TMB of 88 non-small cell lung cancer patients raw The progression-free survival KM curve;

[0056] Figure 5b For 88 patients with non-small cell lung cancer, TMB was calculated based on combined mutation type and immune escape correction. adj The progression-free survival KM curve;

[0057] Figure 5c TMB after mutation type correction in 88 non-small cell lung cancer patients madj The progression-free survival KM curve;

[0058] Figure 6a Original TMB of 88 patients raw Intergroup differences in response to immunotherapy;

[0059] Figure 6b For 88 patients, combined mutation type and TMB after immune escape correction adj Intergroup differences in response to immunotherapy;

[0060] Figure 6c TMB after mutation type correction in 88 patients madj Intergroup differences in response to immunotherapy. Detailed Implementation

[0061] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0062] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0063] like Figure 1 As shown, the present invention provides a method for optimizing molecular biomarker algorithms based on machine learning and immune escape mechanisms, comprising the following steps:

[0064] S1. Obtain NGS sequencing data of tumor tissue and normal tissue, including whole genome sequencing data or whole exon capture sequencing data based on next-generation sequencing.

[0065] S2. Perform quality control and alignment on the NGS sequencing data of tumor tissue and paired normal tissue to obtain the original alignment file. Preprocess the original alignment file to obtain the final alignment file.

[0066] S3. The software is used to detect tumor somatic single nucleotide variants and small fragment insertions / deletions in the final comparison file, and the original variant list obtained by the detection is filtered to obtain somatic variant information.

[0067] S4. Perform HLA typing on NGS sequencing data of normal tissues from tumor patients to obtain the patient's HLA four-position typing results;

[0068] S5. Based on the somatic cell variation information filtered in S3, the NetMHCpan deep learning algorithm is used to predict the tumor neoantigen peptides of patients and screen the tumor neoantigen peptides.

[0069] S6. Based on the final alignment file obtained in S2, software is used to assess the HLA heterozygosity loss of tumor samples;

[0070] S7. Statistically count the number of SNV and INDEL mutations in the TCGA database and the number of predicted high-quality neoantigen peptides to construct a multiple linear regression model to calculate the neoantigen prediction capabilities of SNV and INDEL sites.

[0071] S8. The TMB calculation method was optimized by combining two factors: the different predictive abilities of neoantigens and the immune escape mechanism of tumor patients.

[0072] Specifically, the method of the present invention includes:

[0073] 1.1 Obtaining raw data

[0074] Obtain NGS sequencing data from tumor tissue and normal tissue (or blood cells), including whole-genome sequencing (WGS) or whole-exome capture sequencing (WES) data based on next-generation sequencing. Sequencing data can be paired-end or single-end sequencing. Sequencing depth requirements: WGS sequencing depth > 30X, WES sequencing depth > 100X. In this embodiment of the invention, the samples are paired-end sequencing data, using WES sequencing data from the PE150 sequencing strategy.

[0075] 1.2 Data Preprocessing

[0076] Quality control was performed on the raw sequencing data of tumor tissue and paired normal tissue. FastP software was primarily used to remove low-quality, low-complexity, and adapter sequences from the raw sequencing data. Then, BWA-MEN software was used to align the filtered clean reads to the human reference genome hg19, obtaining the raw alignment file (BAM). The raw BAM file underwent the following processing: ① Sorting using SAMTools software; ② Marking repetitive sequences using Picard software; ③ Performing local realignment and base quality correction using GATK software, resulting in the final version of the BAM file for downstream data analysis.

[0077] 1.3 Gene mutation detection and filtering

[0078] Based on BAM files from tumor samples and paired normal tissues, MuTect2 software was used to detect single nucleotide variants (SNVs) and small insertion / deletion (INDEL) mutations in tumor somatic cells. The raw variant list was filtered using the following steps to obtain highly reliable somatic variant information: ① Non-silent variants were retained, including missense mutations, nonsense mutations, frameshift mutations, and splice site variants; ② High-reliability variant sites were retained, requiring a sequence alignment quality greater than 30, a base sequencing quality greater than 20, no bias in mutation positive / negative strand readings, a variant site depth greater than 50X, and a mutant allele frequency greater than 0.03; ③ High-frequency mutation sites in the human population were removed, retaining sites with a frequency less than 0.01 in human population databases (ExAC, gnomAD).

[0079] 1.4 Human Leukocyte Antigen (HLA) Typing

[0080] HLA typing data of cancer patients is obtained, and HLA typing of sequencing data from normal tissues (or blood cells) of cancer patients is performed using the software OptiType to obtain the patient's four-position HLA typing results. The main steps are as follows:

[0081] Using the HLA DNA reference sequence provided by OptiType as a benchmark, the filtered cleanreads were aligned to the genomic region where HLA is located using Razers3 software to obtain sequencing data within the HLA region. Then, the alignment file (BAM) obtained was converted into a sequencing sequence file (fastq) using samtools software. The fastq file obtained was then analyzed for HLA sequence using OptiType analysis software to obtain the patient's HLA typing results.

[0082] 1.5 Screening for high-quality neoantigen peptides in tumors

[0083] Based on the above-mentioned filtered and annotated mutations, the deep learning algorithm NetMHCpan (https: / / services.healthtech.dtu.dk / services / NetMHCpan-4.1 / ) was used to predict potential tumor neoantigen peptides for each patient. The binding affinity inhibitory concentrations (IC50) between the mutant and paired wild-type at each mutation site were calculated. Only predicted peptides with lengths between 8 and 11 meters were considered. The IC50 of the binding affinity inhibitory concentrations between the mutant and paired wild-type and the predicted peptides were calculated separately. Based on these IC50 concentrations, the differential affinity fold (DAI) for each site was calculated. The predicted peptides were then screened in two ways: If a neoantigen peptide was predicted in the wild-type, the screening criteria were: ① mutant IC50 <= 500 nM and ② DAI > 1; if no neoantigen peptide was predicted in the wild-type, the screening criterion was: mutant IC50 <= 500 nM. Finally, high-quality tumor neoantigen peptides from 6562 patients were obtained for downstream analysis.

[0084]

[0085] 1.6 Determination of HLA heterozygosity

[0086] Based on the BAM files of tumor tissue and paired normal tissue obtained in section 1.2, the HLA loss of heterozygosity in tumor samples was assessed using HLALOH software. Heterozygosity was determined for each patient's HLAA, HLAB, and HLA C gene alleles. If homozygous genes were present in all three HLA genes, loss of heterozygosity was not considered. The calculation method involved analyzing the mismatch sites and genomic coverage of paired alleles in tumor and control tissues. A paired t-test was used to detect the probability of allele imbalance for each HLA gene. A p-value less than 0.1 indicated that an HLA gene allele had experienced a loss of heterozygosity event.

[0087] 1.7 Differential Prediction of Tumor Neoantigens by SNV and INDEL

[0088] 1.7.1 Obtain TCGA mutation data. Following TCGA rules for mutation classification, filter the mutation data, retaining only mutations marked "PASS" in the FILTER column. Then, use the vcf2maf software to convert the MAF file to a standard VCF file, and sort the mutation data in the VCF file according to chromosome location.

[0089] 1.7.2 Obtain HLAI typing data of 7300 TCGA patients from the TCIA database.

[0090] 1.7.3 After matching the mutation information and HLA typing data of the sample, predict and filter the tumor neoantigens of the sample according to the steps in 1.5 above.

[0091] 1.7.4 The number of SNV and INDEL mutations and the number of predicted high-quality neoantigen peptides were statistically analyzed in 6562 patients in the TCGA database. Considering that the ability of tumor mutations to generate neoantigen peptides is mainly related to the binding of mutations to peptides, and that neoantigens do not differ significantly across cancer types, with the limited cancer-specificity actually resulting from cancer-specific mutations, a multiple linear regression model was constructed in 6562 pan-cancer patients to calculate the predictive power of neoantigens at SNV and INDEL sites.

[0092] Count Neo =β+α1*Count SNV +α2*Count INDEL

[0093] Where, Count Neo The number of neoantigens predicted by the patient samples included in the model through the steps described in 1.7.3 above is used to count... SNV To include patient samples in the model for predicting the number of SNV sites for neoantigens, Count INDEL To include patient samples in the model for predicting the number of INDEL sites for neoantigens, α1 represents the weight of the ability of SNVs to generate neoantigens, α2 represents the weight of the ability of INDELs to generate neoantigens, and β represents the intercept of the multiple linear regression model, i.e., the value after controlling for all other independent variables (Count). SNV and Count INDEL After the influence of ), the dependent variable Count Neo The average baseline value.

[0094] 1.7.5 The determination coefficient R2 of the multiple linear regression model established above is 0.92 (R2 is the proportion of the variance explained by the model in the total variance, ranging between 0 and 1, with a larger value indicating a higher degree of explanation of the variance by the model), where α1 is 0.2939 (P<2e-16), and the coefficient / weight α2 of the INDEL variable is 0.5093 (P<2e-16).

[0095] 1.8 Optimization Algorithm for TMB, a Prognostic Molecular Marker in Immunotherapy Patients

[0096] 1.8.1 Original TMBraw Calculation

[0097] Based on the high-quality mutation sites obtained in section 1.3, and considering the coverage of the detection kit, the original TMB was calculated using currently accepted methods. The calculation formula is as follows:

[0098]

[0099] Where, Count snv+indel For the mutation sites ultimately detected in the patient, including SNV and INDEL, Count bed The size of the genome covered by the detection kit.

[0100] 1.8.2 Combined immune escape and mutation type optimization TMBadj

[0101] To optimize the TMB calculation method in section 1.8.1, considering both the different predictive abilities of mutation types and the immune escape mechanisms in cancer patients, the optimization method is as follows: ① Combining the neoantigen prediction list for cancer patients obtained in step 1.5 and the HLA heterozygosity loss status obtained in step 1.6, the proportion of neoantigens bound by the missing HLA types in the total number of neoantigens bound by all HLA alleles (HLAA, HLAB, and HLA C) in the patients is calculated. This proportion is used as the neoantigen ineffectiveness status caused by HLA heterozygosity loss in cancer patients; ② Based on the neoantigen predictive ability of SNVs and INDELs obtained from linear regression in section 1.7, SNVs and INDELs in cancer patients are calculated with different generation weights to obtain the mutation type-corrected TMB value. Combining the above two steps, the neoantigen assessment bias caused by immune escape and mutation type is corrected, and the calculation formula is as follows:

[0102]

[0103] Where, r hla This represents the correction coefficient for TMB based on the immune escape mechanism (using HLA heterozygosity loss status as the criterion for judging immune escape). TMB adj Neo represents the corrected tumor mutation burden value. g* This represents the number of high-quality neoantigen peptides from the screened patients that have bound to the g* type HLA alleles. Here, hhlaloh is the list of HLA alleles identified as lacking heterozygosity in step S6, and gi is the allele in the hhlaloh list; total is the list of all HLA alleles from the patient obtained in step S4, and gj is the allele in the total list. Count snv The number of single nucleotide mutations ultimately detected in the patient, Count indel Count represents the number of insertion / deletion mutations ultimately detected in the patient. bed The size of the genome covered by the detection kit.

[0104] 1.8.3 Optimization of TMBmadj based on mutation type

[0105] The method described in 1.8.2 above requires calculating the number of neoantigens, which necessitates providing the patient's whole genome sequencing (WGS) or whole exon capture sequencing (WES) data. Considering that in actual clinical testing, some patients may not be able to provide WGS or WES next-generation sequencing data or HLA typing results, this patent provides an alternative optimization scheme.

[0106] The optimization method is as follows: Based on the differential weights of SNV and INDEL calculated in step 1.7, the mutation type-corrected TMB value is obtained. Although this scheme does not consider the impact of the patient's immune escape mechanism on neoantigens, it corrects for different mutation types by using different prediction weights. Compared to the method in 1.8.1, it still better reflects the immunotherapy status of cancer patients. The specific calculation formula is as follows:

[0107]

[0108] Where, Count snv The number of single nucleotide mutations ultimately detected in the patient, Count indel Count represents the number of insertion / deletion mutations ultimately detected in the patient. bed The size of the genome covered by the detection kit.

[0109] In a specific embodiment, real queue data is used to calculate the original TMB. raw Combined mutation types and TMB after immune escape correction adj and TMB after mutation type correction madj To illustrate the use of TMB by comparing its clinical predictive capabilities. adj and TMB madj The advantages.

[0110] In embodiments of the present invention, whole-exome capture sequencing (WES) data from next-generation sequencing is included. Paired-end sequencing is preferred in these embodiments, using the PE150 sequencing strategy. The NGS sequencing data is preprocessed to obtain intermediate data files.

[0111] 2.1 Sample Information

[0112] Tumor tissues from 88 non-small cell lung cancer patients who underwent immunotherapy were collected and analyzed by WES, with an average depth greater than 100x.

[0113] 2.2 Specific Implementation Process

[0114] The sequencing data were preprocessed in step 1.2, and then the mutation detection method described in step 1.3 was used to detect the high-confidence mutation list for each sample. Step 1.4 determined the patient's HLA allele type, step 1.5 determined the predicted high-quality neoantigen peptides for the patient, and step 1.6 obtained the patient's HLA loss of heterozygosity status.

[0115] The original TMB is calculated using steps 1.8.1, 1.8.2, and 1.8.3 respectively. raw Combined mutation types and TMB after immune escape correction adj and TMB after mutation type correction madj The study aims to predict the therapeutic efficacy, progression-free survival, and overall survival of immunotherapy samples.

[0116] 2.3 Verification Results

[0117] 2.3.1 Predictive effect of three algorithms on overall survival (OS) of patients

[0118] For the original TMB raw Combined mutation types and TMB after immune escape correction adj and TMB after mutation type correction madj Each group was divided into high and low groups based on its median value as a threshold. Univariate Cox regression analysis was performed using R software, and Kaplan-Meier survival curves for overall survival (OS) were plotted.

[0119] Figure 2 OS forest graphs of three algorithms show TMB adj and TMB madj Compared to the original TMB raw It is better at predicting patients' overall survival status. When grouping patients using the median as the threshold, TMB... adj and TMB madj Both indicators were significantly correlated with the patient's overall survival status. Figure 3a The results show that: with the original TMB raw When the median value is the threshold, the original TMB raw There was no significant difference in overall survival between the high and low groups (HR = 1.9, p = 0.081), and the original TMB... raw It does not reflect the patient's overall survival status;

[0120] Figure 3b The results showed that, based on combined mutation type and immune escape-corrected TMB adj The median was used as a threshold to divide the data into high and low groups, and TMB was used for this purpose. adj Patients in the low-trimester group had significantly worse overall survival than those in the TMB group. adjPatients in the high-risk group had a 2.5 times higher risk of death than those in the low-risk group (HR = 2.5; p = 0.02), indicating that TMB... adj It has predictive value for the overall survival of patients after immunotherapy;

[0121] Figure 3c The results showed that TMB corrected for mutation type madj The median was used as a threshold to divide the data into high and low groups, and TMB was used for this purpose. madj Patients in the low-trimester group had significantly worse overall survival than those in the TMB group. madj Patients in the high-risk group (HR=2.3; p=0.033) indicate that the TMB after mutation type correction is significantly higher. madj It also has predictive value for overall survival in cancer patients after immunotherapy;

[0122] 2.3.2 Predictive value of three algorithms for progression-free survival (PFS)

[0123] The three biomarkers were grouped using the same method as in step 2.3.1 above, and univariate Cox regression analysis was performed using R software to plot the Kaplan-Meier survival curves for PFS.

[0124] Figure 4 The study uses forest graphs to illustrate the prognostic results of three algorithms for progression-free survival (TMB), revealing that TMB... adj and TMB madj Compared to the original TMB raw It is better able to predict the patient's progression-free survival status.

[0125] Figure 5a The results show that: with the original TMB raw When the median value is the threshold, the original TMB raw There was no significant difference in progression-free survival between the high and low groups (HR = 1.7, p = 0.138), and the original TMB... raw It is impossible to predict whether a patient will relapse after immunotherapy;

[0126] Figure 5b The results showed that, based on combined mutation type and immune escape-corrected TMB adj The median was used as a threshold to divide the data into high and low groups, and TMB was used for this purpose. adj Patients in the lower group had significantly worse progression-free survival than those with TMB. adj Patients in the high-risk group (HR = 2.1; p = 0.027) indicate that TMB... adj It has indicative significance for the relapse rate in patients receiving immunotherapy;

[0127] Figure 5c The results showed that TMB corrected for mutation type madjThe median was used as a threshold to divide the data into high and low groups, and TMB was used for this purpose. madj Patients in the low-trimester group had significantly worse overall survival than those in the TMB group. madj The high-risk group (HR=2; p=0.048) indicates that the TMB after mutation type correction is also indicative of relapse in patients receiving immunotherapy.

[0128] 2.3.3 Predictive role of three biomarkers in the efficacy of immunotherapy

[0129] The association between the above three biomarkers and the efficacy of immunotherapy was analyzed, with 56 people in the immunotherapy response group (PRorCR) and 32 people in the non-response group (PD or SD).

[0130] Figure 6a Demonstrates the original TMB of non-small cell lung cancer patients after immunotherapy. raw No statistically significant difference was found between the immunotherapy response and non-response groups (p = 0.053); Figure 6b TMB demonstrates immune escape combined with mutation type correction mdj A significant difference was observed between the immunotherapy response and non-response groups (p = 0.0096), with the TMB of immunotherapy-responsive patients being significantly higher. adj Significantly higher than the non-responder group; Figure 6c TMB demonstrates mutation type correction madj A significant difference was observed between the immunotherapy response and non-response groups (p = 0.038), with the TMB of immunotherapy-responsive patients being significantly higher. madj The response rate was significantly higher than that of the non-responder group.

[0131] The above results fully demonstrate that the TMB corrected by the two correction methods is a better predictor of the efficacy of immunotherapy and the patient's survival outcome compared to the original TMB.

[0132] This invention provides two methods for TMB optimization:

[0133] (1) A more accurate TMB correction method is provided based on the patient's availability of WES or WGS sequencing data; (2) Considering the cost of testing, an alternative TMB correction method is provided when the patient only has mutation data and cannot provide raw sequencing data such as WES or WGS. This invention is mainly based on the basic principle of TMB as a biomarker for immunotherapy and the immune mechanisms involved, and optimizes the shortcomings of existing TMB methods. On the basis of the optimized method being widely applicable, it aims to reflect the patient's true immune status as accurately as possible.

[0134] Through the above description of the embodiments, those skilled in the art can clearly understand that each embodiment can be implemented by means of software plus necessary general-purpose hardware platforms, and of course, it can also be implemented by hardware. Based on this understanding, the above technical solutions, in essence or the part that contributes to the prior art, can be embodied in the form of a software product. This computer software product can be stored in a computer-readable storage medium, such as ROM / RAM, magnetic disk, optical disk, etc., and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute the methods described in the various embodiments or some parts of the embodiments.

[0135] This document uses specific examples to illustrate the principles and implementation methods of the present invention. The descriptions of the above embodiments are only for the purpose of helping to understand the method and core ideas of the present invention. Furthermore, those skilled in the art will recognize that, based on the ideas of the present invention, there will be changes in the specific implementation methods and application scope. Therefore, the content of this specification should not be construed as a limitation of the present invention.

Claims

1. A method for optimizing molecular biomarker algorithms based on machine learning and immune escape mechanisms, characterized in that, Includes the following steps: S1. Obtain NGS sequencing data of tumor tissue and normal tissue, including whole genome sequencing data or whole exon capture sequencing data based on next-generation sequencing. S2. Perform quality control and alignment on the NGS sequencing data of tumor tissue and paired normal tissue to obtain the original alignment file, and preprocess the original alignment file to obtain the final alignment file; S3. Perform tumor somatic cell single nucleotide variant and small fragment insertion / deletion detection on the final comparison file, and filter the original variant list obtained from the detection to obtain somatic cell variant information; S4. Perform HLA typing on the NGS sequencing data of normal tissues from tumor patients to obtain the patient's HLA four-position typing results; S5. Based on the somatic cell variation information filtered in S3, the NetMHCpan deep learning algorithm is used to predict the tumor neoantigen peptides of the patient, and the tumor neoantigen peptides are screened. S6. Based on the final alignment file obtained in S2, assess the HLA heterozygosity loss of tumor samples; S7. Statistically count the number of SNV and INDEL mutations in the TCGA database and the number of predicted neoantigen peptides to construct a multiple linear regression model and calculate the neoantigen prediction capabilities of SNV and INDEL sites. S8. The TMB calculation method was optimized by combining two factors: the different predictive abilities of neoantigens and the immune escape mechanism of tumor patients. S7 uses the number of SNV and INDEL mutations in the TCGA database and the number of predicted neoantigen peptides to construct a multiple linear regression model. The ability to predict neoantigens at SNV and INDEL sites specifically includes: (1) Obtain TCGA mutation data, classify mutations according to TCGA rules, filter mutation data, and retain mutations listed as PASS in the FILTER column; then use vcf2maf software to convert the maf file in the mutation data into a standard vcf file, and sort the mutation data in the vcf file according to chromosome position. (2) Obtain the patient's HLA type I classification data from the TCGA database; (3) After matching the mutation information and HLA typing data of the sample, the tumor neoantigen peptides of the sample are predicted and filtered; (4) Count the number of SNV and INDEL mutations in patients in the TCGA database and the number of predicted tumor neoantigen peptides; construct a multiple linear regression model to calculate the ability of SNV sites and INDEL sites to predict neoantigens; in, The number of neoantigens predicted for patient samples included in the model. To include patient samples in the model for predicting the number of SNV sites for neoantigens, To include patient samples in the model for predicting the number of INDEL sites for neoantigens, Weights for the ability of SNVs to generate neoantigens. Weights for the ability of INDEL to generate neoantigens. This represents the intercept value of the multiple linear regression model; S8 optimizes the TMB calculation method by considering two factors: differences in neoantigen prediction capabilities and the immune escape mechanism in tumor patients. Specifically, this includes: The original TMB calculation formula is as follows: in, The mutation sites ultimately detected in the patient include SNV and INDEL. The size of the genome covered by the detection kit; Considering both the varying predictive abilities of neoantigens and the immune escape mechanisms in cancer patients, the following optimization method for TMB calculation is proposed: (1) Combining the list of neoantigen predicted peptides obtained in S5 and the HLA heterozygosity loss status obtained in S6, the proportion of the number of neoantigens bound by the missing HLA types in the total number of neoantigens bound by all alleles in the patient was calculated, and the ratio was used as the neoantigen failure status caused by HLA heterozygosity loss in the patient. (2) Based on the neoantigen prediction ability of the SNV and INDEL obtained by linear regression, the SNV and INDEL of tumor patients are calculated with different production weights to obtain the TMB value after mutation type correction; (3) Combining the above two steps, the bias in neoantigen assessment caused by immune escape and mutation type is corrected, and the calculation formula is as follows: in, This is the correction factor for TMB based on the immune escape mechanism. This is the corrected tumor mutation burden value. High-quality neoantigen peptides representing patients after screening g* The number of HLA alleles bound, among which hlaloh This is a list of HLA alleles identified as having a loss of heterozygosity in step S6. gi for hlaloh Alleles in the list; total This is the list of all HLA alleles of the patient obtained in step S4. gj for total Alleles in the list, The number of single nucleotide mutations ultimately detected in the patient. The final number of insertion / deletion mutations detected in the patient. The size of the genome covered by the detection kit; S8 also includes optimizations to the TMB calculation method based on the different neoantigen prediction capabilities of SNV and INDEL: Based on the neoantigen predictive ability of SNVs and INDELs obtained from linear regression, the TMB value after mutation type correction was calculated for SNVs and INDELs in tumor patients with different production weights. The calculation formula is as follows: in, The number of single nucleotide mutations ultimately detected in the patient. The final number of insertion / deletion mutations detected in the patient. The size of the genome covered by the detection kit.

2. The method for optimizing molecular marker algorithms based on machine learning and immune escape mechanisms according to claim 1, characterized in that, In S2, the NGS sequencing data of tumor tissue and paired normal tissue are quality controlled and compared to obtain the raw alignment file. The raw alignment file is then preprocessed to obtain the final alignment file, specifically including: The NGS sequencing data of tumor tissue and paired normal tissue were subjected to quality control, and low-quality sequences, low-complexity sequences and adapter sequences in the NGS sequencing data were removed using FastP software. Then, the filtered clean reads were aligned to the human reference genome hg19 using BWA-MEN software to obtain the original alignment file. Preprocessing of the original alignment file: (1) Sorting with samtools software; (2) Marking duplicate sequences with Picard software; (3) Performing local re-alignment and base quality correction with GATK software to obtain the final alignment file.

3. The method for optimizing molecular marker algorithms based on machine learning and immune escape mechanisms according to claim 1, characterized in that, In S3, the final alignment file is subjected to tumor somatic single nucleotide variant and small fragment insertion / deletion detection, and the original variant list obtained is filtered to obtain somatic variant information, specifically including: Based on the final alignment file, MuTect2 software was used to detect single nucleotide variants and small fragment insertions / deletions in tumor somatic cells; The original variant list obtained from the detection is filtered by the following steps to obtain somatic variant information. The filtering steps are as follows: (1) Retain non-silent variants, including missense mutations, nonsense mutations, frameshift mutations, and splice site variants; (2) Retain variant sites with high reliability: the sequence alignment quality of the site is greater than 30, the base sequencing quality is greater than 20, the mutation positive and negative strand readings are unbiased, the variant site depth is greater than 50X, and the mutant allele frequency is greater than 0.03; (3) Remove mutation sites with high mutation frequency and retain sites with a frequency of less than 0.01 in the TCGA database.

4. The method for optimizing molecular marker algorithms based on machine learning and immune escape mechanisms according to claim 1, characterized in that, In S5, based on the somatic cell mutation information filtered in S3, the NetMHCpan deep learning algorithm is used to predict the patient's tumor neoantigen peptides, and the screening of the tumor neoantigen peptides specifically includes: Based on the somatic cell variation information filtered in S3, the NetMHCpan deep learning algorithm is used to predict the tumor neoantigen peptides of patients. The binding affinity inhibition concentrations of site mutants and paired wild-types with the predicted peptides are calculated respectively. Based on the affinity inhibition concentrations, the differential affinity fold for each site is calculated, and then the predicted tumor neoantigen peptides are screened.

Images