Identification of CXCL13 as an Immune-Related Biomarker Associated with Tumorigenesis and Prognosis in Cutaneous Melanoma Patients.

BACKGROUND Melanoma is one of the most lethal tumors and its treatment is still challenging. It is urgent to detect novel therapy targets in melanoma. MATERIAL AND METHODS The GEO dataset was used to obtain a list of DEGS (differentially-expressed genes). Integrative bioinformatics analyses, including HPRD database, TCGA data, and TIMER, were performed to determine the role of CXCL13 in SKCM (skin cutaneous melanoma) progression and the immune environment. Furthermore, Pearson correlation coefficient analysis was used to measure correlations between CXCL13 and its co-expressed genes. Survival analysis, GO, and KEGG enrichment analysis were performed to investigate the role of CXCL13 in SKCM. RESULTS A total of 41 DEGs were identified in 3 GEO datasets, and 4 out of 41 DEGs are hub genes. Among the 4 hub genes, CXCL13 is involved in the most KEGG terms. CXCL13 is co-expressed with well-known immune checkpoint blockade targets, and it was associated with better overall survival. In addition, CXCL13 levels in infiltrating immune cells (neutrophil and myeloid dendritic cells) affect prognosis and survival in SKCM. Functional enrichment analysis clarified that CXCL13-co-expressed top 30 genes were associated with immune signaling pathways. Network analysis identified CXCL13 as a hub gene that interacts with CXCR5 to participate in immune-related biological process. CONCLUSIONS This study found that CXCL13 is associated with SKCM tumorigenesis and prognosis and immune infiltrations. Our result suggests that CXCL13 has great potential in development of novel immunotherapy targets in melanoma.

Background

Melanoma is a malignant tumor that causes approximately three-quarters of skin cancer-related deaths [1]. Melanoma treatment has been revolutionized with the approval of BRAFV600E inhibitors, which can significantly improve response and overall survival for patients with the BRAFV600E mutation [2]. However, patients develop acquired resistance to BRAFV600E inhibitors treatment [3]. Subsequent development of targeted immune checkpoints therapy has prolonged the survival of patients, such as with targeted CTLA-4 (ipilimumab, approved in 2011), and PD-1 (nivolumab, pembrolizumab, approved in 2014) has extended recurrence-free survival in high-risk resected melanoma patients [4,5]. However, melanoma recurrence was approximately 25–30% within 1 year in patients who were targeted by PD1, and PD1 monotherapy was not effective in those patients [5]. To overcome the drug resistance, it is urgent to explore novel therapeutic targets and treatment strategies for melanoma.

In our study, based on 3 datasets from the GEO database, we found that 41 differentially-expressed genes (DEGs) were overexpressed in melanoma compared to normal samples. Through network analysis, 4 (CXCL13, MMP1, SPP1, GZMB) out of the 41 DEGs had a high degree. KEGG enrichment analysis for the 41 DEGs indicated that CXCL1 is involved in 3 pathways. Moreover, Kaplan-Meier survival analysis showed that high CXCL13 expression is associated with longer survival time for tumor patients. Correction analysis further demonstrated that CXCL13 is co-expressed with 3 well-known immune checkpoint blockade targets (PD1(PDCD1), PDL1(CD274), and CTLA4) in SKCM. Functional enrichment analysis showed that CXCL13 and its highly co-expressed genes are involved in some important biological processes, such as immune pathways. Immune analysis showed that CXCL13 expression levels in infiltration of immune cells (CD4 T cell, neutrophil, macrophage, and myeloid dendritic cell) were associated with SKCM survival. These immune infiltration cells contribute to the tumor microenvironment, which can directly or indirectly modulate tumor immunity and have anti-tumor effects [6]. Thus, our results suggest CXCL13 as a new immune target and strategy for SKCM diagnosis and treatment.

Material and Methods

Identification of DEGs

We accessed GEO database to obtain the expression data. GSE15605 (58 melanoma samples and 16 normal samples), GSE46517 (104 melanoma samples and 12 normal samples), and GSE114445 (16 melanoma samples and 18 normal samples) were adopted to perform the differentially expression. The DEGs between melanoma and normal samples were collected using GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r) online service. DEGS were defined as the threshold of |logFC| >1.5 and p value ≤0.01. Next, Venn website (https://bioinfogp.cnb.csic.es/tools/venny/index.html) was applied to detect the overlap DEGs among these 3 datasets.

Network Construction and Functional Gene Sets

The network is an undirected graph G (V, E) where V stands for the genes and edges (i, j) ∈ E are weighted by PCC. Protein interaction networks were sourced from HPRD (http://www.hprd.org) database [7], the version is release 9. Protein self-interactions were removed, resulting in 39240 interactions among 9616 proteins. Cytoscape can create network interaction maps based on input genes and gene interaction [8]. The HPRD was used to extract the protein interaction of DEGs. 41DEGs were mapped to the HPRD, and the interaction was plotted by Cytoscape software (version 3.8.2, https://cytoscape.org/). The degree of each DEG was defined as the number of each DEG connected directly to its neighbor in HPRD. Only the nodes with degree greater than 5 were regarded as high-degree nodes.

Survival Analysis

The gene expression RNAseq and the survival time of TCGA-SKCM were downloaded from UCSC Xena (https://xenabrowser.net/datapages/). Gene expression data were log2(x+1) transformed, and include 20 531 identifiers and 474 samples (473 melanoma samples and 1 normal samples). The survival data include 479 samples. There were 457 common samples overlapped with 473 melanoma samples and 479 samples. The samples with gene expression greater than the 50% quantile of the gene expression were regarded as high expression, and otherwise were regarded as low expression. For CXCL13 gene survival analysis, we first classified samples into CXCL13 high-expression and CXCL13 low-expression groups. Then, we used the overall survival and event information from survival data to evaluate the survival difference between these 2 groups. Overall survival (OS) was defined as the time from the date of SKCM diagnosis to the date of patient death (from any cause) or last contact (whichever occurred first). Survival and survminer packages were utilized to perform survival analysis in R (version 3.6, https://www.r-project.org/).

Survival curves were estimated using the Kaplan-Meier method and compared across groups with the log-rank test. The t test was used to evaluate the statical significance. A p value <0.05 was statistically significant.

Co-Expression Analysis

Gene and gene co-expression were analyzed statistically using Pearson’s correlation coefficient (PCC). The sample size of the gene expression was 473 melanoma samples. The scatter plots of 3 known immune genes (PD1, PDL1 and CTLA4) correlated with CXCL13 were plotted using ggplot2, ggpubr, and ggpmisc packages in R (version 3.6). The heat map of the top 30 highly co-expressed genes with CXCL13 was plotted using the pheatmap package in R (version 3.6, https://www.r-project.org/).

Enrichment Analysis

We then performed GO and KEGG enrichment analyses of the top 30 high-correction-associated genes using the “clusterProfiler” package in R (version 3.6, https://www.r-project.org/).

TIMER

TIMER (http://timer.cistrome.org) is a comprehensive resource for systematic analysis (immune association, cancer exploration, and immune estimation) of immune infiltrates across diverse cancer types [9]. The association between immune infiltrates (CD4 T cells, neutrophils, macrophages, and myeloid dendritic cells) and clinical outcome for CXCL13 was explored using the TIMER online service.

Results

Identification of Differential Genes

To find the potential and reliable target genes, we first extracted genes which were differentially expressed in melanoma and normal samples. Through differential expression analysis, each of these datasets (GSE15065, GSE46517, and GSE114445) produced the corresponding list of DEGs. The 41 overlapped genes (Supplementary Table 1), as shown in the Venn diagram, are common DEGs in the 3 datasets (Figure 1A)

Figure 1

Venn diagram and the most significant KEGG pathways of 41DEGs. (A) The Venn diagram displays the number of DEGs in 3 datasets from the GEO database. (B) Four significant KEGG pathways of 41 DEGs. Adjusted p value <0.01. Figure 1A was produced using the Venny website (https://bioinfogp.cnb.csic.es/tools/venny/index.html, version 2.1). Figure 1B was produced using clusterProfiler packages in R (version 3.6, https://www.r-project.org/).

KEGG Enrichment Analysis of the DEGs

KEGG enrichment of 41 DEGs was performed using the clusterProfiler package in R software, then displayed in bubble charts (P<0.05, Figure 1B). The size of the bubble represents the number of genes in the KEGG terms, and the color of bubbles represents the adjusted P value in the KEGG terms. KEGG enrichment analysis indicated that changes in biologic pathways were significantly enriched in cytokine-cytokine receptor interaction, chemokine signaling pathway, and Toll-like receptor signaling pathway. These results revealed that immune response and inflammatory response play important roles in melanoma tumorigenesis.

PPI Network Construction and Hub Genes Identification

We mapped 41 DEGs to the HPRD database and extracted their directly connected proteins, including 113 nodes (gene) and 92 edges (gene and gene interaction). Then, we input the information of 92 edges into Cytoscape. In Figure 2A, we display the high-degree genes-associated network. Figure 2A shows 4 (MMP1, SPP1, CXCL13, and GZMB) hub genes (the node with yellow color) with a degree greater than 5. Combined with the KEGG enrichment information of 41 DEGs, and among the 4 hub genes, CXC13 was found to be involved in 3 KEGG pathway terms (Supplementary Table 2).

Figure 2

The network and survival plot. (A) The PPI (Protein-Protein Interaction) network of 4 high-degree DEGs (differentially-expressed genes). (B) Kaplan-Meier survival curves show that the expression of CXCL13 is consistently associated with better overall survival (OS) in SKCM. Figure 2A was produced using Cytoscape (version 3.8.2, https://cytoscape.org/). Figure 2B was produced using survival and survminer packages in R (version 3.6, https://www.r-project.org/).

Figure 2A shows that CXCL13 is connected with CXCR5. To investigate whether CXCL13 is co-expressed with CXCR5, we calculated the Pearson correlation coefficient. There was a positive correlation between CXCL13 expression and CXCR5 expression (PCC=0.74). High CXCL13 and CXCR5 expression were significantly related to longer survival time in SKCM patients (P=0.001). CXCL13 can interact with its ligand CXCR5 and play vital roles in apoptosis, proliferation, and differentiation of immune cells. It has been revealed that CXCL13 makes pivotal contributions to the biological process of multiple cancers via CXCL13/CXCR5 signaling [10]. For SKCM, in a previous study, Huang et al found that CXCL13 acts as an immune-related biomarker [11], which is in agreement with our results. The prognostic value and biological function of CXCL13 in SKCM are unclear and need to be better characterized. Here, we screened out available datasets associated with SKCM from public databases to systematically investigate the mechanism of CXCL13 on the progression, prognosis, and microenvironment of SKCM.

CXCL13 Level in Patients with SKCM

Early detection and treatment of SKCM is still a challenge that perplexes clinicians. Therefore, new therapy targets are needed to improve treatment efficacy. We first detected and compared the mRNA level of CXCL13 between SKCM samples and normal samples. Expression information was based on TCGA database, showing that the mRNA expression of CXCL13 was significantly lower in SKCM (the mean of CXCL13=−1.99) than in normal samples (the mean of CXCL13=6.78). We suggest that low CXCL13 expression occurs in SKCM and deserves further clinical validation as a potential therapy marker.

The Relationship Between CXCL13 Expression and Clinical Outcomes of SKCM Patients

We also evaluated the prognostic value of CXCL13 in SKCM from TCGA. Kaplan-Meier survival curves of OS were produced utilizing CXCL13 expression data and survival information of SKCM. The sample with CXCL13 expression greater than the 50th percentile of CXCL13 expression was treated as high CXCL13 expression, and otherwise was treated as low CXCL13 expression. The survival analysis showed that low expression of CXCL13 (P =0.00037) predicted poor OS (Figure 2B). In the above analysis, CXCL13 had lower expression in SKCM than in normal samples. This supports that CXCL13 low expression is related to poor OS. Therefore, it is conceivable that high CXCL13 expression is an independent risk factor and leads to a better prognosis in melanoma patients.

CXCL13 Co-Expression Genes Are Involved in Immune Pathways

We then focused on the function of the top 30 highly co-expressed genes with CXCL13 using GO enrichment analysis, showing that co-expression genes were correlated with CXCL13 in SKCM. PCC was utilized to analyze CXCL13 co-expressed genes from TCGA data. As shown in Figure 3A, the heat map displays the genes that showed significant positive correlations with CXCL13. The X axis represents SKCM samples, and the Y axis represents the rank of genes which are positive correlations with CXCL13. To explore the function of the top 30 highly co-expression genes, GO enrichment analyses were performed. LCK is a kinase that is among the top 30 highly co-expressed genes. Kinases play a crucial role in cancer formation and progress by modulating cancer cell migration, progression, and apoptosis [12]. The results of our study indicated that LCK is the vital regulator of CXCL13 in SKCM.

Figure 3

CXCL13 associated top 30 co-expressed genes and functional enrichment analysis in SKCM. (A) PCC was used to calculate correlations between CXCL13 and expressed genes in SKCM. Heat map showing top 30 genes positively correlated with CXCL13 in SKCM. (B) The significantly enriched GO annotations of the top 30 high CXCL13 co-expression genes in SKCM were analyzed. Figure 3A was produced using pheatmap packages in R (version 3.6, https://www.r-project.org/) and Figure 3B was produced using clusterProfiler packages in R (version 3.6, https://www.r-project.org/).

Figure 3B displays the top 10 most statistically significant GO terms. Biological processes (BP) were significantly enriched in T cell activation, lymphocyte differentiation, regulation of lymphocyte activation, regulation of T cell action, and regulation of cell-cell adhesion, which suggests CXCL13 has an effect on the transcriptome. Therefore, we found that the biological functions of these co-expressed genes are primarily related to immune-related pathways. Previous studies have demonstrated that the immune microenvironment plays key roles in immune evasion and metastasis of various cancers [13]. CXCL13 interacts with other genes to regulate the tumor immune environment, and CXCL13 co-expression genes make a crucial contribution to immune-related pathways. In summary, CXCL13 makes an important contribution to immune biological process.

CXCL13 is Co-Expressed with 3 Well-Known Immune Checkpoint Blockade Targets

To assess the potential of CXCL13 as an immune target, the correction of CXCL13 with 3 well-known checkpoint genes was measured. The most well-known checkpoint blockade used as cancer immunotherapies are anti-CTLA4 and -PD-1 and -PD-L1 antibodies.

To explore the biological effects of CXCL13 in SKCM, we examined CXCL13 co-expression in SKCM. CXCL13 is co-expressed with CTLA4, PD-1 (PDCD1), and PD-L1 (CD274) in SKCM and the relationship is statically significant (Figure 4). Immune checkpoint therapies, including PD-1 or PD-L1, and CTLA4, have extended patient survival times in multiple cancer types. In Figure 4A, X represents the CXCL13 expression data across 473 SKCM samples and the Y represents the PDCD1 expression data across 473 SKCM samples. We found that CXCL13 is co-expressed with these immune checkpoint blockade genes. CXCL13 expression showed a strong positive association with expression of PDCD1 (PCC=0.87, P<0.01), PDL1 (PCC=0.66, P<0.01), and CTLA4 (PCC=0.57, P<0.01). This suggests CXCL13 interacts with immune checkpoint blockade genes and plays a vital in immune function. The positive correlation between the expression of these proteins and CXCL13 expression indicates that CXCL13 could be a predictor for active response to anti-PD-1/PD-L1 immunotherapy in SKCM since PD-L1 expression is a biomarker for predicting the anti-PD-1/PD-L1 immunotherapy response in cancer. The cooperative anticancer effects of anti-CD137 with anti-CTLA4 or/and -PD-1 antibodies in several cancer types exceed that of a single antibody, which shows the advantage of antibody combinations in the immune checkpoint therapies in cancer treatment. PD-L1 expression can act as a indictor for the anti-PD-1/PD-L1 immunotherapy response in cancer [14]. The positive correlation between the expression of PD-1/PD-L1/CTLA4 and CXCL13 expression indicates that CXCL13 upregulation could be a biomarker for the response to anti-PD-1/PD-L1/CTLA4 immunotherapy in SKCM.

Figure 4

Gene expression correlation analysis. The scatter plot shows Pearson correlation of CXCL13 expression with expression of PDCD1 (A), CD274 (B), and CTLA4 (C). These figures were produced using ggplot2, ggpubr, and ggpmisc packages in R (version 3.6, https://www.r-project.org/).

To investigate whether the CXCL13 and the immune checkpoint gene have an effect on survival time, we performed survival analysis. We used the 50th quantile as a threshold to create high-expression and low-expression groups. Figure 5A shows the correlation between CXCL13 expression in SKCM expression data, CTLA4 expression in SKCM expression data, and survival time. Patients with high CXCL13 gene expression and high PD-1 gene expression had a longer survival time in SKCM (p=0.0011). In Figure 5B, patients with high CXCL13 gene expression and high PD-L1 gene expression had a longer survival time in SKCM (P<0.0001). In Figure 5C, patients with high CXCL13 gene expression and high CTLA4 gene expression had a longer survival time in SKCM (P<0.0001). This suggests that the efficacy of the combination of CXCL13 and PD-1/PD-L1/CTLA4 antibodies to treat SKCM. Thus, more investigations are urgently needed on the combination of CXCL13 and other antibodies to advance treatment effects and determine the underlying mechanisms.

Figure 5

Kaplan-Meier survival curves show that CXCL13 expression and PDCD1/CD274/CTLA4 expression are significantly associated with OS in TCGA-SKCM cohorts (A) PDCD1, (B) CD274, (C) CTLA4. These figures were produced using survival and survminer packages in R (version 3.6, https://www.r-project.org/).

Modules Exploring the Association Between Immune Infiltrates and CXCL13 Expression or Clinical Uutcome in TCGA

CXCL13 mediates the migration and localization of immune cells. Increasing information shows that immune cell infiltration can accelerate tumor progression and recurrence and affect immunotherapy and clinical outcome. Correlation between CXCL13 in SKCM expression, abundance of immune infiltrates (neutrophil cells, macrophages, myeloid dendritic cells and CD4+ T cells), and survival time are shown in Figure 6. Patients with high CXCL13 gene expression and high neutrophil infiltrates had a longer survival time than patients with high gene expression and low neutrophil infiltrates (P<0.01). Patients with high CXCL13 gene expression and low macrophage infiltrates tended to have a longer survival time than those with high CXCL13 gene expression and high macrophage infiltrates (P<0.01). Patients with high CXCL13 gene expression and high myeloid dendritic cell infiltrates tended to a have longer survival times (P=0.003). Patients with low CXCL13 gene expression and low macrophage infiltrates tended to have a longer survival time than those with low CXCL13 gene expression and high macrophage infiltrates (P=0.02). CD4+ T cells can inhibit cancer growth via activating M1 macrophages. Patients with low CXCL13 gene expression and low T cell CD4+ infiltrates had a longer survival time than patients with low gene expression and high T cell CD4+ infiltrates (P=0.04). Thus, cumulative curve analysis showed that immune infiltrates were significantly associated with CXCL13 in SKCM, indicating that immune infiltrate cells significantly affect the prognosis. In this study, we found a significant correlation between the survival time and the expression of CXCL13 in the infiltration of the 4 immune cell types – neutrophils, CD4+ T cells, macrophages, and myeloid dendritic cells – indicating that CXCL13 expression in combination with immune cell status can predict prognosis. Further research is needed to explore the potential immune therapy using CXCL13.

Figure 6

Correlation of CXCL13 expression with immune infiltration level in the TIMER database. (A) Neutrophil cells, (B) macrophages, (C) myeloid dendritic cells and (D) CD4+ T cells. These figures were produced using the TIMMER website (http://timer.cistrome.org, version 2.0).

Discussion

Melanoma is a malignant tumor. Although PD1 monotherapy and BRAF inhibitors were widely used in patients with melanoma, these therapies were not active in those patients with drug resistance. Therefore, it is urgent to detect more therapeutic targets and prognostic biomarkers.

In our study, 41 common DEGs were found in 3 datasets. KEGG enrichment analysis indicated that these DEGs play a crucial role in immune response and inflammatory response. After inputting 41 DEGs and their PPI into Cytoscape software, we obtained 4 (MMP1, SPP1, CXCL13, and GZMB) hub genes with a high degree. KEGG enrichment revealed that CXCL13 was enriched in the most KEGG pathway terms. Recent studies have reported that chemokine receptors expressed on melanoma and immune cells are closely associated with the prognosis and the efficacy of melanoma immunotherapy [15,16]. CXCL13 is a chemokine with a pivotal role in the immune system [17]. Through literature searching, we know that CXCL13 is differentially expressed in the skin cutaneous melanoma (SKCM) [18]. Kazanietz et al suggest that how CXCL13 regulates the relationship between tumor cell and cancer microenvironment is vital to detect novel therapeutic targets for cancer treatment [16]. Jiao et al revealed the associations of CXCL13 and immune cell infiltration signature in renal cancer and the potential of CXCL13 as an immune response biomarker in patients with renal cancer [19]. In melanoma, no studies have systematically investigated the interaction between CXCL13 and the microenvironment. Herein, we performed comprehensive analyses to investigate the mechanisms of CXCL13 in the SKCM microenvironment and immunotherapy efficacy.

Our analysis suggests that CXCL13 is differentially expressed in TCGA-SKCM samples. Through network analysis, CXCL13 was shown to have a high degree and can interact with CXCR5. CXCL13/CXCR5 signaling has been reported in many cancers. High CXCL13 expression is an independent risk factor and leads to a better prognosis in TCGA-SKCM samples. We then focused on the function of the top 30 highly co-expressed genes with CXCL13 using GO enrichment analysis. Our results suggested that CXCL13 interacts with other genes to regulate the tumor immune environment, which makes a crucial contribution to immune-related pathways. Further, the positive correlation between the expression of PD-1/PD-L1/CTLA4 and CXCL13 expression indicates that CXCL13 could be a biomarker for the response to anti-PD-1/PD-L1/CTLA4 immunotherapy in SKCM.

In addition, the expressions of CXCL13 were positively correlated with infiltration of 4 immune cell (neutrophils, CD4+ T cells, macrophages, and myeloid dendritic cells). Previous research reported that the degree of immune cell infiltration was related to favorable outcomes [20]. This is similar to our results that increasing levels of neutrophils and myeloid dendritic cells have a relationship with longer survival time in TCGA-SKCM patients.

Our results provide insights for future research on immunotherapy of SKCM. CXCL13 and associated pathways involved in SKCM immunity could be a target to design more effective therapeutic strategies. We anticipate our study results can help design new immunotherapeutic drugs, and to guide clinicians in using appropriate drugs for patients.

It also could be highly useful in the identification of comorbidity patterns between diseases. Bioinformatics tools are novel alternatives for conventional approaches to get these data and integrate them in a biologically meaningful manner. Exploiting the bioinformatics tools, we have introduced CXCL13-related genes and pathways obtained from our network-based analysis. These genes and pathways could be subjected to further empirical investigations to attain additional insights about the biology and immunology of CXCL13 as well as its relationship to other diseases. Our study has some limitations. Analysis on the transcriptional level can reflect some aspects of immune status, but not global changes. Moreover, other independent cohort and in vitro or in vivo studies should be performed to validate our results.

Conclusions

In summary, we systematically analyzed CXCL13 expression, prognosis, biological function, immune infiltration, and its relationship with 3 well-known immune checkpoint blockade targets. CXCL13 expression was closely related to survival rates of SKCM patients. CXCL13 can interact with other genes to regulate the immune pathways. It also has a correction with the expression of PD-1/PD-L1/CTLA4. In addition, CXCL13 is associated with infiltration of multiple immune cells. Research shows that CXCL13 may serve as potential therapeutic target in cutaneous melanoma, and it could be a biomarker for the response to anti-PD-1/PD-L1/CTLA4 immunotherapy in SKCM.

Supplementary Table 1

The 41 common DEGs in GSE15605, GSE46517, and GSE114445.

41 differential expressed genes
LAMB4, CHP2, EXPH5,CDHR1, SCGB1D2, CCL27, GATA3, PRAME, CYP3A5, RORA, POU2F3, SCEL, PDZD2, DSC3, SPP1, KLF5, WIF1, GZMB, PHACTR1, TFAP2B, GDF15, RGS1, IGF2BP3, SCGB2A2, CCL3L3, C1QB, MAGEA6, CYP39A1, IL37, MAGEA6, CITED1, MAGEA12, FCGR2A, SLAMF7, MMP1, CXCL13, MAGEA2B, CXCL9, CCL5, UBD, LCE2B

Supplementary Table 2

The KEGG enrichment of 41 DEGs.

ID	Description	p value	p. adjust	Gene ID	Count
hsa04061	Viral protein interaction with cytokine and cytokine receptor	7.23E-08	4.55E-06	10850/414062/27178/10563/4283/6352	6
hsa04060	Cytokine-cytokine receptor interaction	2.73E-06	8.58E-05	10850/9518/414062/27178/10563/4283/6352	7
hsa04062	Chemokine signaling pathway	6.28E-05	0.00132	10850/414062/10563/4283/6352	5
hsa04620	Toll-like receptor signaling pathway	8.55E-05	0.001347	6696/414062/4283/6352	4