Multivariate Risk Analysis of RAS, BRAF and EGFR Mutations Allelic Frequency and Coexistence as Colorectal Cancer Predictive Biomarkers

Simple Summary The colorectal cancer (CRC) stage and evolution should be described by biomarker profiles. In 60 CRC cases, KRAS, NRAS, BRAF, and EGFR mutations were analyzed by droplet digital PCR (ddPCR). KRAS G12/G13 mutation was present in all patients with variable allelic frequencies. KRAS Q61 mutation was correlated with tumor invasion beyond the subserosa and poor differentiation, both associated with worst prognosis. Tumors with NRAS and BRAF mutations were prevalently localized on the right segment colon. The discovery of the KRAS Q61 role in tumor phenotypes provides the foundation for new therapeutic strategies and perspectives on molecular subtypes classification of CRC. Abstract Background: Biomarker profiles should represent a coherent description of the colorectal cancer (CRC) stage and its predicted evolution. Methods: Using droplet digital PCR, we detected the allelic frequencies (AF) of KRAS, NRAS, BRAF, and EGFR mutations from 60 tumors. We employed a pair-wise association approach to estimate the risk involving AF mutations as outcome variables for clinical data and as predicting variables for tumor-staging. We evaluated correlations between mutations of AFs and also between the mutations and histopathology features (tumor staging, inflammation, differentiation, and invasiveness). Results: KRAS G12/G13 mutations were present in all patients. KRAS Q61 was significantly associated with poor differentiation, high desmoplastic reaction, invasiveness (ypT4), and metastasis (ypM1). NRAS and BRAF were associated with the right-side localization of tumors. Diabetic patients had a higher risk to exhibit NRAS G12/G13 mutations. BRAF and NRAS G12/G13 mutations co-existed in tumors with invasiveness limited to the submucosa. Conclusions: The associations we found and the mutational AF we reported may help to understand disease processes and may be considered as potential CCR biomarker candidates. In addition, we propose representative mutation panels associated with specific clinical and histopathological features of CRC, as a unique opportunity to refine the degree of personalization of CRC treatment.


Introduction
Evidence-based medicine has paved the way for a paradigm in which biomarker profiles are used to create a coherent description of the health status within a personalized medicine [1], with predictive results. As a result, prevention is beginning to play an increasingly important role.
Many metabolic, immunologic, or therapeutic factors regulate the tumor progression by influencing the development of a microenvironment containing cells with heterogenic genetic phenotypes and behaviors [2]. Cellular heterogeneity refers to distinct populations of cells in the same tumor microenvironment, displaying various phenotypical features. On the other hand, genetic heterogeneity is determined in many cases by genomic instability, leading to a high mutation frequency in several genes. This heterogeneity is frequent in cancer patients and is crucial for acquired resistance to therapy, a common cause of relapse [3]. Many patients harbor mutations in different genes and in different exons within the same gene [4]. Thereby, the intratumor heterogeneity raises significant challenges in using molecular prognostic markers to select the patients that might benefit from specific therapies [3]. The accurate characterization of the genomic landscape of colorectal cancers (CRC) could identify the distinctive metastasis signature and increase the life expectancy of CRC patients [5].
Colorectal cancers may develop through three main mechanisms: chromosomal instability (CIN), CpG islands methylation phenotype (CIMP), and microsatellite instability (MSI). According to the Consensus molecular subtypes classification, there are four CRC subtypes and one mixed subtype without any clear designation [6,7]. CMS subtype 1 tumors display a higher percentage of MSI (74%) and CIMP (67%); high BRAF (42%) and lower KRAS (23%) gene mutation status; low somatic copy number alteration (SCNA), immune infiltration, and activation; and worse prognosis, representing 14% of total CRC [6,7]. CMS subtype 2 contributes to 37% of all tumor subtypes. These tumors have high SCNA, are microsatellite stable (MSS); have WNT-and Myc-activated pathways and elevated gene expression of EGFR; and carry mutated p53 gene and KRAS mutations (28%). CMS subtype 3 represents 13% of all tumors, with epithelial characteristics and moderately activated WNT and Myc signaling pathways along with the overexpression of IGBP3. About 68% of tumors harbor mutations in the KRAS gene and only 7% in the BRAF gene, with moderate or low MSI and intermediate CIMP status. CMS subtype 4 exhibits upregulation of genes involved in EMT transition, intense stromal infiltration, and low KRAS (28%) and BRAF (7%) mutation frequency [6,7].
The main predictive biomarkers for the response to anti-EGFR monoclonal antibodiestargeted therapy in metastatic CRC are mutations in the KRAS, NRAS, and BRAF genes [8,9] that explain the low percentage of CRC patients (10-20%) responsive to anti-EGFR monoclonal antibodies single treatment [4]. Kristen-RAS (KRAS) and neuroblastoma-RAS (NRAS) belong to the G protein type called the RAS superfamily. In normal cells, RAS protein is inactive (linked to GDP) and can become activated (linked to GTP) by many cellular receptors (tyrosine kinase receptors, G protein-coupled receptors, and integrin receptors). One such activator of RAS is EGFR, which initiates a signaling cascade [10,11]. As a result, CRC development-and progression-related signaling pathways such as MAPK, PI3K-AKT/mTOR, or Wnt/β catenin are activated [12]. Half of all CRCs harbor KRAS and NRAS activating mutations, often located in codons 12, 13, 59, and 61, that affect the metabolism of cancer cells and drive resistance to commonly used drugs [13][14][15]. Mutations in the BRAF gene are localized in exon 15 and are represented by valine amino acid substitution (V600E, V600K, and V600R). They have been described as a prognostic marker or a predictive factor for resistance to anti-EGFR monoclonal antibodies [4]. A recent study showed that proximal colon tumor localization exhibited a significant correlation with mutations in KRAS and BRAF [8]. Another study found an association between mucinous adenocarcinoma and KRAS mutations, but not with NRAS or BRAF mutations [16]. Moreover, cohort studies that sought out to link demographic data and KRAS mutational status reported contradictory results [17,18] All restricted gDNA samples were analyzed in triplicate, and each ddPCR reaction mixture contained the 1X screening kit assay reagent, which contains a primers-probes optimized mix; wild-type probes were labelled with HEX dye, and the mutant probes were labelled with FAM dye, 1X ddPCR Supermix for probes, and 6 µL of gDNA template (8 ng/µL), adjusted to a final volume of 20 µL with DEPC-treated water. The ddPCR reaction mixture samples were mixed with 70 µL of droplet generator oil for probes (Bio-Rad, Hercules, CA, USA) and partitioned into up to 20,000 droplets using QX200 droplet generator (Bio-Rad, Hercules, CA, USA). Emulsified samples were transferred on 96-well plates, and PCR was performed on a C1000 Touch Thermal Cycler (Bio-Rad, Hercules, CA, USA). The thermal cycling conditions were: 95°C for 10 min, and 40 cycles of 94 • C for 30 s, 55 • C for 1 min, and 98 • C for 10 min with a ramp rate of 2 • C/s, according to the manufacturer's recommendations. After that, the fluorescence of the samples was read using the QX200 droplet reader (Bio-Rad, Hercules, CA, USA), selecting FAM and HEX channels. Every ddPCR run included a negative template control (Wild Type Reference Standard, for each analyzed mutation, used at a concentration similar to the samples, Horizon Discovery, Cambridge, UK) and positive template control (Quantitative Multiplex Reference Standard gDNA, covering the BRAF, KIT, EGFR, KRAS, NRAS, and PIK3CA genes with allelic frequencies between 1-24.5%, 50 ng/µL used at a concentration similar to the samples analyzed, Horizon Discovery, Cambridge, UK) for computing fluorescence thresholds. The data were analyzed with the QuantaSoft Analysis Pro Software v.1.2.10.0 (Bio-Rad, Hercules, CA, USA), providing absolute quantification of target DNA (target copies/µL of reaction). Wells with less than 10,000 accepted droplets were excluded from the analyses. The mutation allele frequency (AF‰, the number of mutant haploid genomes in a total of 1000 haploid genomes) was calculated using the mutant allele concentration in copies/µL (Mut) and the wild-type allele concentration in copies/µL (WT) using the equation: The AF‰ descriptive statistics for all studied mutations are presented in Table S2.
Considering that the copy-number of a fresh tumor sample varies with the amount of background wild-type gDNA by comparing with FFPE tumor tissues that have >85% tumor cells, at least two positive droplets for each investigated mutation in a triplicate had to be present for calling a sample positive for a given mutation [25][26][27]. The threshold was manually set based on positive control samples for each channel, and the threshold for positivity was ≥0.1 mutant copies for 10 3 haploid genomes for all assays. Based on the "Rule of Three (3 positive mutant copies out of 3000 copies)" [28,29], the mutant positive samples were divided into two categories, one that has AF‰ ≥ 1 and the second with AF‰ ≥ 0.1 ranging from 0.1-0.99‰. Due to its ability to accurately quantify mutations with low allelic levels, the ddPCR technique has promising potential to be integrated into medical practice as a sensitive prognostic tool.

Classes of Variables Used in This Study
For this study, we classified statistical variables into the following categories: (I) clinical variables (e.g., weight; BMI; risk comorbidities-diabetes, hypertension; risk behaviors such as smoking); (II) pathology variables: (a) macroscopic (e.g., tumor location, tumor volume, number of examined lymph nodes, invasion of other organs, etc.) and (b) microscopic (e.g., tumor differentiation grade, histopathological phenotype, tumor invasion, histopathological staging); and (III) genetic variables (e.g., number of mutant copies detected per 1000 copies of haploid genomes; etc.).

Data Transformations
Among the variables mentioned above, there were some variables defined over untransformed values, derived from "raw data" such as age, gender, weight, etc. or data such as histopathological phenotype, tumor differentiation grade, and number of mutant copies reported per 1000 copies of haploid genomes. The values of the latter variables were obtained after preliminary processing performed according to the experimental protocols. Another category of variables was those defined by data transformation through different procedures such as the logarithmic transformation of numerical variables or the transformation of some quantitative variables into qualitative variables by defining categories using characteristic position parameters (mean; median; percentiles 25%, 33%, 50%, 67%, and 75%; or limit values observed by visual inspection of numerical data that were grouped into particular categories).

Risk Estimation
The contingency of the ordered qualitative variables was described by calculating the odds ratio (OR) and relative risk (RR), considering as significant correlations for which at least one of these two parameters had values over 1.3. The formulas for these parameters are presented in the first section of the Supplementary Materials, Equations (S3)-(S7).

Mutations Coexistence
Starting from the hypothesis that the studied mutations can influence each other, we estimated the risk that these mutations exist or are absent simultaneously, two by two. In addition, we calculated the risk that the presence of one mutation would expressly exclude the presence of another mutation. The two ways in which we analyzed the coexistence of mutations was defined by Equations (S1) and (S2) in the first section of Supplementary Materials. Table S3 shows the risk estimation and interpretation when applying the model defined by the Equation (S2), and Figure S1 shows the results corresponding to this model.

Predictor and Outcome Variables
We considered all possible correlations between (a) clinical and genetic variables as independent variables on the one hand and (b) histopathological variables and genetic variables as dependent variables, on the other hand. All the values of the histopathological variables were listed in Table S4, and the complete risk estimation of their association with mutations studied was detailed in Table S5.

Statistical Analysis Steps
Statistical analysis using the IBM SPSS Statistics 26 statistical analysis package was performed: normality checks, log transformations for data normalization, comparisons, and correlations. Our objective was to assess the statistically significant correlations between clinical, demographic, and histopathology data and all mutations.

Mutations' Prevalence and the Coexistence of Mutations
All CRC tumors had at least one mutation with an AF‰ ≥ 0.1, and 96.7% of them presented at least one mutation with AF‰ ≥ 1 in any of the four genes investigated; the KRAS G12/G13 was the most prevalent mutation detected in our cohort, followed by NRAS G12/G13 (25%, AF‰ ≥ 1) and KRAS Q61 (21.7%, AF‰ ≥ 1) ( Table 1). In addition, we identified BRAF mutations in 25% of tumors (AF‰ ≥ 0.1), but only 11.7% of them had AF‰ ≥ 1. The EGFR exon 19 deletions were present in only 3 tumor samples in our cohort of 60 patients. Table S2 provides detailed descriptive statistics of AF‰ mutations. When considering AF‰ ≥ 1, we identified 28 tumors with two types of mutations and only 10 tumors with three different mutations (Table 2). Moreover, 73.3% of tumors presented both KRAS G12/G13 and KRAS Q61 mutations (AF‰ ≥ 0.1), and only 33.3% had a coexistence for NRAS G12/G13 and NRAS Q61 mutations (AF‰ ≥ 0.1) (Table 3). Interestingly, at AF‰ ≥ 1, only 13 tumors (21.66%) carried mutations in exon 2 and exon 3 of the KRAS gene. In the case of NRAS, only two tumors (3.33%) presented both NRAS G12/G13 and NRAS Q61. The KRAS G12/G13 and BRAFV600 mutations co-occurred in only eight tumors (13.3%). In addition, we found 19 tumors (31.66%) with simultaneous mutations in two genes and five tumors with concomitant mutations in KRAS, NRAS, and BRAF genes, and one tumor with concurrent mutations in KRAS, NRAS, and EGFR genes with AF‰ ≥ 1. To report the mutations coexistence, we formulated risk estimation parameters OR, RR1, RR2, and RR3, as defined by the Equations (S3)-(S7) in Supplementary Materials and interpreted according to Table S2.
First, we analyzed the coexistence of mutations by calculating the OR, RR1, RR2, and RR3 in terms of the presence (AF‰ ≥ 0.1) or absence (AF‰ = 0) of the mutations. Secondly, we explored the possibility of associations between different levels of the mutations. Thus, for each of the mutations, we segregated the group of patients into subgroups defined by AF‰, as follows: group A (cases with AF‰ = 0 and cases with 0.1 ≤ AF‰ < 1, hereinafter referred to AF‰ < 1), group B (cases with AF‰ = 0 and cases AF‰ ≥ 1), and group C (cases with AF‰ = 0, with AF‰ < 1 and AF‰ ≥ 1). The KRAS G12/G13 mutations were ubiquitous in our cohort and for this mutation we evaluated only the cases included in group C. However, for the KRAS Q61, NRAS Q61, NRAS G12/G13, and BRAF mutations, we distinguished six situations relative to the A, B, and C groups: patients with AF‰ < 1 from group A, patients with AF‰ = 0 from group A, patients with AF‰ ≥ 1 included in group B, patients with AF‰ = 0 from group B, patients with AF‰ < 1 included in group C, and respectively patients with AF‰ ≥ 1 from group C. Consequently, 26 categories defined by the type of mutation and the level of the mutation (AF‰) were obtained. After this stratification of cases, we estimated all possible associations between AF‰ < 1, AF‰ ≥ 1, and AF‰ = 0 for any combination of two mutations. In other words, the number of possible associations was 325 (combinations of 26 categories taken in pairs), presented in Figure S1. Statistically significant associations are shown in Figure 1. Table 3. The contingency table presents the coexistence of genetic mutations: 0 represents the absence of mutation, and 1 the presence of the mutation (AF‰ ≥ 0.1).
The BRAF mutation with an AF‰ ≥ 1 was found to be associated with NRAS Q61 mutation, with an AF‰ < 1 (RR2 = 3.968), and in the absence of these, NRAS Q61 mutations frequencies favored the absence of BRAF mutation with an AF‰ ≥ 1 (RR3 = 2.365) ( Our results challenge the exclusion of the BRAF mutation or any other mutations by the presence of KRAS G12/G13 mutations [17,30]. Moreover, we propose that more consideration should be given to the AF of analyzed genes in the context of mutation co-existence studies. The coexistence of mutations with different AF levels also opens up perspectives for further refining the risk stratification in the target population and personalized therapeutic decisions in patients with AF levels of mutations indicating the disease progression.

Association between Clinical Data and Mutational Status
We calculated OR, RR1, RR2, and RR3 to associate different clinical data i.e., gender, diabetes, BMI groups (normal weight patients /overweight and obese patients), smoking, and age groups as risk factors for mutations (Figures 2 and 3). We estimated these parameters in the following situations: Firstly, we considered the total absence of the mutation (AF‰ = 0) versus the cases corresponding to the mutations with 0.1 ≤ AF‰ < 1.0, and secondly, we considered the absence of the mutation (AF = 0) versus the cases corresponding to the mutations with AF‰ ≥ 1.0. In both situations, the risk estimates remained relatively constant. Thus, risk estimation was performed considering that AF‰ ≥ 0.1 corresponds to the presence of mutation. relatively constant. Thus, risk estimation was performed considering that AF‰ ≥ 0.1 corresponds to the presence of mutation.   . Association between age and mutational status. On the horizontal axis are presented the various threshold ages (a t ) by which the cohort was divided into two groups: 1-patients younger than the threshold ages, and 2-patients older than the threshold ages. These values represent the percentiles 25%, 33%, 50% (or median), 67%, and 75% of the "age" variable. The association was illustrated by the calculated OR (a), RR1 (b), RR2 (c), and RR3 (d).

Morphopathological Association with Mutations Presence
Systematic analysis in terms of the association of genetic mutations with histopathological phenotypes may be the premise of creating individualized treatments for CRC patients, more refined than current therapeutic solutions based on consensus on molecular subtypes.  Table S4.
Thus, histopathological analyses in conjecture with genetic ones can contribute to the revelation of disease current status and the identification of diagnostic and prognostic biomarkers that allow indications of personalized treatment.

Discussion
For the first time in a Romanian cohort, the CRC clinicopathological variables association with the mutational status of KRAS, NRAS, BRAF, and EGFR genes was explored. As these mutations are involved in the CRC development, and are currently considered as parameters to guide post-resection treatment decisions, we believe the findings we report here could contribute to refining post-resection therapy planning. We employed ddPCR, a very sensitive technique that allows absolute target quantification, and enabled us to assess the concomitant presence of low-frequency mutations.
In our cohort, 96.7 % of CRC tumors harbored KRAS G12/G13 mutations and 21.1 % had KRAS Q61 mutations, both with AF‰ ≥ 1. All tumors that bore KRAS Q61 mutations had a KRAS G12/G13 mutations with an AF‰ ≥ 1. NRAS Q61 mutations were present in 15% of CRC tumors, while NRAS G12/G13 mutations were present in 20% with an AF‰ ≥ 1. Notably, only two tumors exhibited both of these mutations at the mentioned AFs. BRAF V600 mutations (AF‰ ≥ 1) were present in 11.7% of the tumors, and all these tumors exhibited KRAS G12/G13 mutations with an AF‰ averaging at 3.2. The percentage of BRAF V600 mutations in our cohort did not differ from that described in previous clinical studies, being associated with the right-side colon localization of the tumor [4,8,31,32]. Other authors reported incidences of KRAS G12/G13 mutations of approximately 30-50% [4,16,31], a lower value than reported by us. The high percentage of mutations in exon 2 of the KRAS gene we reported in our cohort could be explained by the quality of the isolated gDNA and the use of highly sensitive ddPCR. The KRAS G12/G13 kit we used allowed for the detection of all the frequent mutations encountered in exon 2 of the gene (G12A, G12C, G12D, G12R, G12S, G12V, and G13D), and thus any of these mutations in the same tumor sample, can contribute towards classification of the sample as being mutation positive. In addition, cellular and genetic tumor heterogeneity could be another relevant factor [3]. Our results demonstrate that BRAF and KRAS exon 2 mutations were not mutually exclusive and that their concomitant occurrence is not a rare event as claimed by other studies [4,17,33]. Additionally, several studies using paraffinembedded tumor samples and detection technologies less sensitive than ddPCR reported the concomitant presence of KRAS and BRAF mutations, but in a lower percentage (below 4%) [34,35] than that found in our study. Moreover, in another study, BRAF mutational status was only screened in the KRAS wild-type tumor specimens, and for this reason, they found similar relationships with the clinicopathologic features for both mutations [17]. Our data reveals the co-occurrence of KRAS, NRAS, and BRAF mutations (Tables 2 and 3) in CRC, which could explain why only a reduced number of patients (10-20%) respond to single-agent treatments based on anti-EGFR antibodies [36]. The risk of metastasis or progression in CRC requires a rigorous early assessment based on the correlated analyses of clinical factors (gender, BMI, age), genetic mutations, and histopathological characteristics. A recent study of KRAS mutations in exons 2, 3, and 4 in stage I-IV CRC patients concluded that exon 3 mutations predict the worst prognosis, and suggested that mutations of different KRAS exons should be analyzed separately [37]. Our results are in accordance with this study, as we highlight the strong association of the KRAS Q61 mutations with tumors, demonstrating histopathological features with an adverse impact on the disease prognosis ( Figure 4). Thus, KRAS Q61 mutations were the only ones (amongst the mutations studied) associated with the ypT4 stage that were also associated with the absence of any NRAS and BRAF mutations (statistically significant). In addition, the association of KRAS Q61 mutation with poorly differentiated tumors strengthens this aspect (Figure 4, Table S4).
Along with the EGFR mutations, the KRAS Q61 mutations were also linked with invasiveness beyond the serosa, lymphovascular invasiveness, and with the invasion of neighboring structures. The invasion mechanism was strongly correlated to the desmoplastic reaction, which, in our samples, was strongly associated with KRAS Q61 mutations and also with the NRAS Q 61 mutations. Moreover, the KRAS Q61 mutations were associated with moderate or high inflammatory infiltrate in the absence of any of the NRAS mutations. Its absence and that of the NRAS type mutations were associated with ypT2-graded tumors. The presence of KRAS Q61 mutation and the absence of the NRAS type mutations were associated with the ypN2 and ypM1 stages. Thus, our results indicate that stage III-IV tumors, according to the AJCC classification, are strongly associated with the presence of KRAS Q61 mutations.
Our results also reveal unprecedented correlations between clinical or histological data on the one hand, and the NRAS mutations on the other hand. The NRAS Q61 mutations were associated with limited invasiveness up to the submucosal level. Its absence was associated with the mucoid and glandular phenotype ( Figure 4, Table S4). In contrast, another study found an association between mucinous adenocarcinoma and KRAS mutations, but not with NRAS or BRAF mutations [16].
The NRAS type mutations seem to have an impact on tumor localization, with a higher risk towards the right colon. The absence of these mutations was, however, linked to a higher risk probability for left colon tumors. Both these associations are significant. A recent study showed that only tumor location in the right colon exhibited a significant correlation with KRAS and BRAF mutational status [8]. The KRAS G12/G13 mutation was ubiquitous and the NRAS G12/G13 mutation coexisted in the tumor samples we analyzed with the presence of BRAF mutation (AF‰ ≥ 1).
In addition, both mutations were associated with many clinical variables. The NRAS G12/G13 mutations were present in CRC tumors from diabetic patients and were absent in those of non-diabetics, normal-weight patients, and female patients, and also absent in CRC tumors resected from patients under 75 years of age (Figures 2 and 3, Table S4).
As in the case of the NRAS Q61 mutation, the BRAF mutation correlates with tumor invasiveness limited to the submucosa and the right-side colon localization, the latter association being weaker than in the case of NRAS mutations (Figure 4). Our results show that the BRAF mutation is associated with lymphocytes-rich high-grade inflammatory infiltrates that correspond with CMS 1, characterized by immune infiltration and high BRAF mutation percentage [6]. The inflammatory microenvironment is an essential contributor in tumor progression [38]; thus, the association between the levels of specific biomolecules and AF of genetic mutations requires further investigation. In our group, the BRAF mutation was more likely to be present in patients younger than 75 years but older than 60 years and with a BMI > 25.0. In contrast with other studies [39,40], we did not find any association between the BRAF mutation and the presence of mucinous features and poor tumor differentiation, the latter being revealed to strongly associate with KRAS Q61 mutations ( Figure 4). Such quantitative genetic analysis could identify a constellation of specific biomarkers allowing risk stratification of CRC patients, more precise diagnosis, and prognosis prediction through the correlation with relevant histopathological elements. . The filled circles signify the association with the mutation's presence (AF‰ ≥ 0.1), while the empty circles indicate the association with the mutation's absence (AF‰ = 0). In the upper register, we considered that the mutation's presence or absence is correlated with the absence or presence of a risk factor. The significant associations between the KRAS Q61 and NRAS Q61 mutations, on the one hand, and between the NRAS G12/G13 and BRAF mutations, on the other hand, are represented by the black arrows. Thus, the first black arrow shows that the absence of the KRAS Q61 mutation is significantly associated with the presence of the NRAS Q61 mutation. Also, the second black arrow shows that the coexistence of NRAS G12/G13 and BRAF mutations is statistically significant. In the lower register, we resumed the important associations between mutational status and the histopathological characteristics and CRC staging. Integrative diagram of the most significant findings. These associations were evaluated considering the separate presence and absence of each of the studied mutations (labelled in different colors). The filled circles signify the association with the mutation's presence (AF‰ ≥ 0.1), while the empty circles indicate the association with the mutation's absence (AF‰ = 0). In the upper register, we considered that the mutation's presence or absence is correlated with the absence or presence of a risk factor. The significant associations between the KRAS Q61 and NRAS Q61 mutations, on the one hand, and between the NRAS G12/G13 and BRAF mutations, on the other hand, are represented by the black arrows. Thus, the first black arrow shows that the absence of the KRAS Q61 mutation is significantly associated with the presence of the NRAS Q61 mutation. Also, the second black arrow shows that the coexistence of NRAS G12/G13 and BRAF mutations is statistically significant. In the lower register, we resumed the important associations between mutational status and the histopathological characteristics and CRC staging.
A recent study on a Moroccan colon cancer patients cohort showed that KRAS-mutated colon cancers were significantly associated with the female gender, vascular invasion, classical adenocarcinoma phenotype, moderately differentiated tumors, advanced TNM stage (III-IV), left colon tumor localization, and a higher incidence of distant metastases at the time of diagnostic [41]. This study agrees with our findings in the case of KRAS Q61 mutations and their association with advanced TNM stage (III-IV), with KRAS G12/G13 mutation being ubiquitous. In addition, the same study reported a connection between the NRAS-type mutations and less-extensive invasiveness, which is in agreement with our data (Figure 4). In our cohort, the presence of KRAS Q61 and the absence of BRAF-and NRAS-type mutation were associated with T4, albeit, other reports claim that concomitant KRAS-and BRAF-positive mutational status are more prevalent in T3 and T4 tumors [34]. By multivariate non-aprioristic approaches, Isnaldi et al. [4] identified two distinct clinicalmutational profiles. The first profile groups included older patients bearing a BRAF mutation with right-side tumors localization, agreeing with our results and previous works [8,39,42], and the second profile group consisted of younger female patients positive for KRAS and PIK3CA mutations. Our data do not support this latter profile since the only statistically significant association found in our cohort with the female gender was the absence of NRAS G12/13 mutations. In addition, since the KRAS G12/G13 mutations were ubiquitous in our samples, we cannot relate these mutations to any gender. Additionally, we did not find a significant correlation between this mutation and gender in the case of KRAS Q61 mutations either.
In a recent study on a Chinese CRC cohort, the authors concluded that an NRAS mutation is an independent prognostic marker for distant metastasis in stage I-III patients, with shorter metastasis-free intervals than NRAS wild-type patients [33]. In contrast, our data showed that stage III-IV tumors were correlated with the absence of NRAStype mutation (Figure 4 and Table S4). Thus, given the contradictory studies, a range of validated biomarkers, particularly prognostic and predictive markers, are required for the advancement towards personalized cancer treatment [43].

Conclusions
We employed a pair-wise association approach to assess the correlation between several mutations (KRAS Q61, KRAS G12 / G13, NRAS Q61, NRAS G12 / G13, BRAF, and EGFR) and also the associations between the mutations and histopathology features (tumor staging, inflammation, differentiation, and invasiveness). The strongest associations we found and the mutational AF we reported may help to understand disease processes and may be considered as potential CCR biomarker candidates. In addition, we described representative mutation panels associated with specific clinical and histopathological features of CRC.
The KRAS Q61 mutations were associated with most of the invasive features of CRC described by histopathological variables (poor differentiation, microscopic and macroscopic invasiveness, and staging) with consequences on the disease prognosis (ypT4M1N2). The absence of NRAS types of mutation was associated with the same or with other histopathological features with different levels of impact on the aggressiveness of the disease. Our analysis revealed that KRAS Q61 and NRAS mutations have distinct clinical-pathological features, and KRAS G12/G13 mutation with different AFs was ubiquitous in this cohort, probably being essential for the CRC initiation and development.
Thus, our findings suggest refining the CRC consensus molecular subtypes classification by including other mutations such as KRAS Q61 and NRAS-type mutations and the AF levels of CRC-related mutations. Furthermore, based on the AF of the studied mutations, patient cohorts might be organized into different risk groups as per histopathological features. Such risk stratification opens up significant prospects for sensitive technologies such as ddPCR to be used in CCR screening and preventive, personalized treatment. As fresh tissue samples can be easily obtained by routine endoscopic investigations or during resection surgeries, a quantitative mutation analysis offers enormous potential to promote the future development of screening methods. This genetic analysis approach corroborated with histological observations could have a significant potential to indicate progression risk, thus guiding therapeutic indications for more effective treatments and increasing the cancer-free period and overall survival of CRC patients.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/cancers14112792/s1, Figure S1: Risk analysis via OR, RR1, RR2, and RR3 of possible two-by-two associations between the studied mutations, as defined by Table S2;  Table S1: Clinical and demographic data of the studied cohort; Table S2: Descriptive statistics of mutation levels (allelic frequencies) of all studied mutations; Table S3: Expressing associations between mutation levels via risk ratios RR1, RR2, and RR3; Table S4: Models used for OR, RR1, RR2, and RR3 calculation related to mutations status and tumor pathological features; Table S5: Risk estimation calculations using the models presented in Table S3.