Complement C7 (C7), a Potential Tumor Suppressor, Is an Immune-Related Prognostic Biomarker in Prostate Cancer (PC)

Objectives: Prostate cancer (PC) is the second most frequent tumor in men, which has a high recurrence rate and poor prognosis. Therefore, this study aimed to identify novel prognostic biomarkers and therapeutic targets for immunotherapy and small molecule drugs for PC treatment. Materials and Methods: The Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data (ESTIMATE) algorithm was applied to calculate immune scores and stromal scores of TCGA-PRAD data. Differentially expressed genes (DEGs) were identified using R package “limma.” GO, KEGG, and DO analyses were performed to analyze DEGs. Overall survival and disease-free survival analyses were conducted for hub gene identification. To validate the hub gene at the mRNA and protein expression levels, genetic alterations were measured, and CCLE and Cox regression analyses were performed. Connectivity map (CMap) analysis and GSEA were performed for drug exploration and function analysis, respectively. Results: Immune scores ranged from −1795.98 to 2339.39, and stomal scores ranged from −1877.60 to 1659.96. In total, 45 tumor microenvironment (TME)-related DEGs were identified, of which Complement C7 (C7) was selected and validated as a hub gene. CMap analysis identified six small molecule drugs as potential agents for PC treatment. Further analysis demonstrated that C7 expression was significantly correlated with clinical T, pathological N, and immune infiltration level. Conclusions: In conclusion, of the 45 TME-related DEGs, C7 was shown to correlate with PC prognosis in patients, indicating it as a novel prognostic biomarker and immunotherapy target in PC. Additionally, six small molecule drugs showed strong therapeutic potential for PC treatment.


INTRODUCTION
Prostate cancer (PC) is the second most frequent tumor in men, with a high recurrence rate and age dependence (1). In 2018, 1,276,106 new cases of PC and 358,989 associated deaths occurred worldwide (1). Although serum PSA, CT, and magnetic resonance examination have been widely used for PC diagnosis, these methods do not show high accuracy and specificity (2). Thus, novel biomarkers to accurately diagnose PC are needed.
Recently, more studies have shown that immunotherapy could treat cancers effectively and safely (3)(4)(5). With the development of immunotherapy for cancers, more and more researchers focus on finding out more accurate therapeutic targets for immune treatment (6, 7). Thus, in this study, we aimed to screen some novel diagnostic and prognostic biomarkers for PC and therapeutic targets for immunotherapy.
The Cancer Genome Atlas (TCGA), a landmark cancer genomics program, molecularly characterized over 20,000 primary cancer and matched normal samples spanning 33 cancer types. Thus, we retrieved PC samples from TCGA-PRAD (prostate adenocarcinoma) data, which had more data, abundant information and contents characteristics. To the best of our knowledge, this study is the first to apply the Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data (ESTIMATE) algorithm for immune score and stromal score calculation of each PC case from TCGA-PRAD data (8). Then, based on immune scores and stromal scores for PC samples, we attempted to identify and validate some novel immune-related prognostic biomarkers in PC.
In conclusion, our finding indicated that C7 might be a novel prognostic biomarker and therapeutic target for immunotherapy in PC. In addition, six molecule drugs were identified, which showed strong therapeutic potential for PC treatment.

PC Studies Collection and Data Preprocessing
First, gene expression profiles of mRNA of PC samples (TCGA-PRAD data) were retrieved from The Cancer Genome Atlas (TCGA) database (https://genomecancer.ucsc.edu/). Briefly, 495 PC cases with complete clinical information were included in this study. Using the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/), GSE116918 (of GPL25318) (9), GSE16560 (10) of GPL5474, and GSE70770 (11,12) of GPL6801 were downloaded for validation. In total, 248 PCs from GSE116918 with complete survival information, 281 PCs from GSE16560 with survival information, and 204 PCs from GSE70770 with complete clinical information were included in this study. Figure 1 provided a step-by-step flowchart of this study. In brief, all data were preprocessed before analysis. For the TCGA-PRAD data displayed as count number, normalized, and log2 transformation were performed using R package "DEseq.2" (13). Moreover, normalization and transformation of GSE116918 and GSE70770 were conducted using package "affy" (14) in R software, which were shown as raw expression data.

Calculation of Immune and Stromal Scores and Correlation With Clinical Phenotype
Using the ESTIMATE algorithm, immune, and stromal scores of each case from the TCGA database were calculated by R package "estimate" (8). To explore the correlation between immune score or stromal score and clinical phenotype, the immune score levels (or stromal score levels) by gender (female vs. male), neoplasm histological grade (Gx, G1, G2, G3, G4), laterality (left vs. right), pathological stage (I, II, III, IV), and neoplasm cancer status (tumor-free vs. tumor) were assessed using boxplots. Subsequently, 495 PCs were divided into high-and low-immune score groups based on the median immune score. Similarly, cases were divided into high-and low-stromal score groups. Overall survival (OS) and disease-free survival (RFS) analyses were conducted to explore the relationship between immune score (or stromal score) and survival using R package "survival" (15). Unpaired t-test [or One-way Analysis of Variance (ANOVA) test] was performed to measure statistical significance.

Differentially Expressed Gene (DEG) Identification and Hub Gene Identification
Immune-and stromal-related genes were screened to identify DEGs in this study. To detect immune-related DEGs, we first divided 495 PC cases into high-and low-immune groups based on the median immune score, followed by DEG screening using "limma" (16) package in R. Similarly, for stromal-related DEGs, 495 samples were divided into two groups (high-and lowstromal groups), and DEGs were screened by "limma." Genes reaching the standard (Adjust P < 0.05 and | log2 fold change (FC) | ≥ 1.0) were considered as DEGs. Further, overlapping immune-and stromal-related DEGs were selected for subsequent analysis. Gene Ontology (GO) (17) enrichment analysis, Kyoto encyclopedia of Genes and Genomes (KEGG) (18) pathway analysis, and Disease Ontology (DO) (19) were performed through "clusterProfiler" (20) in R software to identify potential functions of DEGs and highly correlated diseases. Gene sets were considered significant when P < 0.05. In order to screen out hub genes among the DEGs, two survival analyses, including OS and RFS, were performed for each DEG by Gene Expression Profiling Interactive Analysis (GEPIA) (http://gepia.cancer-pku.cn/) (21). Only DEGs showing significant P-values in both survival analyses were considered to be correlated with survival and prognosis of patients with PC, which were regarded as hub genes in this study.

Exploration of Hub Gene Expression in Normal Tissues and Tumors at Different Stages
First, we compared expression levels of hub genes in tumors and normal tissues by GEPIA. Based on GSE116918 and GSE70770 with complete T stage information, T stage (T1, T2, T3, and T4) boxplots were generated (GSE116918). Additionally, clinical T stage (T1, T2, and T3) and pathology T stage (T1, T2, and T3) boxplots were plotted. Unpaired t-test and One-way Analysis of Variance (ANOVA) test were conducted to measure statistical significance.

Hub Gene Protein Expression and Survival Validation and Genetic Alteration
Following selection of hub genes, their mRNA expression levels in PC and normal tissues were compared using the Oncomine database (https://www.oncomine.org/) (22). mRNA expression and copy number variation (CNV) levels of hub genes in 40 tumor types were explored using the CCLE database (https://portals.broadinstitute.org/ccle/). Based on GSE116918, biochemical failure-free survival (BFFS) and metastasis free survival (MFS) analyses for hub genes were performed to validate the prognostic value. Similarly, based on GSE16560, OS of patients with PC was assessed. Further, protein expression levels of hub genes were validated using The Human Protein Atlas database (https://www.proteinatlas.org/) (23). Student t-test was used to measure statistical significance. By using the CBio Cancer Genomics Portal (http://www.cbioportal.org/) (24,25), genetic alteration of hub genes was examined. Combined with relative mRNA expression levels, the correlation with copy number variations (CNVs) of hub genes was explored. Following, cases were divided into with and without alterations groups, and survival analyses (OS and RFS) were conducted to explore the association between mutations and survival of PC patients.

Associations Between Hub Genes and Clinical Features of PC Patients
Based on TCGA-PRAD data, 495 PCs were divided into high-(n = 248) and low-(n = 247) expression groups by evaluating the median value of hub gene expression. Subsequently, χ2-test or ANOVA was performed to measure the relationships between hub gene expressions and clinical features (age, laterality, clinical M, clinical T, pathological N, pathological T, and neoplasm cancer status). Expression levels of hub genes and clinical features (age, laterality, Gleason score, clinical M, clinical T, pathological N, pathological T) were also included for univariable Cox analysis of overall survival (OS) using TCGA-PRAD data. Features with P < 0.05 were included for multivariate Cox analysis.

Correlation Between Hub Genes and Immune Microenvironment
Using TIMER (https://cistrome.shinyapps.io/timer/) (26), the correlation between hub gene expression and the abundance of immune infiltrates was explored. Six types of tumor-infiltrating immune cells (B cells, CD8+ T cells, CD4+ T cells, macrophages, neutrophils, and dendritic cells) were included in this analysis. Hub genes with |correlation coefficient (cor) | ≥ 0.3 and P < 0.05 were considered to be significantly related to the immunocyteinfiltrating level. Immune infiltrating levels of tumor-infiltrating immune cells with different CNVs were compared. Moreover, survival analysis was conducted based on high-and lowinfiltrating levels of immunocytes.

Gene Set Enrichment Analysis (GSEA)
GSEA (27) was performed using TCGA-PRAD data to explore the potential functions of hub genes. Briefly, 495 PCs were divided into two groups (high-and low-expression groups) based on the median of hub gene expression levels. Further, "c2.cp.kegg.v7.0.symbols.gmt" was chosen as the reference gene sets. In this study, only KEGG signaling pathways reaching P < 0.05, |ES| > 0.6, gene size ≥100, and FDR < 25% were considered to be significantly enriched.

Identification of Candidate Molecule Drugs
Connectivity map (CMap) (28) is an online webtool to identify molecule drugs showing high correlation to diseases (https:// portals.broadinstitute.org/cmap/), which was applied to the overlapped immune-and stromal-related DEGs in this study to explore molecule drugs with strong correlations to PC. Molecule drugs with the number of instances (n) >5 and P < 0.05 were removed. Moreover, small molecule drugs with |mean| ≥ 0.40 were considered to possess significant therapeutic potential to treat PC.  (Table S1). Further analysis demonstrated that the immune score was significantly associated with pathological T (F = 3.570, P = 0.027, Figure 2A), laterality (F = 6.020, P = 0.005, Figure 2B), pathological N (t = −2.000, P = 0.048, Figure 2C) and neoplasm cancer status (t = −2.660, P = 0.009, Figure 2E). Moreover, stromal score was significantly related to pathological T (F = 3.710, P = 0.024, Figure 2F) and neoplasm cancer status (t = −3.510, P = 0.001, Figure 2J). Unfortunately, stromal score did not show significant relationship with laterality ( Figure 2G), pathologic N ( Figure 2H), and clinical M ( Figure 2I). RFS survival analyses showed a trend of worse RFS in PC patients with high immune scores (P = 0.110, Figure S1B) or high stromal scores (P = 0.190, Figure S1D). No association between immune scores or stromal scores was found with OS in PC patients (Figures S1A,C).

Hub Gene Validation
Using PC data from the Oncomine database (including 12 PC datasets), C7 mRNA expression tended to be lower in PCs than normal tissues ( Figure S3). As shown in Figure 6A, no significant difference was found in C7 protein expression between PCs and normal tissues. In evaluating hub gene expression in various cancer cell lines by GEPIA and CCLE, C7 expression was found to differ in these tumor types. C7 expression was significantly lower in 15 cancer types compared with corresponding normal tissues ( Figure S4A). Moreover, C7 mRNA expression in the prostate ranked as the bottom fifth across all the types of cancer cell lines ( Figure S4B). Further, the CNV level of C7 in the prostate was the fifteenth highest of all cancer types ( Figure S4C). We further validated the prognostic value of C7 by performing survival analysis. PC patients with high expression of C7 had better OS (P = 0.037, Figure 6B). Similarly, low expression of C7 caused worse BFFS (P = 0.048, Figure 6C) and MFS (P = 0.049, Figure 6D), indicating that C7 was significantly associated with survival and prognosis of PC patients.

C7 Genetic Alteration
As shown in Figure 7B, C7 was altered in 26 (19%) of 136 PC patients ( Figure 7A). The main type of alteration was missense mutation (Figure 7B). Figure 7C showed the network containing nine nodes, including C7 and 8, consisted of the most altered neighbor genes. Currently, no drugs target C7, suggesting that C7 might be a novel therapeutic target for PC treatment.

C7 Enrichment in Nine KEGG Signaling Pathways
To explore potential functions of C7, GSEA was conducted, showing that C7 was significantly correlated to nine KEGG     Table 2).

Strong Therapeutic Potential Shown by Six Small Molecule Drugs
Highly associated molecule drugs were identified by CMap. In total, 12 molecule drugs were screened ( Table 3)

DISCUSSION
Tumor microenvironment (TME), the internal environment of the tumor for growth and proliferation, consists of noncancerous cells (29), including fibroblasts, immune cells, cells that comprise the blood vessels, and secreted proteins produced by all cells present in the tumor (30). Two major cell types (immune and stromal cells) in the TME have been reported to have strong correlations with tumor prognosis (30). With the development of precision medicine, immunotherapy has been one of the most effective treatment methods for malignancies, which mainly stimulates autoimmune or alloimmune cells in patients to improve symptoms, prolong survival, and improve prognosis (31,32). Prostate cancer (PC) is one of the most common malignancies in the urinary system and the second most common tumor in men. Currently, surgical resection plus radiotherapy and hormone therapy are the standard treatment methods. Although this method is effective for early and localized PC treatment, PC develops quickly, leading to metastasis. Also, nowadays there are some treatment available for advanced PC, and nearly all have resulted in prolongation of survival. Unfortunately, a recent study has shown that the 5-year survival rate of PC was <30% (33). Because of the poor prognosis of PC, at the same time taking into account the effectiveness and safety of immunotherapy, exploring novel prognostic biomarkers and therapeutic targets for immunotherapy is needed.
In this study, we applied the ESTIMATE algorithm to calculate immune and stromal scores of PC cases from TCGA-PRAD data. Following, 45 overlapping immune-and stromalrelated DEGs were chosen for further analysis. Among them, only one gene C7 was significantly associated with OS and RFS of PC patients, which was considered as a hub gene in this study. Follow-up analysis confirmed that C7 expression was lower in PCs compared to normal tissues. Moreover, Oncomine analysis demonstrated that C7 mRNA expression was higher in normal tissues compared to that in PCs, which was consistent with the TCGA-PRAD data. Furthermore, based on GSE29609, results suggested that C7 expression was significantly associated with OS and BFFS (and MFS), demonstrating that C7 was a powerful prognostic biomarker in PC.
Because immune infiltration level showed a strong correlation with survival in tumors (34), we explored the association between C7 expression level and immune infiltration level in PC by TIMER. C7 expression was found to be significantly negatively related to tumor purity while positively associated with macrophages and dendritic cells neutrophils. Thus, C7 might be a potential therapeutic target for immunotherapy in PC.
C7 can form a membrane attack complex together with complement components C5b, C6, C8, and C9 (35). This complex worked as part of the terminal complement pathway of the innate immune system (35). Moreover, C7 expressed on the cell membrane acted as a regulator of the excessive proinflammatory reaction as reported (36). In addition, Li et al. demonstrated that low expression of C7 in non-small cell lung cancer (NSCLC) might act as a tumor inhibitor, which was associated with tumor progression and prognosis (36). In detail, they firstly found the C7 mRNA content descended gradually in normal, benign, borderline and malignant ovarian tissues based on the fluorescent in situ hybridization assay (Consistent downward tendency of C7 mRNA expression was also found in the course of lung tumor progression). Then they attempted to explore the prognostic potential of C7 by using survival analysis. Unfortunately, they only found that NSCLC patients with decreased expression of C7 had a worse outcome. Finally, some experiments were done to demonstrate that C7 overexpression inhibited proliferation of NSCLC cells in vitro. Although our study was not based on experiments like theirs, similarly, in our study, we found low expression of C7 in PC compared to normal tissues by using some authoritative public databases including GEO, Oncomine, TCGA. Furthermore, we demonstrated that PC patients with decreased expression of C7 had a worse outcome by the public databases we could use. Therefore, we hypothesized that C7 might be a potential tumor suppressor associated with survival and prognosis of PC and might be a new therapeutic target for PC immunotherapy.
Despite the study design and strict filtering criteria at each step, there were some limitations in this study. First, although we validated the hub gene systematically and comprehensively, C7 did not show significant differences at the protein level. Thus, protein expression will be assessed by western blotting in future studies. Second, although we tried our best to validate our findings by using the public databases we could use, this study was lack of experiment. Thus, we will conduct C7 overexpression and colony formation assay in PC and investigate the cell derivation of C7 (local or by infiltrating neutrophils) in the shore future. Third, no clinical trials have been conducted to assess the prognostic value of C7 in this study, which is the focus of upcoming studies.

CONCLUSION
In summary, we identified 45 TME-related DEGs (associated with immune and stromal scores) by applying the ESTIMATE algorithm. A potential immune-related prognostic biomarker named C7 was further screened and validated, which might function as a tumor suppressor in prostate cancer. Moreover, six small molecule drugs were identified that showed strong therapeutic potential for PC treatment.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are openly available in Gene Expression Omnibus (GEO) database at http:// www.ncbi.nlm.nih.gov/geo/ and The Cancer Genome Atlas (TCGA) database at https://genome-cancer.ucsc.edu/.

ACKNOWLEDGMENTS
We would like to acknowledge the TCGA and GEO database developed by the National Institutes of Health (NIH).

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc. 2020.01532/full#supplementary-material     Table S1 | Stromal score and immune score of each patient with PRAD (based on TCGA-PRAD). Table S2 | Differentially expressed genes (DEGs) between low-immune samples and high-immune samples. Table S3 | Differentially expressed genes (DEGs) between low-stromal samples and high-stromal samples.