Identification of driving genes of familial adenomatous polyposis by differential gene expression analysis and weighted gene co-expression network analysis

2.1Gene expression datasets and clinical pathological data

GSE111156, a microarray-based dataset, consists of 48 samples (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE111156). In GSE111156 research, a genome-wide transcriptional analysis of the duodenal samples from patients with FAP was conducted to characterize the alterations that occur in the duodenal adenoma-carcinoma sequence in FAP. Duodenal samples from 12 patients with FAP who did not have duodenal cancer (controls) and 12 patients with FAP who were diagnosed with duodenal cancer from 1988 to 2014 (cases) were analyzed for gene expression. Extraction of RNA from archived cancer, adenoma as well as normal tissue from cases and from archived adenoma tissue from controls. Profiling and analysis of gene expression applying the Affymetrix Human Transcriptome Array 2.0. Three pairwise comparisons were implemented to characterize the adenoma-carcinoma sequence: normal-cancer (No-Ca), adenoma-cancer (Ad-Ca), and normal-adenoma (No-Ad). A total of 48 samples were evaluated, 12 from controls (adenoma) and 36 from cases (cancer, adenoma, normal).

2.2Identification of DEGs

DEGs was determined utilizing open-source software R language (Version 3.3.3, https://www.r-project.org/) together with R package of Bioconductor (http://www.bioconductor.org/) limma package [34, 35]. Screened DEGs at p< 0.05 and logFC = 1.

2.3Construction of WGCNA

Co-expression networks of filtered genes were constructed through the usage of “WGCNA” R package [36, 37]. After assessment of the expression matrix assessed via the average method with the “hclust” function, the clustered genes of gene chips comprising GSM3024204, GSM3024195 and GSM3024196 were defined as biased and were therefore excluded from the further analysis (Fig. 1A). Other samples were applied for calculating the Pearson’s Correlation Matrices. Weighted adjacency matrices were calculated via the formula amn =|cmn|β (amn: the adjacency between gene m and n, cmn: Pearson’ correlation, β: soft power threshold). Subsequently, data on clinical features were loaded and average connectivity and scale independence were evaluated. Moreover, topological overlap measure (TOM) matrices transformed via adjacency matrices were employed for the estimation of their connectivity. Furthermore, the corresponding genetic information was extracted for each module for the further analysis.

2.4Identification of clinically significant module

The traits matrix contains Adenocarcinoma, Adenoma, Normal (case cancer, CaseC), and Control_ademoma 4 clinical traits. After the identification of modules, the module eigengene (ME) was concluded through the first principal component of the expression level of module. Module-trait correlations were estimated utilizing association between clinical traits and ME, which enabled efficient recognition of related modules. To assess the strength of correlation, we measured the module significance (MS), which is determined as the mean absolute gene significance (GS) of all genes that participated in the module [38]. GS was calculated as a log10 transformation of the P-value (lgP) in a linear regression between clinical information and gene expression. The key module was defined as the modules that exhibited the highest correlation (positively and negatively) coefficients among all modules [39].

2.5Tumor immune cell infiltration analysis

The calculation functions for TIP were created in R [40, 41, 42]. The JavaScript, Struts2 and Bootstrap framework were employed to communicate between the web interfaces and R, which were subsequently implemented via JavaScript, and the table results were produced from DataTables. Additionally, bar, heat, scatter and pie plots were produced with HighCharts, while D3 powered box, line, circle and radar plots. The source code for https://github.com/dengchunyu/TIP. In general analysis, TIP receives four categories of expression data from Case tissue as input. One is RNA-seq expression data (TPM or raw counts), and the other is microarray expression data (non-log or log transformed). Phenotyping of tumor cell infiltration was carried out online with the usage of TIP domain (http://biocc.hrbmu.edu.cn/TIP) in accordance with the derived data GSE111156. The samples were divided into CtrlA: Control adenoma tissue; CaseN: case normal tissue; CaseA: case adenoma tissue and CaseC: case cancer tissue group. The analysis followed default parameters.

2.6Identification of TFs in the key module

Enrichr (http://amp.pharm.mssm.edu/Enrichr/) is a comprehensive tool based on the Web [43, 44]. Genetic information of key module was imported into Enrichr to capture the interaction between transcription factors (TFs) and the associated target genes. To minimize the possibility of detecting false positives, we only retrieved TFs from the ENCODE and ChEA gene-set libraries with consensus targets and identified P< 0.05 through Fisher exact test. Later, Cytoscape 3.4.0 software (Cytoscape Consortium, SanDiego, CA, USA) was applied for visualizing the regulatory networks of TF-target gene.

2.7Functional enrichment analysis

Enrichment analysis of GO [45] and KEGG [46] pathway utilizing clusterprofiler package in the R software (Version 3.3.3) [47] for investigating the underlying mechanism of DEGs, the genes in modules of interest, as well as a group of hub genes. FDR < 0.1 was established as the cutoff value. The clusterProfiler is a R package based on ontology that can automate the bio-term classification and the gene cluster enrichment analysis, simultaneously delivers a module for visualization to display the analysis results [48].

2.8Quantitative Real-Time PCR (qRT-PCR)

Extraction of total RNA was performed from 12 pairs of FAP TR and IT tissues through TRIzol reagent (Invitrogen) and reverse transcribed to cDNA. The qRT-PCR was implemented on the ABI PRISM 7500 Real-Time PCR System (Applied Biosystems, Foster City, CA). Primers for CD2, CD3G, CCL5, KLF2 and TLR3 are presented in Table 2. GAPDH was served as the internal control. mRNAs relative expression levels were computed employing the 2-Δ⁢Δ⁢CT approach. All of the experiments were replicated a minimum of 3 times. Our research was authorized by the ethics committee of the First Affiliated Hospital of Kunming Medical University. All the patients provided written informed consent.

3.1Identification of DEGs in FAP and their functional

The differential analysis was performed using GSE111156 data. By comparing case normal tissue (Case N) with control adenoma tissue (Ctrl A), it was possible to determine 63 differentially expressed genes (DEGs) in the sample tissues, of which 3 were raised and 60 were decreased (Table 1). By comparing case adenoma tissue (Case A) with Ctrl A, it was possible to determine the 23 DEGs in the sample tissues, of which 1 gene was raised and 22 were decreased. By comparing case cancer tissue (Case C) with Ctrl A, it was possible to determine 220 DEGs in the sample tissues, of which 44 genes were raised and 176 were decreased. By comparing Case A with Case N, it was possible to determine 10 DEGs in the sample tissues, of which 2 genes were raised and 8 were decreased. In contrast to Case N, 196 DEGs were identified in the samples from Case C, 136 of which were up-regulated and 60 were down-regulated. Then, the number of 20 DEGs were screened in Case C compared with Case A, of which 12 were raised and 8 were decreased.

Figure 1.

Sample clustering was conducted to detect outliers.

Figure 2.

The results of DEGs among different groups. (A) The volcano plot of DEGs in Case N vs Ctrl A. (B) The volcano plot of DEGs in Case A vs Ctrl A. (C) The volcano plot of DEGs in Case C vs Ctrl A. (D) The volcano plot of DEGs in Case A vs Case N. (E) The volcano plot of DEGs in Case C vs Case N. (F) The volcano plot of DEGs in Case C vs Case A. Red points indicate up-regulation genes, blue points indicate down-regulation genes.

Figure 3.

Identification of hub genes associated with FAP. (A) The results of clustering analysis based on genes in Case N vs Ctrl A. (B) T The results of clustering analysis based on genes in Case A vs Ctrl A. (C) The results of clustering analysis based on genes in Case C vs Ctrl A. (D) The results of clustering analysis based on genes in Case A vs Case N. (E) The results of clustering analysis based on genes in Case C vs Case N. (F) The results of clustering analysis based on genes in Case C vs Case A.

Figure 4.

The gene ontology enrichment and KEGG pathway classification of DEGs. It contains three categories including biological process, cellular component and molecular function. The abscissa axis indicates the number of enriched genes. The abscissa represents the annotated genes in KEGG database; the ordinate represents Categories in KEGG database.

Figure 5.

Global visualization of the immunophenotype across samples from TCGA. (A) The heatmap of 23 normalized immune activity scores in 48 samples. (B) The relative proportion of tumor-infiltrating immune cells. (C) Left, the expression pattern of signature genes from the seven-step cancer-immunity cycle. Each row represents a single gene, and each column represents one sample. Right, the principal component analysis (PCA) of signature genes expression for all samples.

Figure 6.

Weighted gene co-expression network analysis. (A) Sample dendrogram and trait heatmap. (B) Scale independence and mean connectivity of various soft-thresholding values (β). (C) Dendrogram of all filtered genes enriched according to a dissimilarity measure (1-TOM) and the cluster module colors. (D) Heatmap of the correlation between the clinical traits and MEs of FAP. The darker the module color, the more significant their relationship. (E) Eigengene adjacency heatmap.

Figure 7.

(A) The scatter plot between the pink module membership and the gene significance for normal. (B) Functional enrichment of genes in the pink module, including Biological process GO terms, Cellular component GO terms, Molecular function GO terms and KEGG analysis.

Figure 8.

The identification and further analysis of hub genes in pink module. (A) Genes network with 381 edges and 79 nodes, red points represented hub genes [20]; (B) Boxplots for the expression of hub genes in normal vs. FAP; (C) Functional annotation for hub genes associated with control ademona, including Biological process GO terms, Cellular component GO terms, Molecular function and KEGG analysis.

Figure 9.

Identification and functional enrichment analysis of transcriptional factors in pink module. (A) The transcriptional factor (TF) network of the 80 genes generated by Webgestalt. The diamond and red color indicate the TFs while the green and eclipse shape indicate the target gene. (B) The GO and KEGG enrichment analysis for TF in pink module.

Figure 10.

(A) The scatter plot between the midnightblue module membership and the gene significance for control ademona. (B) Functional enrichment of genes in the pink module, including Biological process GO terms, Molecular function GO terms and KEGG analysis.

Figure 11.

The identification and further analysis of hub genes in midnight blue module. (A) Genes network with 136 edges and 17 nodes, red points represented hub genes [17]; (B) Functional annotation for hub genes associated with control ademona, including Biological process GO terms and KEGG analysis; (C) Boxplots for the expression of hub genes in control ademona vs. FAP.

Figure 12.

Identification and functional enrichment analysis of TFs in midnightblue module. (A) The transcriptional factor (TF) network of the 48 genes generated by Webgestalt. The diamond and red color indicate the TFs while the green and eclipse shape indicate the target gene. (B) The GO and KEGG enrichment of analysis in midnightblue module.

Figure 13.

The most significant correlation module merge Immune infiltration signature genes, pink (CCL5, CD3G, CD2, TLR3). The GO and KEGG enrichment of four identified genes. The GO item contains (A) biological process, (B) cellular component and (C) molecular function. (D) KEGG pathway analysis.

Figure 14.

The most significant correlation module merge Immune infiltration signature genes: midnightblue (KLF2). The GO and KEGG enrichment of the KLF2. The GO item including (A) biological process, (B) cellular component. (C) KEGG pathway analysis.

Figure 15.

Validation of genes in FAP TR and IT tissues, including CCL5, CD3G, CD2, TLR3 and KLF2.

Figure 2A displays the volcano charts of DEGs expression variations in Case N and Ctrl A. The signal distribution of the plot indicates specific up-regulated (PCK1, SLC15A1, MIR2054) and down-regulated (CEACAM6, POF1B, RNF 43, NQO1, CD44, ANXA10, SNORD123, SLC12A2, CERS6, KIAA1324) transcripts. As shown in Fig. 2B, the differentially expressed genes between Case A with Ctrl A tissue, the specific up-regulated gene contained SNORA76, and down-regulated genes included FABP1, LOC399753, CPS1, DMBT1, AGAP7, ADH1C, SI, MIR548O2, LOC101060495, MIR548X). Figure 2C suggested that the expression of some genes was increased (GREM1, VTRNA2-1, MIR2054, MIR4759, THBS1, VTRNA1-3) whereas that of others was downregulated (ADH1C, CPS1, HMGCS2, PBLD, KIAA1324, CPS1-IT1, ADH4, ADH1A, EPHX2, SLC4A4) between Case C with Ctrl A. Compared Case A with Case N, the DGEs contained up-regulation (CEACAM6 and REG4) and down-regulation (RBP2, GSTA1, GSTA2, ALDOB, SLC15A1, ACE2, MEP1B and AKR1B10) genes (Fig. 2D). Then, compared Case C with Case N, the signal distribution of the plot indicates specific up-regulated (RBP2, ALDOB, APOB, GIP, GSTA1, GSTA2, FOLH1, FOLH1B, CREB3L3, SST) and down-regulated (CEACAM5, SLC2A1, MYOF, CEACAM6, PKM, SNORA71C, SNORA74A, GDF15, SPP1, SULF1) transcripts (Fig. 2E). The differentially expressed genes between Case C with Case A tissue, the specific up-regulated genes contained SULF1, VCAN-AS1, SPP1, POSTN, SNORD114-12, COL3A1, COL1A1, COL1A2, SPARC and GREM1, and down-regulated genes included DEFA5, DEFA6, ADH1C, HMGCS2, KIAA1324, MT1H, MT1G, CPS1 (Fig. 2F).

Hierarchical clustering heatmap was applied to reveal the DEGs profiles among two group, displayed differences in the mode of gene expression between both groups (Fig. 3). Figure 3 presents that these genes could distinctly separate all samples into two groups.

In order to understand their functionality, all of the DEGs were plotted against the terms of the GO database and subsequently compared to the entire reference database. The GO analysis results for the DEGs are presented in Fig. 4A. They were categorized into three groups, namely, the Cellular Component (CC), the Molecular Function (MF), together with the Biological Process (BP). It was revealed that these genes were linked to mRNA catabolic processes, neutrophil degranulation, as well as the neutrophil activation involved in the immune response. For further insight into the bio-functions of the genes determined and to categorize their functional annotations, remarkably enriched pathways were defined through comparison with KEGG database. When comparing case N with Ctrl A, the DEGs were enriched to Ribosome and PPAR signaling pathway. As shown in Fig. 4B, the findings presented that the genes were correlated with catabolic and metabolic process, annotated to digestion and absorption pathways. When comparing Case C vs Ctrl A, DEGs were associated with immune response, annotated to chemical carcinogenesis and metabolism (Fig. 4C). When comparing Case A vs Case N, DEGs were associated with catabolic and metabolic and annotated to metabolism pathways (Fig. 4D). Notably, when comparing Case C vs Case N, DEGs were associated with leukocyte migration (Fig. 4E), while the genes were associated with extracellular organization when comparing Case C vs Case A (Fig. 4F). Some pathways that are strongly associated with immunity were markedly enriched, for instance the “leukocyte migration and acute inflammatory”, which warranted to be further studied.

3.2Immunophenotype profiling

To reveal the status of immune infiltration in FAP, we tracked, analyzed and later validated the anti-cancer immune status and the ratio of tumor-infiltrating immune cells in the 7-step Cancer-Immunity Cycle through the Tracking Tumor Immunophenotype (TIP) website utilizing samples from patients with FAP. Figure 5A exhibited that in contrast to other groups, the scores of the first to third Cancer-Immunity Cycle steps corresponding to cancer cell antigen release (Step 1), the presentation of cancer antigens (Step 2), as well as excitation and activation (Step 3) were increased in Case C. The step 4, which corresponds to the trafficking of immune cells to the tumors showed significant alterations between the samples in each group. Specifically, immune activity scores corresponding to dendritic cell recruiting, Th1 cell recruiting, T cell recruiting, Treg cell recruiting, Th2 cell recruiting, and B cell recruiting raised, while scores corresponding to Th22 cell recruiting, CD8 T cell recruiting, CD4 T cell recruiting, neutrophil recruiting, monocyte recruiting, macrophage recruiting, eosinophil acidophil recruiting, NK cell recruiting, Th17 cell recruiting and basophil recruiting correspondingly declined significantly in both Case A and Case C. In FAP samples, immune cell infiltration of the tumors (Step 5) was diminished, while T cell recognition of cancer cells (Step 6) and killing of cancer cells (Step 7) revealed elevated immune activity scores. Analysis of the infiltration of immune cells in multiple samples (Fig. 5B) suggested that FAP samples mainly infiltrated by CD8 T cells, memory B cells, monocyte and naive B cells. The expression of signature gene on samples was depicted in the heatmap (Fig. 5C). With 48 samples, we identified signature genes participating in the cancer-immune cycle that exhibited diverse patterns of expression among Case A, Case C, Case N and Ctrl A samples (Fig. 5D)

3.3WGCNA analysis for identification of gene modules associated with FAP

After applying the “cluster” function to the expression matrix assessed via the averaging approach [20], a total of four samples including GSM3024207, GSM3024217, GSM3024218 and GSM3024211 whose cluster height above 0.05 were considered as deviations and were therefore excluded from the further analysis (Fig. 6A). R package “WGCNA” was employed to group similarly expressed genes into distinct modules. In this research, β= 8 was taken to ensure low mean connectivity (close to 0) and high scale independence (close to 0.9). The phase anisotropy of modules was set to 0.2, and in total, 10 modules were produced (Fig. 6C). Figure 6D presents the eigengene neighbor-joining heatmap, suggesting that the pink and black modules, as well as some other modules were adjacent to each other. The module-trait association is illustrated in (Fig. 6E). The pink module related to FAP was the darkest (cor = 0.55, P= 6.6e-29) and this module was selected for the further analysis (Fig. 7A). As illustrated in Fig. 7B, GO analysis revealed that these genes were linked to cofactor binding, fatty acid metabolic processes, and organic anion transport. KEGG pathway enrichment analysis indicated that chemical carcinogenesis, cytochrome P450 metabolism of xenobiotics and drug metabolism-cytochrome P450 signaling pathways were the most predominantly enriched.

3.4Construction of co-expression network and identification of hub genes in pink module

In pink module, the edges marking the correlations were filtered with the weight value > 0.2, and in total, 381 edges were derived. 79 nodes were determined after entering them into Cytoscape (Fig. 8A). Therefore, the boxplots proving the associations between the hub genes and the FAP were displayed in Fig. 8B, including SLC2A2, ACE2, ABCG2, ASAH2, MGAM, MS4A10, MEP1A, ASAH2C, TMPRSS15, SLC15A1, ME-AS1, MEP1B, CDHR5, CYP3A, SLC6A19, ADA, GBA3, APOB, MME and SLC7A9, 20 hub gene were significant in distinguishing FAP and normal tissue. The 20 hub genes were performed functional enrichment analysis (Fig. 8C). GO analysis presented that hub genes were related to lipid catabolic process and digestion. KEGG pathway enrichment analysis indicated that protein digestion and absorption signaling pathways was the most predominantly enriched.

Furthermore, we entered these 161 genes into Enrich with a P< 0.01 and identified four TFs, namely NFE2L2, HNF4A, BSR1 and FOXA1. Networks corresponding to these genes were displayed in Fig. 9A. Then performed functional enrichment analysis in Fig. 9B. GO analysis revealed that these TFs were most obviously linked to organic anion transport. KEGG pathway enrichment analysis demonstrated that platinum resistance together with drug metabolism-other enzyme signaling pathways were the most significant enrichment sites.

3.5Function annotation of the midnight blue module of interest and the hub genes

Midnight blue module related to FAP was the darkest (cor = 0.57, p= 1.6e-05) and was selected for the further analysis (Fig. 10A). Figure 10B exhibited the KEGG pathway and GO analysis of the 50 significant genes in midnight blue module. GO analysis indicated that these hub genes were linked to responses to steroid, oxidation, peptide and toxic substance. Enrichment analysis of the KEGG pathway suggested that the MAPK signaling pathway was the primary enrichment pathway.

3.6Construction of co-expression network and identification of hub genes in midnight blue module

In midnight blue module, the edges marking the correlations were filtered with the weight value > 0.2, and in total, 136 edges were derived. 17 nodes were determined after entering them into Cytoscape (Fig. 11A). The 17 nodes were performed GO analysis and KEGG, which was associated with p38MAPK cascade and participated in MAPK signaling pathway (Fig. 11B). Figure 11C showed the mRNA level of 17 nodes (MYADM, DUSP1, KLF2, NR4A1, PLK3, MMP19, JUND, PPP1R15A, LMNA, OTTHUMG00000060830, MIDN, NFIL3, EMP1, MAP2K3, EPHA2, PFKFB3, AFP36) were significantly upregulated. Furthermore, we analyzed TFs in midnight blue module (Fig. 12A), 43 TFs were identified and used for performing functional enrichment analysis. The findings revealed that TFs were also correlated with the MAPK signaling pathway (Fig. 12B).

At last, we merged the genes of the most significant correlation module with immune infiltration signature genes to obtain the potential biomarkers for FAP pathogenesis. For the most positive module (pink module), in total, four genes were eventually acquired which including CCL5, CD3G, CD2 and TLR3 (Fig. 13). GO analysis demonstrated that these four genes were linked to protein self-binding and T cell activation. Enrichment analysis of KEGG pathway suggested that the Toll-like receptor signaling pathway was the predominant enrichment pathway. CD3G is also a gene in pink module correlated to HNF4A, which is one of the types of TFs.

For the most negative module (midnight blue module), we identified only one gene (KLF2), which is also a hub gene in midnight blue correlated to Ctrl A (Fig. 14). KEGG pathway enrichment analysis indicated that FoxO signaling pathway was the most predominantly enriched.

The outcomes of qRT-PCR exhibited that among the genes of cancer-associated pathway, CD2, CD3G, CCL5, KLF2 and TLR3 were all dramatically upregulated in the FAP IT tissues versus tumor rim (TR) tissues (Fig. 15), pointing to their effects in the FAP progression.