Immune-Related Biomarkers Associated with Lung Metastasis from the Colorectal Cancer Microenvironment

Abstract

Immune-associated biomarkers can predict lung metastases from colorectal cancer. Differentially expressed genes (DEGs) were screened from sample data extracted from gene expression omnibus (GEO) database. The DEGs were screened from the lung metastasis (LM) and primary cancer (PC) groups of the Moffitt Cancer Center cohort dataset. Then, the tumor immune microenvironment and abundance of immune cell infiltration analyses were performed, and the immune-related DEGs were retrieved. In addition, the transcription factor (TF)-miRNA-mRNA network was constructed and enrichment analyses of the immune-related DEGs and upregulated and downregulated DEGs were carried out. Then, the protein–protein interaction (PPI) network was conducted and the drug–gene interaction was predicted. A total of 268 DEGs were screened. The Immune_Score of samples in the LM group was significantly higher compared with the PC group. The infiltration ratio of M0 macrophages and M2 macrophages of samples was higher than others. A total of 54 immune-related DEGs in M0 macrophages were screened. Moreover, the TF-miRNA-mRNA network was constructed among 8 miRNA-mRNA and 50 TF-mRNA, and the secreted phosphoprotein 1 was regulated by 12 TFs, and the oxidized low-density lipoprotein receptor 1 was regulated by 3 miRNAs and 3 TFs. The TF SAM pointed domain containing ETS TF was also a downregulated DEG. The Kyoto Encyclopedia of Genes and Genomes pathway analysis showed that the DEGs in the TF-miRNA-mRNA network were mainly involved in the interleukin-7 signaling pathway and cell adhesion molecules. In total, 23 protein interactions in this PPI network of M0 macrophage cells were involved in 27 mRNAs. There were 38 drug–gene interactions of immune-related DEGs of M0 macrophage cells predicted to contain 34 small molecule drugs and 8 mRNAs. Finally, the CON cohort dataset verified that the infiltration ratio of M0 and M2 macrophages of the samples was higher.

Introduction

Colorectal cancer (CRC) is the third most common malignancy in the world (Hadad and others 2019; Sung and others 2021). In China, the incidence of CRC is increasing with changes in people's living and eating habits. The treatment methods for CRC include surgical resection, chemotherapy, radiotherapy, and biotherapy, among which surgical resection is the primary treatment method (Iplik and others 2018). Despite the surgical resection rate of 50%–70%, the 5-year survival rate for colon cancer is only ∼50% due to recurrence or metastasis (Sargent 2005).

Recurrence and metastasis of CRC generally occur within 3 years of diagnosis, and the most common organs of metastasis are the liver and lungs. Studies have shown that the occurrence of liver metastasis has been more than 20% in CRC patients (Kamal and others 2019). In addition, liver metastasis has also been proved to be the worst prognostic factor for CRC (Raoux and others 2020).

Therefore, there are many studies on liver metastasis of CRC. However, unlike liver metastasis, lung metastasis (LM) of CRC has not gained significant attention in clinical studies (Hwang and others 2010). In fact, LM is also an important factor influencing the prognosis of CRC patients. Previous literature has reported that about 5%–15% of CRC patients will have LM. The incidence of LM of CRC was 10%–20%. Moreover, lung metastases have a poorer response to treatment and prognosis than liver metastases. Therefore, early identification of lung metastases is crucial to improve the prognosis of CRC.

Several studies have investigated the possible molecular mechanisms of CRC metastasis initiation and development to identify novel therapeutic targets for CRC metastasis. For instance, Tan and others (2020) illustrated that the silencing of brain-expressed X-linked 2 (BEX2) promotes CRC metastasis through the Hedgehog signaling pathway. In addition, Tamjidifar and others (2021) suggested that the relative expression of miR-506 and the SPON 1 gene could be considered diagnostic or predictive biomarkers for CRC. According to the competing endogenous RNA hypothesis, different RNA transcripts participate in pathological processes primarily by competing for the binding sites of shared miRNAs. Lei and others (2020) uncovered that ELFN1-AS1 accelerates the proliferation and migration of CRC by regulation of the miR-4644/TRIM44 axis.

In recent years, a large number of studies have shown that disorders of immune mediators, such as immune cells and cytokines, play an important role in tumor genesis, development, and metastasis. For instance, Huang and others (2020a) showed that CDC42 and CMTM6 are immune-related genes exhibiting strong predictive power in melanomas. Yuan and others analyzed changes in the tumor microenvironment caused by CRC lymph node metastasis, which is of great significance for understanding tumor immune escape, and also provides assistance for patients with personalized treatment (Huang and others 2020b). Sun and others (2017) found that CRC cells suppress CD4 T cell immunity through canonical Wnt signaling.

GSE131418 was analyzed by Kamal and others (2019), and their results showed transcriptomic differences between primary colorectal adenocarcinomas and distant metastases reveal metastatic CRC subtypes. However, immune-related biomarkers associated with LM from the CRC microenvironment have not been reported. In this study, GSE131418 from the Gene Expression Omnibus (GEO) database was extracted. Then, the sample data of the Moffitt Cancer Center (MCC) cohort were extracted as the discovery group, and the sample data of the CON cohort were extracted as the verify group.

After differentially expressed analysis, differentially expressed genes (DEGs) were obtained for LM and primary cancer (PC) groups of the MCC cohort dataset. Then tumor immune microenvironment and abundance of immune cell infiltration analyses were performed, and the immune-related DEGs were retrieved.

In addition, the transcription factor (TF)-miRNA-mRNA network was constructed, and the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of the immune-related and upregulated and downregulated DEGs were carried out. Then, the protein–protein interaction (PPI) network was conducted and the drug–gene interaction was predicted. Finally, the results of this study were verified through the Consortium (CON) cohort dataset. A flow chart depicting the schematic of this study is shown in Fig. 1. We expected to explore changes in the tumor microenvironment caused by LM from CRC and investigate the immune-related biomarkers associated with LM from the CRC microenvironment.

FIG. 1.

The flow chart of this study.

Materials and Methods

Data sources

The microarray dataset GSE131418 (Species: Homo sapiens) was downloaded from the NCBI GEO database (Barrett and others 2007) (https://www.ncbi.nlm.nih.gov/gds/?term=). A total of 1,135 CRC tumor samples contained 878 primary tumor tissues and 257 colorectal liver metastasis and LM tumor samples, among which 60 LM tumor samples were included. All samples were detected through the Rosetta/Merck Human RSTA Custom Affymetrix 2.0 microarray (HuRSTA_2a520709.CDF) platform.

Data preprocessing

The GEO query package of R language (Davis and Meltzer 2007) was applied to download the GSE131418_RAW.tar data of the GSE131418 dataset from the GEO database (www.ncbi.nlm.nih.gov/geo/). The make.cdf.env function configuration chip of the makecdfenv package of R language was used to annotate the environment package GPL15048_HuRSTA_2a520709.CDF.gz. The expression data of the chip were processed through the background correcting of the affy package of R language (Gautier and others 2004).

Then, the sample data of the MCC cohort were extracted as the discovery group according to Series Matrix File(s), and the sample data of the CON cohort were extracted as the verify group. Then, a box plot was constructed to observe the distribution of the expression of samples in each group. The probe ID was converted to gene symbol using the annotation file GSE131418_cleaned_probes_merck_annotation_Entrez_IDs.csv. The probes that did not correspond to gene symbol were removed. The average value of the diverse probes was considered the eventual expression level of the gene if they corresponded to the identical gene. Finally, the FactoMineR package of R language (Jan and others 2008) was used to draw the principal component analysis diagram to observe the grouping of samples.

Differentially expressed analysis

The DEGs of the LM and PC groups were analyzed using the classical Bayesian modified T-test provided by the limma package of R language (Smyth 2005). The significant differentially expressed cutoff was set as |logFC| > 1 and P < 0.05. Finally, a heat map was drawn utilizing ggscatter functions of the ggpubr package of R language to observe the clustering of samples, and a volcano plot was drawn to show the differences.

Tumor immune microenvironment analysis

Immune score calculation

The Stromal_score, Immune_Score and ESTIMATE_score of all samples were calculated with the ESTIMATE arithmetic (Yoshihara and others 2013) to evaluate the score of stromal and immune cells and analyze the differences of scores between tumor and normal tissues. In addition, a box plot was constructed through the ggpubr package of R language.

Immunocytolysis activity score

The immunocytolysis activity score of all samples was calculated (Rooney and others 2015), and a box plot was constructed through the ggpubr package of R language.

The abundance of immune cell infiltration analysis

The abundance matrix of immune cells of the samples was evaluated through the CIBERSORT deconvolution algorithm (Newman and others 2015) to analyze the abundance of infiltration of immune cells in the samples, and the parameter was set as perm = 100 and QN = TRUE. In addition, the proportion differences of immune cell subsets of the LM and PC groups were calculated. Finally, landscape, heat map, and violin plots were constructed through R language.

Identification of immune-related DEGs

The DEGs and infiltrating abundance of the differentially immune cell subsets between the LM and PC groups were analyzed using a spearman correlation test utilizing the corrplot package of R language. The DEGs with a threshold value of |r| > 0.3 and P < 0.05 were considered immune related. Then, a box plot of the expression quantity of the immune-related DEGs was constructed through the ggpubr package of R language.

Analysis of TF TF-miRNA-mRNA network

Prediction of miRNA-mRNA interaction pairs

The miRNAs of the immune-related DEGs were predicted through the miRWalk3.0 (Dweep and others 2011), TargetScan (Agarwal and others 2015), MiRDB (Wong and Wang 2014), and MirTarBase (Chou and others 2018) databases with a threshold of score >0.95. Based on the results from each database, the miRNAs verified through the MirTarBase database and predicted in at least 1 other database were selected as the final miRNA-mRNA interaction pairs.

Prediction of TF-mRNA interaction pairs

The TF-mRNA interaction pairs of the immune-related DEGs were predicted through the TRRUST tool of the online database (Han and others 2017) (https://www.grnpedia.org/trrust/). To analyze the TF-miRNA-mRNA network, the miRNAs of the immune-related DEGs were predicted through the miRWalk3.0, TargetScan, MiRDB, and MirTarBase databases with a threshold of score >0.95. Based on the results from each database, the miRNAs that were validated through the MirTarBase database and were predicted in at least 1 other database were selected as the final miRNA-mRNA interaction pairs.

The TF-mRNA interaction pairs of the immune-related DEGs were predicted through the TRRUST tool of the online database. Based on the miRNA-mRNA interaction pairs and TF-mRNA interaction pairs, the TF-miRNA-mRNA network was constructed utilizing the Cytoscape software. According to the upregulated and downregulated DEGs, the clusterProfiler package of R language was used to perform GO and KEGG pathway enrichment analyses with the cutoff of P < 0.05 and count ≥2. According to the interaction pairs of 8 miRNA-mRNAs and 50 TF-mRNAs, the TF-miRNA-mRNA network was constructed.

Construction of the TF-miRNA-mRNA network

Based on the miRNA-mRNA and TF-mRNA interaction pairs, the TF-miRNA-mRNA network was constructed utilizing Cytoscape software (Shannon and others 2003).

Enrichment analyses

According to the immune-related DEGs and upregulated and downregulated DEGs of the TF-miRNA-mRNA network, the clusterProfiler package of R language (Yu and others 2012) (Version 3.12.0) was used to perform GO (Ashburner and others 2000) and KEGG pathway (Ogata and others 1999) enrichment analyses with a cutoff of P < 0.05 and count ≥2. The GO functional analysis included biological processes (BPs), cellular components (CCs), and molecular functions (MFs).

PPI network analysis

The Search Tool for the Retrieval of Interacting Genes (STRING) (Szklarczyk and others 2017) (Version 10.0; www.string-db.org/) database was used to analyze the interactions between protein and protein encoded by genes. The PPI score was set as 0.4 (referred as median confidence). Afterward, the PPI network was constructed using Cytoscape software. In addition, significant network modules were analyzed by the Molecular Complex Detection (MCODE) (Bader and Hogue 2003) plugin in Cytoscape software with the cutoff of score ≥9.

Drug–gene interaction prediction

According to the immune-related DEGs, drug–gene interactions were predicted through the Drug–Gene Interaction database (DGIdb) (Cotto and others 2017) (Version: 3.0.2). The drug–gene interaction network was constructed using Cytoscape, and drug information was retrieved from the DrugBank database (Law and others 2014).

Verification analysis

The differentially immune cells screened in this study were verified through the CON cohort dataset.

Statistical analysis

SPSS 22.0 statistical software was used for statistical processing. Continuous variables were expressed as mean ± standard deviation. The enumeration data consistent with normal distribution between the LM and PC groups were compared using the T-test. Statistical differences of P < 0.05 were considered significant.

Results

Data preprocessing and basic information statistics

The sample basic statistics information of the MCC and CON cohort datasets in the GSE131418 dataset is shown in Table 1. The samples were involved in 60,607 probes. In addition, based on the MCC cohort dataset, the number of samples in the PC group was 333, and 43 samples in the LM group. According to the CON cohort dataset, there were 545 samples in the PC group and 17 in the LM group. The HuRSTA_2a520709.CDF probe annotation information included 47,408 gene-probe interactions involving 24,495 genes, and 22,913 gene-probe interactions belonged to multiple probes corresponding to the identical gene.

Table 1.

Sample Basic Information Statistic

Cancer type: CRC	Dataset: GSE131418
Total number of probes	60,607
Number of probes with annotated information	47,408
Gene number	24,495
PC/LM (MCC+CON) sample number (total)	333 + 545/43 + 17 (878/60)

CON, Consortium; CRC, colorectal cancer; LM, lung metastasis; MCC, Moffitt Cancer Center; PC, primary cancer.

The box plot of samples of the MCC and CON cohort datasets is shown in Fig. 2A and B, and the principal component analysis diagram of the samples of the MCC and CON cohort datasets is shown in Fig. 2C and D. Finally, the principal component analysis diagram of the total samples is shown in Fig. 2E. The results showed that the overall expression pattern of each sample was similar, and there were some differences among the samples of each group. Data preprocessing was consistent with the prediction.

FIG. 2.

The DEGs between the LM and PC groups of the MCC cohort dataset. GSE131418 from the GEO database was extracted. The sample data of the MCC cohort were extracted as the discovery group and the sample data of the CON cohort were extracted as the verify group. To observe the distribution of the expression and grouping of samples, a box plot and principal component analysis diagram were constructed through the FactoMineR package of R language. The classical Bayesian modified T-test and provided limma package of R language were applied to screen the DEGs between the LM and PC groups, and the screening criteria were |logFC| > 1 and P < 0.05. A total of 268 DEGs were obtained, containing 173 upregulated and 95 downregulated DEGs. (A) The box plot of samples of the MCC cohort dataset; (B) the box plot of samples of the CON cohort dataset; (C) the principal component analysis diagram of the samples of the MCC cohort dataset; (D) the principal component analysis diagram of the samples of the CON cohort dataset; (E) the principal component analysis diagram of the total samples. (F) Heat map of the DEGs; (G) volcano plot of the DEGs. CON, Consortium; DEGs, differentially expressed genes; GEO, Gene Expression Omnibus; LM, lung metastasis; MCC, Moffitt Cancer Center; PC, primary cancer.

The DEGs between the PC and LM groups of the MCC cohort dataset

Based on the screening criteria set as P < 0.05 and |logFC| > 1, a total of 268 DEGs between the PC and LM groups of the MCC cohort dataset were obtained, containing 173 upregulated and 95 downregulated DEGs (Supplementary File S1). The heat map and volcano plot of the DEGs are shown in Fig. 2F and G. The 2 sample groups could be significantly separated, indicating that the difference analysis results were reliable.

Immune score calculation and immunocytolysis activity score of the MCC cohort dataset

The dataset of the ESTIMATE arithmetic included 10,007 matched genes and 405 nonmatched genes. The box plots of the Stromal_score, Immune_Score, and ESTIMATE_score among the LM and PC groups are shown in Fig. 3A–C. The results illustrated that the Immune_Score of the samples in the LM group was significantly higher compared with the PC group (P = 0.018), and there was no difference between the Stromal_score and ESTIMATE_score of samples in the LM and PC groups (P > 0.05). The box plot of the cytolysis activity score is shown in Fig. 3D, and the results uncovered that there was no difference between the cytolysis activity score of samples in the LM and PC groups (P = 0.2).

FIG. 3.

Immune score calculation and immune cytolysis activity score of the MCC cohort dataset. The Stromal_score, Immune_Score, and ESTIMATE_score of all samples were calculated with the ESTIMATE arithmetic, and a box plot was drawn through the ggpubr package of R language. The DEGs and infiltrating abundance of the differentially immune cell subsets of the LM and PC groups were subjected to a Spearman's correlation test through the corrplot package of R language. DEGs with a threshold value of |r| > 0.3 and P < 0.05 were considered immune related. According to the immune-related DEGs, the clusterProfiler package of R language was used to perform GO and KEGG pathway enrichment analyses with the cutoff of P < 0.05 and count ≥2. (A) The box plot of the Stromal_score among the LM and PC groups; (B) the box plot of the Immune_Score among the LM and PC groups; (C) the box plot of the ESTIMATE_score among LM and PC groups; (D) the box plot of the cytolysis activity score. (E) Landscape of the immune cell infiltration of the 57 PCs and 43 LM cancers samples; (F) heat map of the 22 types of immune cell infiltration; (G) violin plot of the 22 types of immune cell infiltration; (H) the top 8 terms of GO analysis of immune-related DEGs; (I) the top 8 terms of KEGG analysis of immune-related DEGs. BP, biological processes; CCs, cellular components; GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function.

The abundance of immune cell infiltration of the MCC cohort dataset

In total, 515 genes were included in the 547 related genes in the LM22 signal matrix, accounting for 94.15%. The CIBERSORT algorithm estimated the abundance matrix of 22 types of immune cell infiltration in all samples. The results showed that 325 of the 376 samples were effective (P < 0.05), including 282 PCs and 43 LM cancers. The landscape of the immune cell infiltration of 57 PC and 43 LM cancer samples is shown in Fig. 3E, and the results showed that the infiltration ratio of M0 and M2 macrophages of samples was higher than others. In addition, the heat map of the 22 types of immune cell infiltration is shown in Fig. 3F, and the results of the infiltration ratio of M0 and M2 macrophages of samples were higher than others.

Moreover, the violin plot of the 22 types of immune cell infiltration is shown in Fig. 3G, and the results showed that the proportion of infiltration was different among 12 types of immune cells (naive B cells, memory B cells, CD8 T cells, resting memory CD4 T cells, helper follicular T cells, monocytes, M0 macrophages, M2 macrophages, activated dendritic cells, activated mast cells, eosinophils, and neutrophils), and the proportion of infiltration of memory B cells, CD8 T cells, and M0 and M2 macrophages in the LM group was higher compared with the PC group.

The immune-related DEGs of the MCC cohort dataset

According to the 7 immune infiltration cells (naive B cells, resting memory CD4 T-cells, M0 macrophages, activated dendritic cells, activated mast cells, eosinophils, and neutrophils) and DEGs, a correlation analysis was conducted. Then, a total of 138 immune-related DEGs were screened (0 in naive B cells, 26 in resting memory CD4 T cells, 54 in M0 macrophages (Supplementary File S2), 16 in activated dendritic cells, 12 in activated mast cells, 17 in eosinophils, and 13 in neutrophils). As the infiltration rate of M0 macrophage cells was higher compared with other cells, and more immune-related DEGs were obtained, the follow-up analysis focused on the M0 macrophage cells.

The enrichment analysis results showed that the 54 immune-related DEGs in M0 macrophage cells were significantly enriched in 229 GO-BP, 24 GO-CC, and 33 GO-MF functions, as well as 4 KEGG pathways, and the top 8 terms of GO and KEGG analyses of immune-related DEGs are shown in Fig. 3H and I. The GO enrichment analysis of 54 immune-related DEGs revealed that the significantly enriched GO-BP terms included protein-containing complex remodeling, protein-lipid complex remodeling, and plasma lipoprotein particle remodeling. The mainly enriched GO-CC terms contained zymogen granule membrane, chylomicron, and zymogen granules.

Furthermore, it was also found that these genes played essential roles in receptor ligand activity, receptor regulator activity, and carbohydrate binding according to the GO-MF analysis. Meanwhile, KEGG pathway analysis suggested that these genes are involved in numerous pathways, including cholesterol metabolism, arachidonic acid metabolism, interleukin (IL)-7 signaling pathway, and pancreatic secretion.

Construction of the TF-miRNA-mRNA network

According to the interaction pairs of 8 miRNA-mRNAs and 50 TF-mRNAs, the TF-miRNA-mRNA network was constructed (Fig. 4A). The results showed that secreted phosphoprotein 1 (SPP1) was regulated by 12 TFs, including TFCP2, NR3C1, POU2F1, POU2F2, POU5F1, ERG, SP1, CEBPA, FOXD3, HDAC1, HTAIP2, and ING4. Meanwhile, oxidized low-density lipoprotein receptor 1 (OLR1) was regulated by 3 miRNAs (has-miR-6734-5p, has-miR-6867-5p, and has-miR-4733-5p) and 3 TFs (HR, NFKB1, and RELA).

FIG. 4.

Construction of the TF-miRNA-mRNA and PPI network network of immune-related DEGs of M0 macrophage cells. To analyze the TF-miRNA-mRNA network, the miRNAs of the immune-related DEGs were predicted. To analyze the interactions between proteins and proteins encoded by DEGs, the STRING database was utilized to analyze their interactions. (A) The TF-miRNA-mRNA network; (B) the top 8 terms of the GO analysis of upregulated DEGs; (C) the top 8 terms of the KEGG pathway analysis of upregulated DEGs; (D) the top 8 terms of the GO analysis of downregulated DEGs; (E) the top 10 terms of the KEGG pathway analysis of downregulated DEGs; (F) the top 8 terms of the GO analysis of all DEGs; (G) the top 10 terms of the KEGG pathway analysis of all DEGs; (H) the PPI network of the upregulated and downregulated DEGs; (I) the PPI network of the 54 immune-related DEGs of M0 macrophage cells.; MCODE: Molecular Complex Detection; PPI, protein–protein interaction; STRING, Search Tool for the Retrieval of Interacting Genes; TF, transcription factor.

The results also showed that the TF SAM pointed domain containing ETS transcription factor (SPDEF) was also a downregulated DEG. The upregulated DEGs were significantly enriched in 539 GO-BP, 29 GO-CC, and 60 GO-MF functions as well as 8 KEGG pathways, and the top 8 terms of GO analysis and KEGG pathway analysis are shown in Fig. 4B and C. Meanwhile, the downregulated DEGs were significantly enriched in 391 GO-BP, 14 GO-CC, and 37 GO-MF functions, as well as 29 KEGG pathways, and the top 8 terms of GO analysis and the top 10 terms of KEGG pathway analysis are shown in Fig. 4D and E.

In addition, all DEGs were enriched in 26 KEGG pathways and 897 GO functions. In this study, we present all DEGs of the Top 8 terms of the GO-BP functions and top 10 terms of the KEGG pathways (Fig. 4F, G). The GO enrichment analysis of the DEGs revealed that the significantly enriched GO-BP terms included humoral immune response, antimicrobial humoral response, and kidney morphogenesis. The mainly enriched GO-CC terms contained collagen-containing extracellular matrix, blood microparticles, and lamellar bodies.

Furthermore, it was found that these genes play essential roles in receptor regulator activity, glycosaminoglycan binding, and sulfur compound binding according to the GO-MF analysis. Meanwhile, KEGG pathway analysis suggested that DEGs are involved in numerous pathways, including the IL-7 signaling pathway, cell adhesion molecules (CAMs), and arachidonic acid metabolism.

PPI network of immune-related DEGs of M0 macrophage cells

To analyze the interactions between proteins and proteins encoded by DEGs, the PPI network was constructed. The network showed that a total of 107 nodes (Fig. 4H), including 65 upregulated genes and 42 downregulated genes, corresponded to protein. In addition, the PPI network of the 54 immune-related DEGs of M0 macrophage cells is shown in Fig. 4I. Overall, 23 protein interactions in this PPI network were found to be involved in 27 mRNAs, including 10 upregulated and 17 downregulated genes.

Drug–gene interaction prediction of immune-related DEGs of M0 macrophage cells

In total, 38 drug–gene interactions of immune-related DEGs of M0 macrophage cells were predicted to contain 34 small molecule drugs and 8 mRNAs [eg, SPP1, somatostatin (SST), C-C motif chemokine ligand 11 (CCL11), vasoactive intestinal peptide (VIP), alcohol dehydrogenase 1C (class I), gamma polypeptide (ADH1C), phospholipase A2 group IIA (PLA2G2A), transient receptor potential cation channel subfamily A member 1 (TRPA1), and apolipoprotein E (APOE)]. In addition, drug information was retrieved from the DrugBank database (Table 2).

Table 2.

The Drug Information of the Drug–Gene Interaction Prediction of Immune-Related Differentially Expressed Genes of M0 Macrophages Cells

Gene	Drug name
SPP1	Tretinoin
SPP1	Calcitonin
SPP1	Tacrolimus
SPP1	Gentamicin
SST	Cysteamine
SST	Estradiol benzoate
SST	Amphetamine
SST	Captopril
SST	Valproic acid
SST	Streptozotocin
SST	Aspirin
SST	Ganciclovir
CCL11	Mometasone furoate
CCL11	Roxithromycin
VIP	Omeprazole
VIP	Lisinopril
VIP	Ribavirin
VIP	Quazepam
VIP	Flutamide
VIP	Digoxin
ADH1C	Fomepizole
ADH1C	Tretinoin
ADH1C	Cholic acid
PLA2G2A	Etidronate disodium
PLA2G2A	Ascorbate
PLA2G2A	Gentamicin
PLA2G2A	Piroxicam
TRPA1	Menthol
APOE	Prednisone
APOE	Tretinoin
APOE	Ganciclovir
APOE	Triamcinolone
APOE	Albumin human
APOE	Irbesartan
APOE	Vitamin E
APOE	Lorazepam
APOE	Soybean oil
APOE	Gonadotropin, Chorionic

ADH1C, alcohol dehydrogenase 1C (class I), gamma polypeptide; APOE, apolipoprotein E; CCL11, C-C motif chemokine ligand 11; PLA2G2A, phospholipase A2 group IIA; SPP, secreted phosphoprotein 1; SST, somatostatin; SPP1, secreted phosphoprotein 1; TRPA1, transient receptor potential cation channel subfamily A member 1; VIP, vasoactive intestinal peptide.

Verification analysis through the CON dataset

The differentially expressed immune cells screened in this study were verified through the CON cohort dataset. The results of the CIBERSORT calculations showed that 434 of the 562 samples were valid (P < 0.05), including 419 PCs and 15 lung metastatic cancers, among which the landscape of 35 patients with PC and 15 with lung metastatic cancer is shown in Fig. 5A. In addition, the infiltration ratio of the M0 and M2 macrophages of the samples was higher than others, which was consistent with the MCC dataset (Fig. 5B).

FIG. 5.

Verification analysis through the CON dataset. The differentially immune cells screened were verified through the CON dataset. The results of the CIBERSORT calculations showed that 434 of the 562 samples were valid (P < 0.05), including 419 PCs and 15 lung metastatic cancers. (A) The landscape of the 35 patients with PC and 15 with lung metastatic cancer; (B) heat map of immune cell infiltration; (C) violin plot of immune cell infiltration; (D) venn diagram of differential immune cells.

Moreover, the violin plot of immune cell infiltration is shown in Fig. 5C, and the results showed that the proportion of infiltration between the LM and PC group samples of the CON cohort dataset was different among the 9 types of immune cells (naive B cells, plasma cells, resting memory CD4 T cells, delta gamma T cells, M0 macrophages, activated dendritic cells, activated mast cells, eosinophils, and neutrophils). The different immune cells screened together with the MCC cohort dataset included naive B cells, resting memory CD4 T cells, M0 macrophages, activated dendritic cells, activated mast cells, eosinophils, and neutrophils (Fig. 5D).

Discussion

CRC is the most harmful malignant tumor to human health worldwide. In this study, we analyzed changes in the tumor microenvironment caused by LM from CRC and investigated the immune-related biomarkers associated with LM from the CRC microenvironment. A total of 268 DEGs were obtained from the LM and PC groups of the MCC cohort dataset. The Immune_Score of the samples in the LM group was significantly higher compared with the PC group, and there was no difference between the Stromal_score and ESTIMATE_score of samples in the LM and PC groups. In addition, there was no difference in the cytolysis activity score of samples in the LM and PC groups.

The results of immune cell infiltration showed that the infiltration ratio of M0 and M2 macrophages of samples was higher than others. Then, a total of 138 immune-related DEGs were screened, and 54 immune-related DEGs in M0 macrophages were identified. Moreover, the TF-miRNA-mRNA network was constructed among 8 miRNA-mRNAs and 50 TF-mRNAs, SPP1 was found to be regulated by 12 TFs and OLR1 by 3 miRNAs and 3 TFs.

The results also showed that TF SPDEF is a downregulated DEG. KEGG pathway analysis showed that the DEGs in the TF-miRNA-mRNA network were primarily involved in the IL-7 signaling pathway and CAMs. In total, 23 protein interactions in this PPI network of M0 macrophage cells were found to be involved in 27 mRNAs. Overall, 38 drug–gene interactions of the immune-related DEGs of M0 macrophage cells were predicted to contain 34 small molecule drugs and 8 mRNAs (eg, SPP1, SST, CCL11, VIP, ADH1C, PLA2G2A, TRPA1, and APOE). Finally, the CON cohort dataset verified that the infiltration ratio of M0 and M2 macrophages of samples was also higher than others, which was consistent with the results of the MCC cohort dataset.

Macrophages constitute a major component of tumor stroma and are one of the most important players in the tumor microenvironment. A previous study found that the polarization status of tumor-associated macrophages into a proinflammatory type M1 or anti-inflammatory type M2 may influence cancer progression and patient survival (Popēna and others 2018). Keirsse and others (2018) found that the role of hepatic macrophages, both recruited monocyte-derived and tissue-resident Kupffer cells, is more versatile than initially thought.

In this study, the CON cohort dataset verified that the infiltration ratio of M0 and M2 macrophages of samples was higher than other cells, which was consistent with the results of the MCC cohort dataset. In addition, many studies have shown that macrophages play a significant role in the occurrence of CRC. For instance, van der Bij and others (2008) suggested that the inhibition of monocyte migration promotes tumor growth, supporting that not only resident but also newly recruited macrophages limit peritoneal CRC metastasis development, which further implies that macrophages are likely related to the progression of LM from CRC.

SPP1, a gene in the adipokine family, has been reported to have close connection with CRC. For instance, Choe and others (2018) found that SPP1 might be associated with poor survival in CRC. Xu and others (2017) showed that the expression of SPP1 was significantly upregulated in CRC and promoted its metastasis by activating the EMT pathway. OLR1, also known as LOX-1, plays a key role in the occurrence and development of CRC. For instance, Murdocca and others (2016) showed that LOX-1 is a new potential molecular target in CRC.

Meanwhile, Wang and others (2014) found that OLR1 could play crucial roles in the regulation of CRC metastasis. SPDEF has been found to be closely linked with CRC. Lo and others (2017) highlighted that β-catenin activity determines the proliferation or quiescence of CRC cells based on the absence or presence of SPDEF. Interestingly, in this study, the TF-miRNA-mRNA network showed that SPP1 is regulated by 12 TFs, OLR1 is regulated by 3 miRNAs and 3 TFs, and TF SPDEF is also a downregulated DEG. The KEGG pathway analysis showed that the DEGs in the TF-miRNA-mRNA network were primarily involved in the IL-7 signaling pathway and CAMs.

Amatya and others uncovered that IL-17 is the founding member of a novel family of inflammatory cytokines. While the proinflammatory properties of IL-17 are key to its host-protective capacity, unrestrained IL-17 signaling is associated with immunopathology, autoimmune disease, and cancer progression (Amatya and others 2017). Moreover, a previous study showed the abnormal immunoexpression of CAMs in cervical cancer (de Méndez and Bosch 2011). Taken together, we speculate that SPP1, OLR1, and SPDEF might contribute to LM in CRC by activating the IL-7 signaling and CAMs pathways.

SST has been used as a biomarker for the diagnosis of CRC. Liu and others (2017) showed that SST may serve as aberrantly methylation-based biomarkers for the precise diagnosis and treatment of CRC in the future. Liu and others (2016) also found that methylation of the serum SST gene is an independent prognostic marker in CRC. Cho and others (2016) illustrated that the discrepancy in CCL11 and CCL24 expression between glandular and stromal cells may shed light on how CRC evades the immune system, which would enable further development of immunotherapies that target these chemokines.

Zyqulska and others (2019) compared the plasma levels of VIP in patients with gastrointestinal malignancies and healthy controls, and the results showed that the plasma levels of VIP observed in the gastric cancer group were lower than in patients with CRC. Offermans and others (2018) also found that ADH1B rs3811802 and ADH1C rs4147542 significantly modified the alcohol-colon cancer association in women. Meanwhile, He and others (2015) uncovered that PLA2G2A overexpression is associated with poor therapeutic response and inferior outcome in patients with rectal cancer receiving neoadjuvant concurrent chemoradiotherapy. de Almeida and others (2020) highlighted that TRPA1 could be investigated as a target for breast carcinoma pain treatment.

Nevertheless, the role of TRPA1 in CRC development has not yet been explored. It has also been found that APOE is associated with tumor progression and poor survival in CRC (Zhao and others 2018). In this study, there were 38 drug–gene interactions of immune-related DEGs of M0 macrophage cells predicted to contain 34 small molecule drugs and 8 mRNAs (SPP1, SST, CCL11, VIP, ADH1C, PLA2G2A, TRPA1, and APOE).Taken together, we speculate that SPP1, SST, CCL11, VIP, ADH1C, PLA2G2A, TRPA1, and APOE contribute to CRC progression.

However, there are some limitations in this study. First, the small LM sample size of the CON dataset limited the differentially expressed analysis. In addition, relevant experiments, including cell biology assays and animal and clinical studies, need to be performed to verify the multiple candidate targets and signaling pathways identified from our bioinformatics analyses.

Conclusion

In summary, M0 macrophage cells SPP1, OLR1, SPDEF, SST, CCL11, VIP, ADH1C, PLA2G2A, TRPA1, and APOE are immune-related biomarkers associated with LM from the CRC microenvironment, which might provide assistance for patients with personalized treatment.

Footnotes

Authors' Contributions

All authors participated in the conception and design of the study. Conceived and drafted the article: P.Y. and S.Y. Wrote the article: W.D. and Z.J. Analyzed the data: P.Y. and W.Y. Drew figures: X.J. and G.X. All authors read and approved the article.

Acknowledgments

The authors gratefully acknowledge the database available to us for this study.

Availability of Data and Materials

The datasets generated during this study are not publicly available, but were obtained from corresponding authors on reasonable request.

Author Disclosure Statement

No competing financial interests exist.

Funding Information

This research was supported by Zhejiang Provincial Natural Science Foundation of China (No. LQ17H160008) and National Key R&D Program of China (No. 2016YFC1302803).

Supplementary Material

Supplementary File S1

Supplementary File S2

References

Agarwal

, Bell

, Nam

, Bartel

. 2015. Predicting effective microRNA target sites in mammalian mRNAs. eLife 4:e05005. DOI: 10.7554/eLife.05005.

Amatya

, Garg

, Gaffen

. 2017. IL-17 Signaling: the Yin and the Yang. Trends Immunol, 38(5):310–322.

Ashburner

, Ball

, Blake

, Botstein

, Butler

, Cherry

, Davis

, Dolinski

, Dwight

, Eppig

, Harris

, Hill

, Issel

, Kasarskis

, Lewis

, Matese

, Richardson

, Ringwald

, Rubin

, Sherlock

. 2000. Gene Ontology: tool for the unification of biology. Nucleic Acids Res, 27(1):29–34.

Bader

, Hogue

. 2003. An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinformatics, 4(1):2.

Barrett

, Troup

, Wilhite

, Ledoux

, Rudnev

, Evangelista

, Kim

, Soboleva

, Tomashevsky

, Edgar

. 2007. NCBI GEO: mining millions of expression profiles—database and tools. Nucleic Acids Res, 33(Database issue):D562–D566.

Cho

, Lim

, Won

, Bae

, Kim

, Min

, Noh

. 2016. Eosinophils in colorectal neoplasms associated with expression of CCL11 and CCL24. J Pathol Transl Med, 50(1):45–51.

Choe

, Yi

, Chai

, Park

. 2018. Upregulation of the adipokine genes ADIPOR1 and SPP1 is related to poor survival outcomes in colorectal cancer. J Surg Oncol, 117(8):1833–1840.

Chou

, Shrestha

, Yang

, Chang

, Lin

, Liao

, Huang

, Sun

, Tu

, Lee

, Chiew

, Tai

, Wei

, Tsai

, Huang

, Wang

, Wu

, Ho

, Chen

, Chuang

, Hsieh

, Wu

, Chen

, Li

, Wu

, Huang

, Ng

, Buddhakosai

, Huang

, Lan

, Huang

, Weng

, Cheng

, Liang

, Hsu

, Huang

. 2018. miRTarBase update 2018: a resource for experimentally validated microRNA-target interactions. Nucleic Acids Res, 46(D1):D296–D302.

Cotto

, Wagner

, Feng

, Kiwala

, Coffman

, Spies

, Wollam

, Spies

, Griffith

. 2017. DGIdb 3.0: a redesign and expansion of the drug–gene in teraction database. Nucleic Acids Res, 46(D1):D1068–D1073.

10.

Davis

, Meltzer

. 2007. GEOquery: a bridge between the Gene Expression Omnibus (GEO) and BioConductor. Bioinformatics, 23(14):1846–1847.

11.

de Almeida

, Rigo

, De Prá

, Milioli

, Pereira

, Lückemeyer

, Antoniazzi

, Kudsi

, Araújo

DMPA

, Oliveira

, Ferreira

, Trevisan

. 2020. Role of transient receptor potential ankyrin 1 (TRPA1) on nociception caused by a murine model of breast carcinoma. Pharmacol Res 152:104576. DOI: 10.1016/j.phrs.2019.104576.

12.

de Méndez

, Bosch

. 2011. Abnormal immunoexpression of cell adhesion molecules (CAMs) in cervical cancer. Int J Surg Pathol, 19(6):733–742.

13.

Dweep

, Sticht

, Pandey

, Gretz

. 2011. miRWalk-Database: prediction of possible miRNA binding sites by “walking”. the genes of three genomes. J Biomed Inform, 44(5):839–847.

14.

Gautier

, Cope

, Bolstad

, Irizarry

. 2004. affy—analysis of Affymetrix GeneChip data at the probe level. Bioinformatics, 20(3):307–315.

15.

Hadad

, Hazmi

, Alhebshi

, Aldahlawi

, Bassam

. 2019. Lactobacillus rhamnosus enhances the immunological antitumor effect of 5-fluorouracil against colon cancer. Pak J Biol Sci, 22(12):597–606.

16.

Han

, Cho

, Lee

, Yun

, Kim

, Bae

, Yang

, Kim

, Lee

, Kim

, Lee

, Kang

, Jeong

, Kim

, Jeon

, Jung

, Nam

, Chung

, Kim

, Lee

. 2017. TRRUST v2: an expanded reference database of human and mouse transcriptional regulatory interactions. Nucleic Acids Res, 46(D1):D380–D386.

17.

, Lee

, Shiue

, Lee

, Lin

, Chen

, Wu

, Li

. 2015. PLA2G2A overexpression is associated with poor therapeutic response and inferior outcome in rectal cancer patients receiving neoadjuvant concurrent chemoradiotherapy. Histopathology, 66(7):991–1002.

18.

Huang

, Mao

, Lu

, Yu

, Liao

. 2020a. A novel immune-related genes prognosis biomarker for melanoma: associated with tumor microenvironment. Aging (Albany NY), 12(8):6966–6980.

19.

Huang

, Cai

, Liu

, Yuan

. 2020b. Low mutation burden and differential tumor-infiltrating immune cells correlate with lymph node metastasis in colorectal cancer. Int J Clin Exp Pathol, 13(9):2259–2269.

20.

Hwang

, Park

, Kim

, Chang

, Kim

, Choi

, Kim

, Zo

, Oh

. 2010. Early intrapulmonary recurrence after pulmonary metastasectomy related to colorectal cancer. Ann Thorac Surg, 90(2):398–404.

21.

Iplik

, Ertugrul

, Kozanoglu

, Baran

, Cakmakoglu

. 2018. An answer to colon cancer treatment by mesenchymal stem cell originated from adipose tissue. Iran J Basic Med Sci, 21(5):465–468.

22.

Jan

, Josse

, Husson

. 2008. FactoMineR: an R Package for multivariate analysis. J Stat Softw 25(1). DOI:10.18637/jss.v025.i01.

23.

Kamal

, Schmit

, Hoehn

, Amos

, Frost

. 2019. Transcriptomic differences between primary colorectal adenocarcinomas and distant metastases reveal metastatic colorectal cancer subtypes. Cancer Res, 79(16):4227–4241.

24.

Keirsse

, Van Damme

, Geeraerts

, Beschin

, Raes

, Van Ginderachter

. 2018. The role of hepatic macrophages in liver metastasis. Cell Immunol, 330:202–215.

25.

Law

, Knox

, Djoumbou

, Jewison

, Guo

, Liu

, Maciejewski

, Arndt

, Wilson

, Neveu

, Tang

, Gabriel

, Ly

, Adamjee

, Dame

, Han

, Zhou

, Wishart

. 2014. DrugBank 4.0: shedding new light on drug metabolism. Nucleic Acids Res, 42(Database Issue):D1091–D1097.

26.

Lei

, Feng

, Hong

. 2020. ELFN1-AS1 accelerates the proliferation and migration of colorectal cancer via regulation of miR-4644/TRIM44 axis. Cancer Biomark, 27(4):433–443.

27.

Liu

, Li

, Sun

, Wang

, Xing

, Yuan

. 2017. Aberrantly methylated-differentially expressed genes and pathways in colorectal cancer. Cancer Cell Int, 17(1):75.

28.

Liu

, Chew

, Tham

, Tang

, Ong

, Zhao

. 2016. Methylation of serum SST gene is an independent prognostic marker in colorectal cancer. Am J Cancer Res, 6(9):2098–2108.

29.

, Noah

, Chen

, Zou

, Borras

, Vilar

, Shroyer

. 2017. SPDEF induces quiescence of colorectal cancer cells by changing the transcriptional targets of β-catenin. Gastroenterology, 153(1):205–218.

30.

Murdocca

, Mango

, Pucci

, Biocca

, Testa

, Capuano

, Paolesse

, Sanchez

, Orlandi

, di Natale

, Novelli

, Sangiuolo

. 2016. The lectin-like oxidized LDL receptor-1: a new potential molecular target in colorectal cancer. Oncotarget, 7(12):14765–14780.

31.

Newman

, Liu

, Green

, Gentles

, Feng

, Xu

, Hoang

, Diehn

, Alizadeh

. 2015. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods, 12(5):453–457.

32.

Offermans

NSM

, Ketcham

, van den Brandt

, Weijenberg

, Simons

CCJM

. 2018. Alcohol intake, ADH1B and ADH1C genotypes, and the risk of colorectal cancer by sex and subsite in the Netherlands Cohort Study. Carcinogenesis, 39(3):375–388.

33.

Ogata

, Goto

, Sato

, Fujibuchi

, Bono

, Kanehisa

. 1999. KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res, 28(1):27–30.

34.

Popēna

, Ābols

, Saulīte

, Pleiko

, Zandberga

, Jēkabsons

, Endzeliņš

, Llorente

, Linē

, Riekstiņa

. 2018. Effect of colorectal cancer-derived extracellular vesicles on the immunophenotype and cytokine secretion profile of monocytes and macrophages. Cell Commun Signal, 16(1):17.

35.

Raoux

, Maulat

, Mokrane

, Fares

, Suc

, Muscari

. 2020. Impact of the strategy for curative treatment of synchronous colorectal cancer liver metastases. J Visc Surg, 157(4):289–299.

36.

Rooney

, Shukla

, Wu

, Getz

, Hacohen

. 2015. Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell, 160(1–2):48–61.

37.

Sargent

. 2005. Endpoints for Colon Adjuvant Clinical Trials (CACT): recommendations based on individual patient data (IPD) from 20898 patients (pts) and 18 randomized trials. J Clin Oncol, 23(16):3512.

38.

Shannon

, Markiel

, Ozier

, Baliga

, Wang

, Ramage

, Amin

, Schwikowski

, Ideker

. 2003. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res, 13(11):2498–2504.

39.

Smyth

. 2005. LIMMA: linear models for microarray data. In: Bioinformatics and computational biology solutions using R and bioconductor. Springer, pp. 397–420.

40.

Sun

, Liu

, Wang

, Zhang

, Li

, Guo

, Zhang

, Suo

. 2017. Colorectal cancer cells suppress CD4+ T cells immunity through canonical Wnt signaling. Oncotarget, 8(9):15168–15181.

41.

Sung

, Ferlay

, Siegel

, Laversanne

, Soerjomataram

, Jemal

, Bray

. 2021. Global Cancer Statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin, 71(3):209–249.

42.

Szklarczyk

, Morris

, Cook

, Kuhn

, Wyder

, Simonovic

, Santos

, Doncheva

, Roth

, Bork

, Jensen

, von Mering

. 2017. The STRING database in 2017: quality-controlled protein–protein association networks, made broadly accessible. Nucleic Acids Res, 45(D1):D362–D368.

43.

Tamjidifar

, Akbari

, Tarzi

, Sadeghzadeh

, Abolghasemi

, Poursaei

, Shomali

, Mahdavi

. 2021. Prognostic and diagnostic values of miR-506 and SPON 1 in colorectal cancer with clinicopathological considerations. J Gastrointest Cancer, 52(1):125–129.

44.

Tan

, Hu

, Xiao

, Tang

, Chen

, He

, Chen

, Jiang

, Wang

, Yuan

, Ding

. 2020. Silencing of brain-expressed X-linked 2 (BEX2) promotes colorectal cancer metastasis through the Hedgehog signaling pathway. Int J Biol Sci, 16(2):228–238.

45.

van der Bij

, Bögels

, Oosterling

, Kroon

, Schuckmann

, de Vries

, Meijer

, Beelen

, van Egmond

. 2008. Tumor infiltrating macrophages reduce development of peritoneal colorectal carcinoma metastases. Cancer Lett, 262(1):77–86.

46.

Wang

, Xu

, Tu

, Tan

, Zhu

, Han

, Dou

, Zeng

, Tan

, Du

, Jiao

, Cai

. 2014. Comparison analysis in synchronous and metachronous metastatic colorectal cancer based on microarray expression profile. Hepatogastroenterology, 61(136):2215–2218.

47.

Wong

, Wang

. 2014. miRDB: an online resource for microRNA target prediction and functional annotations. Nucleic Acids Res, 43(D1):D146–D152.

48.

, Sun

, Jiang

, Zhou

, Gu

, Liu

, Xu

. 2017. SPP1, analyzed by bioinformatics methods, promotes the metastasis in colorectal cancer by activating EMT pathway. Biomed Pharmacother, 91:1167–1177.

49.

Yoshihara

, Shahmoradgoli

, Martínez

, Vegesna

, Kim

, Torres-Garcia

, Treviño

, Shen

, Laird

, Levine

, Carter

, Getz

, Stemke

, Mills

, Verhaak

. 2013. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun, 4(1):2612.

50.

, Wang

, Han

, He

. 2012. clusterProfiler: an R Package for comparing biological themes among gene clusters. OMICS, 16(5):284–287.

51.

Zhao

, Zou

, Guan

, Wang

, Jiang

, Liu

, Li

, Lin

, Liu

, Yang

, Gao

, Wang

. 2018. Apolipoprotein E overexpression is associated with tumor progression and poor survival in colorectal cancer. Front Genet 9:650. DOI: 10.3389/fgene.2018.00650.

52.

Zygulska

, Furgala

, Kaszuba-Zwoińska

, Krzemieniecki

, Gil

. 2019. Changes in plasma levels of cholecystokinin, neurotensin, VIP and PYY in gastric and colorectal cancer—preliminary results. Peptides 122:170148. DOI: 10.1016/j.peptides.2019.170148.