In Silico Analysis of Ion Channels and Their Correlation with Epithelial to Mesenchymal Transition in Breast Cancer

Abstract

Uncontrolled growth of breast cells due to altered gene expression is a key feature of breast cancer. Alterations in the expression of ion channels lead to variations in cellular activities, thus contributing to attributes of cancer hallmarks. Changes in the expression levels of ion channels were observed as a consequence of EMT. Additionally, ion channels were reported in the activation of EMT and maintenance of a mesenchymal phenotype. Here, to identify altered ion channels in breast cancer patients, differential gene expression and weighted gene co-expression network analyses were performed using transcriptomic data. Protein-protein interactions network analysis was carried out to determine the ion channels interacting with hub EMT-related genes in breast cancer. Thirty-two ion channels were found interacting with twenty-six hub EMT-related genes. The identified ion channels were further correlated with EMT scores, indicating mesenchymal phenotype. Further, the pathway map was generated to represent a snapshot of deregulated cellular processes by altered ion channels and EMT-related genes. Kaplan-Meier five-year survival analysis and Cox regressions indicated the expression of CACNA1B, ANO6, TRPV3, VDAC1 and VDAC2 to be potentially associated with poor survival. Deregulated ion channels correlate with EMT-related genes and have a crucial role in breast cancer-associated tumorigenesis. Most likely, they are potential candidates for the determination of prognosis in patients with breast cancer.

Keywords: RNA-Seq; bioinformatics; interaction networks; membrane proteins; microarray; prognosis.

Figure 1

Workflow for the identification of potential ion channels and its interaction with EMT-related genes in tumor and metastatic samples of patients with breast cancer.

Figure 2

Overlap of differentially expressed genes among HM, HT and TM. (A) Venn diagram depicting total number of differentially expressed ion channels in the RNA-Seq and two microarray datasets (GSE42568 and GSE52604) in HM (blue), HT (yellow) and TM (green). (B) Venn diagram showing the total number of differentially expressed EMT-related genes in the RNA-Seq and two microarray (GSE42568 and GSE52604) datasets in HM (blue), HT (yellow) and TM (green).

Figure 3

Co-expressed ion channel modules based on the non-redundant HT DEGs. (A) Soft-thresholding power to ensure scale-free network model: (i) A soft-thresholding power of 10 was chosen in the microarray dataset (GSE42568) corresponding to non-redundant HT DEGs; (ii) A soft-thresholding power of 7 was chosen in the RNA-Seq dataset corresponding to non-redundant HT DEGs (B) Hierarchical clustering of genes into modules. Modules are assigned different colors as depicted in the horizontal bar below the tree diagram: (i) 4 modules (yellow, blue, turquoise and brown) were obtained in the microarray dataset (GSE42568); (ii) 4 modules (turquoise, brown, yellow and blue) were obtained in RNA-Seq. (C) Correlation between module eigengenes and binary traits—normal, tumor and metastatic. Rows correspond to modules depicted as different colors and columns are the binary traits. Numbers in each cell are the correlation coefficient between module eigengenes and the binary traits and the corresponding p-value. (i) The yellow module was chosen as a significant module from the microarray dataset (GSE42568). (ii) The turquoise module was chosen as a significant module from RNA-Seq. (D) Scatter plot of gene significance (GS) for the binary trait vs. the module membership (MM) in the selected module. (i) GS-MM plot of the yellow module in the microarray (GSE42568) dataset. (ii) GS-MM plot of the turquoise module in the RNA-Seq dataset.

Figure 4

Representation of protein–protein interaction networks (PPINs) of combined gene set. Nodes in the shade of blue represent EMT-related genes and yellow nodes represent ion channels. The nodes are arranged based on the degree centrality measure. Larger nodes represent nodes with high degree centrality. (A) PPINs for microarray dataset (GSE42568) based on the list of non-redundant HT DEGs. (B) PPINs for microarray dataset (GSE42568) based on the list of non-redundant HM DEGs.

Figure 5

Correlation between GS76, MLR and KS scoring methods for samples from different platforms. (A) Correlation between GS76, MLR and KS methods for samples from RNA-Seq TCGA data. (B) Correlation between GS76, MLR and KS methods for samples from GSE42568 microarray data. (C) Correlation between GS76, MLR and KS methods for samples from GSE52604 microarray data.

Figure 6

Correlation between ion channels (ICs) and EMT scores obtained using GS76, MLR and KS methods. In the bubble plot, the x-axis consists of correlation values of ion channels with the KS method, and the y-axis consists of correlation values with the MLR method. Each bubble corresponds to a particular ion channel represented by different colors. The size of the bubble corresponds to the correlation values of ion channels with GS76. TGCA data were used as RNA-Seq data and GSE42568 and GSE52604 were microarray datasets. (A–C) Correlation of ICs identified as interacting with EMT-related genes with GS76, MLR and KS scores; (D–F) Correlation of ICs identified in tumor state with GS76, MLR and KS scores; (G–I) Correlation of ICs identified in metastatic state with GS76, MLR and KS scores.

Figure 7

Depiction of events that may occur in breast cancer patients upon dysregulation of ion channels with EMT-related genes. The pathway map was generated using PathVisio (v3.3.0). Cancer phenotypes may appear due to alterations in important cellular processes, such as calcium signaling, insulin secretion, adipocyte metabolism, nitric oxide signaling and glutamatergic signaling.

Figure 8

Kaplan–Meier 5-year survival curves representing the prognostic relationship between high and low expression of ion channels identified in breast cancer with survival probability. (A) CACNA1B, (B) ANO6, (C) TRPV3, (D) VDAC1, and (E) VDAC2.