Colorectal cancer (CRC) ranks as the third most prevalent cancer around the world, with estimations of approximately 1.9 million new cases and around 935,000 deaths in the year 2022 (Bray et al., 2024). Further estimations suggest an increase to 3.2 million new cases by the year 2040 (Morgan et al., 2023). Even as most of the cases occur among those older than 50 years, there is a significant rise in the rates of younger-onset colorectal cancers (Patel et al., 2022). In Saudi Arabia, CRC ranks first among males at 15.3% and third among females at 9.8% (Elwali et al., 2023). Most patients are present in advanced stages with symptoms of altered bowel habits, abdominal pain, rectal bleeding, and anemia (Dekker et al., 2019). CRC results from molecular changes in the colon or rectum, where there is the formation of a benign polyp that can take 5–15 years to advance into cancer (Huck and Bohl, 2016). The cause of this disease is multifactorial, with interactions between genes and the environment. About 70% of cases are considered sporadic, while 30% are inherited, with 5% being due to syndromes such as Lynch and familial adenomatous polyposis (Jasperson et al., 2010). The genetic basis of CRC has been posed to be consistent with that of the Fearon and Vogelstein model, which describes the progression from a benign adenoma to carcinoma (Fearon and Vogelstein, 1990). This model is thought to account for 70%–90% of CRC cases (Li et al., 2021). A less common way through which CRC arises is the microsatellite instability pathway, due to defects in DNA mismatch repair genes, accounting for 15%–20% of cases (Gupta et al., 2021). Lifestyle factors associated with the development of CRC include obesity, low physical activity, tobacco smoking, consumption of alcohol, and a diet high in processed red meat (Keum and Giovannucci, 2019). Endoscopic testing is considered the gold standard for diagnosing CRC, and surgical removal remains the main form of treatment for the disease (Montminy et al., 2020). In metastatic CRC (mCRC), genetic testing is used to personalize treatment protocols based on actionable somatic mutations on specific genes, namely, BRAF, KRAS, and NRAS, among others (Van Cutsem et al., 2016; Biller and Schrag, 2021). For patients who are not positively testing for such mutations, chemotherapy is administered with combined EGFR inhibitors. The KRAS and NRAS mutation-positive patients are treated with chemotherapy only (FOLFOX or FOLFIRI protocols) (Van Cutsem et al., 2016). Patients with the BRAF V600E mutation slightly benefit from chemotherapy in combination with VEGF inhibitors or a mixture of BRAF and EGFR inhibitors (Santarpia et al., 2012). The determination of somatic mutations in actionable genes (such as BRAF, KRAS, and NRAS) is vital with regards to prognosis and in offering personalized therapies. BRAF, KRAS, and NRAS encode core proteins in the RAS-RAF-MEK-MAPK pathway, which, if altered, can induce aberrant proliferation through dysregulating apoptosis and cell cycle progression from G1 to S phase (Seruca et al., 2009). Importantly, some detection of rare or novel variants in these actionable genes may help in the identification of new biomarkers for better therapeutic outcomes (Mikolajcik et al., 2021). In Saudi Arabia, a retrospective study from two hospitals reported a somatic mutational rate of 2.5% in the BRAF gene and 28.6% in KRAS (Siraj et al., 2014). Moreover, a similar study from another hospital revealed a 2.4% rate in BRAF and 30.1% in KRAS (Beg et al., 2015). Furthermore, a study for patients from the Gulf countries who were treated at two hospitals in the United States showed a mutational rate of 4% in BRAF, 44% in KRAS, and 4% in NRAS (Al-Shamsi et al., 2016). A more recent study from one Saudi hospital has shown slightly different findings for BRAF (0.4%), KRAS (49.6%), and NRAS (2%) (Alharbi et al., 2021). Additionally, a recent review has combined mutational rates from several other Saudi hospitals, showing 0.4%–2.5% rate for BRAF, 28.6%–56% in KRAS, and 2%–2.2% in NRAS (Alfahed, 2023) Our understanding of the mutational landscape of CRC in the Saudi population is lacking, and the percentage of common somatic mutations reported from Saudi hospitals is sparse. Moreover, diagnostic labs may sometimes ignore the rare mutations in CRC tumors, which could affect treatment outcomes. Overall, the genetic characterization of CRC is largely incomplete due to its complex nature. We analyzed CRC samples for somatic mutations in the BRAF, KRAS, and NRAS genes for the TruSight Tumor 15 gene panel (TST 15). The 15 genes covered by the TST 15 panel include AKT1, BRAF, EGFR, ERBB2, FOXL2, GNA11, GNAQ, KIT, KRAS, MET, NRAS, PDGFRA, PIK3CA, RET, and TP53 genes, which are frequently altered in solid tumors (Kopetz et al., 2015; Yan et al., 2015; Van Cutsem et al., 2016; Voskoboynik et al., 2016). Therefore, to understand the pathogenicity and clinical relevance of both known and novel variants, an extensive computational analysis was performed. We have combined raw DNA sequence analysis, variant clinical interpretation, pathogenicity prediction, functional domain mapping, secondary structure and stability analyses, and 3D structure superimposition. In addition, molecular docking was performed on Encorafenib against both the BRAF V600E variant protein and its wild-type counterpart. This integrated approach aims to provide insight not only into the novel biomarkers but also to extend the knowledge of the mutational spectrum of CRC in the Saudi population, helping to further refine precision oncology interventions.

CRC is a challenging disease. Twenty percent of patients present at a metastatic stage and 50% of localized cases progress to mCRC (Dekker et al., 2019). Despite the genetic heterogeneity of CRC tumors, targeted treatments remain limited. BRAF, KRAS, and NRAS express vital cascade signaling proteins that act in the RAS-RAF-MEK-MAPK pathway, ultimately influencing the cell cycle G1 checkpoint. While cell commitment to DNA synthesis is largely controlled by retinoblastoma tumor suppressor protein (Rb), mutations in these signaling proteins can promote cell cycle progression from G1 to S phase by increasing expression of cyclins and cyclin-dependent kinases that are known to nullify the effects of Rb. Such mutations in these actionable genes disrupt normal cell cycle regulation, leading to uncontrolled proliferation and carcinogenesis (Sherr, 1995; Smalley et al., 2008; Puyol et al., 2010). In this study, we collected 86 samples from KAUH hospital: 40 (46.5%) were positive and the rest were negative for mutations in these actionable genes. Fourteen variants were identified through molecular screening and raw DNA sequence analysis: 1 in BRAF, 11 in KRAS, and 2 in NRAS, including NRAS R102Q, previously unreported in CRC data. BRAF mutations, particularly the V600E variant, occur in 66% of melanomas and less frequently in other cancers (OMIM:164757). Mutations in BRAF codon 600, KRAS codons 12 and 13, and NRAS codon 61 confer resistance to anti-EGFR treatments, while BRAF V600E mutations indicate aggressive disease but potential benefit from VEGF inhibitors (Ciardiello et al., 2022). This variant disrupts kinase activity, increasing signaling independently of RAS, with a 138-fold rise in oncogenic activity (Oikonomou et al., 2014). It is a key biomarker in thyroid (Koh et al., 2013), melanoma (Cheng et al., 2018), and CRC (Ciardiello et al., 2022). BRAF mutations are present in 8%–12% of European mCRC patients and associated with poor survival (Van Cutsem et al., 2016; Tabernero et al., 2022). In this study, CancerVar identified this variant as Tier I (Table 1), consistent with ClinVar (RCV001030023.12). The frequency of BRAF mutations in our data is 2.32% (Table 1), aligning with Saudi data (0.4%–4%) (Siraj et al., 2014; Beg et al., 2015; Al-Shamsi et al., 2016; Alharbi et al., 2021; Alfahed, 2023).

The BRAF gene, located at chromosome 7q34, consists of 18 exons encoding a 766-amino acid protein (Huang et al., 2013). BRAF is a serine/threonine kinase that performs several functions in the RAS-RAF-MEK-MAPK pathway (Rauch et al., 2011). Under physiological conditions, growth factors activate RAS proteins at the cell membrane, which then activate RAF proteins, thus amplifying the signal. Balanced regulatory signals autoinhibit RAF as a means of turning off growth signals when not needed. BRAF V600E, a type 1 mutation, functions as a monomer, increasing signaling without RAS dimerization and resisting regulated inhibition (Ciombor et al., 2022). BRAF protein has three conserved domains: CR1, CR2, and CR3. CR1 (residues 120–280), near the N-terminus, includes cysteine-rich and RAS-binding (residues 155–227) subdomains. The cysteine-rich subdomain controls kinase activity by autoinhibition (Cutler et al., 1998; Zaman et al., 2019), while the RAS-binding domain upregulates BRAF at the plasma membrane (Freeman et al., 2013). CR2, a serine/threonine-rich region, links CR1 and CR3 and contains a 14-3-3 scaffold binding site. CR3 (residues 457–717) (Zaman et al., 2019), near the C-terminus, is a kinase-binding domain that activates upon phosphorylation, comprising a glycine-rich ATP-phosphate binding lobe (residues 462–469) and a longer C-terminal lobe (residues 593–623) (Śmiech et al., 2020; Maloney et al., 2021). The most common mutation, V600E, is situated in the CR3 region of BRAF, disrupting its autoinhibition mechanism, hence contributing to oncogenic activity due to destabilization within the inactive form of BRAF (Park et al., 2019).

In our secondary structure analysis, we predicted that BRAF V600E falls into a loop structure identified by Park et al. (2019), as the inhibitory loop (Park et al., 2019). Loops are one important determinant in the stability and functions of proteins. Superimposing BRAF V600E on its wild-type counterpart gave a high RMSD value of 1.7216; DUET thermodynamic stability analysis also predicted that this mutation would highly be destabilizing, with a value of −2.118 kcal/mol. The results predict that this missense variant will have a profound impact on the structure of the protein and severely disrupt normal function. Computational binding assays for Encorafenib with BRAF wild-type and its V600E variant in Figure 3 showed differential binding behavior. The binding energy of the wild-type-Encorafenib complex was −5.915 kcal/mol, displaying two hydrogen interactions at D742 and A762. By contrast, V600E-Encorafenib had only one hydrogen bond interaction at Y746 shown with a binding energy of −3.509 kcal/mol for the conformational pose. With this low binding energy score with two hydrogen bond interactions in the case of a wild-type-Encorafenib complex, it indicates that the wild-type protein binds with the inhibitor very tightly and, hence, is effective under normal conditions. The V600E-Encorafenib has a higher binding energy score and, therefore, a weaker binding affinity, though it is still predicted to be able to inhibit its target protein. Both complexes interacted in proximity to the residues expected to interact with the inhibitors near the C-tail inhibition binding sites, F743 and A749, respectively (Martinez Fiesco et al., 2022), and both complexes interacted nearby. Our results indicate that Encorafenib may have a reduced efficacy in V600E variant proteins. The lower inhibition ability of Encorafenib in patients with the V600E variant could be a factor explaining why these patients may not respond to monotherapy, and often require a combination with other drugs such as Cetuximab (an EGFR inhibitor), and Binimetinib a (MEK inhibitor) to gain benefit from the treatment. It is because blocking the BRAF activity causes an over activation of EGFR, the latter must also be inhibited to gain a better response (Kopetz et al., 2019). Therefore, Encorafenib might be better at targeting wildtype BRAF proteins whether blocking this member of the pathway is needed in designing a new drug regimen for patients without the V600E variant. KRAS and NRAS are gene homologs belonging to the RAS superfamily, which play crucial roles in the RAS-RAF-MEK-MAPK signaling pathways (Fernández-Medarde and Santos, 2011) and PI3K-AKT-mTOR pathways (Castellano and Downward, 2011). Approximately 19% of cancer patients develop a RAS mutation, including approximately 3.4 million new cases yearly that are positive for RAS. Moreover, 75% of mutated RAS isoforms are KRAS (Prior et al., 2020). KRAS mutations occur in 11.6% of all cancer types, as per TCGA data. They are mostly prevalent in pancreatic ductal adenocarcinoma (PDAC) (81.72%), CRC (37.97%), and non-small cell lung cancer (NSCLC) (21.20%) (Yang et al., 2023). In mCRC patients, ∼60% have KRAS mutations (Peeters et al., 2015), primarily in codons 12, 13, 61, 117, and 146 (Serebriiskii et al., 2019). Over 90% of KRAS mutations occur at codons 12 and 13, part of the p-loop of the GTP binding domain, with codon 12 being more associated with oncogenesis and higher mortality than codon 13 (Colussi et al., 2013). The most prevalent KRAS mutations in CRC are G12D (28.04%), followed by G12V (18.50%) and G13D (18.10%), according to http://cBioPortal.org (Yang et al., 2023). Among the 86 samples in our study, 18.60% carried the G12D mutation, 5.81% carried the G12V mutation, and 3.48% carried the G13D, G12A, and G138E mutations. Additionally, 2.32% of samples had K117N and A146T mutations, while 1.16% of samples had G12R, G12S, E49K, or Q61H mutations. Except for E49K and G138E, all identified KRAS mutations are within hotspots. KRAS mutations have been reported in Saudi Arabia in the past at a frequency of 28.6%–56%, with G12D and G12V being the most prevalent variants (Siraj et al., 2014; Beg et al., 2015; Al-Shamsi et al., 2016; Alfahed, 2023). In comparison, our data demonstrated ∼10% lower incidences of G12D than international and regional studies, which may be related to the limited sample size. Surprisingly, G138E, a rare variant that was only reported in COSMIC twice from large intestine samples, was detected in three cases and thus could be more common in Saudi or Arab CRC patients. One patient carried a rare KRAS mutation (E49K), which was concomitantly found with another KRAS mutation (G138E). It was reported only twice so far from large intestine samples in COSMIC. CancerVar has classified G12D as Tier I, strong clinical significance, and nine variants as Tier II, potential clinical significance: G12A; G12R; G12S; G12V; G13D; Q61H; K117N; G138E; and A146T. E49K was classified as Tier III with unknown clinical significance. The clinical significance of these variants aligns with expectations, except for G138E and E49K, which lack sufficient clinical data.

In the cell, KRAS is expressed in two isoforms (KRAS-4B and KRAS-4A) due to alternative splicing of the exon 4. The term KRAS refers to the KRAS-4B transcript, as it is more highly expressed than KRAS-4A (Nuevo-Tapioles and Philips, 2022). The KRAS gene, located at chromosome 12p11.1-12p12.1, encodes a protein of 188 amino acids, divided into five exons (Huang et al., 2021). KRAS is a cascade signaling protein that is activated by GTP and deactivated by GDP, a process that needs to be balanced for normal protein function (Zhu et al., 2021). KRAS and NRAS share two main conserved domains. The GTP-binding (G) catalytic domain (residues 1–165) is highly conserved and identical in both RAS forms. It comprises two lobes: the effector lobe (residues 1–86) contains switch 1 (S1) (residues 30–38), switch 2 (S2) (residues 60–76), and the phosphate-binding region (P-loop) (residues 10–17), all essential for GTP binding. The allosteric lobe (residues 87–166) is less conserved, with 70%–80% similarity. The hypervariable region (HVR) (residues 167–188) at the C-terminus is the least conserved and crucial for post-translational modifications for membrane anchoring (Jancík et al., 2010; Zeitouni et al., 2016; Pantsar, 2020; Han et al., 2021; Nuevo-Tapioles and Philips, 2022). Our identified mutations in codons 12 and 13 (G12A, G12D, G12R, G12S, G12V, and G13D) are in the P-loop region of the G domain. These mutations disturb the active site’s geometry, preventing GTP hydrolysis and locking the RAS protein in its active state (Malumbres and Barbacid, 2003). The rare E49K and the Q61H mutations, located in exon 3, affect the G domain; Q61H is in the S2 region, potentially disrupting GTPase activity by blocking GTP transition to GDP. Mutations K117N, G138E, and A146T are in exon 5, near the HVR. K117N and A146T, found in hotspot codons, are known to predict resistance to Cetuximab plus Irinotecan (Loupakis et al., 2009). Codon 146 mutations occur in ∼4% of CRCs and are more common than codon 61 mutations (Edkins et al., 2006). Functional studies on Q61H and K117N show increased proliferation and high levels of GTP-bound KRAS, indicating they affect GTPase activity by maintaining the active form (Stolze et al., 2015). G138E is a rare mutation located near the HVR with no current functional data. Secondary structure analysis of KRAS variants revealed all variants were in loop regions except E49K, which is on a ß-pleated sheet. No obvious changes were seen in secondary structure elements, except an increased loop size in the G138E mutant (Supplementary Figure 1). Loop regions are important for protein stability, and mutations can disturb this stability. E49K showed no obvious effects. Superimposing KRAS mutant and wild-type forms revealed significant alterations at affected residues for all variants (Figure 3), with deviations ranging from 1.1634 to 1.5958, except for Q61H and K117N, which showed higher deviations (2.3574 and 2.789, respectively). These deviations likely affect the binding site geometry, particularly at codon 12, where sidechain addition disturbs GTP hydrolysis complex binding. Q61H and K117N significantly disrupt the amino acid structure, altering overall protein function and GTPase activity. Protein stability analysis indicated seven KRAS variants to destabilize the structure (G12D, G12R, G12S, G13D, K117N, G138E, A146T). Other four variants further stabilize the structure (G12A, G12V, E49K, and Q61H). The increased or decreased free energy values demonstrate the influence of these mutations on protein stability. A146T and K117N are highly destabilizing, predicting decreased catalytic ability, while Q61H is highly stabilizing. E49K weakly increases stability, demonstrating low impact on function.

NRAS is a homolog to KRAS, with similar protein structure and function (Colicelli, 2004). NRAS mutations occur in 3.03% of all cancers (Consortium et al., 2017), predominantly in skin (10%–25%) (Jakob et al., 2012), acute myeloid leukemia (9.76%) (Consortium et al., 2017), and colon (∼4%) (Cicenas et al., 2017). In CRC, the majority of NRAS mutations occur at codons 12, 13, and 61, with codon 61 being the most frequently mutated (Vaughn et al., 2011). In this study, two NRAS variants were identified, Q61K and R102Q (2.32%). Previous studies from Saudi Arabia have reported NRAS variants occurring at a rate of 2%–4% (Al-Shamsi et al., 2016; Alharbi et al., 2021; Alfahed, 2023). Our findings align with this but are 2% lower than the international CRC prevalence. The R102Q variant, although present in low VAF, is considered novel and has not been previously reported in CRC. The low VAF suggests it could be an artifact, present in a small subset of cancer cells, or with potential biological relevance (i.e., low- frequency driver mutations arising in later-stage subclones). CancerVar classifies Q61K as Tier II—potential clinical significance, as expected. R102Q is predicted to be benign by CancerVar pathogenicity prediction tools and is classified as a Tier III variant of unknown significance, though our study’s selected tools suggest it is pathogenic. Therefore, its true pathogenicity remains to be functionally tested for confirmation. NRAS is located at chromosome 1p13.2 (Mitchell et al., 2008) and encodes a protein of 189 amino acids across seven exons (Eisfeld et al., 2014). The Q61K variant is mapped to the S2 region of exon 3, within the G catalytic domain. Due to this location, it is expected to affect NRAS proteins by disrupting GTP hydrolysis, thus activating the protein. R102Q is located in exon 4, also a part of the catalytic domain. Therefore, it is predicted to similarly influence the catalytic site’s conformation. 3D structural modeling showed significant deviations between the variant and wild-type structures of Q61K (4.263) and R102Q (2.4139). Both variants increased protein stability (Q61K: 0.486; R102Q: 0.012). Q61K likely makes the G domain more rigid, altering the binding site, while R102Q’s impact on G domain function cannot be predicted based on our data.

We compared our selected tools with better algorithms against CancerVar tools to make accurate predictions for identified somatic variants. Our pathogenicity prediction tools showed a significant difference against CancerVar tools. Variants such as KRAS G12S, E49K, and NRAS R102Q were predicted as benign by CancerVar but pathogenic by our tools. The other variants were consistently identified as pathogenic (Tables 2, 3). KRAS G12S, located at codon 12, is expected to confer drug resistance due to increased activation. CancerVar’s tools mostly predicted it as benign, while our tools as pathogenic. ROC analysis demonstrated the superior performance of 6/7 (85%) of our selected tools (PrimateAI, BayesDel_addAF, REVEL, FATHMM-MKL, CADD_Phred). Significant z-scores further supported these findings, highlighting tool performance variability. Overall, study-selected tools outperformed CancerVar tools, probably due to the advanced machine learning algorithms these tools use to integrate diverse data sources. The best-performing tools in this study rely on the more complex algorithms such as deep learning, multiple kernel learning, support vector machines, deep neural networks, and advanced Bayesian statistical frameworks (PrimateAI, REVEL, FATHMM-MKL, CADD_Phred, DANN, and BayesDel_addAF, respectively). They also incorporate additional features, such as evolutionary conservation, functional annotations, and population allele frequencies, providing different angles of prediction. By comparison, the CancerVar tools SIFT, FATHMM, and PolyPhen2_HDIV are older and employ simpler algorithms and most likely out-of-date training data. Being foundational models, their much lower complexity may be affecting their accuracy relative to newer tools. The newer tools also have the largest and most diverse training datasets, covering more variant types and contexts. This study emphasizes that advanced and diverse prediction tools, when possible, must be used in combination for comprehensive variant assessment. It is, however, paramount to complement such predictions with clinical or functional evidence for the correct determination of pathogenicity for a somatic variant.

In conclusion, our study, one of the few from Saudi Arabia, supports pathogenicity for some known variants and predicts the pathogenicity of other rare variants. Functional studies, however, are needed for validation. We reported on the prevalence of somatic variants from CRC patients in a Saudi hospital, identifying 14 variants in actionable genes, including one novel variant (NRAS R102Q) and two rare variants (KRAS E49K and G138E). KRAS G138E, present in three patients, was previously reported only twice from large intestine samples and might be more common in Saudi or Middle Eastern populations. The KRAS E49K co-occurred with KRAS G138E, while the NRAS R102Q co-occurred with KRAS G12D. Wild-type BRAF proteins showed higher binding affinity to Encorafenib, suggesting more effective inhibition compared to the variant protein V600E. Our genotype-protein-phenotype analyses emphasized the pathogenicity of identified variants, highlighting the necessity of effective targeted therapies. Whether the identified rare or novel variants were functionally shown to be activating their respective proteins, clinical decisions can be made to tailor treatment for patients harboring these variants. The prevalence of somatic variants in Saudi Arabia within the BRAF, KRAS, and NRAS genes is roughly 47.7%, which is consistent with global estimates. As cell cycle modulators, these genes continue to be attractive therapeutic targets. Overall, this study highlights the significance of comprehensive molecular screening and bioinformatics in understanding the mutational landscape of CRC in the Saudi population, with the goal of improving targeted drug therapy.