- 1Department of General Surgery, Shijiazhuang People's Hospital, Shijiazhuang, China
- 2Department of Medical Oncology, The Fourth Hospital of Hebei Medical University, Hebei Tumor Hospital, Shijiazhuang, China
- 3Translational Medicine, YuceBio Technology Co., Ltd., Shenzhen, China
- 4Internal Medical, University of Occupational and Environmental Health, Fukuoka, Japan
Pancreatic cancer is one of the most challenging cancer types in clinical treatment worldwide. This study aimed to understand the tumorigenesis mechanism and explore potential therapeutic targets for patients with pancreatic cancer. Single-cell data and expression profiles of pancreatic cancer samples and normal tissues from multiple databases were included. Comprehensive bioinformatics analyses were applied to clarify tumor microenvironment and identify key genes involved in cancer development. Immense difference of cell types was shown between tumor and normal samples. Four cell types (B cell_1, B cell_2, cancer cell_3, and CD1C+_B dendritic cell_3) were screened to be significantly associated with prognosis. Three ligand–receptor pairs, including CD74-MIF, CD74-COPA, and CD74-APP, greatly contributed to tumorigenesis. High expression of BUB1 (BUB1 Mitotic Checkpoint Serine/Threonine Kinase) was closely correlated with worse prognosis. CD1C+_B dendritic cell_3 played a key role in tumorigenesis and cancer progression possibly through CD74-MIF. BUB1 can serve as a prognostic biomarker and a therapeutic target for patients with pancreatic cancer. The study provided a novel insight into studying the molecular mechanism of pancreatic cancer development and proposed a potential strategy for exploiting new drugs.
Introduction
In 2020, 495,773 new cases of pancreatic cancer were diagnosed all over the world, showing an increase of more than 1/3 cases compared with the data in 2015 (1, 2). High-income countries, especially Europe and Northern America, have a high incidence of pancreatic cancer, which was partially attributed to the improved diagnosis and an aging population (1). Smoking is the most studied risk factor and has been demonstrated to be highly associated with pancreatic cancer (3). Evidence showed that smokers have two- to three-fold risk much higher than non-smokers, and an obvious dosage–risk relationship is observed (3). Other risk factors, such as dietary factors (4), low physical activity (5), and obesity (6), are also associated with pancreatic cancer.
Pancreatic cancer is one of the most lethal cancers with high mortality, ranking the seventh in cancer-related mortality. Patients with pancreatic cancer were already in an advanced stage when diagnosed, which increases difficulty in treatment. In addition, pancreatic cancer is not sensitive to conventional therapies such as chemotherapy, radiotherapy, or molecular-targeted therapy. Surgery is still the main strategy for non-metastatic or local advanced patients. The 5-year survival rate is about 15–25% of patients accepting surgery compared with <7% of an averaged 5-year survival in all individuals (7). However, not more than 20% of patients with pancreatic cancer could be treated by surgery, whereas over 60% of surgery-managed patients will undergo relapse in the first year (8). Therefore, effective biomarkers for early screening pancreatic cancer or molecular targets for personalized therapy are urgently needed.
CA19-9 is a well-known biomarker of pancreatic cancer with a sensitivity of about 80% for diagnosis. It can also serve as a monitor to evaluate the response systematic treatment in the neoadjuvant or surgery (9). However, CA19-9 level may be affected by other diseases such as biliary obstruction (10). Other biomarkers related to transforming growth factor-beta (TGF-β), angiogenesis, inflammation, and immune response have been explored (11). For example, a meta-analysis demonstrated that loss of SMAD4 (SMAD Family Member 4) expression was associated with worse prognosis of patients with pancreatic cancer (12). Although a number of biomarkers for predicting prognosis have been developed, the mechanisms and treatments of pancreatic cancer should be improved. Molecular-targeted therapy is promising strategy for managing pancreatic cancer.
The BUB1 gene is located on human chromosome 2q14. Its encoded protein is a platform protein for spindle physical examination. It is the basis for other cell components to be accurately located in the spindle. BUB1 gene plays an indispensable role in maintaining correct chromosome separation and reducing aneuploid formation during mitosis. Increasing studies have shown that BUB1 plays an important role in the occurrence and development of tumors. Qi et al. (13) identified the expression of BUB1 in liver hepatocellular carcinoma (LIHC) tissue by immunohistochemistry. BUB1 expression is accompanied by immune cell infiltration into LIHC tissue. Yun et al. (14) found that BUB1 is a key gene in colorectal cancer. Gao et al. (15) observed that BUB1 can be used as a prognostic marker of gastric cancer. Alam et al. (16) found that BUB1 may be used as a biomarker in breast cancer. Jiang et al. (17) reported that BUB1 mediates STAT3 signaling pathway to drive the occurrence and development of bladder cancer. These studies have shown that BUB1 is an important gene in cell cycle control and DNA damage repair, but it is rarely reported in pancreatic cancer. The development of single-cell sequencing technology allows an in-depth understanding of the occurrence and development mechanism of pancreatic cancer, the heterogeneity of tumor microenvironment, and the formation mechanism of drug resistance, so as to find new therapeutic targets. Earlier, Zhao et al. (18) identified the metabolic reprogramming of human colon immune cells in different locations and disease states by using single-cell RNA sequencing (scRNA-seq). Lai et al. (19) constructed a prognostic model for predicting the survival rate of human glioblastoma by comprehensive analysis of scRNA-seq dataset and RNA-seq dataset. Wang et al. (20) identified the heterogeneity of CD8+ T cells and novel biomarker genes in hepatocellular carcinoma by single-cell sequencing. These studies have shown that it is effective to use single-cell technology to find tumor heterogeneity and key biomarkers. In the recent years, single-cell technology has been greatly improved and applied in studying mechanisms from different aspects in pancreatic cancer (21–23). Single-cell data provide an expanded space for exploring tumor microenvironment, heterogeneity, and other aspects in cancer. Therefore, in this study, we applied the scRNA-seq data of pancreatic cancer to distinguish differential cell types between normal and tumor samples. Workflow is as shown in Figure 1A. Together with other independent data, we discovered hub genes highly associated with prognosis. The study provides a new insight into understanding the mechanisms of pancreatic cancer development.
Figure 1. Work flowchart and clustering and dimensionality reduction of one normal sample and two tumor samples by UMAP. (A) Work flowchart. (B) The distribution of three samples labeled by three colors after dimensionality reduction. (C) The distribution of one normal sample and tumor samples. (D) 20 subgroups from C0 to C19 of three samples. (E) The proportion of different subgroups in each sample. IPMN, PASC, and the normal sample located from inside to outside, respectively. IPMN, intraductal papillary mucinous neoplasm; PASC, pancreatic adenosquamous carcinoma.
Materials and Methods
Data Source
The scRNA-seq data (GSE165399) (24) of an intraductal papillary mucinous neoplasm (IPMN) sample, a pancreatic adenosquamous carcinoma (PASC) sample, and a normal sample were downloaded from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The Cancer Genome Atlas (TCGA) dataset, including 177 tumor samples, was downloaded from TCGA database (https://portal.gdc.cancer.gov/). Other normal samples of pancreatic tissue were downloaded from UCSC Xena database (http://xena.ucsc.edu/). Other pancreatic adenocarcinoma (PAAD) samples, including GSE21501 (25), GSE28735 (26), GSE57495 (27), GSE62452 (28), and GSE85916, were downloaded from GEO database. The clinical characteristics of each dataset are shown in Table 1.
Data Preprocessing
For scRNA-seq data, quality control was performed, where each gene expressed at least three cells and each cell expressed at least 250 genes. “PercentageFeatureSet” in SEURAT R package (29) was conducted to calculate the proportion of mitochondria and rRNA. Mitochondria fraction <30% and each cell expressing at least 500 genes were set to screen data. Samples from TCGA and UCSC Xena were combined and grouped by normal and tumor samples. “NormalizeBetweenArrays” in limma R package (30) was performed to normalize the data. Finally, 177 tumor samples and 167 normal samples remained with a total of 24,210 genes, and they were defined as TCGA-PAAD dataset. For GSE21501, GSE28735, GSE57495, GSE62452, and GSE85916 datasets, samples without clinical information were excluded, and only tumor samples remained. Probes were converted to gene symbol. “RemoveBatchEffect” function was conducted to remove batch effects and “NormalizeBetweenArrays” function was used to normalize the data.
Dimensionality Reduction by UMAP
Three single-cell samples were merged and normalized by log-normalization. “FindVariableFeatures” function was conducted to identify highly variable genes. “ScaleData” function was used to scale data and principal component analysis (PCA) was applied to reduce dimensionality. “FindNeighbors” and “FindClusters” functions were performed to cluster cells under the conditions of dim = 40 and resolution = 0.5. The top 40 components were selected, and uniform manifold approximation and projection (UMAP) was used to further reduce dimensionality. UMAP is a nonlinear dimensionality-reduction technique for visualizing single-cell data and process high-dimensional data (31).
Definition of Cell Types
Specific marker genes of different tissues were obtained from CellMarker database (http://biocc.hrbmu.edu.cn/CellMarker/). Markers of pancreas, pancreatic acinar tissue, fetal pancreas, peripheral blood, and blood were selected. By using enricher function in clusterProfiler R package (32), the enrichment score of these marker genes of each subgroup was calculated. According to the high or low enrichment of these marker genes in different subgroups, different cell types were defined. One cell type having multiple subgroups was classified into different cell types. “Minimum.spanning.tree” function was used to calculate the minimum distance between two cell types, and a minimum-cost spanning tree (MST) was constructed. MST and UMAP plots were used to estimate whether two subgroups could be combined.
Analysis of Cell Trajectory and Identification of key Regulator Genes
Cell trajectory reflected the development time and differentiation degree of different cell types. Monocle 2 toolkit was applied to generate the trajectory of cell development (33). Monocle is an unsupervised algorithm for processing high-dimensionality data and dynamically presents single-cell data. Branched expression analysis modeling (BEAM) measurement in monocle 2 was employed to identify regulator genes (33). BEAM is a regression model for detecting critical genes involved in cell development.
ReactomeGSA for Functional Analysis
When analyzing functional pathways of subgroups, ReactomeGSA R package (34) was implemented to calculate the enrichment score of each pathway in each subgroup. The top 20 differentially enriched pathways were selected. ReactomeGSA linking to reactome database can allow to analyze functional pathways from multi-omics and derive novel biomedical insights (34).
Ligand–Receptor Interactions Analyzed by CellPhoneDB
CellPhoneDB, a public tool with curated receptors and ligands, was applied to assess cell–cell communication (35). Databases, including UniProt, Ensembl, PDB, the IMEx consortium, and IUPHAR, were utilized by CellPhoneDB for comprehensively assessing ligand–receptor interactions. CellPhoneDB is convenient to screen important interactions between different cell types for single-cell data. Mean expression of each ligand–receptor and their p-values were calculated. Ligand–receptor network was used to manifest the overall interactions among cell types. Thick and thin interactions between two cell types represent strong and weak interactions between them, respectively. Dot plots presenting specific ligand–receptor pairs were used to display interactions between specific cell types under a condition of mean interaction >1.
Constructing Protein–Protein Interaction Network
Univariate Cox regression analysis was performed to screen marker genes of C4, C10, C16, and C18 associated with prognosis (p < 0.001). For these screened genes, rcorr function in Hmisc R package (https://cran.r-project.org/web/packages/Hmisc/index.html) was applied to analyze the correlation among each genes. Correlation coefficient| > 0.9 and p < 0.001 were selected to screen gene pairs. Based on the screened marker genes associated with prognosis, STRING database (https://string-db.org/) and Cytoscape software (version 3.6.1) (36) were introduced to construct protein–protein interaction (PPI) network for identifying hub genes involved in cancer development.
Statistical Analysis
All statistical analyses were performed in R platform (version 3.4.2). Student's t-test was conducted to test the significance of expression between two groups. ANOVA was performed to test the distribution of cell types between normal and tumor samples. Log-rank test was conducted in the Kaplan–Meier survival analysis and univariate Cox regression analysis. p < 0.05 was considered as significant. All parameters of tools not specifically shown were default.
Results
Clustering Single-Cell Data of Pancreatic Cancer Based on Dimensionality Reduction
Single-cell data (GSE165399), including three samples (one normal samples and two tumor samples), was preprocessed by SEURAT R package. Parameters of one gene expressing at least in three cells and one cell expressing at least 250 genes were set to screen gene expression data. “PercentageFeatureSet” function was conducted to calculate the proportion of mitochondria and rRNA in cells. One cell expressing at least 500 genes with not more than 30% percentile mitochondria was set to ensure the quality. The quality control and total cell counts before and after filtration were displayed (Supplementary Figures S1–S3). Then, data of three samples were merged and normalized by log-normalization. High-variable genes were identified by “FindVariableFeatures” function, and the top 20 high-variable genes were shown (Supplementary Figure S4).
The PCA was performed to construct a two-dimensional distribution of gene expression (Supplementary Figure S5). Then, function of “FindNeighbors” and “FindClusters” was performed to cluster cells, and 20 subgroups were generated. The dimensionality of expression data was further reduced by UMAP function based on the top 40 components (Figures 1B,C). In total, 20 subgroups from C0 to C19 were identified, and their proportions in each sample were shown (Figures 1D,E). We observed that three samples had significantly different distribution of 20 subgroups even within two tumor samples, indicating that there was a difference in cancer development within difference cancer types. Using “FindAllMarkers” function, we screened the marker genes of each subgroup significantly differential expressed among the subgroups (p < 0.05), and only the top five marker genes were exhibited (Figure 2).
Figure 2. The top 5 markers of 20 subgroups. Yellow represents relatively high expression and plum represents relatively low expression.
Defining Cell Types for 20 Subgroups
To accurately classify 20 subgroups into different cell types, we downloaded cell markers of specific tissues, including pancreas, pancreatic acinar tissue, fetal pancreas, peripheral blood, and blood, from CellMarker database. By implementing “enricher” function in clusterProfiler R package, 20 subgroups were defined into 12 different cell types, and 6 out of 12 cell types could be further classified into different cell subtypes (Table 2).
To further classify the 6 cell types, MST was used to calculate the minimum distance between the two subgroups (Figure 3A). C4 and C10 could not be combined, as the closest distance to C4 was C19. C4 and C10 were defined as B cell_1 and B cell_2, respectively. Cancer cells had three subgroups (C9, C15, and C18), where C18 was the closest to C17, and C9 was close to both C12 and C15. Therefore, three subgroups of cancer cells were defined into three cell types. CD1C-CD141-dendritic cells and CD1C+_B dendritic cells (DCs) had multiple subgroups, but their distance was not close. Fibroblasts (C2 and C11) were close to C13, and natural killer cells (C5 and C14) were close to C4. The multiple branches of C4 and C13 indicated that they may develop into different subgroups. The development of one cell type was affected by the cancer microenvironment; therefore, we did not merge different subgroups in both fibroblasts and natural killer cells.
Figure 3. Confirmation of 20 subgroups. (A) MST of 20 subgroups. Yellow points labeled with different numbers represents different subgroups. (B) Violin plots showing marker genes expressed in different subgroups.
Analysis of marker genes in these cell types demonstrated that different cell types expressed specific markers, supporting that these subgroups should be classified into different cell types (Figure 3B). For example, CD1C+_B dendritic cell_1 specifically expressed LAMA3 and S100A16, CD1C+_B dendritic cell_2 specifically expressed IL1R2 and LGALS2, and CD1C+_B dendritic cell_3 specifically expressed S100B and CXCL10 (Figure 3B).
Cell Development and Key Regulatory Genes of 20 Subgroups
To delineate the cell development of different cell types, we applied monocle to predict their cell trajectory. Pseudotime was used to evaluate the degree of their cell division, and 20 subgroups were divided into three major branches (Figure 4). An obvious difference of cell distribution was observed between normal and tumor samples. Tumor samples had a significantly higher proportion of cells locating in the early pseudotime than the normal sample, suggesting that tumor samples had a large number of undifferentiated or incompletely differentiated cells, which probably led to high potential of cell proliferation. Three branches were defined as three different status (states 1, 2, and 3).
Figure 4. Cell trajectory of 20 subgroups grouped by normal and tumoral. (A) Pseudotime of different subgroups. Colors from red to blue indicates pseudotime from early to late. (B) The location of 20 subgroups in the trajectory. (C) Three states with different colors defined by three branches.
By using branched expression analysis modeling, we then screened a series of regulatory genes of three branches. The result showed a total of 127 key regulatory genes highly associated with cell development (corrected p < 0.0001). Gene expression of 127 regulatory genes were shown in a heatmap grouped by three status (Supplementary Figure S6). We found that the expression of 127 genes significantly varied by the pseudotime, indicating that these 127 genes may serve as key regulators in the development of different cell types. There was also an obvious difference in the expression pattern of 127 genes among three branches, with cluster 1 mostly expressed in state 2, cluster 2 mostly expressed in states 1 and 3, and cluster 3 mostly expressed in state 3 (Supplementary Figure S7). The distribution of different 20 cell types in the three branches is shown in Table 3. State 1 included the majority of CD1C-CD141- DCs, CD1C+_B dendritic cell_2, myeloid DCs, CD1C+_B dendritic cell_3, and beta cells. Plasmacytoid DCs, B cells, and natural killer cells were mostly accumulated in state 2. Fibroblasts, CD1C+_B dendritic cell_1, cancer cells, endothelial cells, and AXL+SIGLEC6+ DCs were mostly enriched in state 3. The results suggested that these regulatory genes may play different roles in the different status of cancer cell development.
ReactomeGSA for Identifying Differential Enriched Pathways
To evaluate functional pathways of 20 cell types, we applied ReactomeGSA linking to reactome database to analyze the enrichment score of each pathway in each cell type. The top 20 differentially enriched pathways of 20 cell types were visualized in a heatmap (Figure 5). We observed that different cell types in one cell type manifested different enrichment patterns of these pathways. For instance, CD1C-CD141- dendritic cell_2 had two more enriched pathways (COX reactions and CYP2E1 reactions) than CD1C-CD141- dendritic cell_1. Three subgroups of cancer cells exhibited obviously different enrichment in six pathways, for example, beta Klotho-mediated ligand binding was enriched in cancer cell_1, alanine metabolism and degradation of GABA were enriched in cancer cell_2, and COX reactions were enriched in cancer cell_3. The specific enrichment score of all 20 pathways in each subgroup was also displayed (Supplementary Figure S8).
Figure 5. A heatmap of functional pathways differentially enriched in 20 subgroups. Enrichment score from high to low was labeled as colors from red to blue. The right line indicates 20 subgroups and the bottom line indicates 20 functional pathways. ES, enrichment score.
Differential Distribution of Subgroups Between Normal and Tumor Samples
Next, we used CIBERSORT measurement to analyze the proportion of 20 subgroups in normal cells and tumor cells in an independent dataset (TCGA-PAAD) according to the screened marker genes in single-cell data. The results showed that 16 of 20 subgroups were differentially enriched between normal cells and tumor cells (p < 0.01, Figure 6). B cells, cancer cell_2, cancer cell_3, and natural killer cell_1 were more enriched in normal cells, whereas beta cells, cancer cell_1, DCs, endothelial cells, and fibroblasts were more enriched in tumor cells.
Figure 6. The fraction of 20 cell subgroups in TCGA-PAAD dataset. Red indicates tumor samples and green indicates normal samples.
In the independent comparison of the proportion of 16 subgroups between normal and tumor cells, we found that except for fibroblast_1, 15 out of 16 subgroups were significantly distributed between them (p < 0.05, Figure 7). Interestingly, normal cells had limited proportion of beta cells, cancer cell_1, CD1C-CD141- dendritic cell_1, CD1C+_B dendritic cell_1, endothelial cells, fibroblast_2, and natural killer cell_2. The proportion of these cell types were all markedly elevated in tumor samples. Furthermore, we analyzed the relation between 16 subgroups and survival in tumor samples, and observed that only 4 subgroups were associated with prognosis (p < 0.05, Supplementary Figure S9). B cell_1, B cell_2, and cancer cell_3 were positively associated with prognosis, whereas low enrichment of CD1C+_B dendritic cell_3 had more favorable prognosis. We considered that these four cell types possibly played critical roles in cancer development.
Figure 7. The independent proportion of 16 subgroups in normal and tumor samples. High and low indicates the enrichment of cell types. *p < 0.05.
Ligand–Receptor Interactions Among Different Subgroups
As B cell_1 (C4), B cell_2 (C10), cancer cell_3 (C18), and CD1C+_B dendritic cell_3 (C16) were identified to be significantly associated with prognosis, we therefore analyzed their interactions among them and with other subgroups based on CellPhoneDB. The ligand–receptor interactions among 20 cell types were displayed in a network (Figure 8A). The interactions of four cell types (C4, C10, C18, and C16) with other types were independently displayed (Figures 8B–E). Within these four cell types, CD1C+_B dendritic cell_3 (C16) had the most interactions with other types, especially CD1C-CD141- dendritic cell_1 (C0), CD1C+_B dendritic cell_1 (C6), and myeloid DC (C8) (Figure 8E). Both B cell_1 (C4) and B cell_2 (C10) had close interactions with CD1C-CD141- dendritic cell_2 (C1) and CD1C+_B dendritic cell_3 (C16) (Figures 8B,C). Cancer cell_3 (C18) strongly interacted with CD1C+_B dendritic cell_1 (C6) and CD1C+_B dendritic cell_3 (C16) (Figure 8D).
Figure 8. Ligand–receptor interactions among 20 subgroups analyzed by CellPhoneDB. (A) A network of ligand–receptor interactions among 20 subgroups. (B–E) Ligand–receptor interactions of B cell_1 (C4), B cell_2 (C10), cancer cell_3 (C18), and CD1C+_B dendritic cell_3 (C16) to other subgroups, respectively. Circles with different colors indicate different subgroups. Thick lines and fine lines between subgroups indicate strong and weak interactions between them.
In addition, we screened the specific ligand–receptor interactions of four cell types to others under the condition of averaged interaction value >1 (Figure 9). Three ligand–receptor pairs, including CD74_MIF, CD74_COPA, and CD74_APP, were the most active among these interactions. Overall, these four cell types had close interactions with DCs, indicating that different types of DCs probably played important roles in constructing microenvironment beneficial to tumor growth.
Figure 9. Dot plots of specific ligand–receptor interactions of C4 (A), C10 (B), C18 (C), and C16 (D) to other subgroups.
Screening the hub Gene BUB1 of Pancreatic Cancer Through PPI Networks
To identify the hub genes responsible for cancer development, we focused on four cell types (C4, C10, C16, and C18) and conducted univariate Cox regression analysis for their marker genes. A total of 73 genes were strongly correlated with prognosis, with 2 marker genes in C4, 15 genes in C10, 54 genes in C16, and 2 genes in C18. Through conducting rcorr function in Hmisc package, the correlation and significance among these marker genes were calculated. In total, 14 pairs of marker genes were screened under |correlation coefficient > 0.9 and p < 0.001 (Supplementary Table S1). The PPI network revealed that BUB1 gene was in the center, and it was the bridge to link CD1C+_B DCs and B cells (Figure 10A). Survival analysis of the relation between BUB1 expression and survival showed a significant association, with low BUB1 expression showing more favorable prognosis (p = 0.00066, Figure 10B). Moreover, we determined BUB1 expression in two tumor samples and one normal sample. It almost had no expression in the normal sample but was obviously expressed in tumor samples, specially expressed in natural killer cell_2 and CD1C+_B dendritic cell_3 (Figures 10C–F). Integrating the data of normal samples and tumor tissues of pancreas from GTEX and TCGA showed that BUB1 in tumor samples was significantly higher than that in normal samples (Figure 10G). In addition, in the transcriptome samples, we also analyzed the relationship between BUB1 expression and immunity. By evaluating the immune cell infiltration score of patients with ciberport and mcpcounter, we can observe that the monocytes, CD8 T cells, endothelial cells, cytotoxic lymphocytes, and T-cell infiltration of patients with low expression of BUB1 were significantly higher than that of patients with high expression of BUB1 (Figures 11A,B). Similarly, estimate was used to evaluate the immune microenvironment infiltration of patients, it can be observed that the immune microenvironment infiltration in patients with low BUB1 expression was significantly higher than that in patients with high BUB1 expression (Figure 11C). Furthermore, the enrichment score of each patient in Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway was evaluated by single-sample gene set enrichment analysis. We calculated the correlation between BUB1 expression and pathway, selected KEGG pathway with significant correlation and correlation coefficient >0.4, and identified a total of 33 key pathways, which included cell_ Cycle, DNA_REPLICATION and MISMATCH_REPAIR and other important pathways related to cell proliferation, such as P53_SIGNALING_PATHWAY, SMALL_CELL_LUNG_CANCER and BLADDER_CANCER, and other pathways closely related to tumors (Figure 11D). These results showed that the abnormality of BUB1 was closely related to the occurrence and development of tumors. To validate the above results, we introduced other expression data, including GSE21501, GSE28735, GSE57495, GSE62452, and GSE85916, to evaluate BUB1 expression and its relation to prognosis. Five datasets were preprocessed and normalized for removing batch effects. The expression profiles of five datasets were combined without bias, and their distribution before and after preprocessing was shown in PCA plots (Supplementary Figures S10A,B). The Kaplan–Meier survival analysis demonstrated that low BUB1 expression was still significantly correlated with overall survival (p = 0.023, Supplementary Figure S10C). In addition, using SangerBox (http://vip.sangerbox.com) to analyze the relationship between BUB1 expression and clinical features, we observed that although there was no significant expression difference of BUB1 in different clinical stages, but there were significant expression differences in patients with different depth of invasion (Supplementary Figure S10D). We also evaluated the relationship between the expression of BUB1 in Pan-cancer and prognosis. It can be observed that BUB1 had a significant prognostic correlation in many other tumors, such as renal cancer, glioma, liver cancer, leukemia, lung cancer, and so on (Supplementary Figure S10E), indicating that BUB1 acted as an essential role in pancreatic cancer.
Figure 10. Screening of hub genes highly contributing to cancer development. (A) A PPI network of marker genes of C4, C10, C16 and C18. Red triangle indicates cell types and blue circle indicates marker genes. (B) Kaplan–Meier survival plot of BUB1 grouped by high and low expression. (C–E) UMAP plots showing distribution and expression of BUB1 in IPMN, PASC and normal samples. Colors from gray to blue indicate expression from low to high. (F) Violin plot of BUB1 expression in 20 subtypes. (G) Difference in expression and distribution of BUB1 between cancer and normal samples.
Figure 11. BUB1 expression and immune and functional analysis. (A) Distribution difference of 22 kinds of immune cell infiltration in patients with high and low expression of BUB1. (B) Distribution difference of immune cell infiltration in 10 kinds of patients with high and low expression of BUB1. (C) Differences in the distribution of immune infiltration in patients with high and low expression of BUB1. (D) BUB1 expression was significantly related to KEGG pathway.
Discussion
The application of single-cell sequencing technology improves the past understanding of tumor microenvironment (TME), pathological process, molecular characteristics in different subtypes of pancreatic cancer (21, 23). In the current study, we characterized the distribution of different cell types between normal samples and tumor samples, and aimed to identify a molecular target highly associated with pancreatic cancer development for potentially used in drug exploration.
We included three samples (IPMN, PASC, and normal samples) with scRNA-seq data and identified 20 subgroups (C0 to C19). Obvious differences of subgroup distribution among three samples were observed, indicating the possible distinct mechanisms of cancer development between IPMN and PASC and great alternation of TME. In the normal pancreatic tissue, myeloid DCs (C8), cancer cell_1 (C9), fibroblast_2 (C11), and endothelial cells (C12) contributed to an over 3/4 proportion of all cells. On the contrary, these cell types were less enriched or almost disappeared in two tumor samples. In IPMN, CD1C-CD141- dendritic cell_1 (C0), fibroblast_1 (C2), plasmacytoid DCs (C3), B cell_1 (C4), and natural killer cell_1 (C5) were the most enriched. In PASC, two cell types of CD1C-CD141- dendritic cell_2 (C1) and CD1C+_B dendritic cell_1 (C6) consisted of about 3/4 in all cell types. Notably, (DCs comprise a majority number in all cell types, and it largely varied in these three samples, suggesting that DCs underwent immense transformation or abnormal division during tumorigenesis.
A line of studies have demonstrated that DCs play key roles in activating immune response in tumors, such as capturing antigens, processing, and presenting them as antigenic peptides to T cells (37). However, its function can be impaired by tumors through suppressing DC accumulation, activation, and antigen presentation (38). With the exploring study on DCs, DCs have been considered as promising vaccines for integrating tumor peptides and presenting antigens against tumors (37). In the present study, different types of DCs in TME may play different roles in tumor development.
The distinct distribution of cell types between normal and tumor samples was also validated in the independent dataset. Importantly, we discovered that four cell types (B cell_1, B cell_2, cancer cell_3, and CD1C+_B dendritic cell_3) were significantly associated with prognosis. We assumed that these four cell types may be essential in tumor progression. Therefore, ligand–receptor interactions among different cell types were analyzed. Within the above four cell types, CD1C+_B dendritic cell_3 was found to have strong cell–cell communications with other cell types. Especially, three ligand–receptor pairs, including CD74-MIF, CD74-COPA, and CD74-APP, were much more enriched in the interactions of CD1C+_B dendritic cell_3 to other cell types.
CD74-MIF has been reported to regulate immune activity (39). CD74 is an essential receptor to DCs for their migration and mediating immune response and regulates the development of T and B cells (40). Macrophage migration inhibitory factor (MIF) is a key cytokine involving in inflammatory diseases such as pulmonary inflammation through CD74 signaling (41). Tanese et al. (42) found that CD74-MIF serves as mediators to activate PI3K/AKT signaling pathway in melanoma. Moreover, a study illustrated that blockading CD74-MIF on macrophages and DCs can recover antitumor activity of immune system in metastatic melanoma (43). Apart from melanoma, CD74 and MIF have been considered as therapeutic targets in other cancer types such as prostate cancer and gastric cancer (44, 45). In pancreatic cancer, CD74-MIF is possibly a promising target for molecular therapy, but further experimental study is needed. CD74-COPA and CD74-APP were seldom reported, so their high enrichment also indicated the important function in tumor progression in pancreatic cancer, but this needs further verification.
As four cell types were identified to be associated with prognosis, we considered that there were key marker genes in them regulating tumorigenesis in pancreatic cancer. Therefore, the marker genes of B cell_1, B cell_2, cancer cell_3, and CD1C+_B dendritic cell_3 were screened by regression analysis. Through correlation analysis and constructing the PPI network, we discovered that BUB1 gene displayed as a bridge to communicate CD1C+_B dendritic cell_3 and B cells.
BUB1 is a multitask protein kinase for chromosome segregation in eukaryotes. BUB1 impairment or dysregulation leads to chromosomal instability and results in tumorigenesis (46, 47). BUB1 has been widely reported to be associated with tumorigenesis in various cancer types including gastric cancer (48), breast cancer (49), and pancreatic ductal adenocarcinoma (50). Piao et al. (50) proved that high expression of BUB1 was associated with worse overall survival of pancreatic ductal adenocarcinoma, which was consistent with our observation. BUB1 could serve as a biomarker for predicting prognosis of patients with pancreatic cancer, as it also presented robust performance in the independent dataset. Importantly, we discovered that BUB1 was specifically expressed in natural killer cell_2 and CD1C+_B dendritic cell_3, suggesting that BUB1 probably played a critical role in modulating the expression of MIF that had close relation to tumorigenesis. To date, no study has reported this novel possible mechanism in cancer development. Therefore, BUB1 had great potential to be a therapeutic target for molecular therapy, but this needs more evidence in cell or animal experiments.
Immunotherapy is recognized as a promising strategy in many cancer types. According to our observations, DCs play a critical role against tumor cells, and they were impaired in TME, probably resulted from the dysregulation of cytokines or chemokines in TME. Currently, surgery is the only treatment with curative possibility for patients with pancreatic cancer. Considering the insensitivity of pancreatic cancer to chemotherapy and radiotherapy, other molecular-targeted drugs are particularly needed. In this study, CD74-MIF and BUB1 were found to be potential therapeutic targets for treating pancreatic cancer.
In conclusion, this study identified four cell types that are significantly associated with pancreatic cancer prognosis based on single-cell data and integrated bioinformatics analysis. Cell types largely varied in different cancer types and between normal and tumor samples. DCs, especially CD1C+_B dendritic cell_3, played a critical role in cancer progression probably through CD74 signaling. Three pairs of ligand–receptor interactions (CD74-MIF, CD74-COPA, and CD74-APP) were considered to be closely involved in tumorigenesis. Importantly, BUB1 could serve as a biomarker for predicting prognosis and a therapeutic target for treating pancreatic cancer, but this needs further experiments.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Author Contributions
DW and JS designed the study, reviewed, and edited the manuscript. MY contributed to the literature research. ZZ, XC, and JS analyzed and interpreted the data. ML, XD, and YX wrote the initial draft of the manuscript. All the authors read and approved the manuscript and its publication. All authors contributed to the completion of this study.
Conflict of Interest
YX, ZZ, XC, and DW were employed by YuceBio Technology Co., Ltd.
The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher's Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpubh.2022.900853/full#supplementary-material
References
1. Kleeff J, Korc M, Apte M, La Vecchia C, Johnson CD, Biankin AV, et al. Pancreatic cancer. Nat Rev Dis Primers. (2016) 2:16022. doi: 10.1038/nrdp.2016.22
2. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin. (2021) 71:209–49. doi: 10.3322/caac.21660
3. Bosetti C, Lucenteforte E, Silverman DT, Petersen G, Bracci PM, Ji BT, et al. Cigarette smoking and pancreatic cancer: an analysis from the International Pancreatic Cancer Case-Control Consortium (Panc4). Ann Oncol. (2012) 23:1880–8. doi: 10.1093/annonc/mds491
4. Larsson SC, Wolk A. Red and processed meat consumption and risk of pancreatic cancer: meta-analysis of prospective studies. Br J Cancer. (2012) 106:603–7. doi: 10.1038/bjc.2011.585
5. Behrens G, Jochem C, Schmid D, Keimling M, Ricci C, Leitzmann MF. Physical activity and risk of pancreatic cancer: a systematic review and meta-analysis. Eur J Epidemiol. (2015) 30:279–98. doi: 10.1007/s10654-015-0014-9
6. Genkinger JM, Kitahara CM, Bernstein L, Berrington de Gonzalez A, Brotzman M, Elena JW, et al. Central adiposity, obesity during early adulthood, and pancreatic cancer mortality in a pooled analysis of cohort studies. Ann Oncol. (2015) 26, 2257–66. doi: 10.1093/annonc/mdv355
7. He J, Ahuja N, Makary MA, Cameron JL, Eckhauser FE, Choti MA, et al. 2564 resected periampullary adenocarcinomas at a single institution: trends over three decades. HPB. (2014) 16:83–90. doi: 10.1111/hpb.12078
8. Varadhachary GR, Tamm EP, Abbruzzese JL, Xiong HQ, Crane CH, Wang H, et al. Borderline resectable pancreatic cancer: definitions, management, and role of preoperative therapy. Ann Surg Oncol. (2006) 13:1035–46. doi: 10.1245/ASO.2006.08.011
9. Bauer TM, El-Rayes BF, Li X, Hammad N, Philip PA, Shields AF, et al. Carbohydrate antigen 19-9 is a prognostic and predictive biomarker in patients with advanced pancreatic cancer who receive gemcitabine-containing chemotherapy: a pooled analysis of 6 prospective trials. Cancer. (2013) 119:285–92. doi: 10.1002/cncr.27734
10. Mann DV, Edwards R, Ho S, Lau WY, Glazer G. Elevated tumour marker CA19-9: clinical interpretation and influence of obstructive jaundice. Eur J Surg Oncol. (2000) 26:474–9. doi: 10.1053/ejso.1999.0925
11. Hasan S, Jacob R, Manne U, Paluri R. Advances in pancreatic cancer biomarkers. Oncol Rev. (2019) 13:410. doi: 10.4081/oncol.2019.410
12. Shugang X, Hongfa Y, Jianpeng L, Xu Z, Jingqi F, Xiangxiang L, et al. Prognostic Value of SMAD4 in Pancreatic Cancer: A Meta-Analysis. Transl Oncol. (2016) 9:1–7. doi: 10.1016/j.tranon.2015.11.007
13. Qi W, Bai Y, Wang Y, Liu L, Zhang Y, Yu Y, et al. BUB1 Predicts Poor Prognosis and Immune Status in Liver Hepatocellular Carcinoma. APMIS. (2022). doi: 10.1111/apm.13219. [Epub ahead of print].
14. Yun ZJ, Wang HJ, Yu YX, Sun ZY, Yao SK. [Screening of differentially expressed genes for colorectal cancer and prediction of potential traditional Chinese medicine: based on bioinformatics]. Zhongguo Zhong Yao Za Zhi. (2022) 47:1666–76. doi: 10.19540/j.cnki.cjcmm.20211108.402
15. Gao F, Wang J, Li C, Xie C, Su M, Zou C, et al. Risk-related genes and associated signaling pathways of gastrointestinal stromal tumors. Int J Gen Med. (2022) 15:3839–49. doi: 10.2147/IJGM.S357224
16. Alam MS, Rahaman MM, Sultana A, Wang G, Mollah MNH. Statistics and network-based approaches to identify molecular mechanisms that drive the progression of breast cancer. Comput Biol Med. (2022) 145:105508. doi: 10.1016/j.compbiomed.2022.105508
17. Jiang N, Liao Y, Wang M, Wang Y, Wang K, Guo J, et al. BUB1 drives the occurrence and development of bladder cancer by mediating the STAT3 signaling pathway. J Exp Clin Cancer Res. (2021) 40:378. doi: 10.1186/s13046-021-02179-z
18. Zhao Q, Zhang T, Yang H. ScRNA-seq identified the metabolic reprogramming of human colonic immune cells in different locations and disease states. Biochem Biophys Res Commun. (2022) 604:96–103. doi: 10.1016/j.bbrc.2022.03.034
19. Lai W, Li D, Kuang J, Deng L, Lu Q. Integrated analysis of single-cell RNA-seq dataset and bulk RNA-seq dataset constructs a prognostic model for predicting survival in human glioblastoma. Brain Behav. (2022) 12:e2575. doi: 10.1002/brb3.2575
20. Wang H, Fu Y, Da BB, Xiong G. Single-Cell Sequencing Identifies the Heterogeneity of CD8+ T Cells and Novel Biomarker Genes in Hepatocellular Carcinoma. J Healthc Eng. (2022) 2022:8256314. doi: 10.1155/2022/8256314
21. Bernard V, Semaan A, Huang J, San Lucas FA, Mulu FC, Stephens BM, et al. Single-cell transcriptomics of pancreatic cancer precursors demonstrates epithelial and microenvironmental heterogeneity as an early event in neoplastic progression. Clin Cancer Res. (2019) 25:2194–205. doi: 10.1158/1078-0432.CCR-18-1955
22. Elyada E, Bolisetty M, Laise P, Flynn WF, Courtois ET, Burkhart RA, et al. Cross-species single-cell analysis of pancreatic ductal adenocarcinoma reveals antigen-presenting cancer-associated fibroblasts. Cancer Discov. (2019) 9:1102–23. doi: 10.1158/2159-8290.CD-19-0094
23. Peng J, Sun BF, Chen CY, Zhou JY, Chen YS, Chen H, et al. Single-cell RNA-seq highlights intra-tumoral heterogeneity and malignant progression in pancreatic ductal adenocarcinoma. Cell Res. (2019) 29:725–38. doi: 10.1038/s41422-019-0195-y
24. Zhao X, Li H, Lyu S, Zhai J, Ji Z, Zhang Z, et al. Single-cell transcriptomics reveals heterogeneous progression and EGFR activation in pancreatic adenosquamous carcinoma. Int J Biol Sci. (2021) 17:2590–605. doi: 10.7150/ijbs.58886
25. Stratford JK, Yan F, Hill RA, Major MB, Graves LM, Der CJ, et al. Genetic and pharmacological inhibition of TTK impairs pancreatic cancer cell line growth by inducing lethal chromosomal instability. PLoS ONE. (2017) 12:e0174863. doi: 10.1371/journal.pone.0174863
26. Zhang G, He P, Tan H, Budhu A, Gaedcke J, Ghadimi BM, et al. Integration of metabolomics and transcriptomics revealed a fatty acid network exerting growth inhibitory effects in human pancreatic cancer. Clin Cancer Res. (2013) 19:4983–93. doi: 10.1158/1078-0432.CCR-13-0209
27. Chen DT, Davis-Yadley AH, Huang PY, Husain K, Centeno BA, Permuth-Wey J, et al. Prognostic fifteen-gene signature for early stage pancreatic ductal adenocarcinoma. PLoS ONE. (2015) 10:e0133562. doi: 10.1371/journal.pone.0133562
28. Yang S, He P, Wang J, Schetter A, Tang W, Funamizu N, et al. A Novel MIF Signaling Pathway Drives the Malignant Character of Pancreatic Cancer by Targeting NR3C2. Cancer Res. (2016) 76:3838–50. doi: 10.1158/0008-5472.CAN-15-2841
29. Gribov A, Sill M, Lück S, Rücker F, Döhner K, Bullinger L, et al. SEURAT: visual analytics for the integrated analysis of microarray data. BMC Med Genomics. (2010) 3:21. doi: 10.1186/1755-8794-3-21
30. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. (2015) 43:e47. doi: 10.1093/nar/gkv007
31. Becht E, McInnes L, Healy J, Dutertre CA, Kwok IWH, Ng LG, et al. Dimensionality reduction for visualizing single-cell data using UMAP. Nat Biotechnol. (2018). doi: 10.1038/nbt.4314
32. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics. (2012) 16:284–7. doi: 10.1089/omi.2011.0118
33. Qiu X, Mao Q, Tang Y, Wang L, Chawla R, Pliner HA, et al. Reversed graph embedding resolves complex single-cell trajectories. Nat Methods. (2017) 14:979–82. doi: 10.1038/nmeth.4402
34. Griss J, Viteri G, Sidiropoulos K, Nguyen V, Fabregat A, Hermjakob H. ReactomeGSA - Efficient Multi-Omics Comparative Pathway Analysis. Mol Cell Proteom. (2020) 19:2115–25. doi: 10.1074/mcp.TIR120.002155
35. Efremova M, Vento-Tormo M, Teichmann SA, Vento-Tormo R. CellPhoneDB: inferring cell-cell communication from combined expression of multi-subunit ligand-receptor complexes. Nat Protoc. (2020) 15:1484–506. doi: 10.1038/s41596-020-0292-x
36. Su G, Morris JH, Demchak B, Bader GD. Biological network exploration with Cytoscape 3. Curr Protocols Bioinform. (2014) 47:11-24. doi: 10.1002/0471250953.bi0813s47
37. Yang J, Shangguan J, Eresen A, Li Y, Wang J, Zhang Z. Dendritic cells in pancreatic cancer immunotherapy: Vaccines and combination immunotherapies. Pathol Res Pract. (2019) 215:152691. doi: 10.1016/j.prp.2019.152691
38. Wculek SK, Cueto FJ, Mujal AM, Melero I, Krummel MF, Sancho D. Dendritic cells in cancer immunology and immunotherapy. Nat Rev Immunol. (2020) 20:7–24. doi: 10.1038/s41577-019-0210-z
39. Leng L, Metz CN, Fang Y, Xu J, Donnelly S, Baugh J, et al. MIF signal transduction initiated by binding to CD74. J Exp Med. (2003) 197:1467–76. doi: 10.1084/jem.20030286
40. Su H, Na N, Zhang X, Zhao Y. The biological function and significance of CD74 in immune diseases. Inflamm Res. (2017) 66:209–16. doi: 10.1007/s00011-016-0995-1
41. Takahashi K, Koga K, Linge HM, Zhang Y, Lin X, Metz CN, et al. Macrophage CD74 contributes to MIF-induced pulmonary inflammation. Respir Res. (2009) 10:33. doi: 10.1186/1465-9921-10-33
42. Tanese K, Hashimoto Y, Berkova Z, Wang Y, Samaniego F, Lee JE, et al. Cell Surface CD74-MIF interactions drive melanoma survival in response to interferon-γ. J Invest Dermatol. (2015) 135:2775–84. doi: 10.1038/jid.2015.204
43. Figueiredo CR, Azevedo RA, Mousdell S, Resende-Lara PT, Ireland L, Santos A, et al. Blockade of MIF-CD74 Signalling on macrophages and dendritic cells restores the antitumour immune response against metastatic melanoma. Front Immunol. (2018) 9:1132. doi: 10.3389/fimmu.2018.01132
44. Meyer-Siegler KL, Iczkowski KA, Leng L, Bucala R, Vera PL. Inhibition of macrophage migration inhibitory factor or its receptor (CD74) attenuates growth and invasion of DU-145 prostate cancer cells. J Immunol (Baltimore, Md: 1950). (2006) 177:8730–9. doi: 10.4049/jimmunol.177.12.8730
45. Zheng YX, Yang M, Rong TT, Yuan XL, Ma YH, Wang ZH, et al. CD74 and macrophage migration inhibitory factor as therapeutic targets in gastric cancer. World J Gastroenterol. (2012) 18:2253–61. doi: 10.3748/wjg.v18.i18.2253
46. Baker DJ, Jin F, Jeganathan KB, van Deursen JM. Whole chromosome instability caused by Bub1 insufficiency drives tumorigenesis through tumor suppressor gene loss of heterozygosity. Cancer Cell. (2009) 16:475–86. doi: 10.1016/j.ccr.2009.10.023
47. Ricke RM, Jeganathan KB, van Deursen JM. Bub1 overexpression induces aneuploidy and tumor formation through Aurora B kinase hyperactivation. J Cell Biol. (2011) 193:1049–64. doi: 10.1083/jcb.201012035
48. Grabsch H, Takeno S, Parsons WJ, Pomjanski N, Boecking A, Gabbert HE, et al. Overexpression of the mitotic checkpoint genes BUB1, BUBR1, and BUB3 in gastric cancer–association with tumour cell proliferation. J Pathol. (2003) 200:16–22. doi: 10.1002/path.1324
49. Wang Z, Katsaros D, Shen Y, Fu Y, Canuto EM, Benedetto C, et al. Biological and Clinical Significance of MAD2L1 and BUB1, Genes Frequently Appearing in Expression Signatures for Breast Cancer Prognosis. PLoS ONE. (2015) 10:e0136246. doi: 10.1371/journal.pone.0136246
Keywords: pancreatic cancer, single-cell data, dendritic cells, ligand–receptor interactions, CD74, Bub1, therapeutic targets, bioinformatics analysis
Citation: Li M, Duan X, Xiao Y, Yuan M, Zhao Z, Cui X, Wu D and Shi J (2022) BUB1 Is Identified as a Potential Therapeutic Target for Pancreatic Cancer Treatment. Front. Public Health 10:900853. doi: 10.3389/fpubh.2022.900853
Received: 21 March 2022; Accepted: 03 May 2022;
Published: 13 June 2022.
Edited by:
Yuanpeng Zhang, Nantong University, ChinaReviewed by:
Shuxia Ma, Jiamusi University, ChinaXin Chen, Second Hospital of Hebei Medical University, China
Alvarado Anabell, University of Guadalajara, Mexico
Haitao Xu, Anqing Hospital Affiliated to Anhui Medical University, China
Copyright © 2022 Li, Duan, Xiao, Yuan, Zhao, Cui, Wu and Shi. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Dongfang Wu, d3VkZiYjeDAwMDQwO3l1Y2ViaW8uY29t; Jian Shi, amppYW5zaGkwMCYjeDAwMDQwOzE2My5jb20=
†These authors have contributed equally to this work