Differential expression of FZD7 between tumor and normal tissues
In the TCGA databases of mRNA expression, FZD7 gene expression was significantly upregulated in tumor tissues of glioblastoma multiforme (GBM), brain lower-grade glioma (LGG), lung squamous cell carcinoma (LUSC), liver hepatocellular carcinoma (LIHC), stomach adenocarcinoma (STAD), Wilms tumor (WT), pancreatic adenocarcinoma (PAAD), testicular germ cell tumors (TGCT), acute lymphoblastic leukemia (ALL), acute myeloid leukemia (LAML), adrenocortical carcinoma (ACC), and cholangiocarcinoma (CHOL). Contrary to this, FZD7 was downregulated in uterine corpus endometrial carcinoma (UCEC), breast invasive carcinoma (BRCA), cervical and endocervical cancer (CESC), lung adenocarcinoma (LUAD), kidney renal papillary cell carcinoma (KIRP), colon adenocarcinoma (COAD), prostate adenocarcinoma (PRAD), kidney renal clear cell carcinoma (KIRC), skin cutaneous melanoma (SKCM), bladder urothelial carcinoma (BLCA), thyroid carcinoma (THCA), rectum adenocarcinoma (READ), ovarian serous cystadenocarcinoma (OV), and kidney chromophobe (KICH) tumor tissue (Figure 1A). In addition, we performed a correlation analysis between FZD7 expression and the pathological stage of tumors. According to Figure 1B-1J and Figure S1A-S1B, expression of the FZD7 gene increased with tumor pathological progression in LIHC (P = 0.00033), STAD (P = 0.049), COAD (P = 0.029), BLCA (P = 0.023), KICH (P = 0.04), KIRP (P = 3.6e-05), TGCT (P = 0.0025), and UVM (P = 0.033), while lack of correlation was observed in KICH (P = 0.14), indicating that the FZD7 potential dispensable for tumor proliferation and metastasis in several cancers.
Effect of FZD7 on prognosis
In TCGA pan-cancer, we performed Cox regression analysis to elaborate on the relationship between FZD7 expression and overall survival and disease-specific survival. As shown in Figure 2A, upregulation of the FZD7 gene correlated with poor OS in multiple tumors, including ACC, COAD, GBM, KIRP, LGG, LIHC, PAAD, PCPG, STAD, THCA, UCEC, and UVM (all P-values < 0.05). Kaplan-Meier curves indicated that upregulated FZD7 expression was associated with shorter overall survival in LIHC (P = 0.0036), BLCA (P = 0.015), KIRP (P = 0.0015), LGG (P = 0.0033), STAD (P = 0.038), UCEC (P = 0.006), and UVM (P = 0.00018) (Figure 2B-H).
In addition to this, we explored the impact of the FZD7 gene on disease-specific survival. Forest plots showed that FZD7 was a risk factor for DSS in ACC (P = 0.009), COAD (P = 0.018), GBM (P = 0.031), KIRP (P < 0.001), LGG (P = 0.001), LIHC (P = 0.001), STAD (P = 0.002), THCA (P = 0.001), UCEC (P < 0.001), UVM (P < 0.001) (Figure 3A). Elevated FZD7 expression predicted shorter DSS in LIHC (P = 0.018), COAD (P = 0.035), KIRP (P = 0.001), LGG (P = 0.002), STAD (P = 0.019), UCEC (P = 0.005), and UVM (P < 0.001) (Figure 3B-D, F-I). In LUSC (P = 0.022) (Figure 3E), the FZD7-high expression group had a better prognosis than the low expression group. These findings indicate that upregulated FZD7 expression has a profound impact on the prognosis of some tumor patients.
Mutation landscape of FZD7 in pan-cancer
Mutations in specific genes in tumors have a crucial impact on their treatment and prognosis. In this section, we examined genomic alterations and mutations in TCGA pan-cancers according to FZD7 in cBioPortal. The “mutation” was the most frequent alteration in UCEC of FZD7 gene alteration frequency, while copy number “amplification” and “deep deletion” were the most frequent in LUSC and CESC, respectively (Figure 4A). In SKCM, DLBC, COAD / READ, KIRP, MESO, ACC, THYM, and LGG, “mutations” were the only classes of genomic alteration. Multiple copy number alterations were only present in LIHC. As shown in Figure 4B, we found a total of 110 mutation sites for FZD7, which included 101 missense mutations, 8 truncating, and 1 inframe. As shown in Figure 4C, FZD7 mutation information was exhibited for 9349 samples out of 9896 tumor samples.
Tumor mutational burden (TMB) refers to the total number of somatic mutations in defined regions of the tumor genome and can also be defined as the total number of nonsynonymous point mutations per coding region of the tumor genome. Tumor mutational burden and microsatellite instability can respond to anti-tumor immune and predict tumor responsiveness to immunotherapy[26]. We investigated the association between the FZD7 gene and TMB and MSI. We found that FZD7 expression was positively correlated with TMB in THYM, KICH, and ACC, whereas conversely, it was negatively correlated in LIHC, PAAD, THCA, STAD, UCEC, COAD, UCS, KIRP, BRCA, PRAD, BLCA, and ESCA. In MSI, FZD7 expression was positively correlated in TCGT, LUSC, and KICH, but negatively correlated in UCEC, STAD, SARC, PRAD, KIRP, and COAD (Figure 4D).
FZD7 expression and immune infiltration
Immune infiltrating cells in tumors have an indispensable impact on tumor growth, proliferation, metastasis, and even patient survival. We used the ESTIMATE algorithm to calculate the immune score and stromal score in pan-cancer, after which the correlation of FZD7 gene expression with immune score and the stromal score was examined separately. Immune score, stromal score, and FZD7 expression were strongly correlated in LIHC (P = 1.63e-10 in the ImmuneScore, P = 4.33e-10 in StromalScore), PRAD (P < 2.2e-16 in ImmuneScore, P < 2.2e-16 in stromal score), and PAAD (P = 1.9e-9 in ImmuneScore, P < 2.2e-16 in stromal score) (Figure 5A). After that, correlation analysis was then performed between FZD7 expression and immune checkpoint genes including BTLA, CD161, CD27, CD274, CD276, CD28, CD40, CD70, CD80, CD86, CTLA4, HAVCR2, HHLA2, ICOS, ICODLG, IDO, IDO2, LAG3, PDCD1, TIGIT, TNFRSF9 and TNFSF9 (Figure 5B). Nearly all immune checkpoint genes were positively correlated with FZD7 expression in LIHC, PAAD, PCPG, PRAD, and UVM, implicating that FZD7 gene potentially upregulates immune checkpoint genes while also influencing the efficacy of immune checkpoint inhibitors.
The EPIC, TIMER, XCELL, TIDE, CIBERSORT, CIBERSPRT-ABS, MCPCOUNTER, and QUANTISEQTIDE were used to examine the correlation of FZD7 expression with myeloid-derived suppressor cells, Macrophage cells (M1 and M2), T cell regulatory cells, T cell CD8+, and CD4+ Th1, and CD4+ Th2. Intriguingly, FZD7 expression was also positively correlated with myeloid-derived suppressor cells but negatively correlated with CD4+ Th1 and CD4+ Th2 (Figure 5C). Lee et al. demonstrate that M2 macrophages differentiated from tumor-associated macrophages (TAMs) in HCC promote tumor progression by supporting angiogenesis. In addition, the degree of Treg infiltration in the tumor was positively correlated with the poor prognosis of HCC[27]. Our study found a relatively strong positive correlation between the degree of infiltration by MDSCs and FZD7 expression in both LIHC (R = 0.41, P = 6.95e−04) and ACC (R = 0.33, P = 4.32e−03) (Figure S2A). Moreover, in many tumors, increased infiltration levels of TAMs cells are associated with upregulated FZD7 expression, such as LIHC (R = 0.41, p = 1.97e−15), HNSC-HPV– (R = 0.354, p = 2.93e−13), PRAD (R = 0.383, p = 5.59e−16), PAAD (R = 0.533, p = 6.13e−14), TGCT (R = 0.699, p = 7.86e−23), THYM (R = 0.433, p = 1.31e−06), UVM (R = 0.537, p = 4.89e−07), and READ (R = 0.401, p = 9.03e−05) (Figure S2B). Whereas in CHOL (R = 0.42, p = 1.21e−02), DLBC (R = 0.359, p = 2.11e−02), LIHC (R = 0.3, p = 1.39e−08), GBM (R = 0.317, p = 1.64e−04) and UVM (R = 0.34, p = 2.49e−03), the degree of regulatory T cell infiltration in the tumor was likewise positively correlated with FZD7 expression(Figure S3A).
Enrichment analysis associated with the FZD7 gene
To further dissect the involvement of the molecular signaling pathways associated with the FZD7 in tumorigenesis, we performed a protein-protein interaction network of the FZD7 gene through the STRING database with filtering criteria of confidence level > 0.4. A total of 10 genes that met the criteria were included in the PPI network (Figure 6A). In addition, we screened for top 100 genes (Supplemental Table S1) with a most significant correlation of FZD7 in TCGA and GTEx through the GEPIA2 database. The correlation heatmap obtained the top 6 genes that fit the criteria, LRP6 (R = 0.42, P < 2.2e-16), PTPN14 (R = 0.5, P < 2.2e-16), TCF7L (R = 0.44, P < 2.2e-16), ROR2 (R = 0.43, P < 2.2e-16), SOCS5 (R = 0.53, P < 2.2e-16), and DVL3 (R = 0.46, P < 2.2e-16) (Figure 6B). Nearly all tumors in pan cancers showed a positive correlation with FZD7-associated genes (Figure 6C). The Venn diagram revealed three genes common to FZD7 correlated genes and FZD7 interacted proteins, LRP6, DVL3, and ROR2 (Figure 6D). The FZD7 interacted proteins and FZD7 associated genes were integrated for GO and KEGG signaling analysis. As shown in Figure 6E, the more concentrated signaling pathways in KEGG were “Wnt signaling pathway”, “mTOR signaling pathway”, “Hippo signaling pathway”, and “hepatocellular carcinoma”. In GO, the more enriched signaling pathways were “G protein-coupled receptor binding”, “Wnt protein binding”, “beta-catenin binding” in biological process (BP), “cell−substrate junction”, “cell−substrate adherens junction” and “focal adhesion” in cellular component (CC), “regulation of Wnt signaling pathway” and “neural tube development” in molecular function (MF) (Figure S4A). This suggests that FZD7 plays an inexorable role in tumor initiation, proliferation, and metastasis. CancerSEA was used to perform single-cell analysis in several tumors, and as shown in Figure S4B, differentiation, DNA repair, inflammation, and metastasis can be activated by the FZD7 gene to participate in tumorigenesis. Further, we performed gene set enrichment analysis (GSEA) using RNA sequencing data of LIHC from TCGA based on the FZD7 gene. GSEA indicated denser enrichment of FZD7 genes in immune response regulatory signaling pathways, immunoglobulin complexes, lymphocyte-mediated immune responses, and T cell receptor complexes (Figure S4C). Enrichment plots of transcription factor targets indicated that the FZD7 gene could be involved in HCC by targeting HOXB7, NFE2L3, NKX6, and SOX11(Figure S4C).
Somatic mutation landscape and function analysis in FZD7-high and low groups
We performed a mutational landscape analysis of the FZD7-high and FZD7-low groups, which differed in mutational frequency (Figure 7A, B). Although CTNNB1, TTN, TP53, MUC16, and ALB were the most common mutations, their relative mutation frequencies differed between the FZD7-high and FZD7-low groups. FZD7-low showed a higher frequency of CTNNB1, TTN, and TP53 with respective frequencies of 37%, 29%, and 24%. In FZD7-high, followed by 36% TP53, 27% TTN, and 17% CTNNB1. The differences in mutation frequencies among different FZD7 expression groups may potentially facilitate the development of individualized HCC targeted agents. Consistent with the results of Figure 4D, here, the expression of FZD7 was negatively correlated with TMB of HCC (Figure 7C). Gene set variation analysis (GSVA) showed that the main enriched pathways in FZD7 high group were “Notch signaling”, “Apical surface”, “Apical junction”, and “TGF beta signaling”. In the FZD7 low group, the main enriched signaling pathways were “Cholesterol homeostasis” and “Coagulation” (Figure 7D). After that, we performed consensus clustering analysis of HCC samples according to the mutated genes (a total of 30 mutated genes were identified) in the FZD-7 high and low expression groups. After k-mean clustering, we finally identified two subtypes in TCGA HCC samples (Figure 7E, S5A). Kaplan-Meier curves indicated that overall survival was superior in the C1 subtype compared with the C2 subtype (Figure S5B). After that, immune infiltration analysis found distinct differences in the immune microenvironment of C1 and C2. Immune microenvironment score, macrophages, NK T cells, and CD4+ Th1 cells were significantly higher in C1 than C2 (Figure S5C).
Differentially expressed genes (DEGs) analysis and construction the prognostic model
We identified the DEGs in different FZD7 groups to understand the molecular mechanisms that regulate tumor microenvironment and prognosis. We found a total of 580 dysregulated genes, including 465 upregulated and 115 downregulated genes (Figure 8A) (Supplemental Table S2), and the heatmap of differential expressed genes in different FZD7 group was displayed in Figure 8B. A total of 25 potentially prognostic related genes were selected based on differential expression genes between FZD7-high and FZD7-low groups by univariate Cox regression analysis (Figure 8C). After that, we included these 25 genes in the LASSO analysis and further obtained 8 genes that determine the prognosis of HCC for risk model construction (Figure 8D)(Table 1). We filtered to eight genes used to construct the HCC prognostic risk model with the following formula: RiskScore=0.0371* LPCAT1+0.0646*MACIR+0.0948*G6PD+0.0046*PTK7+0.0024*NT5DC2+0.0581*S100A9+0.0944*CDC20-0.0333*PON1. HCC patients in TCGA and ICGC-JP were divided into high-risk and low-risk groups based on the median values of the risk score formula. In the LIHC data from TCGA, Kaplan-Meier curves indicated that the low-risk group had a significantly superior prognosis compared with the high-risk group of HCC patients (P=3.44e-07) (Figure S6B). In the TCGA HCC sample, as the risk score for patients with HCC increased, the proportion at high risk increased, as did the proportion of patients with HCC who died (Figure S6A). In TCGA HCC, the values of area under the curve (AUC) in receiver operating characteristic (ROC) at 1, 2 and 3 years were 0.766, 0.705 and 0.702, respectively (Figure S6C). In ICGC-JP, Kaplan-Meier curves indicated that the prognosis of HCC patients in the low-risk group was significantly better than that in the high-risk group (P<0.001) (Figure S6E). Similarly, the proportion of patients with HCC in the high-risk group and the proportion of patients who died of HCC both increased as the risk score increased (Figure S6D). In ICGC-JP, the AUC in ROC at 1, 2 and 3 years were 0.780, 0.698 and 0.699, respectively (Figure S7A). After that, we included clinical variables and riskScore in the univariate and multivariate Cox regression analyses in the TCGA HCC sets. In TCGA HCC, univariate Cox regression analysis indicated that tumor pathological stage (HR=1.879, P<0.001) and riskScore (HR=6.296, P<0.001) were independent risk factors for HCC prognosis. Multivariate Cox regression analysis indicated that stage (HR=1.675, P=0.046) and riskScore (HR=4.930, P<0.001) were independent risk factors for HCC (Figure S7B-C).
Clinicopathological relevance of FZD7 in HCC
To gain further support of the role of the FZD7 in the development and progression of HCC, and possibly other types of cancer, and to extend our observations to a clinicopathologically relevant context, we first examined the protein level of FZD7 by immunohistochemical staining of human tissue arrays. A series of carcinoma samples from esophagus, stomach, colon, rectum, liver, pancreas cancer patients were included, with each type of carcinoma had at least six tumor samples and two adjacent normal tissues. Tissue microarray analysis by immunohistochemical staining showed that FZD7 protein was mainly detected in the cytoplasm of adjacent cells, and that its expression was upregulated in carcinomas of liver, esophagus, stomach, colon, rectum, pancreas, compared with their respective adjacent normal cells (Figure 9A).
To address the clinicopathological relevance of FZD7, we collected 20 samples of HCC paired with adjacent normal liver tissues and analyzed by Real-time PCR for the expression of FZD7. The results showed that the levels of FZD7 were higher in tumor samples than in adjacent tissues (Figure 9B).