Co-expression network analysis identified atypical chemokine receptor 1 (ACKR1) association with lymph node metastasis and prognosis in cervical cancer

Abstract

Cervical cancer (CC) is one kind of female cancer. With the development of bioinformatics, targeted specific biomarkers therapy has become much more valuable. GSE26511 was obtained from gene expression omnibus (GEO). We utilized a package called “WGCNA” to build co-expression network and choose the hub module. Search Tool for the Retrieval of Interacting Genes Database (STRING) was used to analyze protein-protein interaction (PPI) information of those genes in the hub module. A Plug-in called MCODE was utilized to choose hub clusters of PPI network, which was visualized in Cytoscape. Clusterprofiler was used to do functional analysis. Univariate and multivariate cox proportional hazards regression analysis were both conducted to predict the risk score of CC patients. Kaplan-Meier curve analysis was done to show the overall survival. Receiver operating characteristic (ROC) curve analysis was utilized to evaluate the predictive value of the patient outcome. Validation of the hub gene in databases, Gene set enrichment analysis (GSEA) and GEPIA were completed. We built co-expression network based on GSE26511 and one CC-related module was identified. Functional analysis of this module showed that extracellular space and Signaling pathways regulating pluripotency of stem cells were most related pathways. PPI network screened GNG11 as the most valuable protein. Cox analysis showed that ACKR1 was negatively correlated with CC progression, which was validated in Gene Expression Profiling Interactive Analysis (GEPIA) and datasets. Survival analysis was performed and showed the consistent result. GSEA set enrichment analysis was also completed. This study showed hub functional terms and gene participated in CC and then speculated that ACKR1 might be tumor suppressor for CC.

Keywords

Bioinformatics analysis lymph node metastasis atypical chemokine receptor 1 (ACKR1)cervical cancer

Abbreviations

CC	Cervical cancer
GEO	Gene expression omnibus
WGCNA	Weighted gene expression network analysis
DEGs	Differentially expressed genes

ACKR1	Atypical chemokine receptor 1
STRING	Search Tool for the Retrieval of Interacting Genes Database
PPI	Protein-protein interaction
MCODE	Molecular Complex Detection
ROC	Receiver operating characteristic
GSEA	Gene set enrichment analysis
TOM	Topological overlap matrix

Figure 1.

Study design and data preprocessing. Note: (A) Flow diagram of study. (B) Sample clustering after detecting two outliers.

1. Introduction

Cervical cancer is the fourth-most prevalent female cancer worldwide, and contributes to estimated 530,000 new cases annually with 270,000 deaths [1]. The persistent infection of human papilloma virus especially the type 16/18, increasing number of lifetime sexual partners, acquisition of new male partners and non-monogamous male partners have been identified as high-risk factors for cervical cancer [2]. In recent years, bioinformatics has developed rapidly. We can use the database to analyze tumor-related molecules and have guiding significance for targeted therapy [3].

Gene expression profiling [4] and weighted gene expression network analysis (WGCNA) [5] were currently the most widely used analytical method in bioinformatics. In this research, we established co-expression network based on dataset GSE26511 including 19 lymph node-positive CC samples and 20 lymph node-negative CC samples in order to determine the hub module. Then we built PPI network for the genes in hub module. We further performed survival regression analysis to screen out hub genes and finally validate hub genes using databases and website. Combined with the above results, it is concluded that they may be potential targets for CC.

2. Materials and methods

2.1 Research plan and data pre-calculation

The research was designed based on the flow chart (Fig. 1A). Gene expression profile matrix files were obtaied from the Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo/). And the platform of the GSE26511 [6] dataset was the GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array [7]. This dataset contained 19 lymph node-positive CC tissue samples and 20 lymph node-negative CC tissue samples. The series matrix and platform TXT files were both downloaded. A Perl language command was applied to transformed the gene probe IDs to the gene symbols to obtain a matrix files including the international standard gene name. Normalized data from GSE26511 contained a total of 21655 genes. R package was used to screen genes in the top 25% of variance and 5413 genes were identified for further analysis. We next constructed PPI and co-expression networks so as to confirm hub terms and genes of this study.

2.2 Construction and analysis of WGCNA

First, we check the expression data profile of the genes for their rationality. The two obvious outliers (GSM651842, GSM651846) and removed them from the cohort . Obvious outliers means sample clustering according to the distance between different samples in Pearson’s correlation matrices and average linkage. Secondly, we used the R package called “WGCNA” to build the co-expression network based on the genes [8, 9]. The Pearson’s correlation matrices was prepared for all pair-wise genes. Next, a power function amn $=$ $|$ cmn $|\beta$ (cmn $=$ Pearson’s correlation between gene m and n; amn $=$ adjacency between gene m and n) was chosen to erect weighted adjacency matrix. Then, we transformed the adjacency to topological overlap matrix (TOM), which could measure network connectivity of a gene. It could be defined to be the sum of its adjacency with all other genes for network generation. We created average linkage hierarchical clustering due to the TOM-based dissimilarity measure with a minimum size (gene group) of 50 for the genes dendrogram. Therefore genes which owned similar expression profiles might be classified into same module. Eigengenes were the main parts of microarray gene expression data, which represented the oscillatory characteristics of gene expression patterns [10]. We calculated the dissimilarity of module eigengenes, selecting a cut line for module dendrogram and so that we can merge them. We tried to choose modules which could influence CC, thereby functional analysis was done based on the modules.

2.3 PPI network and hub cluster functional analysis

We utilized STRING (Search Tool for the Retrieval of Interacting Genes Database) (http://www.string-db.org/) to evaluate PPI (protein-protein interaction) information for hub modules [11]. Further, we exported the results from STRING database to Cytoscape software for visual display [12]. MCODE plug-in was functioned for searching clustered sub-networks [13]. We set the following default parameters: node score cut-off $=$ 0.2, k-core $=$ 2, degree cut-off $=$ 5, and max. depth $=$ 100. Clusterprofiler was conducted for functional analysis of hub clusters.

2.4 Construction of a prognostic signature

Univariate Cox proportional hazards regression analyses were performed based on the genes in hub module. Genes which had prognostic signature were identified by the cutoff of $p<$ 0.05. Next, a multivariate Cox proportional hazards regression model was built. Cox proportional hazards regression with a cutoff of $p<$ 0.05 was constructed to predict the risk score of CC patient based on the expression of the hub genes. Compared with mean risk score, patients with different risk score were divided into a low-risk group and high-risk group. Kaplan-Meier curve analysis was done to show the overall survival of the low-risk group and high-risk group. $P<$ 0.05 was seen to be significant. Receiver operating characteristic (ROC) curve analysis using 5 years as the predicted time was also completed to evaluate the predictive value of the CC patient outcome. The area under the ROC curve along with sensitive and specificity were identified as evaluation indices.

Figure 2.

Hub module selection. Note: (A) Dendrogram of all genes clustered based on a dissimilarity measure (1-TOM). (B) Dendrogram of consensus module eigengenes obtained by WGCNA on the consensus correlation. The red line was the merging threshold, and groups of eigengenes below the threshold represent modules whose expressions profiles should be merged due to their similarity. (C) Correlation between modules and traits. The upper number in each cell referred to the correlation coefficient of each module in the trait, and the lower number was the corresponding $p$ -value. Among them, the green module was the most relevant modules with cancer traits.

2.5 Validation of hub genes

The hub gene needed to be verified in GEPIA (Gene Expression Profiling Interactive Analysis) [14]. GSE39001 and GSE6791 were also used to validate hub gene expression. GSE39001 [15] included 43 CC tissue samples and 12 normal tissue samples. GSE6791 [16] included 20 CC tissue samples and 8 normal tissue samples.

2.6 Gene set enrichment analysis (GSEA)

Based on hub gene expression level, samples were divided into two groups. In order to study the potential function of this gene, GSEA (http://software.broadinsti tute.org/gsea/index.jsp) [17] was used to research a series of priori defined biological pathway which might enriched in the gene rank derived from hub gene between the two groups. With FDR $<$ 0.05, Terms which were enriched in all real hub genes were exhibited.

Figure 3.

Select hub genes in hub modules. Note: (A) A scatter plot of GS for CC versus the MM in the green module. Intramodular analysis of the genes found in the green module, which contained genes that had a high correlation with CC, with $p<$ 3.9e-38 and correlation $=$ 0.53. (B) A heatmap of all genes. The intensity of the red color indicated the strength of the correlation between pairs of modules on a linear scale. (C) Dendrogram of merged module eigengenes obtained by WGCNA. Heatmap plot of the adjacencies of modules. Red represented high adjacency (positive correlation) and blue represented low adjacency (negative correlation).

2.7 Statistical analysis

All data was processed by R version 3.5.1 (https://www.r-project.org/) with the following packages: “ggstatsplot”, “ggplot2” and WGCNA. All analysis were done three times and showed data from three separate experiments. Two-tailed Student’s $t$ -test was utilized for significance of differences between subgroups. Statistical analysis processed by SPSS 16.0. Statistical significance was set with values $p<$ 0.05.

Figure 4.

Cluster analysis of the PPI network. Note: (A) The genes in green module were filtered into the PPI network that contained 476 nodes and 298 edges. (B) Histogram of key proteins. The y-axis represented the name of genes, the x-axis represented the number of adjacent genes, and height was the number of gene connections.

3. Results

3.1 Weighted co-expression network construction and analysis

Expression data profile of GSE26511 including 19 lymph node-positive CC tissue samples and 18 lymph node-negative CC tissue samples was used to construct co-expression networks and identify hub genes (Fig. 1B). We conducted the R package “WGCNA” and set the power $\beta$ to a soft-thresholding parameter $\beta=$ 8 (scale free $R^{2}=$ 0.84) to ensure a scale-free network (Fig. S1A–D). We pre-calculated the data so as to deeply analyzed the modules which owned mostly related genes, and 8 modules were found (Fig. 2A). MEDissThres was set as 0.25 in order to merge similar modules (Fig. 2B). Eight modules were produced. Black module included 257 genes. Blue module included 1037 genes. brown module included 986 genes. green module included 508 genes. grey module included 342 genes. red module included 408 genes. turquoise module included 1037 genes. yellow module included 839 genes. Genes in gray module did not belong to any of the functional modules, therefore needed not to be analyzed. Besides, an intramodular analysis of Gene significance (GM) and module membership (MM) of the genes in the 7 modules mentioned was conducted. As GS and MM illustrated meaningful correlation, the result intimated that green module was most strongly related to CC lymph node metastases among the 7 modules (Fig. 2C). Moreover, genes belong to module green had the highest negative correlation with lymph node metastases (Fig. 3A). The genes selected were identified in the heatmap (Fig. 3B).

Figure 5.

Survival prognosis model of the hub genes. Note: (A) Multivariate Cox proportional hazards regression analysis further screened 2 hub genes. (B) Survival analysis showed that the patients in the high risk group had statistically significantly worse overall survival than those in low risk group in TCGA cohort. (C) ROC analysis was performed to find out the most optimal cutoff value to divide the CC patients into high risk and low risk group. (D and E) The risk scores for all patients in TCGA cohort are plotted in ascending order and marked as low risk (blue) or high risk (red), as divided by the threshold (vertical black line). (F) Five expression and risk score distribution in TCGA cohort by z-score, with red indicating higher expression and light blue indicating lower expression.

It was interesting that some of these modules had similar expression profiles. In order to research the relationships and interactions among the mentioned 8 co-expressed modules, we analyzed the connectivity of eigengenes and completed a cluster analysis. Finally 7 clusters were divided into 3 clusters. one contained three branches and two contained two branches (Fig. 3C).

Functional analysis were performed on the genes in green module. Go analysis showed that genes were mostly enriched in extracellular spaceï¼Œproteinaceous extracellular matrix and extracellular region (Fig. S2A). KEGG analysis showed that genes were mostly enriched in Signaling pathways regulating pluripotency of stem cells, Melanoma and TGF-beta signaling pathway (Fig. S2B).

3.2 PPI network and analysis of hub clusters

Via the STRING website, 508 genes in green module mentioned were screened into PPI network. The complex owned 476 nodes and 298 edges (Fig. 4A). STRING software identified 10 prominent proteins. Among them GNG11 was seen to be the most important one connecting 20 nodes (Fig. 4B).

Figure 6.

Validation of hub gene. Note: (A) The result of GEPIA website validation. (B) The validation of ACKR1 expression in GSE39001. (C) The validation of ACKR1 expression in GSE6791. (D) Survival analysis of ACKR1.

Then, we applied the plug-in named MCODE for scoring. So that we could find parameters. The parameters were able to show the qualified results for PPI. After that, clusters could be found in the PPI network. Eighteen clusters were produced due to k-core $=$ 2. Among those clusters, cluster 1, owned 16 nodes and 65 edges, got the highest score (Fig. S3A). Cluster 2, included 6 nodes and 15 edges, got the secondly high score (Fig. S3B). Cluster 3, had 5 nodes and 10 edges, got the thirdly high score (Fig. S3C). These results indicated that genes in the top 3 clusters might be related with CC lymph node metastases.

MCODE was applied to analyze the top 3 significant clusters. Functional analysis of the 2 clusters were also conducted using clusterprofiler (Fig. S4A–D).

Figure 7.

GSEA using GSE26511. Note: Only listed the 4 most common functional gene sets enriched in CC samples with ACKR1 highly expressed.

3.3 Construction of a prognostic signature

We set GM (cor.geneModuleMembership) $>$ 0.2 and MM (cor.geneTraitSignificance) $>$ 0.8 as the criteria, 44 genes from the green module were screened. Univariate Cox proportional hazards regression analysis was performed on the 44 genes. We set $p<$ 0.05 as the criteria, five genes were screened including ACKR1, SERPINF1, NAP1L2, EMCN and PRKAR2B (Fig. S5). Multivariate Cox proportional hazards regression analysis was further performed on the 5 genes, which screened SERPINF1 and ACKR1. The risk assessment score for predicting overall survival was calculated: Risk score $=$ $-$ 0.164*SERPINF1 $-$ 0.154* ACKR1. According to the $p<$ 0.05, ACKR1 was chosen as the hub gene, which was identified as tightly linked to CC prognosis and lymph node metastasis (Fig. 5A). Based on the risk score, CC lymph node metastases patients were divided into low-risk and high-risk groups. Kaplan-Meier survival analysis suggested that low-risk patients had much more significantly better overall survival than high-risk patients in TCGA cohort (Fig. 5B). ROC curve analysis was also completed according to the 5-year survival of the area under the receiver operating characteristic curve (AUC) value, the specificity and sensitivity were both highest when the risk score was 0.657 (Fig. 5C). The distribution of risk score, survival status, and the expression of five genes for each patient were also analyzed (Fig. 5D–F).

3.4 Hub gene validation

To further demonstrate ACKR1 was related to CC, we used the GEPIA website to verify it (Fig. 6A). ACKR1 was proved to be positively correlated with CC prognosis and interestingly had higher expression level in normal cervical tissues than in cervical cancer tissues. The expression of ACKR1 in GSE39001 and GSE6791 showed the consistent result (Fig. 6B–C). Survival analysis of ACKR1 showed that patients with higher ACKR1 expression had better overall survival (Fig. 6D).

3.5 Gene set enrichment analysis (GSEA)

To identify the potential function of the ACKR1 in CC, GSEA was conducted to search KEGG pathways enriched in the highly-expressed samples. Four gene sets ( $n=$ 306) “abc transporters”, “cell adhesion molecules cams”, “chemokine signaling pathway” and “cytokine receptor interaction” were enriched (FDR $<$ 0.05) (Fig. 7A–D).

4. Discussion

Cervical cancer is the fourth-most prevalent cancer in women worldwide, and contributes to estimated 530,000 new cases annually with 270,000 deaths. Its development has a complex mechanism, and many genes are involved in its regulation. In recent years, bioinformatics has developed rapidly. We can use the database to analyze tumor-related molecules and have guiding significance for targeted therapy. In this research, we established co-expression network based on dataset GSE26511 including 19 lymph node-positive CC samples and 20 lymph node-negative CC samples in order to determine the hub module. Then we built PPI network for the genes in hub module. We further performed survival regression analysis to screen out hub genes and finally validate hub genes using databases and website. Combined with the above results, it is concluded that they may be potential targets for CC.

In these findings, WGCNA analysis suggested that module green possessed significantly relevant expression pattern. Functional analysis were performed on the genes in green module. Go analysis showed that genes were mostly enriched in extracellular space. KEGG analysis showed that genes were mostly enriched in Signaling pathways regulating pluripotency of stem cells.

Wu et al. found that extracellular vesicles might be emerging targets in cancer [18]. Pickup et al. found that the extracellular matrix modulated the hallmarks of cancer [19]. Song et al. The diversity of characteristics of pancreatic cancer stem cells via regulating sonic hedgehog signaling pathway [20]. The role of these pathways in CC has not been studied and deserved further exploration.

A PPI network was built based on the green module. According to a k-core of 2, we found 7 clusters which might have impact on CC. Ten prominent proteins were identified, among the 10 proteins GNG11 was seen to be the most important one, which suggested that GNG11 might have inhibition effect in CC.

GNG11 (G-protein subunit $\gamma$ 11) forms a complex with the G $\beta$ 1 subunit and shuttles between the plasma membrane and the Golgi apparatus. Takauji et al. found that GNG11 suppressed cell growth with induction of reactive oxygen species and abnormal nuclear morphology [21]. However, the role of GNG11 in tumors, especially cervical cancer, has not yet been elucidated.

We combined survival regression analysis and then identified that ACKR1 was negatively associated with CC. The validation of GEPIA showed that ACKR1 was positively correlated with CC prognosis. ACKR1 was demonstrated to have higher expression level in normal cervical tissues than in cervical cancer tissues. The expression of ACKR1 in GSE39001 and GSE6791 showed the consistent result.

Atypical chemokine receptor 1 (ACKR1) was known to encode a glycoprotein expressing the Duffy blood group antigens [22], which was proved to regulate hematopoiesis on nucleated erythroid cells [23]. Wan W et al. proved that ACKR1 deficiency reduced atherogenesis in ApoE-knockout mice [24]. The functional mechanism between ACKR1 and CC deserved further exploration.

GSEA results also provided some valuable details. For example, ABC transporters were proved to be mediators of drug resistance and contributors to cancer cell biology [25]. Cell adhesion molecules cams have been shown to be involved in tumor metastasis and had a certain value in tumor therapy [26, 27]. Xiao et al. found that CXCL16/CXCR6 chemokine signaling mediated breast cancer progression [28]. Chemokine pathways were also proven to be involved in the process of hepatobiliary cancer [29]. These pathways were closely related to tumors, and their roles in CC were worth further research.

5. Conclusion

In summary, our study illustrated the hub genes and pathways involved in the progression of CC. We finally speculated that ACKR1 should be a tumor suppressor for predicting and treating cervical cancer.

Supplementary data

The supplementary files are available to download from http://dx.doi.org/10.3233/CBM-190533.

sj-docx-1-cbm-10.3233_CBM-190533.docx - Supplemental material

Supplemental material, sj-docx-1-cbm-10.3233_CBM-190533.docx

Footnotes

Acknowledgments

This work was supported by the National Nature Science Foundation of China (81872119 and 81472442) and the Jiangsu province medical innovation team (CXTDA2017008).

Conflict of interest

The authors report no conflicts of interest in this work.

References

Small

, Jr. Bacon

M.A.

Bajaj

Chuang

L.T.

Fisher

B.J.

Harkenrider

M.M.

Jhingran

Kitchener

H.C.

Mileshkin

L.R.

Viswanathan

A.N.

and Gaffney

D.K.

, Cervical cancer: A global health crisis, Cancer 123 (2017), 2404–2412.

Chelimo

Wouldes

T.A.

Cameron

L.D.

and Elwood

J.M.

, Risk factors for and prevention of human papillomaviruses (HPV), genital warts and cervical cancer, J Infect 66 (2013), 207–17.

Tsimberidou

A.M.

, Targeted therapy in cancer, Cancer Chemother Pharmacol 76 (2015), 1113–32.

Xia

Shen

Meng

Yan

Zhang

Chen

Wang

and Xu

, ANLN functions as a key candidate gene in cervical cancer as determined by integrated bioinformatic analysis, Cancer Manag Res 10 (2018), 663–670.

Langfelder

and Horvath

, WGCNA: an R package for weighted correlation network analysis, BMC Bioinformatics 9 (2008), 559.

Zhang

Xiao

Liang

Liao

Xiang

Cao

Chen

and Zhou

, Identification of differentially expressed genes in cervical cancer by bioinformatics analysis, Oncol Lett 16 (2018), 2549–2558.

Harbig

Sprinkle

and Enkemann

S.A.

, A sequence-based identification of the genes detected by probesets on the Affymetrix U133 plus 2.0 array, Nucleic Acids Res 33 (2005), e31.

Mason

M.J.

Fan

Plath

Zhou

and Horvath

, Signed weighted gene co-expression network analysis of transcriptional regulation in murine embryonic stem cells, BMC Genomics 10 (2009), 327.

Horvath

and Dong

, Geometric interpretation of gene coexpression network analysis, PLoS Comput Biol 4 (2008), e1000117.

10.

Alter

Brown

P.O.

and Botstein

, Singular value decomposition for genome-wide expression data processing and modeling, Proc Natl Acad Sci U S A 97 (2000), 10101–6.

11.

Szklarczyk

Franceschini

Wyder

Forslund

Heller

Huerta-Cepas

Simonovic

Roth

Santos

Tsafou

K.P.

Kuhn

Bork

Jensen

L.J.

and von Mering

, STRING v10: protein-protein interaction networks, integrated over the tree of life, Nucleic Acids Res 43 (2015), D447–52.

12.

Shannon

Markiel

Ozier

Baliga

N.S.

Wang

J.T.

Ramage

Amin

Schwikowski

and Ideker

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

13.

Wang

L.G.

Han

and He

Q.Y.

, clusterProfiler: An R package for comparing biological themes among gene clusters, Omics 16 (2012), 284–7.

14.

Tang

Kang

Gao

and Zhang

, GEPIA: a web server for cancer and normal gene expression profiling and interactive analyses, Nucleic Acids Res 45 (2017), W98–W102.

15.

Zhou

Shi

Cao

Dinh

T.L.A.

Wan

and Zhao

, Investigation of differentially-expressed microRNAs and genes in cervical cancer using an integrated bioinformatics analysis, Oncol Lett 13 (2017), 2784–2790.

16.

Zhao

L.P.

R.H.

Han

D.M.

Zhang

X.Q.

Nian

G.X.

M.X.

Feng

Zhang

and Sun

Z.G.

, Independent prognostic Factor of low-expressed LncRNA ZNF667-A Fig. S for cervical cancer and inhibitory function on the proliferation of cervical cancer, Eur Rev Med Pharmacol Sci 21 (2017), 5353–5360.

17.

Subramanian

Kuehn

Gould

Tamayo

and Mesirov

J.P.

, GSEA-P: A desktop application for Gene Set Enrichment Analysis, Bioinformatics 23 (2007), 3251–3.

18.

Xing

S.Y.

and Watabe

, Extracellular vesicles as emerging targets in cancer: Recent development from bench to bedside, Biochim Biophys Acta Rev Cancer 1868 (2017), 538–563.

19.

Pickup

M.W.

Mouw

J.K.

and Weaver

V.M.

, The extracellular matrix modulates the hallmarks of cancer, EMBO Rep 15 (2014), 1243–53.

20.

Song

Chen

Gao

Zhang

Wang

and Liu

, Ski modulate the characteristics of pancreatic cancer stem cells via regulating sonic hedgehog signaling pathway, Tumour Biol (2016).

21.

Takauji

Kudo

Matsuo

Hossain

M.N.

Nakabayashi

Miki

Fujii

and Ayusawa

, GNG11 (G-protein subunit gamma 11) suppresses cell growth with induction of reactive oxygen species and abnormal nuclear morphology in human SUSM-1 cells, Biochem Cell Biol 95 (2017), 517–523.

22.

Yin

Srivastava

Gebremedhin

Makuria

A.T.

and Flegel

W.A.

, Long-range haplotype analysis of the malaria parasite receptor gene ACKR1 in an East-African population, Hum Genome Var 5 (2018), 26.

23.

Duchene

Novitzky-Basso

Thiriot

Casanova-Acebes

Bianchini

Etheridge

S.L.

Hub

Nitz

Artinger

Eller

Caamano

Rulicke

Moss

Megens

R.T.A.

von Andrian

U.H.

Hidalgo

Weber

and Rot

, Atypical chemokine receptor 1 on nucleated erythroid cells regulates hematopoiesis, Nat Immunol 18 (2017), 753–761.

24.

Wan

Liu

Lionakis

M.S.

Marino

A.P.

Anderson

S.A.

Swamydas

and Murphy

P.M.

, Atypical chemokine receptor 1 deficiency reduces atherogenesis in ApoE-knockout mice, Cardiovasc Res 106 (2015), 478–87.

25.

Fletcher

J.I.

Williams

R.T.

Henderson

M.J.

Norris

M.D.

and Haber

, ABC transporters as mediators of drug resistance and contributors to cancer cell biology, Drug Resist Updat 26 (2016), 1–9.

26.

Xin

Dong

X.W.

and Guo

X.L.

, Role of the interaction between galectin-3 and cell adhesion molecules in cancer metastasis, Biomed Pharmacother 69 (2015), 179–85.

27.

Okegawa

Pong

R.C.

and Hsieh

J.T.

, The role of cell adhesion molecule in cancer progression and its application in cancer therapy, Acta Biochim Pol 51 (2004), 445–57.

28.

Xiao

Wang

Cheng

Hao

Sun

Xue

and Wang

, CXCL16/CXCR6 chemokine signaling mediates breast cancer progression by pERK1/2-dependent mechanisms, Oncotarget 6 (2015), 14165–78.

29.

Ehling

and Tacke

, Role of chemokine pathways in hepatobiliary cancer, Cancer Lett 379 (2016), 173–83.