Extraction of mutation signatures in CRC
A total of 192905 mutations were detected in 535 samples with a median of 91, including single nucleotide variants (SNVs) and small insertions and deletions (Indels). SNV was the main mutational type, in which C>T displayed the highest frequency followed by C>A and T>C. Among the top 10 most frequently mutated genes, APC possessed the highest number of delete mutations (236) and TTN possessed the highest number of missense mutations (646) (Fig S1). To gain insights into the potential mutation generation processes operative in patients with CRC, we decomposed the mutation signatures via NMF algorithm (Fig S2A). Subsequently, eight mutation signatures were extracted from the CRC genomic data and annotated them against the (COSMIC) signature nomenclature based on cosine similarity (Fig S2B). Therefore, the extracted mutation signatures were ultimately called as cosmic signature 1, 6, 10, 15, 20, 28, 29, and 30 (Fig 1A). The clocklike signature 1 is thought to be connected with age-related accumulation of spontaneous deamination of 5-methylcytosine. Signatures 6, 15, and 20 are all associated with defective DNA mismatch repair. Signature 10 featured by altered activity of the error-prone polymerase POLE is often found in six cancer types including CRC. Signature 29 exhibits transcriptional strand bias for C>A mutations due to tobacco chewing habit. Signatures 28 and 30 have been observed in a subset of stomach and breast cancers remaining unknown etiology.
Fig. 1 The extraction of mutation signatures and generation of the mutation signature relevant subtypes in CRC. A. Eight mutation signatures were deciphered (mutation signature 1, 6, 10, 15, 20, 28, 29, and 30) based on NMF algorithm and COSMIC signatures. B. The consensus score matrix of all samples when k = 2. A higher consensus score between two samples indicates they are more likely to be grouped into the same cluster in different iterations. C-D. Kaplan–Meier analysis for OS (C) and DFS (D) between MSC-1 and MSC-2. E. Pie charts show the relative proportion of eight categories of mutation patterns in MSC-1 and MSC-2, respectively. F. The difference of APOBEC enrichment score between MSC-1 and MSC-2. G. Mutational oncoplot of 11 APOBEC associated genes in two subtypes.
Generation of the mutation signature relevant subtypes
Based on the mutation signatures deciphered in the CRC genome, consensus clustering analysis revealed the two subtypes was the optimal choice (Fig 1B). The CDF curve, PAC value, and Nbclust results further confirmed the stable and reliable of the cluster results (Fig S2C-SE). We annotated the two subtypes as mutation signature cluster (MSC) 1 and 2 respectively. The Kaplan-Meier survival analysis suggested MSC-1 was significantly associated with favorable prognosis (Log-rank OS: p =0.005; DFS: p =0.070) (Fig 1C-D). Of note, MSC-1 (n=226) had dominant signature 6, 15, and 20, linked to defective DNA mismatch repair. MSC-2 (n=309) was characterized by signature 1 and 29, which was associated with age and tobacco chewing habit (Fig 1E; Fig S2F). We also observed the APOBEC signature enrichment score was significantly higher (p =0.003) in MSC1, which indicated MSC-1 had a higher occurrence of C>T transition in TpCpW motifs (Fig 1F). Previous study demonstrated the mutation of APOBEC family might contribute to high tumor mutation burden (TMB) [26]. Therefore, we further explored the mutation in the APOBEC family and found the rare mutation in the CRC patients. Of note, MSC-1 had the more mutation proportion relative to MSC-2, in line with the APOBEC enrichment score (p =0.003) (Fig 1G), which might give rise to the high mutation rate of MSC-1.
Somatic mutation landscape of two subtypes
The tumor mutation burden (TMB) in MSC-1 was significantly higher than MSC-2 (p <0.001) (Fig 3SA). The higher TMB may tend to occur in patience with defective DNA mismatch repair relative mutation signatures[7]. We further determined 28 mutation driven genes with MutSigCV q-value <0.05 and mutation frequency >10% in CRC (Table S5; Fig 2A). Out of these 28 genes, 18 genes have been reported in at least one CRC associated research, such as APC, TP53, KRAS, SYNE1, PIKSCA, and FBXW7 et.al. Besides, ten novel drivers were identified including DNAH11, USHA2, HMCN1, HYDIN, MDN1, DST, VPS13B, DNAH8, EYS, and NBEA. We also dissected the prognostic role of these genes. The mutation of EYS prolonged DFS, and the mutation of USH2A suggested an unfavorable OS (Fig 2B-C). In two MSC subtypes, 22 out of 28 drivers exhibited significant mutation differences (Fig 2D). Consistent with the high TMB in MSC-1, the mutation frequency of most divers was also superior in MSC-1, such as ATM, SOX-9, and PRKDC. Of note, APC and KRAS, the early mutation event of colon adenoma–carcinoma process[27], dominantly mutated in MSC-2, which implied that familial adenomatous polyposis (FAP) may be one of main causes of MSC-2. Further analyses revealed a predominant mutation co-occurrence relationship between KRAS and SYNE1, TP53 and SYNE1, as well as APC and USH2A et.al (Fig S3B). Interestingly, we found some specific co-occurrence phenomenon such as BRAF-HMCN and DNAH17-MDN1 that appeared only in MSC-1, which suggested the specific co-occurrences of BRAF-HMCN and DNAH17-MDN1 could be promisingly employed to distinguish different subtypes (BRAF-HMCN: p <0.001; DNAH17-MDN1: p <0.001) (Fig 2E-F). In addition, for the first time, we revealed the prognostic value of some co-occurrence relationship, the co-occurrence of APC-TP53 demonstrated a favorable DFS (Fig S3C), and the co-occurrence of APC-KRAS, KRAS-TP53, and KRAS-SYNE1 was significantly associated with worse DFS (Fig 2G-I). Furthermore, the literature had confirmed that CRC patients with a defected mismatched repair (MMR) system could lead to hypermutation and microsatellite instability (MSI)[28]. Hence, we investigated the mutation status of nine known MMR genes, and the results exhibited MSC-1 had most mutations of MRR genes (Fig S3D), and the cases with MMR genes mutation relatively high in MSC-1 (26% vs. 7%; p <0.001) (Fig 2J), which was in line with its specific mutation signatures such as signature 6, 15, and 20.
Fig. 2 The mutation driven genes in CRC. A. Mutational oncoplot of mutation driven genes in MSC-1 and MSC-2. Rows are genes and columns are tumor samples. B-C. Kaplan–Meier survival analysis of EYS (B) and USH2A mutations (C). D. The mutation frequency of mutational drivers in two subtypes, ***P < 0.001. E-F. The relative proportion of BRAF-HMCN (E) and DNAH17-MDN1 (F) co-occurrences in two subtypes. J. The relative proportion of patients with the MMR genes mutations in two subtypes. G-I. Kaplan–Meier survival analysis of APC−KRAS (G), KRAS−TP53 (H), and KRAS−SYNE1 (I) co-occurrence.
SCNA investigation of two subtypes
In arm level, both the gain and loss load were significantly higher in the MSC-2 relative to MSC-1 (p <0.05). Although there was no statistical significance in the focal level load between the two subtypes, the slightly superior trends were also shown in MSC-2 (Fig S4A). Different from MSC-1 which was characterized by higher mutation load, MSC-2 might be dominant in the alteration of copy number. According to the GISTIC algorithm, we eventually identified 39 driven segments encompassing 14 amplification segments and 25 deletion segments (Table S6-S7; Fig S4B). We further compared the alteration frequency of 39 segments between the two subtypes, and found MSC-1 had generally low frequency compared to MSC-2, in accordance with the CNA load (Fig 3A). We also found a multitude of oncogenes and tumor suppressor genes in these driven segments, which might play an essential role in the tumorigenesis and progression of CRC, such as MYC (8q24.21), CCND3 (6p21.1), ERBB2 (17q12), PTEN (10q23.31), SMAD4 (18q21.2), and APC (5q22.2) et.al (Fig 3B). Although MSC-2 had generally frequent SCNA events of these genes, there were still high fraction amplification or deletion occurred in MSC-1, such as MYC, FTK3, MCC as well as NOTCH and TGF-beta pathway associated genes. Interestingly, we found oncogenes with only amplification and tumor suppressor genes with only deletion. Thus, the gene expression difference between gain and no-gain or loss and no-loss was further explored, and we found oncogenes with gain was more prone to overexpress such as ERBB2, MYC, and MLST8, and the expression of tumor suppressor genes with loss was predominantly lower than no-loss group such as APC, SMAD4, and PTEN (p <0.001) (Fig 3C; Fig S4C). These results suggested CNA status played a master regulator role in the aberrant expression of oncogenes and tumor suppressor genes in CRC. Further survival analysis demonstrated the prognosis significance of these genes (Fig S4D-E). For the first time reported, the gain of MLST8 and MAP2K2 could prolong OS (Fig 3C; Fig S4F), the gain of CCND3 indicated the worse DFS (Fig 3D), as well as the loss of CTNN6, DKK1, APC, MCC and SMAD4 was also associated with unfavorable DFS (Fig S4F). Moreover, we also investigated the CNA of MMR genes, and found the fraction of patients with the MMR genes deletion was higher in MSC-2 relative to MSC-1 (62% vs. 53%; p =0.042) (Fig 3F-G). It was nonnegligible that some MMR genes displayed high level of loss frequency such as MLH3, MSH4, MSH3, and MLH1, so that might diminish the expression of MMR genes and give rise to MSI of CRC.
Fig. 3 The driven segments identified from GISTIC algorithm in CRC. A. The amplification (orange) and deletion (purple) frequency of 39 driven segments in two subtypes. B. The distribution of CNA relevant oncogenes and tumor suppressive genes in two subtypes. C. The expression difference of ERRB2 and MYC between the gain and no-gain groups, as well as APC and SMAD4 between the loss and no-loss groups. ***, P < 0.001. D-E. Kaplan–Meier survival analysis of MLST (D) and CCND3 (E) gain. F. The relative proportion of patients with the MMR genes deletions in two subtypes. G. Oncoplot for the deletion of nine MMR-related genes in two subtypes.
Methylation driven genes
To identify methylation driven genes (MDGs) in CRC, the MethyMix package and the Wheeler criterion were employed. The MethyMix algorithm screened 608 genes that their expression was significantly related to methylation events, and the Wheeler criterion filtered out 147 epigenetic silencing genes (Table S8). Eventually, we determined a total of 69 MDGs by the intersection of the two methods. Further univariate Cox regression uncovered the prognosis significance of these MSGs (Table S9). The high methylation of TBX1, GREB1L, and CNNM1 was significantly associated with the unfavorable OS (Fig 4A-C). In addition, we defined the subtype-specific MDG (ssMDG), and 13 ssMDGs were significantly different in their expression and methylation between the two MSC subtypes (Fig S5A-S5B). In terms of these ssMDGs, we observed significant negatively correlation between the expression and methylation level (Fig 4D). MSC-2 dominated in the hypermethylation of AQP5 and ZNF304 compared to MSC-1. To our knowledge, AQP5 was a potential epigenetic driver of tumor development[29]. The other 11 ssMDGs were specific for MSC-1, such as ADAM32, SLC35D3, and TMEM150C. Of note, the MLH1 was also a specific ssMDG of MSC-1. As illustrated, the methylation level of MLH-1 was dominantly higher relative to MSC-2, and the expression level was the opposite (Fig S5A-S5B). Previous report demonstrated the hypermethylation of MLH-1 was a potential mechanism contributing to the MSI of CRC[30]. We thus divided CRC patients into the methylation cases and unmethylation cases based on the threshold of beta =0.3, and found all methylation cases occurred in MSC-1 (22% vs. 0%; p <0.001) (Fig 4E), which explained its specific MSI-associated mutation signatures to some extent.
Fig. 4 The methylation driven genes in CRC. A-C. Kaplan–Meier survival analysis of TBX18 (A), GREB1L, (B) and CNNM1 (C) methylation. D. The correlation analysis between the methylation and mRNA expression levels of13 ssMDGs. E. The relative proportion of patients with the MMR genes methylation events between two subtypes.
Functional status, immune cell infiltration and immunogenicity assessment
We performed the biological process and KEGG pathway enrichment analysis through the GSEA approach. The MSC-1 subtype was tightly associated with immune related pathways such as adaptive immune response, antigen processing and presentation, response to interferon-gamma, and Th1 and Th2 cell differentiation (Fig 5A and 5C). The MSC-2 subtype was significantly enriched in reactive stroma related pathways such as epidermis or mesenchymal morphogenesis, mesenchymal cell proliferation, transform growth factor beta (TGF-β) signaling pathway, and Wnt signaling pathway (Fig 5B and 5D). Further GSVA Hallmark pathways assessment suggested the similar result to the above, also elucidated the MSC-1 was dominant in the immune activation such as canonical T cell excitation pathway: interferon-gamma, and the MSC-2 was master in the stromal activation such as TGF-β process (Fig 5E). In addition, we also evaluated the subpopulations difference of eight immune cells and two nonimmune cells between the two subtypes (Fig 5F). Consistently, the immune killing cells such as T cells, CD8+ T cells, cytotoxic lymphocytes, and nature killer cells were superior in MSC-1, and MSC-2 was characterized by higher fibroblasts. The leukocyte and stomal fraction data retrieved from Thorsson et al. study also demonstrated the dominant role of the two subtypes in immune activation and stromal activation, respectively (Fig 5G-H).
Furthermore, 17 indicators were applied to decode the immunogenicity features of the two subtypes (Table S2; Fig 5I). In line with the specific mutation signatures such as 6, 10, and 15 in MSC-1, the nonsilent mutation rate and MSI score were higher in MSC-1 (Fig 5J-K), meanwhile, SNV and Indels neoantigens were also more prone to occur in MSC-1 relative to MSC-2 (Fig 5L-M), but there was no significance in term of CTA score (Fig S6A). Conversely, the CNV-relevant indicators were slightly high in MSC-2 such as AS, ITH, number or fraction of segments alteration, HRD, and LOH, although most of them did not reach statistical significance (Fig 5N and Fig S6B-S6G). These results implied the immunogenicity of two subtypes might be derived from different genome alterations. In addition, the BCR/TCR diversity and CYT that may reflect a robust antitumor response and cytolytic activity were also higher in MSC-1 (Fig S6H-S6K; Fig 5O). Overall, although there was heterogeneity between the two subtypes in different aspects of immunogenicity, MSC-1 still displayed the stronger immunogenicity compared to MSC-2, which might arise from the predominant mutation pattern. The stronger immunogenicity further conferred the superior immune activation in MSC-1.
Fig. 5 Functional status, immune cell infiltration and immunogenicity assessment. A-B. The biological process significantly enriched in MSC-1 (A) and MSC-2 (B). C-D. The KEGG pathways significantly enriched in MSC-1 (C) and MSC-2 (D). E. The specific Hallmark pathways in MSC-1 and MSC-2. F. The infiltration abundance of eight immune cells and two nonimmune cells populations in MSC-1 and MSC-2. ns, P > 0.05; *, P < 0.05; **, P < 0.01; ***, P < 0.001. G-H. The distribution of leukocyte (G) and stomal (H) fraction in MSC-1 and MSC-2. I. The comparison of 17 immunogenicity associated indicators between two subtypes, the cell represented by the mean value of corresponding cluster divided by the overall mean value. J-O. The distribution of nonsilent mutation rate (J), MSI score (K), SNV neoantigens (L), Indel neoantigens (M), fraction of segments alteration (N), and cytolytic activity (CYT) (O) in MSC-1 and MSC-2.
The expression and regulation of immune checkpoint molecules
We next explored the expression and regulation difference of 75 immune checkpoint molecules (ICMs) between the two MSC subtypes in multi-omics dimensions (Table S3). Obviously, the expression of ICMs was generally high in MSC-1 (Fig 6A; Fig S7A-S7B). For example, the MHC molecules displayed relatively low expression in the MHC-2 (Fig 6B). We further calculated the antigen processing and presenting machinery score (APS) via ssGSEA algorithm, and observed the MSC-2 also presented a lower APS (Fig 6C). That suggested antigen presentation capacity might be impaired. In line with the immune activation status, MSC-1 demonstrated higher expression of stimulatory ICBs such as CCL5, CD40, and ITGB2 (Fig S7A). Meanwhile, the inhibitory ICBs such as IDO1, PDCD1, CTLA4, and CD274 also predominantly expressed in MSC-1, which implied the overexpression of inhibitory ICMs might be the immune escape mechanism of MSC-1 (Fig S7B).
Furthermore, we integrated the mutation, SCNA, and methylation profile to decipher the regulation of ICMs. Notably, although the mutation of ICMs displayed generally rare frequency (Fig 6A), it still presented some effects on the expression of ICMs, for example, the mutation of HLA-B and ITGB2 displayed significant lower expression only in MSC-1, but a slightly high was observed in MSC-2 (Fig 6D; Fig S7C). On the contrary, the SCNA of ICMs demonstrated the relatively prevalent frequency (Fig 6A). CD40 had the highest amplification frequency, but there was no significant expression difference between gain and no-gain groups (Fig S7D). Consistent with the deletion status, the expression of CD276, ICOSLG, TNFRSF9, TNFRSF14, and TNRSF18 was relatively low in loss group compared to no-loss group (Fig 6E-G and Fig S7E-S7F). In addition, the hypermethylation of ICMs also played a critical regulation role in a number of ICMs such as HLA-B, CXCL10, and CD40, we observed their expression was significantly negative correlation with the methylation profile (HLA-B: r =-0.51; CXCL10: r =-0.43; CD40: r =-0.46; all p <0.001) (Fig 6H-J).
Fig. 6 Multi-omics analysis of 75 immunomodulators in the TCGA-CRC cohort. A. From left to right: mRNA expression (z-score), mutation frequency, amplification frequency, deletion frequency, and expression versus methylation (gene expression correlation with DNA-methylation beta value) of 75 immunomodulators in MSC-1 and MSC-2. B. The expression difference of MHC molecules between two subtypes. ns, P > 0.05; *, P < 0.05; **, P < 0.01; ***, P < 0.001. C. The distribution of APS score in MSC-1 and MSC-2. D-G. The expression difference of HLA-B (D) between the mutant and wild groups, as well as CD276 (E), TNFRSF9 (F), and TNFRSF14 (G) between loss and no-loss groups. H-J. The correlation analysis between the methylation and mRNA expression levels of HLA-B (H), CXCL10 (I), and CD40 (J).
Identify the reliable gene pair markers for prognosis and distant metastasis
A total of 108 differentially expression genes (DEGs) were screened between the two MSC subtypes (Fig 7A; Table S11). We further transformed the gene expression matrix into the two gene expression relationships matrix. Based on the pipeline to screen consensus prognosis gene pair signature (CPGPS), we eventually determined three gene pairs with dominant prognostic significance in at least five cohorts, which were FAM83A|IDO1, FABP4|KLK12, and FABP4|GBP5 (Fig 7B-C; Fig S7A). Of note, the gene pairs with a single relationship ratio >90% of cases in corresponding cohort would be expurgated. Ultimately, the FAM83A|IDO1 was removed in the GSE103479 and GSE87211 cohorts, the FABP4|KLK12 was absented in the GSE103479, GSE87211, GSE18105, GSE21510, GSE27854, and GSE71222 cohorts, and the FABP4|GBP5 was deleted in the TCGA-CRC, GSE103479, GSE72970, GSE87211, GSE18105, GSE21510, GSE27854, and GSE71222 cohorts. The expression relationships of FAM83A and IDO1 was significantly associated prognosis in 7/9 of cohorts (Fig 7B and Fig 7D-L), and it was a poor prognostic factor when FAM83A > IDO1 in terms of the mRNA expression level. Although there was no significance in GSE17537 and GSE72970, the FAM83A|high group still indicated the adverse prognosis, which might be due to their relatively small sample sizes (Fig 7F and 7L). The gene pair FABP4|KLK12 was also a prognosis marker which exhibited significance in 6/9 of cohorts. It turned out FABP4 > KLK12 was predominantly associated with unfavorable prognosis (Fig 7C and Fig 7M-U). In addition, the patients with FABP4 > GBP5 were more prone to have a poor prognosis in 5/7 of cohorts (Fig S1A-S1H). Further multivariate Cox analysis revealed FAM83A|IDO1 was an independent prognosis factor in most cohorts (7/9) (Table S11). Conversely, the two gene pairs FABP4|KLK12 and FABP4|GBP5 did not perform well here.
Fig. 7 Identification of gene pairs with the ability to predict prognosis of CRC patients. A. Volcano plot of differentially expressed genes (DEGs) between MSC-1 and MSC-2. The abscissa is log2FC, and the ordinate is −log10 (FDR). The red and blue points in the plot represent DEGs with statistical significance (FDR < 0.05 and |log2FC| > 1). B. Forest plot of IDO1-high versus FAM83A-high groups in nine cohorts. C. Forest plot of KLK12-high versus FABP4-high groups in nine cohorts. D-L. Kaplan–Meier survival analysis for FAM83A|IDO1 in the TCGA-CRC (D), GSE17536 (E), GSE17537 (F), GSE29621 (G), GSE38832 (H), GSE39084 (I), GSE39852 (J), GSE71187 (K), and GSE72970 (L) cohorts. M-U. Kaplan–Meier survival analysis for FABP4|KLK12 in the TCGA-CRC (M), GSE17536 (N), GSE17537 (O), GSE29621 (P), GSE38832 (Q), GSE39084 (R), GSE39852 (S), GSE71187 (T), and GSE72970 (U) cohorts.
We then dissected the predictive role in the CRC metastasis of the three gene pairs. Interestingly, the distant metastatic features of CRC were significantly distinct between FAM83A|high and IDO1|high groups in all cohorts (TCGA-CRC: 29% vs. 15%, p =0.007; GSE29621: 46% vs. 17%, p =0.031; GSE39084: 44% vs. 16%, p =0.028; GSE18105: 49% vs. 23%, p =0.021; GSE21510: 47% vs. 25%, p =0.026; GSE27854: 45% vs. 18%, p =0.006; GSE71222: 24% vs. 8%, p =0.012) (Fig 8A-G). Due to the overrepresented single relationship of gene pair, the FABP4|KLK12 and FABP4|GBP5 were retained in three and two cohorts with metastasis information, respectively. Although the statistical significance was not reached in most cohorts, the proportion of patients with metastasis varied between FABP4|high and KLK12|high groups or FABP4|high and GBP5|high groups. For example, the FABP4|high group had more metastasis ratio compared to KLK12|high group (TCGA-CRC: 24% vs. 12%, p =0.001; GSE29621: 57% vs. 25%, p = 0.173; GSE39084: 31% vs. 31%, p =1.000) (Fig 8H-J), and the dominant fraction of metastasis cases occurred in the FABP4|high group relative to GBP5|high group (GSE29621: 44% vs 23%, p =0.199; GSE39084: 44% vs. 27%, p =0.278) (Fig 8K-L). Therefore, the predictive performance of FABP4|KLK12 and FABP4|GBP5 was much weaker than that of FAM83A|IDO1. Taken together, the expression relationships of FAM83A and IDO1 was a very promising biomarker for prognosis and distant metastasis of CRC patients.
Fig. 8 The predictive ability of three prognostic relevant gene pairs for distant metastasis. A-G. The relative proportion of patients with distant metastasis between FAM83A|high and IDO1|high groups in the TCGA-CRC (A), GSE29621 (B), GSE39084 (C), GSE18105 (D), GSE27854 (F), and GSE71222 (G) cohorts. H-J. The relative proportion of patients with distant metastasis between FABP4|high and KLK12|high groups in the TCGA-CRC (H), GSE29621 (I), and GSE39084 (J) cohorts. K-L. The relative proportion of patients with distant metastasis between FABP4|high and GBP5|high groups in GSE29621 (K) and GSE39084 (L) cohorts. M0, no metastasis; M1, metastasis.
Verified the role of FAM83A|IDO1 in prognosis and metastasis using qRT-PCR
QRT-PCR assay was performed in 30 paired CRC tissues and matched adjacent nontumor tissues (Table S4). We observed FAM83A was overexpressed in tumors relative to adjacent nontumor tissues, and the expression of IDO1 was the opposite (p <0.001) (Fig 9A-B). The role of FAM83A|IDO1 in prognosis and metastasis was further explored in the qRT-PCR assay. The clinical outcome details (including survival status and metastasis status) of 30 CRC patients were shown in Fig 9C. There was no correlation between the expression of FAM83A and IDO1. In line with the previous results, when the expression of FAM83A was higher than IDO1, patients had worse OS and DFS (Log-rank p <0.001) (Fig 9D-E), as well as stronger tendency of distant metastasis (71% vs. 13%, p =0.007) (Fig 9F).
Fig. 9 Verified the role of FAM83A|IDO1 in prognosis and metastasis using qRT-PCR. A-B. The expression difference of FAM83A (A) and IDO1 (B) between two subtypes. C. The mRNA expression of FAM83A and IDO1 as well as the clinical outcomes in our cohort. The abscissa is the expression of FAM83A, and the ordinate is the expression of IDO1. Under the line y = x, FAM83A > IDO1, while above it, FAM83A < IDO1. M0, no metastasis; M1, metastasis. OS-0, alive; OS-1, death or censoring; DFS-0, disease free; DFS-1, disease or censoring. D-E. Kaplan–Meier analysis of OS (D) and DFS (E) for FAM83A|IDO1 in our cohort. F. The relative proportion of patients with distant metastasis between FAM83A|high and IDO1|high groups in our cohort.