3.1 Screening of differentially expressed genes
The differentially expressed genes and lncRNAs in 12 breast cancer bone metastasis samples and 12 normal samples were identified using the R limma package for differential expression. We used breast cancer bone metastasis expression profile GSE66206. A total of 133 differentially expressed genes and 23 differentially expressed lncRNAs (fold change >1.2 or fold change <5/6) were identified. Among them, 27 genes and 3 lncRNAs were significantly differentially expressed (p < 0.05).
3.2 Constructing an interaction network related to bone metastasis in breast cancer
We calculated Spearman correlation coefficients for the co-expression of PCGs and lncRNAs, which we corrected using the rank-sum test. We identified 747,870 PCG–lncRNA, PCG–PCG, and lncRNA–lncRNA interactions, of which 7,451 were differentially expressed PCG–lncRNA interactions (| r | ≥ 0.3, p < 0.05). We then visualized the network (Figure 1A) and analyzed the node degrees. The Cytoscape analysis revealed that the degree of the network nodes ranged from 1 to 905. In order to identify the hub nodes in the network, we set the threshold to 20 and extracted all possible core driver genes. Twenty-eight different core driver genes were identified (Figure 1B and Table 1). We also identified regulatory interactions between some differentially expressed genes and lncRNAs (Figure 1C). We consulted the Human Protein Atlas Interactive Analysis to validate the protein and RNA levels of the core driver genes in breast cancer, including the levels of IDI1, BNIP3, IFRD1, COQ10B, and ZBTB10 (Figure 2). Immunohistochemistry and RNA expression analysis indicated that IDI1, BNIP3, IFRD1, and ZBTB10 were significantly differentially expressed in cancer and healthy tissues.
Table 1. The Gene Ontology annotation of the core driver genes of breast cancer bone metastasis.
Gene symbol
|
Gene ID
|
Biological Process (GO)
|
CSPP1
|
79848
|
GO:0051781 positive regulation of cell division
|
COQ10B
|
80219
|
GO:0006743 ubiquinone metabolic process
|
FLVCR1
|
28982
|
GO:0043249 erythrocyte maturation
|
IDI1
|
3422
|
GO:0050993 dimethylallyl diphosphate metabolic process
|
PTP4A1
|
7803
|
GO:0030335 positive regulation of cell migration
|
ODC1
|
4953
|
GO:0009445 putrescine metabolic process
|
JMJD1C
|
221037
|
GO:0033169 histone H3-K9 demethylation
|
GABRB2
|
2561
|
GO:0071420 cellular response to histamine
|
IRF1
|
3659
|
GO:0034124 regulation of MyD88-dependent toll-like receptor signaling pathway
|
KDM6B
|
23135
|
GO:0071557 histone H3-K27 demethylation
|
IFRD1
|
3475
|
GO:0048671 negative regulation of collateral sprouting
|
RBM25
|
58517
|
GO:0000381 regulation of alternative mRNA splicing, via spliceosome
|
NCOA4
|
8031
|
GO:0006879 cellular iron ion homeostasis
|
SERPINA1
|
5265
|
GO:0048199 vesicle targeting, to, from or within Golgi
|
HSPB1
|
3315
|
GO:0038033 positive regulation of endothelial cell chemotaxis by VEGF-activated vascular endothelial growth factor receptor signaling pathway
|
BNIP3
|
664
|
GO:1902109 negative regulation of mitochondrial membrane permeability involved in apoptotic process
|
LARP4
|
113251
|
GO:0034250 positive regulation of cellular amide metabolic process
|
UXT
|
8409
|
GO:0047497 mitochondrion transport along microtubule
|
ABHD2
|
11057
|
GO:0048240 sperm capacitation
|
LAMTOR1
|
55004
|
GO:0060620 regulation of cholesterol import
|
MRPL51
|
51258
|
GO:0070126 mitochondrial translational termination
|
ELF2
|
1998
|
GO:0050855 regulation of B cell receptor signaling pathway
|
TPD52
|
7163
|
GO:0030183 B cell differentiation
|
PAQR5
|
54852
|
GO:0048477 oogenesis
|
3.1 WGCNA of co-expression
Gene expression profiles for 1,533 PCG and lncRNA genes related to breast cancer bone metastases were constructed with the WGCNA R package to build a co-expression network (merge Cut Height = 0.25, verbose = 3). The optimal threshold for the WGCNA was 12 (Figure 3A). The sample clustering chart is shown in Figure 3B. We excavated 5 modules from the breast cancer bone metastasis gene expression profile, including the blue (112 genes), turquoise (1,274 genes), brown (37 genes), yellow (20 genes), and grey (90 genes) modules (Figure 3C). The results of the co-expression analysis are shown in Figure 3D. The relationships between the various modules and traits related to bone metastases in breast cancer are shown in Figure 3E. The genes in the yellow module are related to the disease characteristics of breast cancer bone metastases and breast cancer. In addition, the occurrence of bone metastases positively correlated with the expression of the genes in the yellow module.
3.2 Functional enrichment analysis and verification
The yellow module comprises 14 genes (Table 2). We used KOBAS to perform functional enrichment analysis using gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) for the genes in the yellow module [6]. With a significance threshold of p < 0.05, we identified 7 significant KEGG functional pathways, which revealed that the core genes in the yellow module are enriched in the prolactin signaling pathway [7] (Figure 4A). GO enrichment analysis revealed 300 related GO functions (p < 0.05) involving a large number of pathways that may be related to cancer development, including regulation of signal transduction (GO: 0009966, p = 0.00017), cellular processes (GO: 0009987, p = 0.00024), cellular communication regulation (GO: 0010646, p = 0.00030), and positive regulation of biological processes (GO: 0048518, p = 0.000826) (Figure 4B).
Enrichr was used to enrich the KEGG and GO functions of the 5 lncRNAs in the yellow module (p < 0.05). Some pathways enriched for the lncRNAs were related to amino acid transport (Figure 4C and D), including alanine transport (GO: 0015808, p = 0.00249) and amino acid transmembrane transport (GO: 0003333, p = 0.00747). We identified the GABAergic synapse (p = 0.02205) in the lncRNA-enriched KEGG pathway, which is closely related to breast cancer metastasis [8, 9]. Therefore, the yellow module is likely to be a pathogenic module.
The lncRNA RP11-317-J19.1 and the gene PTP4A1 act on the gene BNIP3, which plays a key role in cell apoptosis and autophagy (Figure 4E and F) [10-13]. The gene–lncRNA interaction diagram reveals a subtle interaction between the 5 differentially expressed genes and differentially expressed lncRNAs in the yellow module (Figure 5A and B). To validate the gene–lncRNA interactions in the yellow module, we created gene correlation scatter plots for IDI2-AS1 and IDI1, IDI2-AS1 and PTP4A, IDI2-AS1 and BNIP3, PTP4A and IDI1, BNIP3 and IDI1, and PTP4A and BNIP3 (Figure 5C). The results revealed that expression of BNIP3 and PTP4A significantly positive correlated with that of IDI1 (p < 0.05), whereas expression of BNIP3 and PTP4A significantly negatively correlated with that of IDI2-AS1 (p < 0.05). Moreover, Kaplan Meier-plotter analysis showed that 7 genes and 2 lncRNAs in the core driver gene network and co-expression yellow module correlated with overall survival, including IDI1, PTP4A1, BNIP3, IFRD1, ZBTB10, DISP1, COQ10B, IDI2-AS1, and SLC38A3 (Figure 6A–I).
Table 2. The GO annotations of the genes in the yellow module.
Gene symbol
|
Gene ID
|
Biological Process (GO)
|
IDI1
|
3422
|
GO:0050993 dimethylallyl diphosphate metabolic process
|
PTP4A1
|
7803
|
GO:0030335 positive regulation of cell migration
|
BNIP3
|
664
|
GO:1902109 negative regulation of mitochondrial membrane permeability involved in apoptotic process
|
NCKAP5
|
344148
|
GO:0007019 microtubule depolymerization
|
SOCS2
|
8835
|
GO:0060396 growth hormone receptor signaling pathway
|
ADAMTS3
|
9508
|
GO:1900748 positive regulation of vascular endothelial growth factor signaling pathway
|
ANXA11
|
311
|
GO:0032506 cytokinetic process
|
DYRK1A
|
1859
|
GO:0043518 negative regulation of DNA damage response, signal transduction by p53 class mediator
|
EPS8L2
|
64787
|
GO:1900029 positive regulation of ruffle assembly
|
DISP1
|
84976
|
GO:0007225 patched ligand maturation
|
WDFY3
|
23001
|
GO:0035973 aggrephagy
|
SLC38A5
|
92745
|
GO:1904557 L-alanine transmembrane transport
|
NOVA1
|
4857
|
GO:0120163 negative regulation of cold-induced thermogenesi
|
SLC38A3
|
10991
|
GO:2000487 positive regulation of glutamine transport
|