For the literature (TCGA) dataset, we extracted MSI statuses for 1,670 available RNAseq samples from the Broad Firehose webpage. These MSI statuses obtained using IHC or PCR profiling were then considered as the gold standards for the assessment of transcriptomic signatures. As only MSI-high tumors are considered for specific therapeutic options, we pooled MSI-low and MSS (microsatellite stable) samples in a single class for further analyses. Totally, we obtained 1,340 MSI-low/MSS and 330 MSI-high profiles. These samples represented CRC, endometrial carcinoma, gastric cancer, uterine cancer, and esophageal cancer (Table 1). This was higher than the samplings used previously to validate Li, Pacinkova and Popovici, and Danaher signatures in the original studies (a total of 1,302, 626, and 689 samples, respectively; Table 1). We checked RNAseq gene signatures in binary classifier mode.

For the experimental group, we profiled gene expression by RNAseq using FFPE tumor tissue blocks for a total of 23 CRC patients. In addition, we also analyzed a control group of 13 non-CRC tumor blocks to assess MSI signature performance on these samples as well. Among them, five patients had cervical cancer, two had breast cancer, two had gastric cancer, two had glioblastoma, one had ovarian cancer, and one had endometrial carcinosarcoma (Supplementary Table S1). In total, the experimental group (n = 36) represented 27 female and nine male patients. The patient age varied from 31 to 84 years; the mean patient age in the experimental group was 60.36 years. More detailed patient information is given in Supplementary Table S1.

We performed RNAseq for each tumor sample and obtained ∼3.75–78.02 million reads uniquely mapped on known human Ensembl genes (genome version GRCh38 and transcriptome annotation GRCh38.89), on the average ∼15.5 million gene-mapped reads per library.

For these samples, “gold standard” MSI statuses were determined by PCR test for seven marker microsatellite loci: BAT25, BAT26, BAT40, NR21, NR24, NR27, and CAT25 that are included in a routinely used clinical panel that requires no healthy tissue control (Suraweera et al., 2002). When there were ≥2 marker loci with detected unstable microsatellite length, these samples were considered MSI-high. Otherwise, the samples were put to the common MSI-low/MSS group. In the experimental group, there were a total of seven MSI-high and 29 MSI-low/MSS samples (Table 1, Supplementary Table S2).

By performing PubMed and Google Scholar literature search with keywords “gene signature,” “gene expression,” “RNA sequencing,” “microsatellite instability,”and “MSI” in March 2021, we extracted 73 hits that were manually processed and returned three recent original publications. These three unrelated research papers authored by Li et al. (2020), Pačínková and Popovici (2019), and Danaher et al. (2019) communicated different gene signatures of MSI status. All these signatures were deduced and initially validated on TCGA CRC samples available at the date of research (Table 1). For all the signatures identified, the initial bioinformatic validation cohorts were smaller than those extracted from TCGA in the current study (Table 1).

The signatures included 15 genes (Li), 25 genes (Pacincova and Popovici), and 14 genes (Danaher) (Figure 1). We compared gene compositions of different signatures and found that they were largely different and shared only one common gene, MLH1, which encodes for mutL homolog 1 that can heterodimerize with mismatch repair endonuclease PMS2 to form MutL alpha, part of the DNA mismatch repair system (Figure 1). Li signature shared four other genes with Danaher signature: EPM2AIP1, RNLS, SMAP1, and TTC30A. These genes encode for EPM2A interacting protein 1, renalase, small ArfGAP 1, and tetratricopeptide repeat domain 30A, respectively. Pacincova and Popovici signature also had two other common genes with Li signature: RPL22L1 and SHROOM4 encode for ribosomal protein L22 like 1 and shroom family member 4, respectively. Pacincova and Popovici signature had no other common genes with the Danaher signature (Figure 1).

The experimental and literature samples were then used to assess the performances of those three signatures. All signature values were calculated as described in the original papers. We created and made publicly available the code for signature calculation at Gitlab: https://gitlab.com/ef.viktor/msi_signatures.

The signatures were validated using TCGA RNAseq datasets for tumor samples annotated by MSI status: CRC (n = 619), endometrial carcinoma (n = 533), gastric cancer (n = 380), uterine carcinosarcoma (n = 55), and esophageal cancer (n = 83) datasets and on the set of experimental CRC RNAseq samples (n = 23) and control experimental dataset for non-CRC cancer samples (n = 13). To assess signature biomarker quality, we used area under the ROC curve (ROC AUC) value as the measure. AUC reflects biomarker robustness and depends on its sensitivity and specificity (Borisov et al., 2020b). It varies between 0.5 and 1, and the typical discrimination threshold is 0.7, where greater values denote high-quality biomarkers and vice versa (Boyd, 1997). AUC is often used for scoring different types of molecular biomarkers in oncology (Liu et al., 2018; Tanioka et al., 2018; Chen et al., 2019; Sorokin et al., 2020a). AUC and 95% confidence intervals were calculated using DeLong’s method implemented in pROC R-package. The entire experimental dataset contained different cancer types; therefore, AUC was calculated only for the CRC subgroup of the experimental samples.

In our analysis, Li MSI signature (Figure 2A) scored AUC = 1.0 for the experimental CRC dataset, AUC = 0.9462 for the TCGA CRC, AUC = 0.9397 for the TCGA uterine corpus endometrial carcinoma (UCEC), AUC = 0.9664 for the TCGA STAD dataset, and AUC = 0.9981 for the TCGA joint dataset of UCS + ESCA samples. Pacincova and Popovici signature (Figure 2B) performed as high as AUC = 0.9412 for the experimental CRC dataset, AUC = 0.9583 for the TCGA CRC dataset, AUC = 0.6946 for the TCGA UCEC, AUC = 0.8827 for the TCGA STAD dataset, and AUC = 0.9515 for the TCGA joint dataset of UCS + ESCA samples. In turn, Danaher signature (Figure 2C) showed AUC = 0.9902 for the experimental CRC dataset, AUC = 0.9396 for the TCGA CRC dataset, AUC = 0.9442 for the TCGA UCEC, AUC = 0.9589 for the TCGA STAD dataset, and AUC = 1 for the TCGA joint dataset of UCS + ESCA samples dataset.

Similar to variations in AUC metrics for the three signatures tested, their extents related differently to the true-positive or true-negative MSI statuses (Figures 3A–C).

In the experimental CRC group, there were 6 MSI-high and 17 MSI-low samples. However, in the experimental control group that included non-CRC cancers, there was only one MSI-high sample for endometrial carcinosarcoma, whereas all other samples were MSI-low (Supplementary Table S2). All three signatures supported the true MSI status of samples in the control group (Figures 3A–C).

Assessment of MSI signatures is summarized in Table 2. It can be seen that Li signature showed the highest AUC in the experimental CRC group, followed by Danaher and Pacincova and Popovici signatures, respectively (Table 2). Also, all three signatures performed accurately on TCGA CRC, esophageal cancer, and uterine carcinosarcoma samples with AUC 0.94-1 and highly overlapping 95% confidence intervals. However, in the endometrial carcinoma (UCEC) cohort of TCGA data, Pacincova and Popovici signature showed low AUC below 0.7 threshold, whereas two other signatures showed AUC of at least 0.94. The latter also showed lower performance for TCGA gastric cancer samples (AUC = 0.88 vs. 0.96–0.97 in the other two signatures).

Thus, we conclude that in our tests, all three signatures were equally effective for the CRC, esophageal cancer, and uterine carcinosarcoma samples, whereas for the endometrial carcinomas and gastric cancer samples, the Danaher and Li signatures were found more effective.

We also separately analyzed only early-stage (stages I, IA, and IB) cancer patients from TCGA. In this case, statistical analysis could be performed only for CRC and gastric cancer groups because there were no early-stage MSI-high patients in the other groups. There were 16/13 MSI-high and 89/42 MSI-low samples in CRC and gastric cancer groups, respectively (Supplementary Figure S1). All three signatures performed accurately on early-stage TCGA CRC with AUC 0.966–0.997 and highly overlapping 95% confidence intervals (Supplementary Figure S2). AUC for Li signature was the highest for predicting MSI status in gastric cancer (AUC = 0.956), followed by Danaher (AUC = 0.934) and Pacincova and Popovici (AUC = 0.919) signatures (Supplementary Figure S2).

In this study, we investigated MSI status-annotated RNAseq profiles for a total of 1,693 cancer samples (one sample per individual patient). Among them, there were 619 literature CRC samples from TCGA cohort, 533 TCGA UCEC samples, 380 TCGA gastric cancer samples, 55 TCGA uterine carcinosarcoma samples, 83 TCGA esophageal cancer samples, and 36 experimental samples profiled by RNA sequencing in this study. TCGA RNAseq samples were extracted from five source datasets: COAD (colon cancer, n = 389) and READ (rectal cancer, n = 230) for “CRC,” UCEC (endometrial carcinoma, n = 533), STAD (gastric cancer, n = 380), UCS (uterine carcinosarcoma, n = 55), and ESCA (esophageal cancer, n = 83). MSI annotated TCGA data were downloaded from https://gdac.broadinstitute.org/.

The experimental dataset included 23 colon cancer samples, five cervical cancer samples, two breast cancer, two gastric cancer samples, two glioblastoma samples, one ovarian cancer sample, and one endometrial carcinosarcoma sample. All experimental specimens were stored in the form of FFPE tissue blocks.

To isolate RNA, 10-µM thick paraffin slices were trimmed from each FFPE tissue block using a microtome. RNA preps were extracted using QIAGEN RNeasy FFPE Kit. RNA 6000 Nano or Qubit RNA Assay kits were used to measure RNA concentration. RNA integrity number was measured using Agilent 2100 bio-Analyzer. For depletion of ribosomal RNA and library construction, KAPA RNA Hyper with rRNA erase kit (HMR only) was used. Different adaptors were used for multiplexing samples in one sequencing run. Library concentrations and quality were measured using Qubit ds DNA HS Assay kit (Life Technologies) and Agilent Tapestation (Agilent). RNA sequencing was done using Illumina NextSeq 550 equipment for single-end sequencing, 50-bp read length, for approximately 30 million (mln) raw reads per sample. Data quality check was done on Illumina SAV. De-multiplexing was performed with the Illumina Bcl2fastq2 v 2.17 program. Sequencing data were deposited in National Center for Biotechnology Information Sequencing Read Archive under accession ID PRJNA744404.

RNAseq FASTQ files were processed with STAR aligner (Dobin et al., 2013) in “GeneCounts” mode with the Ensembl human transcriptome annotation (Build version GRCh38 and transcript annotation GRCh38.89). Ensembl gene IDs were converted to HUGO Gene Nomenclature Committee (HGNC) gene symbols using the Complete HGNC dataset (https://www.genenames.org/, database version from July 13, 2017). Totally, expression levels were established for 36,596 annotated genes with the corresponding HGNC identifiers. Quantile normalization (qnorm python package) was used to normalize gene expression values.

MSI RNAseq signature described by Li et al. (2020) was calculated according to the original paper. This signature defines LYG1, MSH4, and RPL22L1 genes as “plus”-genes and DDX27, EPM2AIP1, HENMT1, MLH1, NHLRC1, NOL4L, RNLS, RTFDC1, SHROOM4, SMAP1, TTC30A, and ZSWIM3 as “minus”-genes. The final score is a sum of log10-transformed normalized gene expression levels with consideration of each gene sign.

MSI RNAseq signature described by Pacincova and Popovici was calculated according to the original paper (Pačínková and Popovici, 2019). This signature defines AGR2, TNNT1, VNN2, TNFSF9, TRIM7, and RPL22L1 genes as “plus”-genes and ACSL6, ARID3A, ASCL2, AXIN2, EPDR1, GGT7, GNG4, KHDRBS3, KRT23, MLH1, NKD1, PLAGL2, PRR15, RUBCNL, SHROOM2, SHROOM4, TFCP2L1, TNNC2, and VAV3 genes as “minus”-genes. The final score is a sum of log10-transformed gene expression levels with consideration of each gene sign.

MSI RNAseq signature described by Danaher et al. (2019)was calculated according to the original paper. This signature includes MLH1, MSH2, MSH6, and PMS2 genes for calculating MMR loss score (MLS). First, a minimal Z-score (Zmin) of log2-transformed gene expressions was found. The final MLS = (Zmin + 1.03)/0.69, where 1.03 and 0.69 are the theoretical expectation and standard deviation of the minimum of four standard normal random variables, respectively.

Hypermutation predictor score was calculated by multiplying log2-transformed expressions of EPM2AIP1, TTC30A, SMAP1, RNLS, WNT11, SFXN1, SREBF1, TYMS, EIF5AL1, and WDR76 genes by coefficients from the table given in the original article. The final hypermutation predictor score is a Z-score of the calculated value. The resulting MSI predictor score was calculated as follows: min ( MLS , 0 ) 2 + max ( HPS , 0 ) 2

The MSI predictor score is further used as a predictor of MSI-high status.

KEGG and GO analyses were performed using the R clusterProfiler package. EnrichKEGG and enrichGO functions were used to implement enrichment analysis. GSEA analysis was performed using the web service http://www.webgestalt.org. The following non-default parameters were selected: KEGG pathways were used as a functional database, and the minimum number of genes for a category was set to 3. We used Benjamini–Hochberg false discovery rate correction method and applied a p-value threshold of 0.05 as a cutoff value for filtering pathways and GO terms.

Genomic DNA was isolated from FFPE tissue sections using the QIAamp DNA FFPE Tissue Kit (Qiagen, Valencia, CA).

We performed MSI analysis using a set of five so-called “main” mononucleotide repeat markers: BAT25, BAT26, NR21, and NR24 selected from the revised Bethesda panel (Suraweera et al., 2002) and NR27 selected from the modified pentaplex panel (Buhard et al., 2006). Two additional mononucleotide repeat markers were also included: BAT40, as it was shown to improve the sensitivity of MSI testing in both CRC and extra-colonic tumors (Hartmann et al., 2002; Pagin et al., 2013) and CAT-25, which was reported to increase the sensitivity for identifying dMSH6 tumors (Takehara et al., 2018).

The primer sequences were taken from previous reports (Hartmann et al., 2002; Suraweera et al., 2002; Buhard et al., 2006; Takehara et al., 2018). The sequences of fluorescently labeled oligonucleotides are listed in Table 3.

The marker DNA products were PCR amplified using the qPCRmix-HS (Evrogen, Russia). PCR was carried out in a 20-μl final volume containing 1× qPCRmix-HS, 2 pmoles of each primer, and approximately 20 ng of DNA template.

The marker sets (1) BAT25, BAT26, NR21, and NR27 and (2) BAT-40 and CAT-25 were co-amplified in one PCR tube per set. The marker NR-24 was amplified in a separate PCR tube.

PCR conditions for the tetraplex and duplex assays consisted of an initial 2-min denaturation step at 94C, followed by 37 cycles at 94°C for 20 s, 54°C for 10 s, and 72°C for 12 s, with a final extension at 72°C for 2 min. Conditions for monoplex reaction differed in annealing temperature: 53°C.

Amplified PCR products were analyzed by capillary electrophoresis performed on ABI prism 3130 × l System (Applied Biosystems, United States). The microsatellite marker lengths were detected by Sequence Scanner software (Applied Biosystems, United States).

The cutoff for MSI status classification was chosen on the basis of the threshold of approximately 40%, according to Umar A. et al. (2004). Tumors with instability at ⩾2 of the five main mononucleotide markers were defined as MSI-H. Samples with instability at one main marker were further tested with the additional markers. Tumors with at least one unstable additional marker were defined as MSI-high. Otherwise, tumors were classified as MSI-low/MSS.