Identification of DEMs and DEGs.
To identify DEMs and DEGs between primary CRC and liver metastases, we analyzed two miRNA datasets (GSE72199 and GSE54088), and one mRNA dataset (GSE81582). According to the criteria of p < 0.05 and |logFC| ⩾ 1, we screened 842 DEMs including 160 up-regulated DEMs and 563 down-regulated DEMs in GSE72199, and 27 DEMs including 15 up-regulated and 12 down-regulated DEMs in GSE54088(Fig. 1ab). Subsequently, in accordance with above criteria, we identified 325 DEGs, including 219 up-regulated DEGs and 106 down-regulated DEGs (Fig. 1c).
In order to obtain the common up-regulated and down-regulated DEMs in the two miRNA datasets, we took their intersection through the Venn diagram. In the end, we got 4 up-regulated DEMs and 5 down-regulated DEMs (Fig. 1de, Table 1).
Table 1
Common differentially expressed miRNAs obtained from GSE72199 and GSE54088.
Expression change
|
GSE72199
|
GSE54088
|
LogFC
|
P-value
|
LogFC
|
P-value
|
Down-regulation
|
1.07
|
0.0274
|
1.06
|
0.0353
|
Down-regulation
|
1.44
|
0.0425
|
1.23
|
0.0005
|
Down-regulation
|
1.22
|
0.0081
|
1.47
|
0.0045
|
Down-regulation
|
1.79
|
0.0254
|
4.05
|
0.0001
|
Up-regulation
|
-1.49
|
0.0229
|
-1.10
|
0.0054
|
Up-regulation
|
-1.26
|
0.0383
|
-1.06
|
0.0033
|
Up-regulation
|
-1.59
|
0.0315
|
-1.05
|
5.87E-05
|
Up-regulation
|
-1.11
|
0.0292
|
-1.03
|
0.0005
|
Up-regulation
|
-3.52
|
0.0069
|
-1.01
|
0.0045
|
GO functional analysis of DEGs.
For the sake of understanding the biological processes that DEGs participate in during CRC liver metastasis, we respectively uploaded the up-regulated DEGs and down-regulated DEGs to DAVID for GO functional analysis. The results show that up-regulated DEGs are mainly involved in negative regulation of endopeptidase activity, platelet degranulation, blood coagulation, inflammatory response, acute-phase response and innate immune response (Fig. 2a). As for those down-regulated DEGs, they mainly modulate cell adhesion, extracellular matrix organization, muscle contraction, negative regulation of canonical Wnt signaling pathway, collagen fibril organization, negative regulation of BMP signaling pathway (Fig. 2b).
KEGG pathway enrichment analysis of DEGs.
In the KEGG pathway enrichment analysis, we found that up-regulated DEGs are mainly enriched in Metabolic pathways, Complement and coagulation cascades, Chemical carcinogenesis, Retinol metabolism, Drug metabolism - cytochrome P450, Metabolism of xenobiotics by cytochrome P450, Bile secretion, Cholesterol metabolism, Drug metabolism - other enzymes and Steroid hormone biosynthesis (Fig. 3a). In terms of down-regulated DEGs, the pathway enrichment analysis results focused on Focal adhesion, Human papillomavirus infection, ECM-receptor interaction, Protein digestion and absorption, Vascular smooth muscle contraction, PI3K-Akt signaling pathway, PPAR signaling pathway, Wnt signaling pathway, cGMP-PKG signaling pathway and Proteoglycans in cancer (Fig. 3b).
Construction of PPI network and the most significant model.
We construct a PPI network containing 225 nodes and 1376 edges via STRING (Fig. 4a), with importing the results into Cytoscopy software on its heels. The first thing we did was to analyz the top ten genes under each of the 12 algorithms using the cytoHubba plugin, where nuclear receptor subfamily 1 group H member 4 (NR1H4), coagulation factor II, thrombin(F2), serpin family A member 10(SERPINA10) and complement C3(C3) were identified as hub genes that present in at least 6 algorithms. Then the MCODE plugin was used and the most remarkable model was obtained (Fig. 4b). For a deeper insight into the pathways involved in this model, we performed a KEGG pathway enrichment analysis and the results are shown in Table 2.
Table 2
The enrichment pathway of the model for PPI network.
Term
|
Count
|
P-value
|
Complement and coagulation cascades
|
6
|
1.29E-12
|
Metabolic pathways
|
4
|
0.0030
|
Cholesterol metabolism
|
2
|
0.0002
|
Staphylococcus aureus infection
|
2
|
0.0004
|
Systemic lupus erythematosus
|
2
|
0.0015
|
Phagosome
|
2
|
0.002
|
Herpes simplex virus 1 infection
|
2
|
0.0186
|
Potential target genes for DEMs.
For the purpose of predicting the target genes of DEMs, we used miRNet database to analyze the target genes of up-regulated DEMs and down-regulated DEMs, respectively. The network between DEMs and target genes was then demonstrated in miRNet's built-in visualization (Fig. 5ab). Immediately after, we used Venn diagrams to show common genes between target genes of down-regulated DEMs and up-regulated DEGs (Fig. 5c), and common genes between target genes of up-regulated DEMs and down-regulated DEGs were analyzed using the same method (Fig. 5d). 8 up-regulated target genes (Table 3) and 6 down-regulated target genes (Table 4) were obtained eventually.
Table 3
Target genes of upregulated miRNA.
miRNA
|
Targets
|
Gene name
|
hsa-miR-17-3p
|
ENPP2
|
ectonucleotide pyrophosphatase/phosphodiesterase 2
|
BNC2
|
basonuclin 2
|
hsa-miR-33b-5p
|
GAS1
|
growth arrest specific 1
|
ZEB1
|
zinc finger E-box binding homeobox 1
|
hsa-miR-122-5p
|
MYH11
|
myosin heavy chain 11
|
CALD1
|
caldesmon 1
|
Table 4
Target genes of downregulated miRNA.
miRNA
|
Targets
|
Gene name
|
hsa-miR-576-5p
|
CFHR3
|
complement factor H related 3
|
hsa-miR-329-3p
|
KYNU
|
kynureninase
|
|
SLCO1B3
|
solute carrier organic anion transporter family member 1B3
|
|
C6
|
complement C6
|
|
UGT2B4
|
UDP glucuronosyltransferase family 2 member B4
|
hsa-miR-379-5p
|
SERPINA1
|
serpin family A member 1
|
hsa-miR-137
|
ASGR2
|
asialoglycoprotein receptor 2
|
|
SERPINA3
|
serpin family A member 3
|
Validation of target DEGS expression levels.
The GEPIA database was thereafter utilized to validate the expression levels between 14 DEGs in CRC tissues and normal intestinal mucosa, and Fig. 6 shows 9 genes with significant expression differences. Of all the upregulated target DEGs, only the expression level of SERPINA3 in CRC was significantly lower than that in normal intestinal mucosa, which is contrary to what we observed in liver metastases. As for the down-regulated target DEGs, it is visible that their expression levels in CRC are significantly reduced.
Prognostic analysis of DEGs and DEMs.
We divided differentially expressed genes into high and low expression groups based on their median expression level in colorectal cancer tumors and then analyzed the differences in survival time between the two groups. The results showed that only GAS1 had a significant difference in survival outcome between the high and low expression groups. Moreover, the high expression group had shorter overall survival and worse prognosis compared to the low expression group (HR = 1.6, Logrank p = 0.037) (Fig. 7a). Subsequently, we used the online database OncoLnc to analyze the effect of DEMs on survival time in colon cancer (Fig. 7b) and rectal cancer (Fig. 7c) separately. Also divided into high and low expression groups using the median expression level as a limit. Due to the low expression of hsa-miR-122-5p in colorectal cancer and the low expression of hsa-miR-329-3p in colon cancer, we were unable to obtain the survival data of these two miRNAs in corresponding tumors. The results showed that only the expression level of hsa-miR-33b-5p had a significant effect on the survival time of patients with colon cancer.
Infiltration analysis of GAS1-associated immune cells.
Based on the results of DEGs survival analysis, we further investigated the relationship between GAS1 and immune cell infiltration in tumor tissues. Figure 8a and Fig. 8b respectively show the correlation between GAS1 and 6 immune cells in colon and rectal cancers after purity correction. The results showed that the low expression of GAS1 was negatively correlated with the distribution of CD4 + T cells, CD8 + T cells, macrophages, myeloid dendritic cells (DCs) and Neutrophils in both colon and rectal cancer, but positively correlated with B cells.