The essential bioinformatics tools used in this research could be found in Supplementary Table S1. Tumor Immune Estimation Resource (TIMER2.0) is a comprehensive resource for the systematical analysis of immune infiltrates across diverse cancer types (Li et al., 2020). We studied the differential expression of the BRAF gene between tumor samples and adjacent normal tissues across all TCGA tumors by using the “Gene_DE” module of the TIMER2.0 database. The transcripts per million (TPM) values of transcription factors were log2-converted. For certain cancer types without adjacent normal tissues, we further explored Gene Expression Profiling Interactive Analysis 2 (GEPIA2) portal to investigate the BRAF mRNA expression difference between tumor samples and matched TCGA normal and Genotype-Tissue Expression (GTEx) data. Additionally, the GEPIA2 data portal was also used to generate violin plots of BRAF expression across pathological stages for all TCGA cancer types. The log2 (TPM +1) transformed expression data were applied for the violin plots here. GEPIA2 is an interactive web server for analyzing the RNA sequencing expression data from the TCGA and the GTEx projects and provides customizable functions such as tumor/normal differential expression analysis, survival analysis, and so on (Tang et al., 2019). Next, we analyzed the methylation difference of BRAF and its downstream genes of the MAPK signaling pathway between tumor samples and adjacent normal tissues in various TCGA cancer types by using the “TCGA Cancer-Methylation” module in the Gene Set Cancer Analysis (GSCALite) platform. Furthermore, the correlation between methylation and gene expression of BRAF and downstream genes was also visualized by the GSCALite platform. GSCALite is a web-based analysis platform that integrates cancer genomics data to analyze gene methylation, drug sensitivity, and so on (Liu et al., 2018).

The cBioPortal is an open-access portal that provides an interactive investigation of multidimensional cancer genomics and clinical data (Gao et al., 2013). In this study, we selected the “TCGA PanCancer Atlas Studies” covering 10,967 samples across 32 cancer types to further explore BRAF alterations. Data files including copy-number alterations, mutations, mRNA Expression, Log2 copy-number values (CNV), and clinical data were downloaded from cBioportal. The BRAF mRNA expression data were performed based on RSEM (batch normalized from Illumina HiSeq_RNASeqV2) and then log10 transformed. For the BRAF CNV data, the log-ratio value represents: 2 = amplification; 1 = gain; 0 = diploid; -1 = shallow deletion; and −2 = deep deletion. In addition, BRAF fusion data were collected from the TCGA Fusion Gene Database, which enables researchers to query cancer-associated transcript fusions in an interactive manner (Hu et al., 2018).

The Kaplan-Meier plotter is an open-access online database that enables researchers to assess the effect of a candidate gene on survival analysis in pan-cancer (Gyorffy, 2021). The correlations between BRAF mRNA expression and clinical prognosis across various cancer types were analyzed with the “Start KM Plotter for pan-cancer” module. Data including hazard ratio (HR), p-value, and 95% confidence interval (CI) were collected to draw the forest plots to summarize survival analysis.

The statistical analysis was performed with Graphpad PRISM software Version 8.0. Student’s t-test, Cox regression analysis, and linear regression analysis were conducted when appropriate. p < 0.05 was defined as a statistically significant difference.

Aberrant expression of BRAF gene has been demonstrated in various cancer types (Wang et al., 2021; Zhao et al., 2021). Previous studies on BRAF expression in cancer have used inconsistent research methods and have been limited to small sample sizes and/or to single or limited cancer types. In this study, we conducted a more comprehensive analysis of BRAF expression in pan-cancer. At the outset, we explored the mRNA expression pattern of BRAF between tumor samples and adjacent normal tissues in pan-cancer by TIMER2.0 (Figure 1A). Compared with the corresponding adjacent normal tissues or metastatic lesions, significantly differential expression of BRAF was found in 14 cancer types, with 9 tumor types up-regulated (cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), cholangiocarcinoma (CHOL), colon adenocarcinoma (COAD), esophageal carcinoma (ESCA), liver hepatocellular carcinoma (LIHC), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), stomach adenocarcinoma (STAD) and uterine corpus endometrial carcinoma (UCEC)) and 5 tumor types down-regulated (breast invasive carcinoma (BRCA), kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), SKCM, THCA). After adding GTEx normal tissue samples as a control, we further investigated BRAF differential expression by GEPIA2. As shown in Figure 1B, BRAF expression was up-regulated in pheochromocytoma and paraganglioma (PCPG) and CHOL, down-regulated in testicular germ cell tumors (TGCT) and uterine carcinosarcoma (UCS). We further analyzed the correlation between BRAF expression and pathological stages in pan-cancer. As shown in Figure 1C, we found that BRAF expression was correlated with pathological stages in several tumor types, including COAD, KIRC, LUSC, and ovarian serous cystadenocarcinoma (OV) (p < 0.05).

A growing body of evidence has suggested that DNA methylation is strongly correlated with gene alteration in cancers (Wang et al., 2020; Zhang et al., 2020). Therefore, we searched the methylation profiles of BRAF and its downstream genes in TCGA cancers by using the GSCALite database (Figures 1D,E). The results indicated that the methylation of BRAF was down-regulated in LUSC, head and neck squamous cell carcinoma (HNSC), and UCEC. Then, we evaluated the correlation between methylation and BRAF expression in pan-cancer. The results revealed that the expression profiles of BRAF and downstream genes were generally negatively correlated with methylation in various cancers.

The overall somatic mutation frequency of BRAF was 7.7% for all cancer samples (848/10,976) and 7.0% for all patients (767/10,953) across the 32 TCGA cancer types. And the detailed information on 848 BRAF somatic mutations was shown in Supplementary Table S2. The sample size of each tumor type varied from 36 (CHOL) to 1,084 (BRCA), and the cancer types with a small sample size might not reflect the general spectrum of BRAF mutation status (Supplementary Table S3). As shown in Figure 2A, the most frequent cancer types with BRAF mutations were THCA (59.3%), SKCM (53.6%), colon adenocarcinoma/rectum adenocarcinoma (COADREAD) (10.6%), LUAD (7.2%) and UCEC (4.7%). Instead, almost no BRAF mutations were observed in kidney chromophobe (KICH), acute myeloid leukemia (LAML), LIHC, TGCT, thymoma (THYM), and uveal melanoma (UVM).

Based on the Pfam database (http://pfam.xfam.org/protein/braf_human), BRAF harbors 3 functional domains, including the RBD domain (156-227 aa), C1_1 domain (235-280 aa), and Pkinase_Tyr domain (458-712 aa). The 848 BRAF somatic mutations were observed in various cancer types and widely distributed across different functional domains of the BRAF gene. The most common one was the Pkinase_Tyr domain (722 samples), followed by the other domains whose functions were barely known (75 samples), the C1_1 domain (11 samples), and the RBD domain (8 samples). Fusions (32 samples) were also observed in BRAF somatic mutations across all cancer types. The location distribution of BRAF mutations was dramatically different among numerous cancers. Mutations in THCA, SKCM, COADREAD, and LUAD were most frequently distributed in the Pkinase_Tyr domain. However, mutations in UCEC were predominantly located in the other domains amounting to half of the total mutations. Furthermore, fusions were mainly distributed in THCA and SKCM (Figure 2B and Supplementary Table S4).

Fusion genes generated by cleavage and re-splicing at the genome level are often the targets for tumor diagnosis and treatment. We analyzed fusion transcripts of BRAF across various cancer types by using the TCGA Fusion Gene Database (Figure 3). BRAF fusion transcripts were detected in THCA (17), SKCM (9), prostate adenocarcinoma (PRAD) (3), pancreatic adenocarcinoma (PAAD) (2), READ (2), LIHC (1), LUSC (1), STAD (1), KIRP (1), brain lower-grade glioma (LGG) (1) and bladder urothelial carcinoma (BLCA) (1). The highest number of fusion transcripts was found in THCA (three SND1_BRAF, one BRAF_SND1, one AGK_BRAF, one BRAF_AGK, one MACF1_BRAF, one BRAF_MACF1, one FAM114A2_BRAF, one BRAF_ FAM114A2, one CCNY_BRAF, one MKRN1_BRAF, etc.). AGK_BRAF and BRAF_AGK were also detected in SKCM. The vast majority of these BRAF fusion transcripts were classified as in-frame, while three BRAF fusion transcripts (one BRAF_HIBADH and one HIBADH_BRAF in SKCM, one TMPRSS2_BRAF in PRAD) were classified as out-of-frame and three (one BRAF_MRPS33 in BLCA and one in STAD, one BRAF_CUL1 in KIRP) were classified as CDS-5UTR.

According to functional impact on protein coding, these 848 BRAF somatic mutations could be classified into four categories: missense (778 mutations), truncating (32 mutations), fusion (32 samples), and in-frame (6 mutations) (Supplementary Figure S1A). The 600aa in the Pkinase_Tyr domain was the most mutated position, which was observed in 590 samples (545 samples with V600E, 39 with V600M, 4 with V600G, 1 with V600_K601delinsE and 1 with V600_R603del). Mutations at V600E were primarily distributed in THCA samples (284/545) and SKCM samples (193/545) (Supplementary Figure S1B–C). V600E is known to be oncogenic and serves as target for BRAF inhibitors approved by FDA, while other mutation types (V600M/G) are likely oncogenic. Other tumors with mutations at this position were COADREAD (48 samples), LUAD (9 samples), glioblastoma multiforme (GBM) (5 samples), KIRP (2 samples), BLCA (1 sample), LGG (1 sample), CHOL (1 sample) and HNSC (1 sample). Studies on their role in COADREAD and LUAD are underway (Planchard et al., 2017; Kopetz et al., 2019) and its function remains little known in other cancer types.

Based on the oncogenic effect and predictive significance, the 848 BRAF somatic mutations could be classified into four categories. As shown in Figure 4A, 616 (72.6%) BRAF mutations were oncogenic, 98 (11.6%) likely oncogenic, 1 (0.1%) inconclusive and 133 (15.7%) unknown. Although a major portion of BRAF somatic mutations was distributed in the functional categories, there were still some mutations in the unknown class deserving further study to characterize the potential functional significances of these mutations. As displayed in Figure 4B, mutations distributed in the functional categories comprised the majority of BRAF mutations in several cancers such as THCA, SKCM, COADREAD, and LUAD. However, more than two-thirds of mutations belonged to the unknown class in UCEC, LUSC, BRCA, and STAD.

The 848 BRAF mutations could be divided into five levels by the clinical targeted therapy implication, containing level NA (193 mutations), level 4 (22 mutations), level 3B (326 mutations), and level 3A (14 mutations), and level 1 (293 mutations). Only level 1 mutations are indicated for targeted therapy with FDA-approved drugs. All level 1 mutations were found in SKCM (236), COADREAD (48), and LUAD (9). These mutations (including V600E, V600M, and V600G) were concentrated in 600aa of the Pkinase_Tyr domain. Although all mutations in THCA belonged to oncogenic (289/302)/likely oncogenic (13/302), all of them were in level 3B (301/302) and level NA (1/302) without treatment implications (Figures 5A,B).

The overall BRAF CNV frequency was about 42.7% (4,684/10,967 samples). The majority of BRAF CNV types were gain (3,555 samples), followed by shallow deletion (972 samples), application (133 samples), and deep deletion (24 samples). The most common tumor types with BRAF CNVs were GBM (81.3%), TGCT (69.1%), ESCA (63.7%), OV (62.9%), adrenocortical carcinoma (ACC) (62.0%), KIRP (61.1%), UCS (59.6%) and LUSC (59.5%). By contrast, THYM (16.3%), PCPG (15.7%), LAML (12.0%), UVM (8.8%) and THCA (5.2%) harbored a relatively low BRAF CNV frequency (Figure 6A). Among the 848 samples with BRAF mutations described above, 276 also had BRAF CNV changes, of which 236 with gain, 23 with amplification, 14 with shallow deletion, and 3 with deep deletion. SKCM had the highest numbers of amplification or gain among different cancer types. Meanwhile, SKCM and LUSC were the two cancer types with the highest numbers of shallow deletion, and deep deletion only occurred in THCA (Figure 6B, Supplementary Table S2).

BRAF mRNA expression was compared across 32 TCGA cancer types and exhibited a relatively consistent trend, suggesting that there may be a common mechanism to promote BRAF expression. Based on the interquartile range, BRAF expression was widely distributed in COADREAD and SKCM, and narrowly distributed in UVM and mesothelioma (MESO), which may be attributed to the fact that some cancer types had more than one subtype and therefore more genetic diversity (Supplementary Figure 2A). In addition, we explored the correlation between BRAF mRNA expression and CNVs. The results showed that there was a positive correlation between BRAF mRNA expression and CNVs in pan-cancer (r = 0.1716, p < 0.0001) (Supplementary Figure 2B).

The combined BRAF alteration (CNV and mutation) frequency in all cancers was about 8.3% (905 of 10,967 samples). As shown in Figure 7A, BRAF alteration frequency among various cancer types was dramatically different. While KICH, TGCT, THYM, and UVM had neither BRAF mutation nor BRAF CNVs, BRAF alterations were most frequently observed in THCA (59.6%) and SKCM (53.8%), in which mutation was more common than CNV, with BRAF mutation rates of 57.4 and 49.3%, respectively. Other cancers with dominant BRAF mutation but at much lower mutation rate included COADREAD (10.4%), LUAD (7.2%), UCEC (4.7%), BLCA (3.2%), STAD (3.0%), LUSC (2.9%), CHOL (2.8%) and DLBC (2.1%). Amplification was more common in OV (7.7 vs. 0.3%), SARC (1.6 vs. 0.4%) and LGG (1.4 vs. 0.6%). Deep deletion was mainly distributed in LAML (1.5%), SARC (1.2%) and ESCA (1.1%).

BRAF mutation location and its CNVs occurrence appeared to be associated. 201 of 664 (30.3%) mutations in the Pkinase_Tyr domain and 21 of 75 (28%) mutations in the other function-unknown domain were accompanied by copy gain. Amplification was mainly distributed in the Pkinase_Tyr domain. Mutations in the RBD domain and C1_1 domain rarely had concurrent CNVs (Figure 7B).

In order to assess the clinical significance of BRAF expression, we analyzed patient survival in pan-cancer and showed that increased BRAF expression was associated with poor patient overall survival (OS) in LIHC, OV, and UCEC. Interestingly, increased BRAF expression was correlated with better prognosis in BRCA, HNSC, and KIRC (Figure 8A). In addition, survival analysis between BRAF expression and patient relapse-free survival (RFS) across various cancer types exhibited that increased BRAF expression was associated with poor RFS in LIHC, LUSC, and UCEC, while high BRAF expression was correlated with better RFS in BRCA and OV (Figure 8B). The contradictory results in OV may be attributed to insufficient sample size and diverse genetic backgrounds.