Method for optimizing molecular marker algorithm based on machine learning and immune escape mechanism

Abstract

The invention discloses a method for optimizing a molecular marker algorithm based on machine learning and an immune escape mechanism, and the method comprises the steps: obtaining NGS sequencing data of a tumor tissue and a normal tissue, carrying out the quality control and comparison of the NGS sequencing data of the tumor tissue and the NGS sequencing data of the paired normal tissue, obtaining an original comparison file, preprocessing the original comparison file to obtain a final comparison file; performing tumor somatic cell mononucleotide variation and small fragment insertion and deletion detection on the file through software comparison to obtain somatic cell variation information; constructing a multiple linear regression model, and calculating the neoantigen prediction capability of the SNV site and the INDEL site; the TMB calculation method is optimized by combining two factors, namely different neoantigen prediction capabilities and an immune escape mechanism of a tumor patient. According to the invention, the real immune state of the patient is reflected more truly, and the immunotherapy response and prognosis of the patient are prompted better.

Description

Technical Field

[0001] The present invention relates to the field of molecular marker algorithm optimization, and in particular to a method for optimizing a molecular marker algorithm based on machine learning and immune escape mechanism. Background Art

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

[0003] In theory, a higher mutation load can increase the possibility of new antigens, thereby increasing the possibility of immune recognition and killing of tumor cells, thus affecting the effect of immunotherapy. Therefore, the focus of existing TMB-related scientific research or clinical detection products is to accurately detect the patient's mutation status by optimizing sequencing strategies and mutation detection algorithms. The number of high-quality mutations detected is counted, combined with the coverage of the detection product, to evaluate the tumor mutation load of each patient, and then provide appropriate immunotherapy plans for tumor patients.

[0004] The main way for tumor neoantigens produced by mutations to exert immune effects is for HLA class I molecules to transmit neoantigen peptides to CD8+T cells, thereby triggering downstream immune responses. However, in the late stages of cancers such as lung cancer, a large number of tumor patients will have loss of heterozygosity in their HLA genes. The missing HLA type will lose the ability to bind to the original mutant peptides, and not all predicted neoantigen peptides will be presented to CD8+T cells. Therefore, the traditional method only considers the number of neoantigens produced by mutations, and does not consider the patient's HLA heterozygosity loss status, which is not comprehensive enough. Based on the development of existing technologies, we can calibrate biomarkers based on immune escape mechanisms. Summary of the invention

[0005] The purpose of the present invention is to provide a method for optimizing a molecular marker algorithm based on machine learning and immune escape mechanism, which more truly reflects the patient's true immune status and also has a better prompting effect on the patient's immunotherapy response and prognosis.

[0006] To achieve the above object, the present invention provides the following solutions:

[0007] A method for optimizing a molecular marker algorithm based on machine learning and immune escape mechanism comprises the following steps:

[0008] S1. Obtain NGS sequencing data of tumor tissues and normal tissues, including whole genome sequencing data or whole exome capture sequencing data based on second-generation sequencing;

[0009] S2. Perform quality control and comparison on the NGS sequencing data of the tumor tissue and the paired normal tissue to obtain the original comparison file, and pre-process the original comparison file to obtain the final comparison file;

[0010] S3. Detect tumor somatic single nucleotide variations and small fragment insertion and deletion on the final comparison file, and filter the original variation list obtained by the detection to obtain somatic variation information;

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

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

[0013] S6. Based on the final comparison file obtained in S2, assessing the HLA loss of heterozygosity of the tumor sample;

[0014] S7. Count the number of SNV and INDEL mutations in the TCGA database and the number of predicted neoantigen peptides to construct a multivariate linear regression model and calculate the ability to predict neoantigens at SNV and INDEL sites.

[0015] S8. Optimize the TMB calculation method by combining the different predictive abilities of neoantigens and the immune escape mechanisms of tumor patients.

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

[0017] The NGS sequencing data of tumor tissue and paired normal tissue were quality controlled, 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: (1) samtools software was used for sorting; (2) Picard software was used to mark repeated sequences; (3) GATK software was used for local realignment and base quality correction to obtain the final alignment file.

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

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

[0022] The original mutation list obtained by detection was filtered by the following steps to obtain somatic mutation information. The filtering steps are as follows: (1) retain non-silent mutations, including missense mutations, nonsense mutations, frameshift mutations, and splice site mutations; (2) retain highly reliable mutation sites: the sequence alignment quality of the site is greater than 30, the base sequencing quality is greater than 20, there is no preference for positive and negative strand reads of the mutation, the mutation 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 variation information filtered in S3, the NetMHCpan deep learning algorithm is used to predict the patient's tumor neoantigen peptides, and the tumor neoantigen peptides are screened, specifically including:

[0024] Based on the filtered somatic variation information in S3, the NetMHCpan deep learning algorithm was used to predict the patient's tumor neoantigen peptides. The binding affinity inhibition concentrations of the site mutant and paired wild-type with the predicted peptides were calculated respectively. Based on the affinity inhibition concentration, the differential affinity fold of each site was calculated, and then the predicted tumor neoantigen peptides were screened.

[0025] Preferably, in S7, the number of SNV and INDEL mutations in the TCGA database and the number of predicted neoantigen peptides are counted to construct a multivariate linear regression model, and the ability to calculate neoantigen predictions at SNV sites 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 by FILTER; 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 patient's HLA typing data in 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, as well as the number of predicted tumor neoantigen peptides; construct a multivariate linear regression model to calculate the ability of neoantigen prediction at SNV sites and INDEL sites;

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

[0031] Among them, Count Neo The number of neoantigens predicted for the patient samples included in the model, Count SNV Count is the number of SNV sites in the patient samples included in the model for predicting neoantigens. INDEL is the number of INDEL sites in the patient samples included in the model for predicting neoantigens, α1 is the ability weight of SNV to produce neoantigens, α2 is the ability weight of INDEL to produce neoantigens, and β is the intercept value of the multivariate linear regression model.

[0032] Preferably, in S8, the TMB calculation method is optimized by combining the different predictive abilities of neoantigens and the immune escape mechanisms of tumor patients, specifically including:

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

[0034]

[0035] Among them, Count snv+indel The mutation sites finally detected for the patient, including SNV and INDEL, Countbed To provide genome-wide coverage for the test kit;

[0036] Combining the different predictive abilities of neoantigens and the immune escape mechanisms of tumor patients, the TMB calculation method is optimized as follows:

[0037] (1) Combining the predicted peptide list of neoantigens of tumor patients obtained in S5 and the HLA heterozygosity loss status obtained in S6, the ratio of the number of neoantigens bound by the missing HLA type to the number of neoantigens bound by all alleles of the patient is calculated, and the ratio is used as the neoantigen failure status caused by HLA heterozygosity loss of tumor patients;

[0038] (2) Based on the neoantigen prediction ability of SNVs and INDELs obtained by linear regression, the SNVs and INDELs of tumor patients were calculated with different generation weights to obtain the TMB value corrected for the mutation type;

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

[0040]

[0041] Among them, r hla is the correction factor for TMB based on immune escape mechanism, TMB adj is the corrected tumor mutation load value, neo g* Represents the number of high-quality neoantigen peptides after patient screening that bind to g*-type HLA alleles, where hlaloh is the list of HLA alleles determined to be in the heterozygous loss state in step S6, gi is the allele in the hlaloh list; total is the list of all HLA alleles of the patient obtained in step S4, gj is the allele in the total list, Count snv The number of single-base nucleotide mutations finally detected in the patient, Count indel The number of insertion and deletion mutations finally detected in the patient, Count bed The assay kit covers a genome-wide size range.

[0042] Preferably, it is characterized in that S8 also includes optimizing the TMB calculation method based on the different neoantigen prediction capabilities of SNV and INDEL:

[0043] Based on the neoantigen prediction ability of SNV and INDEL obtained by linear regression, the SNV and INDEL of tumor patients were calculated with different generation weights to obtain the TMB value after correction of mutation type. The calculation formula is as follows:

[0044]

[0045] Among them, Count snv The number of single-base nucleotide mutations finally detected in the patient, Count indel The number of insertion and deletion mutations finally detected in the patient, Count bed The assay kit covers a genome-wide size range.

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

[0047] The present invention combines machine learning and other methods to mine tumor neoantigen data from public databases, construct a multivariate linear regression model, calculate the weights of the neoantigen prediction abilities of two different mutation types, SNV and INDEL, and correct the TMB calculation process based on the weights; the present invention evaluates the patient's immune escape status through a cutting-edge algorithm to reflect the patient's true neoantigen status; the algorithm is optimized from two levels, which more realistically reflects the patient's true immune status compared to existing algorithms, and also has a better prompting effect on the patient's immunotherapy response and prognosis. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative labor.

[0049] Figure 1 A flow chart of a method for optimizing a molecular marker algorithm based on machine learning and immune escape mechanism provided by the present invention;

[0050] Figure 2 Forest plots of overall survival for the three algorithms;

[0051] Figure 3a The original TMB of 88 patients with non-small cell lung cancer raw KM curve of overall survival;

[0052] Figure 3b TMB after combined mutation type and immune escape correction in 88 patients with non-small cell lung cancer adj KM curve of overall survival;

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

[0054] Figure 4 Forest plots of progression-free survival for the three algorithms;

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

[0056] Figure 5b TMB after combined mutation type and immune escape correction in 88 patients with non-small cell lung cancer adj KM curve of progression-free survival;

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

[0058] Figure 6a The original TMB of 88 patients raw Analysis of intergroup differences in response to immunotherapy;

[0059] Figure 6b TMB after combined mutation type and immune escape correction for 88 patients adj Analysis of intergroup differences in response to immunotherapy;

[0060] Figure 6c TMB after mutation type correction for 88 patients madj Analysis of inter-group differences in response to immunotherapy. DETAILED DESCRIPTION

[0061] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0062] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is 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 a molecular marker algorithm based on machine learning and immune escape mechanism, comprising the following steps:

[0064] S1. Obtain NGS sequencing data of tumor tissues and normal tissues, including whole genome sequencing data or whole exome capture sequencing data based on second-generation sequencing;

[0065] S2. Perform quality control and comparison on the NGS sequencing data of the tumor tissue and the paired normal tissue to obtain the original comparison file, and pre-process the original comparison file to obtain the final comparison file;

[0066] S3. Use software to detect tumor somatic single nucleotide variations and small fragment insertion and deletion in the final comparison file, and filter the original variation list obtained by the detection to obtain somatic variation information;

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

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

[0069] S6. Based on the final alignment file obtained in S2, the HLA loss of heterozygosity of the tumor sample was assessed using software;

[0070] S7. Count the number of SNV and INDEL mutations in the TCGA database and the number of predicted high-quality neoantigen peptides to construct a multivariate linear regression model and calculate the ability to predict neoantigens at SNV and INDEL sites.

[0071] S8. Optimize the TMB calculation method by combining the different predictive abilities of neoantigens and the immune escape mechanisms of tumor patients.

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

[0073] 1.1 Obtaining raw data

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

[0075] 1.2 Data Preprocessing

[0076] The original sequencing data of tumor tissue and paired normal tissue were quality controlled, mainly using fastp software to remove low-quality sequences, low-complexity sequences and adapter sequences in the original sequencing data; then the filtered clean reads were aligned to the human reference genome hg19 using BWA-MEN software to obtain the original alignment file (BAM). The original BAM file was processed as follows: ① Sorting by samtools software; ② Marking of repeated sequences by Picard software; ③ Local realignment and base quality correction by GATK software to obtain the final version of the BAM file for downstream data analysis.

[0077] 1.3 Gene mutation detection and filtering

[0078] Based on the BAM files of tumor samples and paired normal tissues, MuTect2 software was used to detect tumor somatic single nucleotide variants (SNVs) and small fragment insertion and deletions (INDELs). The original variant list obtained by detection was filtered by the following steps to obtain highly reliable somatic variant information. The filtering steps are as follows: ① Retain non-silent variants, including missense mutations, nonsense mutations, frameshift mutations, and splice site variants; ② Retain highly reliable variant sites, which need to meet the sequence alignment quality of the site greater than 30, the base sequencing quality greater than 20, the mutation positive and negative chain reads have no preference, the variant site depth is greater than 50X, and the mutant allele frequency is greater than 0.03; ③ Remove mutation sites with high population frequency and retain sites with a frequency of less than 0.01 in the population database (ExAC, gnomAD).

[0079] 1.4 Human leukocyte antigen (HLA) typing

[0080] Obtain the HLA typing data of the tumor patient, use the software OptiType to perform HLA typing on the sequencing data of the tumor patient's normal tissue (or blood cells), and obtain the patient's HLA four-digit typing results. The main steps are as follows:

[0081] Based on the HLA DNA reference sequence provided by OptiType, the filtered cleanreads were aligned to the genomic region where HLA was located using razers3 software to obtain sequencing data in the HLA region. The alignment file (BAM) obtained by the alignment was then converted into a sequencing sequence file (fastq) using samtools software. The OptiType analysis software was used to perform HLA sequence analysis on the fastq file obtained above, and then the patient's HLA typing results were obtained.

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

[0083] Based on the above filtered and annotated mutations, the deep learning algorithm NetMHCpan (https: / / services.healthtech.dtu.dk / services / NetMHCpan-4.1 / ) was used to predict the possible tumor neoantigen peptides for each patient, and the mutant type and paired wild type of each mutation site were calculated. Among them, only the predicted peptides with a length of 8 to 11mer were considered, and the binding affinity inhibition concentration (IC50) of the site mutant type and paired wild type with the predicted peptide was calculated respectively. Based on the above affinity inhibition concentration, the differential affinity multiple (DAI) of each site was calculated, and then the predicted peptides were screened, which were divided into the following two cases. When the wild type predicted the neoantigen peptide, ① the mutant affinity inhibition concentration <= 500nM and ② DAI> 1 were used as the screening criteria; when the wild type trait of the site did not predict the neoantigen peptide, the mutant affinity inhibition concentration <= 500nM was used as the screening criteria. Finally, high-quality tumor neoantigen peptides of 6562 patients were obtained for downstream analysis.

[0084]

[0085] 1.6 Determination of HLA loss of heterozygosity

[0086] Based on the BAM files of tumor tissue and paired normal tissue obtained in 1.2, the HLA loss of heterozygosity of tumor samples was evaluated using HLALOH software. The loss of heterozygosity of each patient's HLAA, HLAB, and HLA C gene alleles was determined. If the patient has homozygous genes in the three HLA genes, loss of heterozygosity will not occur. The calculation method is to analyze the mismatch sites of paired alleles in tumor tissue and control tissue and their coverage in the genome, and to detect the probability of imbalance of each HLA gene allele by paired t test. When the P value is less than 0.1, the allele of the HLA gene is considered to have an HLA loss of heterozygosity event.

[0087] 1.7 Differential prediction of tumor neoantigens by SNV and INDEL

[0088] 1.7.1 Obtain TCGA mutation data, classify mutations according to TCGA rules, filter mutation data, and only retain mutations listed as PASS by FILTER. Then use vcf2maf software to convert the maf file to a standard vcf file, and sort the mutation data in the vcf file according to chromosome position.

[0089] 1.7.2 Apply to obtain HLA-I typing data of 7300 patients in TCGA 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 Count the number of SNV and INDEL mutations of 6562 patients in the TCGA database, and the number of predicted high-quality neoantigen peptides. Considering that the ability of tumor mutations to produce neoantigen peptides is mainly related to the combination of mutations and peptides, the neoantigens of different cancer types are not much different, and a small amount of cancer specificity is actually caused by mutations specific to different cancer types. Therefore, a multivariate linear regression model was constructed in 6562 pan-cancer patients to calculate the ability to predict neoantigens at SNV sites and INDEL sites.

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

[0093] Among them, Count Neo Count is the number of neoantigens predicted by the above 1.7.3 step for the patient samples included in the model. SNV Count is the number of SNV sites in the patient samples included in the model for predicting neoantigens. INDEL is the number of INDEL sites in the patient samples included in the model for predicting neoantigens, α1 is the ability weight of SNV to produce neoantigens, α2 is the ability weight of INDEL to produce neoantigens, and β is the intercept value of the multivariate linear regression model, which means that when all other independent variables (Count SNV and Count INDEL ),the dependent variable Count Neo The average baseline value.

[0094] 1.7.5 The determination coefficient R2 of the multivariate linear regression model established above is 0.92 (R2 is the proportion of variation explained by the model in the total variation, ranging between 0 and 1, and a larger value means that the model explains a higher degree of variation), among which α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 for patients receiving immunotherapy

[0096] 1.8.1 Original TMBraw calculation

[0097] Based on the high-quality mutation sites obtained in 1.3 and the coverage of the detection kit, the original TMB is calculated using the current common calculation method. The calculation formula is as follows:

[0098]

[0099] Among them, Count snv+indel The mutation sites finally detected for the patient, including SNV and INDEL, Count bed The assay kit covers a genome-wide size range.

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

[0101] The TMB calculation method in 1.8.1 above is optimized by combining the two factors of different mutation type prediction capabilities and the immune escape mechanism of tumor patients. The optimization method is as follows: ① Combined with the tumor patient neoantigen prediction list obtained in step 1.5 above and the HLA heterozygosity loss status obtained in step 1.6 above, the number of neoantigens bound to the missing HLA type is counted as the proportion of the number of neoantigens bound to all alleles of the patient's HLAA, HLAB and HLA C, and this ratio is used as the neoantigen failure state caused by HLA heterozygosity loss in tumor patients; ② Based on the neoantigen prediction capabilities of SNVs and INDELs obtained by linear regression in 1.7, the SNVs and INDELs of tumor patients are calculated with different generation weights to obtain the TMB value corrected for the mutation type. Combining the above two steps, the neoantigen assessment deviations caused by the two levels of immune escape and mutation type are corrected, and the calculation formula is as follows:

[0102]

[0103] Among them, r hla It is the correction coefficient of TMB based on immune escape mechanism (with HLA heterozygosity loss status as the criterion for immune escape). adj is the corrected tumor mutation load value, neo g* Represents the number of high-quality neoantigen peptides after patient screening that bind to g*-type HLA alleles, where hlaloh is the list of HLA alleles determined to be in the heterozygous loss state in step S6, gi is the allele in the hlaloh list; total is the list of all HLA alleles of the patient obtained in step S4, gj is the allele in the total list, Count snv The number of single-base nucleotide mutations finally detected in the patient, Count indel The number of insertion and deletion mutations finally detected in the patient, Count bed The assay kit covers a genome-wide size range.

[0104] 1.8.3 Optimizing TMBmadj based on mutation type

[0105] In the above 1.8.2 method, since the number of neoantigens needs to be calculated, the patient's whole genome sequencing (WGS) or whole exome capture sequencing (WES) data needs to be provided. Considering that in actual clinical testing, some patients may not be able to provide WGS or WES second-generation sequencing data or HLA typing results, this patent provides an alternative optimization solution.

[0106] The optimization method is as follows: Based on the differentiated weights of SNV and INDEL calculated in step 1.7, the TMB value after mutation type correction is obtained. Although this scheme does not consider the impact of the patient's immune escape mechanism on neoantigens, it corrects the prediction weights of different mutation types from different aspects. Compared with the method in 1.8.1, it can still better reflect the immunotherapy status of tumor patients. The specific calculation formula is as follows:

[0107]

[0108] Among them, Count snv The number of single-base nucleotide mutations finally detected in the patient, Count indel The number of insertion and deletion mutations finally detected in the patient, Count bed The assay kit covers a genome-wide size range.

[0109] In the specific embodiment, real cohort data are used to calculate the original TMB raw , TMB after combined mutation type and immune escape correction adj and TMB corrected by mutation type madj , to illustrate the use of TMB in terms of clinical predictive ability adj and TMB madj Advantages.

[0110] In the embodiment of the present invention, the whole exon capture sequencing (WES) data based on the second generation sequencing is included. In the embodiment of the present invention, double-end sequencing is preferred, and the sequencing strategy PE150 is used. The NGS sequencing data is preprocessed to obtain an intermediate data file.

[0111] 2.1 Sample Information

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

[0113] 2.2 Specific implementation process

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

[0115] Calculate the raw TMB using steps 1.8.1, 1.8.2, and 1.8.3 respectively. raw , TMB after combined mutation type and immune escape correction adj and TMB corrected by mutation type madj , predict the therapeutic efficacy, progression-free survival, and overall survival of immunotherapy samples.

[0116] 2.3 Verification Results

[0117] 2.3.1 Prediction of overall survival (OS) of patients by three algorithms

[0118] Original TMB raw , TMB after combined mutation type and immune escape correction adj and TMB corrected by mutation type madj The high group and the low group were divided by the median value as the threshold. R software was used for univariate Cox regression analysis, and the Kaplan-Meier survival curve of OS was drawn.

[0119] Figure 2 OS forest plots of the three algorithms show TMB adj and TMB madj Compared with the original TMB raw It can better predict the overall survival status of patients. When the median value is used as the threshold for grouping, TMB adj and TMB madj Both indicators were significantly correlated with the patients' overall survival status. Figure 3a The results showed that: in the original TMB raw When the median value is the threshold, the raw TMB raw The overall survival rate of patients in the high group and the low group did not reach a significant difference (HR=1.9, p=0.081). raw It does not reflect the patient's overall survival status;

[0120] Figure 3b The results showed that TMB after combined mutation type and immune escape correction adj The median value is used as the threshold to divide the high group and the low group, TMB adj The overall survival of patients in the low group was significantly worse than that in the TMB adjThe risk of death in the high-risk group was 2.5 times that of the latter group (HR=2.5; p=0.02), indicating that TMB adj It has a predictive effect on the overall survival of patients after immunotherapy;

[0121] Figure 3c The results showed that TMB corrected by mutation type madj The median value is used as the threshold to divide the high group and the low group, TMB madj The overall survival of patients in the low group was significantly worse than that in the TMB madj Patients in the high group (HR=2.3; p=0.033) showed that the TMB after correction of mutation type madj It also has a predictive effect on the overall survival of tumor patients after immunotherapy;

[0122] 2.3.2 Prediction of progression-free survival (PFS) of patients by three algorithms

[0123] The three biomarkers were grouped in the same way as in step 2.3.1 above, and univariate Cox regression analysis was performed using R software, and the Kaplan-Meier survival curve of PFS was drawn.

[0124] Figure 4 The forest plot shows the prognostic results of the three algorithms for progression-free survival. The results show that TMB adj and TMB madj Compared with the original TMB raw It can better predict the patient's progression-free survival status.

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

[0126] Figure 5b The results showed that TMB after combined mutation type and immune escape correction adj The median value is used as the threshold to divide the high group and the low group, TMB adj The progression-free survival of patients in the low group was significantly worse than that in the TMB adj Patients in the high group (HR=2.1; p=0.027) showed that TMB adj It has suggestive significance for the recurrence of patients receiving immunotherapy;

[0127] Figure 5c The results showed that TMB corrected by mutation type madjThe median value is used as the threshold to divide the high group and the low group, TMB madj The overall survival of patients in the low group was significantly worse than that in the TMB madj Patients in the high group (HR=2; p=0.048), indicating that the TMB after correction of mutation type is also indicative of recurrence in patients receiving immunotherapy;

[0128] 2.3.3 Predictive effects of three biomarkers on the efficacy of immunotherapy

[0129] The correlation 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 Demonstrated original TMB in patients with non-small cell lung cancer treated with immunotherapy raw There was no statistically significant difference between the immunotherapy response and non-response groups (p = 0.053); Figure 6b Demonstrated immune escape combined with mutation type correction of TMB mdj There was a significant difference between the immunotherapy response and non-response groups (p = 0.0096). adj Significantly higher than those in the non-responding group; Figure 6c Demonstrated mutation type-corrected TMB madj There was a significant difference between the immunotherapy response and non-response groups (p = 0.038). madj significantly higher than those in the non-responder group.

[0131] The above results fully demonstrate that the TMB corrected by the two correction methods can better predict the patient's immunotherapy efficacy and survival outcomes compared with the original TMB.

[0132] The present invention provides two methods for TMB optimization:

[0133] (1) Based on the patient's WES or WGS sequencing data, a more accurate TMB correction method is provided; (2) Considering the detection cost, when the patient only has mutation data and cannot provide WES or WGS sequencing raw data, a TMB correction method that can replace (1) is provided. The present invention is mainly based on the basic principle of TMB as an immunotherapy biomarker and the immune mechanism involved, optimizes the defects of the existing TMB method, and reflects the patient's true immune status as accurately as possible on the basis that the optimized method can be widely used.

[0134] Through the description of the above implementation methods, those skilled in the art can clearly understand that each implementation method can be implemented by means of software plus a necessary general hardware platform, and of course, it can also be implemented by hardware. Based on this understanding, the above technical solution is essentially or the part that contributes to the prior art can be embodied in the form of a software product, and the computer software product can be stored in a computer-readable storage medium, such as ROM / RAM, a disk, an optical disk, etc., including a number of instructions for a computer device (which can be a personal computer, a server, or a network device, etc.) to execute the methods described in each embodiment or some parts of the embodiments.

[0135] The principles and implementation methods of the present invention are described in this article using specific examples. The description of the above embodiments is only used to help understand the method and core idea of ​​the present invention. At the same time, for those skilled in the art, according to the idea of ​​the present invention, there will be changes in the specific implementation methods and application scope. In summary, the content of this specification should not be understood as limiting the present invention.

Claims

1. A method for optimizing a molecular marker algorithm based on machine learning and immune escape mechanism, characterized in that: The following steps are involved: S1. Obtain NGS sequencing data of tumor tissues and normal tissues, including whole genome sequencing data or whole exome capture sequencing data based on second-generation sequencing; S2. Performing quality control and comparison on the NGS sequencing data of the tumor tissue and the paired normal tissue to obtain an original comparison file, and preprocessing the original comparison file to obtain a final comparison file; S3, performing tumor somatic single nucleotide variation and small fragment insertion and deletion detection on the final comparison file, and filtering the original variation list obtained by the detection to obtain somatic variation information; S4, performing HLA typing on the NGS sequencing data of normal tissues of tumor patients to obtain the patient's HLA four-digit typing results; S5. Based on the somatic variation information filtered in S3, the NetMHCpan deep learning algorithm is used to predict the patient's tumor neoantigen peptides, and the tumor neoantigen peptides are screened; S6. Based on the final comparison file obtained in S2, assessing the HLA loss of heterozygosity of the tumor sample; S7. Count the number of SNV and INDEL mutations in the TCGA database and the number of predicted neoantigen peptides to construct a multivariate linear regression model and calculate the ability to predict neoantigens at SNV and INDEL sites. S8. Optimize the TMB calculation method by combining the different predictive abilities of neoantigens and the immune escape mechanisms of tumor patients.

2. The method for optimizing molecular marker algorithms based on machine learning and immune escape mechanism according to claim 1, characterized in that: In S2, the NGS sequencing data of the tumor tissue and the paired normal tissue are quality controlled and compared to obtain an original comparison file, and the original comparison file is preprocessed to obtain a final comparison file, specifically including: The NGS sequencing data of the tumor tissue and the paired normal tissue were quality controlled, 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; The original alignment file was preprocessed: (1) samtools software was used for sorting; (2) Picard software was used to mark repeated sequences; (3) GATK software was used for local realignment and base quality correction to obtain the final alignment file.

3. The method for optimizing molecular marker algorithms based on machine learning and immune escape mechanism according to claim 1, characterized in that: In S3, the final comparison file is subjected to tumor somatic single nucleotide variation and small fragment insertion and deletion detection, and the original variation list obtained by the detection is filtered to obtain somatic variation information, including: Based on the final alignment file, the MuTect2 software was used to detect tumor somatic single nucleotide variations and small fragment insertions and deletions; The original mutation list obtained by detection was filtered by the following steps to obtain somatic mutation information. The filtering steps are as follows: (1) retain non-silent mutations, including missense mutations, nonsense mutations, frameshift mutations, and splice site mutations; (2) retain highly reliable mutation sites: the sequence alignment quality of the site is greater than 30, the base sequencing quality is greater than 20, there is no preference for positive and negative strand reads of the mutation, the mutation 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 mechanism according to claim 1, characterized in that: 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 tumor neoantigen peptides are screened, specifically including: Based on the somatic variation information filtered in S3, the NetMHCpan deep learning algorithm was used to predict the patient's tumor neoantigen peptides, and the binding affinity inhibition concentrations of the site mutant type and the paired wild type with the predicted peptides were calculated respectively. Based on the affinity inhibition concentration, the differential affinity multiples of each site were calculated, and then the predicted tumor neoantigen peptides were screened.

5. The method for optimizing molecular marker algorithms based on machine learning and immune escape mechanism according to claim 1, characterized in that: S7 counts the number of SNV and INDEL mutations in the TCGA database and the number of predicted neoantigen peptides to build a multivariate linear regression model. The capabilities of calculating neoantigen predictions at SNV and INDEL sites include: (1) Obtain TCGA mutation data, classify mutations according to TCGA rules, filter mutation data, and retain mutations listed as PASS by FILTER; 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 typing data in 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, as well as the number of predicted tumor neoantigen peptides; construct a multivariate linear regression model to calculate the ability of neoantigen prediction at SNV sites and INDEL sites; Count Neo =β+α1*Count SNV +α2*Count INDEL Among them, Count Neo The number of neoantigens predicted for the patient samples included in the model, Count SNV Count is the number of SNV sites in the patient samples included in the model for predicting neoantigens. INDEL is the number of INDEL sites in the patient samples included in the model for predicting neoantigens, α1 is the ability weight of SNV to produce neoantigens, α2 is the ability weight of INDEL to produce neoantigens, and β is the intercept value of the multivariate linear regression model.

6. The method for optimizing molecular marker algorithms based on machine learning and immune escape mechanism according to claim 1, characterized in that: S8 combines the different predictive abilities of neoantigens and the immune escape mechanisms of tumor patients to optimize the TMB calculation method, including: The original TMB calculation formula is as follows: Among them, Count snv+indel The mutation sites finally detected for the patient, including SNV and INDEL, Count bed To provide genome-wide coverage for the test kit; Combining the different predictive abilities of neoantigens and the immune escape mechanisms of tumor patients, the TMB calculation method is optimized as follows: (1) Combining the predicted peptide list of neoantigens of tumor patients obtained in S5 and the HLA heterozygosity loss status obtained in S6, the ratio of the number of neoantigens bound by the missing HLA type to the number of neoantigens bound by all alleles of the patient is calculated, and the ratio is used as the neoantigen failure status caused by HLA heterozygosity loss of tumor patients; (2) Based on the neoantigen prediction ability of the SNV and INDEL obtained by linear regression, the SNV and INDEL of the tumor patient are calculated with different generation weights to obtain the TMB value corrected by the mutation type; (3) Combining the above two steps, the bias in neoantigen assessment caused by immune escape and mutation type is corrected. The calculation formula is as follows: Among them, r hla is the correction factor for TMB based on immune escape mechanism, TMB adj is the corrected tumor mutation load value, neo g* Represents the number of high-quality neoantigen peptides after patient screening that bind to g*-type HLA alleles, where hlaloh is the list of HLA alleles determined to be in the heterozygous loss state in step S6, gi is the allele in the hlaloh list; total is the list of all HLA alleles of the patient obtained in step S4, gj is the allele in the total list, Count snv The number of single-base nucleotide mutations finally detected in the patient, Count indel The number of insertion and deletion mutations finally detected in the patient, Count bed The assay kit covers a genome-wide size range.

7. The method for optimizing molecular marker algorithms based on machine learning and immune escape mechanism according to claim 6, characterized in that: S8 also includes optimization of TMB calculation method based on the different prediction capabilities of neoantigens between SNV and INDEL: Based on the neoantigen prediction ability of the SNV and INDEL obtained by linear regression, the SNV and INDEL of the tumor patient were calculated with different generation weights to obtain the TMB value after correction of the mutation type. The calculation formula is as follows: Among them, Count snv The number of single-base nucleotide mutations finally detected in the patient, Count indel The number of insertion and deletion mutations finally detected in the patient, Count bed The assay kit covers a genome-wide size range.

Images