RU2759205C1 - Calculation of the burden of tumour mutations using tumour rna sequencing data and controlled machine learning - Google Patents

Abstract

FIELD: medicine.

SUBSTANCE: group of inventions relates to medicine, namely to the field of molecular diagnostics, personalised medicine and clinical oncology. Proposed is a method for calculating the tumour mutation burden (TMB) on a million base pairs based on the data on RNA sequencing from a tumour sample as an alternative for sequencing a full exome or a target gene panel. Filtration of mutations based on the data on RNA sequencing, used in software sets for identifying mutations, is supplemented by filtration after identifying and annotating mutations using machine learning with a true set of mutations, produced from tumours and corresponding normal tissues tested by sequencing a full exome.

EFFECT: group of inventions provides an increase in the correlation coefficient for TMB measured by means of RNA sequencing, an increase in the correlation coefficient between the TMBs measured by two different methods from 0.18 to 0.6.

14 cl, 12 dwg, 4 tbl, 2 ex

Description

Scope of invention

The invention relates to the field of molecular diagnostics, personalized medicine and clinical oncology, in particular to the calculation of tumor mutational load (TMB) per million base pairs using tumor RNA sequencing data and supervised machine learning algorithms.

State of the art

As the global cancer threat continues to grow, new treatment strategies such as targeted therapies and immunotherapy are being developed. Many of these treatments are only useful for a small number of patients, highlighting the need to personalize drug prescriptions and find reliable biomarkers of therapy response.

Tumor mutation load (TMB) per million base pairs is a well-known predictive marker of the effectiveness of immune checkpoint inhibitors (Zhu J, et al., Association Between Tumor Mutation Burden (TMB) and Outcomes of Cancer Patients Treated With PD-1 / PD-L1 Inhibitions: A Meta-Analysis. Front Pharmacol. 2019; 10: 673). TMB can be defined as the number of somatic mutations per million base pairs (megabases) in the protein-coding region of the genome, including both substitutions and deletions / insertions, without considering the functional role of mutations (Chalmers ZR, et al., Analysis of 100,000 human cancer genomes reveals the landscape of tumor mutational burden. Genome Med. 2017; 9:34). Strongly mutated tumors are more likely to produce tumor neoantigens and are more "visible" to the immune system, therefore TMB is considered a reliable estimate of tumor burden with neoantigens (Fancello L, Gandini S, Pelicci PG, Mazzarella L. Tumor mutational burden quantification from targeted gene panels: major advancements and challenges. J Immunother Cancer. 2019; 7: 183).

Today, the TMB assessment is commercially available in the form of diagnostic tests for clinical and research purposes (Büttner R, et al. Implementing TMB measurement in clinical practice: considerations on assay requirements. ESMO Open. 2019; 4: e000442). Calculating TMB is part of many genetic test systems. All currently available genetic tests use tumor DNA sequencing data to calculate TMB. Most of them calculate TMB using data from a targeted panel of exome sequencing, less often data from whole exome sequencing (WES) or whole genome sequencing (WGS).

Gene expression profiling using RNA sequencing data is another type of analysis that can aid in therapy prescription (Buzdin A, et al. RNA sequencing for research and diagnostics in clinical oncology. Semin Cancer Biol. 2020 Feb; 60: 311-323). This type of analysis will have the added benefit of including a TMB score. RNA sequencing data contains only expressed genes; therefore, TMB calculated from RNA sequencing data can potentially be a more reliable predictor of the effectiveness of immune checkpoint inhibitors than TMB calculated from WES, WGS, or targeted DNA sequencing of a panel of genes.

Although theoretically possible, in practice, reliable detection of somatic mutations for calculating TMB using tumor-only RNA sequencing data is problematic due to a number of factors, such as tumor heterogeneity, heterogeneous coverage of gene-coding regions, and artifacts that often occur in as a result of fixation of the tumor sample. Machine learning (ML) is currently widely used for various tasks in bioinformatics and, in particular, for its application in clinical oncology (Tkachev V, et al., Flexible Data Trimming Improves Performance of Global Machine Learning Methods in Omics-Based Personalized Oncology Int J Mol Sci. 2020 Jan 22; 21 (3)). There is a need to develop a proven ML algorithm to improve the calculation of TMB based on tumor-only RNA sequencing data.

The essence of the invention

Embodiments of the invention disclosed herein relate to methods, systems, and products for calculating tumor mutational load (TMB) using tumor-only RNA sequencing data from tumor samples. The goal is to develop an efficient approach for calculating TMB per million base pairs (megabase) from RNA sequencing data using somatic mutation filtering with a supervised machine learning (ML) model. This approach is based on a modification of the standard set of programs for identifying mutations in tumor-only RNA sequencing data, with an additional filtering step performed by the ML model. The technical result consists in increasing the accuracy of calculating TMB according to sequencing data of only tumor RNA.

In one embodiment, there is provided a method for computing tumor mutation load (TMB) using sequencing data of tumor-only RNA from a tumor sample, the method comprising the steps of: (a) obtaining RNA sequencing data from a tumor sample, and predicting mutations in RNA using a suite of programs to identify and annotate mutations, whereby a preliminary mutation profile is generated for a tumor sample; (b) applying a pre-trained supervised machine learning (ML) model to calculate the TMB of a tumor sample by adjusting the pre-mutation profile; where pre-training of the ML model includes: I. obtaining a training sample containing a sample of mutations from RNA sequencing data and a second sample of mutations from the corresponding DNA sequencing data taken from the same tumor samples; II. annotating mutations in the first set of mutations as TRUE or FALSE mutations based on the second set of mutations; III. performing supervised ML model training using the first and second mutation samples.

Several modifications or alternative embodiments of the method are possible without departing from the spirit of the invention; some of these alternative implementations are listed below.

In one implementation, the method further comprises the step of validating the ML model after pre-training by testing the pre-trained ML model with a teacher on a second training set. Separation of data samples into training set and test set can be performed using different thresholds in accordance with Y. Xu and R. Goodacre, J Anal Test. 2018; 2 (3): 249-262. A sufficient number of samples should be left in the test sample to reliably evaluate the performance of the model by calculating the correlation between the actual and predicted TMB.

In another embodiment, the tumor sample and training samples are obtained from tumor tissue embedded in paraffin wax and fixed in formalin (FFPE).

In another embodiment, the tumor sample and training kit samples are obtained from freshly frozen tumor tissue.

In another implementation, the pre-mutation profile correction is performed using a gradient boosting algorithm.

In another implementation, the pre-mutation profile correction is performed using a random tree algorithm.

In another embodiment, supervised ML model pre-training further includes the following steps: identifying features that characterize annotated mutations in the first mutation sample; selection of extracted elements and a suitable model for supervised ML; classification of the selected traits as corresponding to a TRUE or FALSE mutation.

In another implementation, supervised ML model pre-training is performed until the metric AUC (area under the ROC curve) of the model parameters is equal to or greater than 90% on cross-validation. Methods for calculating the AUC metric are known in the art, see, for example, Green Dm, S. J. A. (1966) Signal Detection Theory and Psychophysics, New York.

In another embodiment, a computing product is disclosed, wherein the computing product comprises a computer-readable storage medium containing a plurality of instructions for controlling a computing system to perform an operation in one of the aforementioned methods.

In another embodiment of the invention, a system is disclosed for calculating tumor mutational load (TMB) using sequencing data of only tumor RNA from a tumor sample, the system comprising: at least one data carrier configured to store a sequencing data set, including at least measure, RNA sequencing data and DNA sequencing data taken from tumor samples; and at least one processor operably associated with at least one storage medium, the at least one processor configured to (a) obtain RNA sequencing data from a tumor sample and predict mutations in RNA using a set of programs for identification and annotation, resulting in a preliminary mutation profile for a tumor sample; (b) applying a supervised pre-trained machine learning (ML) model to calculate TMB for a tumor sample by adjusting the preliminary mutation profile; where the preliminary training of the ML model includes: I. obtaining a training set containing the first sample of mutations from the RNA sequencing data taken from tumor samples, and a second sample of mutations from the corresponding DNA sequencing data taken from the same tumor samples; II. annotating mutations in the first set of mutations as TRUE or FALSE mutations based on the second set of mutations; III. performing supervised ML model training using the first and second mutation samples.

Elements of any of the disclosed embodiments may be used in combination with each other without limitation. In addition, other functions and advantages of the present invention will become apparent to those skilled in the art upon consideration of the following detailed description and accompanying illustrations.

BRIEF DESCRIPTION OF THE FIGURES

A detailed description is set out with reference to the accompanying figures.

FIG. 1. An example of a set of programs for identifying mutations. Data files are shown in ovals, program steps are shown in rectangles. RNA sequencing (RNAseq) files are on the left panel, total exome sequencing (WES) of tumor and WES of normal tissue are on the right panel.

FIG. 2. TMB calculated using WES data for regions also covered by fresh frozen tissue RNA sequencing data (Y-axis, "isec" represents the intersection of RNA and WES sequencing data) compared to the entire WES dataset (X-axis). The dot shapes correspond to the localizations of the TCGA tissue samples. Disease Codes TCGA - COAD: colon adenocarcinoma; GBM: glioblastoma multiforme; LUAD: lung adenocarcinoma; STAD: gastric adenocarcinoma; SKCM: Skin melanoma. Spearman's correlation coefficient = 0.92 (p-value <2 * 10 ^-16 ), Pearson's correlation coefficient = 0.88 (p-value <2 * 10 ^-16 ).

FIG. 3. TMB calculated from RNA sequencing (RNAseq) data of fresh frozen tissue filtered using standard Variant ID filtering (Y-axis) versus WES data (X-axis). The dot shapes correspond to the localizations of the TCGA tissue samples. TCGA COAD disease codes: colon adenocarcinoma; GBM: glioblastoma multiforme; LUAD: lung adenocarcinoma; STAD: gastric adenocarcinoma; SKCM: Skin melanoma. Spearman's correlation coefficient = 0.06 (p-value = 0.52), Pearson's correlation coefficient = 0.09 (p-value = 0.359).

FIG. 4. TMB calculated from RNA sequencing (RNAseq) data of fresh frozen tissue filtered using standard Variant ID filtering as well as filtering known polymorphisms from the ExAC database (Y-axis) versus WES data (X-axis). The dot shapes correspond to the localizations of the TCGA tissue samples. Disease Codes TCGA - COAD: colon adenocarcinoma; GBM: glioblastoma multiforme; LUAD: lung adenocarcinoma; STAD: gastric adenocarcinoma; SKCM: Skin melanoma. Spearman's correlation coefficient = 0.4 (p-value = 3.73 * 10 ^-5 ), Pearson's correlation coefficient = 0.44 (p-value = 3.09 * 10 ^-6 ).

FIG. 5. TMB calculated from RNA sequencing data (RNAseq) from fresh frozen tissue filtered using XGBoost (Y-axis) versus WES data (X-axis). The dot shapes correspond to the localizations of the TCGA tissue samples. Disease Codes TCGA - COAD: colon adenocarcinoma; GBM: glioblastoma multiforme; LUAD: lung adenocarcinoma; STAD: gastric adenocarcinoma; SKCM: Skin melanoma. Spearman's correlation coefficient = 0.52 (p-value = 0.00709), Pearson's correlation coefficient = 0.67 (p-value = 0.000165).

FIG. 6. TMB calculated from RNA sequencing (RNAseq) data of fresh frozen tissue filtered using random tree method (Y-axis) versus WES data (X-axis). The dot shapes correspond to the localizations of the TCGA tissue samples. Disease Codes TCGA - COAD: colon adenocarcinoma; GBM: glioblastoma multiforme; LUAD: lung adenocarcinoma; STAD: gastric adenocarcinoma; SKCM: Skin melanoma. Spearman's correlation coefficient = 0.82 (p-value <2 * 10 ^-16 ), Pearson's correlation coefficient = 0.82 (p-value <2 * 10 ^-16 ).

FIG. 7. TMB calculated using WES data for regions also covered by waxed tissue RNA sequencing (FFPE) data (Y-axis, "isec" represents the intersection of RNA and WES sequencing data), compared to the entire WES dataset (X-axis). The dot shapes correspond to the localizations of the TCGA tissue samples. Disease Codes TCGA - COAD: colon adenocarcinoma; LUAD: lung adenocarcinoma; BRCA: breast cancer; UCEC: carcinoma of the uterus and uterine corpus; BLCA: bladder cancer; PRAD: prostate adenocarcinoma. Spearman's correlation coefficient = 0.82 (p-value = 3.28 * 10 ^-7 ), Pearson's correlation coefficient = 0.82 (p-value = 3.97 * 10 ^-7 ).

FIG. 8. TMB calculated from waxed tissue RNA sequencing data (FFPE RNAseq) filtered using standard Variant ID filtering (Y-axis) versus WES data (X-axis). The dot shapes correspond to the localizations of the TCGA tissue samples. Disease Codes TCGA - COAD: colon adenocarcinoma; LUAD: lung adenocarcinoma; BRCA: breast cancer; UCEC: carcinoma of the uterus and uterine corpus; BLCA: bladder cancer; PRAD: prostate adenocarcinoma. Spearman's correlation coefficient = 0.18 (p-value = 0.377), Pearson's correlation coefficient = 0.26 (p-value = 0.207).

FIG. 9. TMB calculated from waxed tissue RNA sequencing (FFPE RNAseq) using standard Variant ID filtering as well as filtering known polymorphisms from the ExAC database (Y-axis) versus WES data (X-axis). The dot shapes correspond to the localizations of the TCGA tissue samples. Disease Codes TCGA - COAD: colon adenocarcinoma; LUAD: lung adenocarcinoma; BRCA: breast cancer; UCEC: carcinoma of the uterus and uterine corpus; BLCA: bladder cancer; PRAD: prostate adenocarcinoma. Spearman's correlation coefficient = 0.49 (p-value = 0.0122), Pearson's correlation coefficient = 0.45 (p-value = 0.0211).

FIG. 10. TMB calculated from waxed tissue RNA sequencing data (FFPE RNAseq) filtered using XGBoost (Y-axis) versus WES data (X-axis). The dot shapes correspond to the localizations of the TCGA tissue samples. Disease Codes TCGA - COAD: colon adenocarcinoma; LUAD: lung adenocarcinoma; BRCA: breast cancer; UCEC: carcinoma of the uterus and uterine corpus; BLCA: bladder urothelium cancer; PRAD: prostate adenocarcinoma. Spearman's correlation coefficient = 0.6 (p-value = 0.00141), Pearson's correlation coefficient = 0.74 (p-value = 1.68 * 10 ^-5 ).

FIG. 11. TMB calculated from waxed tissue RNA sequencing data (FFPE RNAseq) filtered using random tree method (Y-axis) versus WES data (X-axis). The dot shapes correspond to the localizations of the TCGA tissue samples. Disease Codes TCGA - COAD: colon adenocarcinoma; LUAD: lung adenocarcinoma; BRCA: breast cancer; UCEC: carcinoma of the uterus and uterine corpus; BLCA: bladder cancer; PRAD: prostate adenocarcinoma. Spearman's correlation coefficient = 0.58 (p-value = 0.00202), Pearson's correlation coefficient = 0.69 (p-value = 8.96 * 10 ^-5 ).

FIG. 12. General diagram of a computing device (120), which provides data processing necessary for the implementation of the claimed methods.

DETAILED DESCRIPTION OF THE INVENTION

As used herein, the terms "includes", "including" and / or "having" as used herein are defined as containing (ie, open language). These terms are not intended to be construed as “consists only of”. The term "other" as used herein is defined as at least a second or more. The term "plurality" as used herein is defined as two or more than two. Throughout this document, references to the words “one implementation”, “specific implementations”, “implementation”, “implementation”, “example” or other similar terms mean that a particular feature, structure, or characteristic described in connection with an embodiment is included in at least one embodiment of the present invention. Thus, the appearance of such phrases in different places of this description does not necessarily refer to the same implementation option. In addition, specific features, structures, or characteristics may be combined in any suitable manner in one or more implementations without limitation. Technical and scientific terms in this description have the standard meanings generally accepted in the scientific and technical literature, unless otherwise indicated.

A somatic mutation is a genetic variation that occurs during oncogenesis. Such mutations are found in tumors, but absent in healthy tissue of the same person.

Germline mutations are variations that have been inherited by a person from his / her parents or have occurred in gametes or early embryonic development.

“Supervised machine learning model” means a machine learning problem that trains a function that maps inputs to outputs based on examples of I / O pairs from labeled training data.

"Mutational load of a tumor per megabase (TMB)" is defined as the number of somatic mutations per megabase of the analyzed coding region, including single nucleotide substitutions, insertions, deletions, excluding the functional effects of mutations.

In some embodiments of the invention, the processing of RNA sequencing data using a set of programs for identifying somatic mutations in a tumor sample analysis mode includes such steps as: 1) aligning the sequencing reads with a reference genome or transcriptome, 2) using a set of programs for identifying mutations without taking into account normal tissue, with or without a panel of normal (healthy) tissue.

In preferred embodiments of the invention, RNA sequencing data from a tumor sample is analyzed without any matched normal (healthy) tissue sample, for example, obtained from blood or adjacent normal tissue from the same person.

To generate a preliminary mutation profile for a tumor sample, in preferred embodiments of the invention, RNA sequencing data is obtained from a tumor sample and the mutation load of the tumor is predicted by algorithmic search and mutation annotation. Algorithmic search should be understood as the use of standard, specialized software tools known to those skilled in the art. Some of them are disclosed below, for example, the GATK-like mutation identification pipeline.

"Normal tissue data panel" means a set of sequenced healthy tissue samples from a group of individuals. This kit is intended to be representative of germline variants prevalent in the general population, and is also used to detect technical artifacts in order to improve the results of mutation identification.

"Corresponding healthy tissue sample" means a healthy tissue sample from the same individual, including, but not limited to, blood or healthy tissue adjacent to a tumor.

Supervised machine learning requires data sampling for training. In some embodiments of the invention, the training data sample contains tumor RNA sequencing samples with annotated mutations. Abstract implies the study of the status of the presence of a mutation in a tumor and corresponding samples of healthy tissue in a separate experiment or a series of experiments to validate true somatic mutations. If a mutation is found in a tumor and not found in a normal specimen, that variant is annotated as a true somatic mutation. If the status of a mutation cannot be estimated, for example, due to low sequencing coverage, such a mutation should be excluded from the training dataset.

As used herein, “pretraining” an ML model means performing training on the ML model prior to actual use to correct the preliminary mutation profile for calculating TMB for a tumor sample.

"Separate experiment" means any experiment or series of experiments necessary to determine the mutational status (somatic or germ line) in a tumor. Such experiments compare tumor data with corresponding healthy tissue samples and may include, but are not limited to, whole genome or exome sequencing, a targeted DNA sequencing panel, and polymorphism detection by microarray DNA hybridization.

Characteristics of mutations identified by an appropriate set of programs for identifying somatic mutations to be used for supervised learning can be obtained from read parameters, alignment parameters, mutation context and position in the genome, as well as from external annotations such as databases of known polymorphisms and functional annotation of mutations, etc.

RNA and DNA sequencing can be performed on nucleic acids isolated from human tissue stored in a variety of conditions: fresh frozen tissue, archived formalin-fixed tissue embedded in paraffin wax blocks (FFPE), tissue stored in RNAlater or other stabilizing solution, for example are the following scientific publications (Suntsova M. et al., Atlas of RNA sequencing profiles for normal human tissues. Sci Data. 2019 Apr 23; 6 (1): 36; Hedegaard J. et al. Next-generation sequencing of RNA and DNA isolated from paired fresh-frozen and formalin-fixed paraffin-embedded samples of human cancer and normal tissue.PLoS One. 2014 May 30; 9 (5): e98187; Choi Y. et al. Optimization of RNA Extraction from Formalin-Fixed Paraffin -Embedded Blocks for Targeted Next-Generation Sequencing. J Breast Cancer. 2017 Dec; 20 (4): 393-399). These methods are included in the description of the invention.

In some embodiments, the supervised learning algorithm parameters may be cross-validated, for example, 5 × or 10 × . After that, the selection of model parameters is performed on the training data set.

In a preferred embodiment, a separate annotated sample set that is not used in the initial training (or fit) of the model can be stored as an additional test sample to establish filtering efficiency. A suitable and validated model can be used to calculate TMB in tumor samples from cancer patients.

It is known that different ML models will base their predictions on the characteristics (signs) of mutations in completely different ways and there are no single solutions. Thus, preliminary supervised learning of the ML model can include various stages. For example, in the case of linear regression, you can find out which features are accepted by the model as important - and the models will assign higher coefficients in absolute value to such features. It is very difficult for neural networks to determine the importance of each characteristic. There are special methods for explaining the behavior of a trained neural network. Typically, pre-training requires only a set of RNA sequencing datasets with correctly labeled mutation states (such as TRUE or FALSE), and then each ML model will identify important features for classification in its own style. In some embodiments of the invention, for all categorical features of the mutation (when the variable takes not numeric values, but, for example, types of mutations A-> T, G-> A, etc.), the one hot encoding method is used to encode the functions.

In some embodiments of the invention, the cross-validation of the parameters of the ML model can be performed as follows. To select the hyperparameters of the ML model, the classical method of dividing the known labeled dataset into three samples is used: training, test, and validation. A smaller portion is usually set aside for use as a test sample. Most are divided into k parts, usually 5 or 10, depending on the amount of data. Then the training data is generated in various ways from k-1 such parts, and the remaining k-th part is transformed into a test sample. Using this cross-validation technique, you can determine the exact hyperparameters. In the last step, the selected model is validated using a test sample.

In some embodiments, methods are provided for calculating tumor mutation load per megabase (TMB) in RNA sequencing samples. The approach is based on the modification of a standard set of programs for identification of mutations on sequencing data of tumor RNA only with an additional stage of filtration after identification. Filtration is performed using an ML model trained on RNA sequencing data, in which the presence of somatic mutations was established using an alternative technology, including the study of the status of mutations in tumors and healthy tissues. To fully understand how an ML model can be trained and implemented, this application will describe the invention in the context of using the XGBoost model (T. Chen, C. Guestrin. XGBoost: a scalable tree boosting system. ArXiv, 1603 (2016). ArXiv: 1603.02754 ) and random tree models (A. Liaw and M. Wiener (2002). Classification and Regression by randomForest. R News 2 (3), 18-22) for filtering RNAseq mutations. Examples of calculating TMB are provided below to illustrate the characteristics of the present invention and should in no way be construed as limiting the scope of the invention. It should be understood that other ML models can be successfully trained and used in various embodiments of the invention.

Examples of such ML models are linear regression, logistic regression, neural networks (Perceptron, Multilayer Perceptron (MLP), backpropagation method, stochastic gradient descent, Hopfield network, radial basis function network (RBFN)), deep learning algorithms (convolutional neural networks (CNN), recurrent neural networks (RNNs), long short term memory networks (LSTMs), auto-encoders, Boltzmann machine (DBM), deep trust networks (DBN)), decision tree methods (classification and trust trees (CART), iterative dichotomizer 3 (ID3), C4.5 and C5.0 (various versions), Chi-square automatic interaction detection (CHAID), Decision Stump, M5 model, conditional decision trees), support vector machine (SVM), vector networks supervised quantization (LVQ), self-organizing Kohonen map (SOM), locally weighted learning (LWL), random forest method, Gaussian processes, stochastic gradient descent, coefficient estimation method linear regression model (LASSO), Least Angle Regression Elastic Network (LARS), ridge regression, Bayesian algorithms (Naive Bayes, Gaussian Naive Bayes, Multinomial Naive Bayes, Averaged Single Dependency Estimators (AODE), Bayesian Trust Networks , Bayesian Network (BN)), k-nearest neighbors, dimensionality reduction algorithms (principal component analysis (PCA), principal component regression (PCR), partial least squares regression (PLSR), Sammon Mapping, multidimensional scaling (MDS), best projection search, linear discriminant analysis (LDA), mixture discriminant analysis (MDA), quadratic discriminant analysis (QDA), flexible discriminant analysis (FDA)), gradient boosting (e.g. GBM, XGBoost, LightGBM, CatBoost), ensemble methods ( eg AdaBoost, Weighted Average, Bootstrapped Aggregation) and others.

Example 1. Calculation of TMB in 119 RNA sequencing samples.

To calculate TMB using RNA sequencing data (RNAseq) and compare them with WES data, 119 cases of five different types of cancer (colon adenocarcinoma, glioblastoma multiforme, lung adenocarcinoma, skin melanoma and gastric adenocarcinoma) from The Cancer Genome Atlas ( TCGA) were analyzed. TCGA is a cancer patient genetic database containing molecular profiling results from over 20,000 cancer patients and corresponding normal specimens, covering 33 types of cancer (Cancer Genome Atlas Research Network, Weinstein JN, Collisson EA, et al., The Cancer Genome Atlas Pan- Cancer analysis project Nat Genet 2013 Oct; 45 (10): 1113-20.). The WES and RNAseq data for each case of identified mutations were analyzed as described below. 19 samples were used to train the ML model, while 100 were used to validate the model.

Table 1. Fresh frozen TCGA biological samples used to train and validate the XGBoost model for filtering RNAseq mutations.

TCGA fabric code Cancer type Localization code Number of patients COAD Colon adenocarcinoma 3L 1 COAD Colon adenocarcinoma 4N 1 COAD Colon adenocarcinoma 4T 1 COAD Colon adenocarcinoma 5M 5 COAD Colon adenocarcinoma A6 23 GBM Glioblastoma 02 2 GBM Glioblastoma 06 17 LUAD Lung adenocarcinoma 05 twenty LUAD Lung adenocarcinoma 44 ten SKCM Melanoma 3N 3 SKCM Melanoma Bf 17 STAD Adenocarcinoma of the stomach B7 4 STAD Adenocarcinoma of the stomach BR 15

Sample data were downloaded from the GDC portal (https://portal.gdc.cancer.gov/). All the resulting bam files were converted to fastq format using the Picard v2.18.17 SamToFastq tool (http://broadinstitute.github.io/picard), followed by alignment and mutation identification.

For the RNAseq data, a GATK-like mutation identification pipeline was used (FIG. 1). GATK refers to the Genome Analysis ToolKit (Van der Auwera GA, Carneiro MO, et al., From FastQ data to high confidence variant calls: the Genome Analysis Toolkit best practices pipeline. Curr Protoc Bioinformatics. 2013; 43: 11.10.1-11.10. 33), developed by the Broad Institute (broadinstitute.org) and contains several modules for processing NGS data. Reads were aligned to version 38 of the human genome using STAR v2.6.1d software in 2-PASS mode (Dobin et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013 Jan; 29 (1): 15-21) ... The non-default parameters were: sjdbOverhang 100, twopass1readsN 10000000, twopassMode Basic. Exon coordinates were taken from Ensembl annotation version 89. Samtools v1.3.1 was used for bam indexing (Li H., Handsaker B., Wysoker A., et al., 1000 Genome Project Data Processing Subgroup. The Sequence alignment / map (SAM) format and SAMtools Bioinformatics 2009 Aug 15; 25 (16): 2078-9). Editing read group information (all reads were assigned to the same group) and marking duplicates was done using Picard (http://broadinstitute.github.io/picard) AddOrReplaceReadGroups and MarkDuplicatesRead, respectively. The GATK v3.8.0 SplitNCigarReads module was used to separate reads that were aligned to the alternative splicing regions, and the ReassignOneMappingQuality was used to ensure compatibility of STAR alignment quality with subsequent analysis. Re-calibration of the base quality index was performed using the GATK v4.beta.1 BaseRecalibrator and ApplyBQSR modules. GATK4 Mutect2 (Cibulskis K, Lawrence MS, et al. Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples. Nat Biotechnol. 2013 Mar; 31 (3): 213-9) in tumor-only ( tumor-only) mode with dbSNP version 146 genetic variants database (Sherry ST, Ward MH, et al., dbSNP: the NCBI database of genetic variation.Nucleic Acids Res. 2001 Jan 1; 29 (1): 308-11) and the Mills database containing standard insertions and deletions found in the 1000 genomes project (Mills RE, Luttig CT, et al., An initial map of insertion and deletion (INDEL) variation in the human genome. Genome Res. 2006 Sep; 16 (9): 1182-90). The norms panel was not used. Variants were identified only in exons (according to the GENCODE annotation (Frankish A, Diekhans M, et al. GENCODE reference annotation for the human and mouse genomes. Nucleic Acids Res. 2019 Jan 8; 47 (D1): D766-D773), only chromosomes 1 -22, X and Y) and the PCR_indel_model parameter was set to "HOSTILE". Variants were filtered using GATK4 FilterMutectCalls (all variants are stored in VCF, only the FILTER field is editable). The effect of mutations on the protein sequence was annotated using the annovar software (Wang K, Li M, Hakonarson H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data.Nucleic Acids Res. 2010 Sep; 38 (16): e164) ... For the ExAC (Lek M. et al. Analysis of protein-coding genetic variation in 60,706 humans. Nature. 2016 Aug 18; 536 (7616): 285-91) annotations, a portion containing no TCGA normal tissues was used, and therefore not including germline variants found in TCGA patients. All three or more allelic sites were excluded from further analysis.

To analyze the WES data, a GATK-like pipeline for the identification of somatic mutations was used (Figure 1). The readings were mapped to version 38 of the human genome using the BWA mem v0.7.17 software. (Li H, Durbin R. Fast and accurate long-read alignment with Burrows-Wheeler transform. Bioinformatics. 2010 Mar 1; 26 (5): 589-95). The non-default options were: -k 15, -r 2. Converting sam files to bam and sorting bam files were done using samtools software (Li H., Handsaker B., Wysoker A., Fennell T., Ruan J., Homer N., Marth G., Abecasis G., Durbin R. and 1000 Genome Project Data Processing Subgroup. The Sequence alignment / map (SAM) format and SAMtools. Bioinformatics. 2009 Aug 15; 25 (16): 2078-9) ... The rest of the preprocessing steps were identical to those for RNAseq data, except for the omitted read quality editing and mapping steps. To identify mutations, GATK4 Mutect2 (Cibulskis K, Lawrence MS et al. Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples. Nat Biotechnol. 2013 Mar; 31 (3): 213-9) was used simultaneously for tumor and corresponding normal samples using the same dbSNP and insertion / deletion databases as for RNAseq. Subsequent post-processing steps included GATK4 FilterMutectCalls filtering and annovar annotation. All three or more allelic sites were excluded from further analysis. To control parallel computing tasks, the GNU parallel software was used (Tange, O. GNU Parallel-the command-line power tool. USENIX. 2011; 36, 42-47).

To process DNA and RNA data, other sets of programs for identification of mutations can also be used, for example: FreeBayes (E. Garrison, Marth G. Haplotype-based variant detection from short-read sequencing (2012) arXiv preprint arXiv: 12073907), SAMtools (H Li. A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics, 27 (21) (2011), pp. 2987-2993), Platypus (A. Rimmer, H. Phan , I. Mathieson, Z. Iqbal, SR Twigg, AO Wilkie, et al., Integrating mapping-, assembly-and haplotype-based approaches for calling variants in clinical sequencing applications Nat Genet, 46 (8) (2014), pp. 912-918), SNVSniffer (Y. Liu, M. Loewer, S. Aluru, Schmidt B. SNVSniffer: an integrated caller for germline and somatic single-nucleotide and indel mutations BMC Syst Biol, 10 (2) (2016), p 47), VarScan2 (DC Koboldt, Q. Zhang, DE Larson, D. Shen, MD McLellan, L. Lin, et al. Varscan 2: somatic mu tation and copy number alteration discovery in cancer by exome sequencing. Genome Res, 22 (3) (2012), pp. 568-576) and others, for example, presented in (C. Xu. A review of somatic single nucleotide variant calling algorithms for next-generation sequencing data. Comput Struct Biotechnol J. 2018 Feb 6; 16: 15-24).

For ML filtering, 32 traits were selected, such as coverage of the reference allele, median base quality, or the number of events in a haplotype. Of these, 24 parameters were taken directly from the VCF Mutect2 output file, and one was obtained from the VCF annotation with the ExAC database. The other 7 parameters were plotted using data from the output VCF Mutect2: 4 boolean parameters indicating whether the mutation is (1) an insertion, (2) a deletion, (3) a C-> T (G-> A) substitution, (4) a substitution C-> A (G-> T) and 3 numerical parameters: (5) coverage, and (6-7) insertion / deletion length.

The training dataset was used to train the XGBoost model (T. Chen, C. Guestrin. XGBoost: a scalable tree boosting system. ArXiv, 1603 (2016). ArXiv: 1603.02754). The model's hyperparameters were selected during a series of randomized grid searches. The following parameters were set: learning_rate, n_estimators, min_child_weight, gamma, subsample, colsample_bytree, max_depth, reg_alpha, reg_lambda; among them, the first two had the greatest weight. Cross-validation was 5-fold and ROC AUC was used as a metric for hyperparameter selection.

The training dataset was also used to train a random forest model (A. Liaw and M. Wiener (2002). Classification and Regression by randomForest. R News 2 (3), 18-22) with 20 trees. Other supervised machine learning (ML) models can be used.

First, the TMB in WES and TMB in the WES regions also covered in RNAseq were compared (Fig. 2). This was done to evaluate the best possible predictive performance of RNAseq. The TMB calculated using WES-based mutations in the areas covered in RNAseq is strongly correlated with the TMB calculated for the full WES dataset (Pearson's correlation coefficient = 0.88, P-value <2e-16).

Filtering performance of RNAseq data was evaluated using Mutect2 software with default settings (Fig. 3, Pearson's correlation coefficient = 0.09, P-value = 0.359) and default settings, combined with filtering germline variants from the non-TCGA subset of the ExAC database (Fig. 4, Pearson correlation coefficient = 0.5, P-value = 3.09e-6). The correlation of TMB (RNAseq with XGBoost filtering) with TMB (WES) was calculated (Fig. 5, Pearson correlation coefficient = 0.67, P-value = 0.000165). The correlation of TMB (RNAseq with random forest filtering) with TMB (WES) was also calculated (Fig. 6, Pearson correlation coefficient = 0.82, P-value = 2e-16). It turned out that filtering using the ML model was significantly superior to naive (standard) filtering in the TMB assessment.

In addition, the predictive efficacy of TMB based on RNAseq was tested using the calculation of the area under the ROC curve (AUC) in the same comparison groups as described above. For this, the most common thresholds for the WES data were used (TMB> 6,> 10 and> 20, respectively (Stenzinger A. et al. Tumor mutational burden standardization initiatives: Recommendations for consistent tumor mutational burden assessment in clinical samples to guide immunotherapy treatment decisions. 2019 Aug; 58 (8): 578-588) and calculated whether TMB from the corresponding RNAseq profile is a predictor of TMB by WES.The ML-based approach was found to be superior to ExAC filtering and naive filtering and can be used for reliable estimation TMB for all of the above thresholds with AUC> 0.86 (Table 2).

Table 2. ROC AUC values for TMB prediction (thresholds:> 6,> 10 and> 20) in WES using RNAseq data.

Filtration method TMB> 6 TMB> 10 TMB> 20 WES mutations covered in RNAseq data 0.977 0.962 0.934 Variant ID Filtering Standard 0.516 0.551 0.660 Standard Variant ID filtering + mutation exclusion from the ExAC exome database (non-TCGA subset) 0.656 0.723 0.79 Oncobox ML filtering using XGBoost method 0.899 0.863 0.959 Oncobox ML filtering using random tree method 0.911 0.885 0.954

Example 2. Calculation of TMB from RNA sequencing data in 50 samples fixed in formalin and embedded in paraffin blocks (FFPE).

To calculate TMB using RNA sequencing data (RNAseq) and compare it with WES data, 50 patients of eight different types of cancer were analyzed (bladder cancer, invasive breast cancer, cervical squamous cell carcinoma and cervical adenocarcinoma, colon adenocarcinoma, renal cell carcinoma, lung adenocarcinoma, prostate adenocarcinoma and uterine endometrial cancer) with RNAseq data available for formalin-fixed paraffin-embedded paraffin-embedded (FFPE) specimens in The Cancer Genome Atlas (TCGA) database (Table 3).

TCGA is a cancer patient genetic database containing molecular profiling results from over 20,000 cancer patients and corresponding normal specimens, covering 33 types of cancer (Cancer Genome Atlas Research Network, Weinstein JN, Collisson EA, et al., The Cancer Genome Atlas Pan- Cancer analysis project Nat Genet 2013 Oct; 45 (10): 1113-20). The WES and RNAseq data for each case of identified somatic mutations for 27 samples were used to train the ML model, while 23 were used to validate the model.

Table 3. TCGA FFPE biological samples used to train and validate the XGBoost model for filtering mutations identified from RNAseq data.

TCGA fabric code Cancer type Localization code Number of patients BLCA Bladder cancer BL 3 BRCA Mammary cancer A7 7 BRCA Mammary cancer AC 2 BRCA Mammary cancer PL 4 CESC Squamous cell carcinoma of the cervix and endocervical adenocarcinoma PN 1 COAD Colon adenocarcinoma A6 13 KIRC Renal cell carcinoma B2 2 LUAD Lung adenocarcinoma 44 12 PRAD Adenocarcinoma of the prostate HC 2 UCEC Endometrial carcinoma of the uterus BK 4

RNAseq and WES analysis as well as ML filtration were performed in the same manner as in Example 1.

In the case of FFPE samples, TMB calculated using WES-based mutations in areas covered by RNAseq strongly correlates with TMB calculated for the full WES dataset (Fig. 7, Pearson correlation coefficient = 0.82, P-value = 3.97e -7).

RNAseq data filtering performance was evaluated using Mutect2 software with default settings. (Cibulskis K, Lawrence MS, Carter SL, Sivachenko A, Jaffe D, Sougnez C, et al. Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples. Nat Biotechnol. 2013 Mar; 31 (3): 213-9) (Fig. 8, Pearson correlation coefficient = 0.26, P-value = 0.207) and with default settings combined with filtering germline variants from the non-TCGA subset of ExAC (Lek M. et al. Analysis of protein-coding genetic variation in 60,706 humans. Nature. 2016 Aug 18; 536 (7616): 285-91) (Fig. 9, Pearson's correlation coefficient = 0.45, P-value = 0.0211). The correlation of TMB (RNAseq with XGBoost filtering) with TMB (WES) was calculated (Fig. 10, Pearson correlation coefficient = 0.74, P-value = 1.68e-5). The correlation of TMB (RNAseq with filtering by the method of random trees) with TMB (WES) was also calculated (Fig. 11, Pearson correlation coefficient = 0.69, P-value = 8.96e-5). Filtering by ML model significantly outperformed naive (standard) filtering in TMB estimation.

The predictive efficacy of RNAseq TMB was then tested using the AUC ROC assay in the same comparison groups as described above. Were used the same thresholds as in example 1 for the WES data (TMB> 6,> 10 and> 20) to assess whether TMB from the corresponding RNAseq profile is a predictor of TMB from WES. The ML approach outperformed ExAC and naive filtering in most cases and could be used to reliably estimate TMB in FFPE samples for all of the above thresholds with AUC ≥ 0.85 (Table 4).

Table 4. ROC AUC values for TMB prediction (thresholds:> 6,> 10 and> 20) in WES using FFPE RNAseq data.

Filtration method TMB> 6 TMB> 10 TMB> 20 WES mutations covered in RNAseq data 0.947 0.983 0.952 Variant ID Filtering Standard 0.601 0.717 0.79 Variant ID Standard Filtering + Variants Exclusion from the Exome Exom Database (non-TCGA subset) 0.707 0.833 0.857 Oncobox ML filtering using XGBoost method 0.88 1 1 Oncobox ML filtering using random tree method 0.85 0.975 0.971

FIG. 12 shows a general diagram of a computing device (120), which provides data processing necessary for the implementation of the claimed methods. Typically, the device (120) contains the following components: one or more processors (121), at least one memory (122), at least one data carrier (123), input / output (I / O) interfaces (124) , data input / output tools (125), network tools (126). The processor (121) performs the basic computational operations required for the operation of the device (120) or the functionality of one or more of its components. The processor (121) executes the necessary machine-readable instructions located in the main memory (122). Memory (122) is usually made in the form of random access memory and contains the necessary program logic that provides the necessary functionality. The data carrier (123) is presented in the form of hard disks or SSD disks, arrays, network storages, flash memory, optical storage devices (CD, DVD, MD, Blue-Ray), etc. The data carrier (123) provides long-term storage of various types of information, for example, the aforementioned RNA and DNA sequencing data, ML algorithms, user IDs, etc.

I / O interfaces (124) are usually standard tools for connecting and operating from the server side, such as USB, RS232, RJ45, LPT, COM, HDMI, PS / 2, Lightning, FireWire, etc. The choice of I / O interfaces (124) depends on the specifics of the device (120), which can be a personal computer, mainframe, server cluster, smartphone, laptop, etc.

It is preferable to use a keyboard as a means for input / output data (125) in any implementation of a system that uses the described methods. The hardware version of the keyboard can be any known version: it can be a built-in keyboard used on a laptop or netbook, or a stand-alone device connected to a desktop computer, server, or other computing device. The connection can be wired, in which the keyboard connecting cable is connected to the PS / 2 or USB port located on the system unit of the desktop computer, or wireless, in which the keyboard communicates, for example, wirelessly; a radio channel with a base station, which, in turn, is directly connected to the system unit, for example, to one of the USB ports. Besides the keyboard, I / O data can also include: joystick, display (touchscreen), projector, touchpad, mouse, trackball, light pen, speakers, microphone, etc.

Network communication tools (126) or network tools include a device that provides network reception and transmission of data, such as an Ethernet card, WLAN / Wi-Fi module, Bluetooth module, BLE module, NFC module, IrDa, RFID module, GSM modem, etc. Using network tools (126), it is possible to organize data exchange via a wired or wireless data channel, for example, WAN, PAN, LAN, LAN, Intranet, Internet, WLAN, WMAN or GSM.

Finally, the components of the device (120) are usually connected via a common data bus (127).

The above description discloses and describes only exemplary embodiments of the present invention. As will be understood by those skilled in the art, the present invention may be embodied in other specific forms without departing from its spirit or essential characteristics. Numerous modifications and variations of the present invention are possible in light of the above. Accordingly, the disclosure of the present invention is intended to illustrate, but not limit, the scope of the invention, which is set forth in the following claims. The disclosure, including any readily discernible variations on the teachings provided herein, partly defines the scope of the claims terminology, so that no subject matter is intended to be publicly disclosed.

All publications, patent applications, patents, and other references mentioned herein are incorporated by reference in their entirety. In addition, the materials, methods, and examples are illustrative only and are not intended to be limiting unless otherwise indicated.

Claims (33)

1. Computer-implemented method for calculating tumor mutation load (TMB) using RNA sequencing data from a tumor sample, including the following steps: (a) obtaining RNA sequencing data from a tumor sample and predicting the mutational load of the tumor using algorithmic search and annotation of mutations, whereby a preliminary mutation profile for the tumor sample is generated; (b) apply a pretrained supervised machine learning (ML) model to calculate TMB for a tumor sample by adjusting the preliminary mutation profile, where pre-training the ML model includes: i. obtaining a training set containing a first set of mutations from RNA sequencing data taken from tumor samples and a second set of mutations from DNA sequencing data taken from the same tumor samples, respectively; ii. annotating mutations in the first set of mutations as TRUE or FALSE mutations based on the second set of mutations; iii. performing ML model training using the first and second sets of mutations. 2. The method of claim 1, further comprising the step of verifying the ML model after pre-training by testing the pre-trained ML model on a second training set. 3. The method according to claim 1, wherein the tumor sample and samples from the training kit are obtained from tumor tissue fixed in formalin and embedded in paraffin. 4. The method of claim 1, wherein the tumor sample and training kit samples are obtained from fresh frozen tumor tissue. 5. The method of claim 1, wherein the correction of the pre-profile of the mutations is performed using a gradient boosting algorithm. 6. The method of claim 1, wherein correcting the preliminary mutation profile is performed using a random forest algorithm. 7. The method of claim 1, wherein pre-training the ML model further comprises the following steps: extracting features that characterize annotated mutations in the first set of mutations; selection of extracted features and training of the ML model; classification of the selected traits as corresponding to a TRUE or FALSE mutation. 8. The method of claim 1, wherein the ML model is pre-trained until the cross-validation AUC of the model parameters is equal to or greater than 90%. 9. A computer-readable medium intended for long-term storage of information, storing computer-readable instructions, which, when executed by a system including at least one computing device, initiate the implementation of a method for calculating tumor mutational load (TMB) using RNA sequencing data from a tumor sample according to claim 1 ... 10. System for calculating tumor mutation load (TMB) using RNA sequencing data from a tumor sample, including: at least one storage medium configured to store a plurality of sequencing data, including at least RNA sequencing data and DNA sequencing data taken from tumor samples; and also at least one processor operatively associated with at least one storage medium, wherein the at least one processor is configured to: (a) obtaining RNA sequencing data from a tumor sample and predicting tumor mutation load by algorithmic search and annotation of mutations, thereby creating a preliminary mutation profile for the tumor sample; (b) applying a pre-trained supervised machine learning (ML) model to calculate TMB for a tumor sample by adjusting the preliminary mutation profile, where pre-training the ML model includes: i. obtaining a training sample containing a first set of mutations from RNA sequencing data taken from tumor samples and a second set of mutations from DNA sequencing data taken from the same tumor samples, respectively; ii. annotating mutations in the first set of mutations as TRUE or FALSE mutations based on the second set of mutations; iii. performing ML model training using the first and second sets of mutations. 11. The system according to claim 10, in which the preliminary training of the ML model additionally contains the following steps: extracting features that characterize annotated mutations in the first set of mutations; selection of extracted features and training of the ML model; classification of the selected traits as corresponding to a TRUE or FALSE mutation. 12. The system of claim 10, wherein correcting the preliminary mutation profile is performed using a gradient boosting algorithm. 13. The system of claim 10, wherein the pre-mutation profile correction is performed using a random forest algorithm. 14. The system of claim 10, wherein the ML model is pre-trained until the cross-validation AUC of the model parameters is equal to or greater than 90%.

Images