The data were collected from the published literature (Table 1). For different sequencing methods of single-cell data, specific analysis procedures were applied. For Drop-seq single-cell data, Cell Ranger software (Freytag et al., 2018) (3.0.1) was adopted to calculate the cell expression counts. For Smart-seq2 single-cell sequencing data, we operated cell expression matrixes provided in the original article. The expression matrix file was then imported into R 3.6.2 for subsequent analysis.

First, we used Seurat (Stuart et al., 2019) (3.1.4) to filter the quality of cells and delete all cells with more than 6000 expressed genes or less than 201 genes. A total of 39,057 UGIC cells and 215,291 other cancer cells were obtained. Next, standardized integration processing was performed on the cell level, sample level, and study level.

After PCA dimensionality reduction was performed on 39,057 UGIC cells, nine cell sets were obtained by T-distributed stochastic neighbor embedding (t-SNE) clustering. In order to identify the cell types, we calculated highly expressed genes on each cluster through the FindMarkers (Butler et al., 2018) function in Seurat. Then, through the artificial gene annotation on the CellMarker (RRID:SCR_018503) database (Zhang X. et al., 2019), the marker genes and the corresponding cell type were finally annotated. We show the statistical graph of cell types identified by EPCAM in the published articles as an example in the CellMarker (RRID:SCR_018503) database (Supplementary Figure S5). Then, we analyzed the subtypes of cancer stem cells, obtained a total of six subclasses, and calculated the differentially expressed genes (DEGs) of each subclass.

A total of 71 single-cell sequencing data (Supplementary Table S2) from six other cancers were collected. We used the same method to process other tumor single-cell data to ensure the consistency of the analysis process. First, the quality control of single-cell data obtained a total of 215,291 cells. After standardization at the cell level, sample level, and study level, we used PCA and t-SNE visualization to reduce the dimension of those single-cell data and obtained 29 cell collections. We calculated the highly expressed genes of 29 cell collections and used the CellMarker database (Zhang X. et al., 2019) to annotate the cell types. Then, we marked cancer stem cells, which are subtypes 4 and 7. The relevant marker annotations are shown in Supplementary Figure S3.

We gathered bulk RNA-seq data of UGICs in the TCGA database (Aldape et al., 2015). We obtained the expression matrix data using the cBioPortal (Cerami et al., 2012; Gao et al., 2013), including 522 HNSCC samples, 185 EC samples, and 415 GC samples. Three types of UGICs were congregated with the data label “hnsc_tcga,” “esca_tcga,” and “stad_tcga”. DESeq2 (RRID:SCR_000154) (Love et al., 2014) (1.26.0) software was used to measure the DEGs in the cancer sample and the corresponding normal sample.

We annotated the function and pathway information of the significantly different genes in the Gene Ontology (GO) (Ashburner et al., 2000; Ashburner, 2021) database and Kyoto Encyclopedia of Genes and Genomes (KEGG) (RRID:SCR_012773) (Kanehisa et al., 2021) database using the clusterProfiler (RRID:SCR_016884) (Yu et al., 2012) (3.14.3) package in R (3.6.2) software. The top 15 terms are presented in Figure 4.

We adjusted the gene set enrichment score between the specific differential genes and the cancer-related gene sets through the Gene Set Enrichment Analysis (GSEA) (Subramanian et al., 2005) software (4.1.0). “C6: oncogenic signatures” was selected as the existing cancer-related gene set in GSEA software. We filtered several parameters to draw gene enrichment results. The normalized enrichment score (NES) was larger than 1. The normalized significance level (NSL) p-value was lower than 0.05.

We collected all human entries in the String database (RRID:SCR_005223) and deleted low-quality and text mining entries (Szklarczyk et al., 2021). After removing the duplicated edges and self-loops, we constructed a human protein–protein interaction network (PPIN) with 19,267 proteins and 1,689,887 edges by Cytoscape software (RRID:SCR_003032) (Shannon et al., 2003) (3.7.1). Then, 174 genes specifically expressed in UGCSC were mapped to the PPIN. After removing outlier proteins, a regulatory sub-network composed of 144 protein nodes and 545 edges was constructed. Next, we appraised the topological attributes of the network and selected the degree and clustering coefficient (CC) to measure the function of the sub-network (Sporns, 2013). The degree represents the number of connections through a particular node, which measures the importance of the node in the network. CC represents the closeness of connections between a node and the surrounding nodes, which demonstrates the network closeness and function similarity. The formula (2) is described as follows: C i = 2 e i d i ( d i − 1 ) , (2) where C i represents the CC of gene i. d i represents the count of adjacent nodes of gene i. e i represents the number of interconnected nodes among all adjacent nodes of gene i.

We manually reviewed the tumor-related literature studies published since 2000 to screen functions and pathways of DEGs. Then, we formulated UGCSC function networks by integrating different genes and the known inflammation and Wnt pathways (Figure 6).

We collected 39,057 tumor single-cell sequencing data from 30 patients including 4762 HNSCC cells, 366 EC cells, and 33,927 GC cells (Figure 1A). It is noteworthy that EC samples applied Smart-seq2 single-cell sequencing technology (Picelli et al., 2014) which is manually sequencing each cell. So the EC group has few cells but higher confidence. After quality filtering (see Methods) and removing the batch effect, more than 70 million transcripts were obtained from 39,057 cells. Subsequently, we classified cells into different clusters by using T-distributed stochastic neighbor embedding (t-SNE) methods in Seurat software (Supplementary Figure S1). Through marker genes, these identified cell clusters could be assigned to known cell lineages: T cells, B cells, epithelial cells, natural killer cells, fibroblasts, plasma cells, cancer stem cells, mast cells, and endothelial cells (Figure 1B). To corroborate these profiles, we showed the high expression gene distribution heatmap of each cell type and the expression abundance of marker genes of each type (Figure 1C; Supplementary Table S1). Each cell type has specific marker genes: CD3D, KRT8, MS4A1, PDGFRA, IGHG3, EPCAM, ECSCR, TPSB2, and CCL3 (Figure 1D). The violin plot of marker genes shows that the expression of most marker genes is specific, which indicates that the classification of cell types is accurate and is very helpful for subsequent analysis (Figure 1E). Taken together, these results indicate that the cell classification was accurate, and most of the cells were classified into the correct cell type. The distribution of samples and cancer types is shown in Supplementary Figures S1, S2. We also counted the number and frequency of all cell types in HNSCC, EC, and GC and provided the results in Supplementary Figure S6.

We focused on cancer stem cell types in order to reveal the pathogenesis and distant metastasis mechanism of UGIC. We collected a total of 1,586 CSCs (Figures 2A,B) including 136 HNSCC cells, 23 EC cells, and 1427 GC cells. Due to the heterogeneity of CSCs, there are differences in the same type of cancer while similarities exist in different types of cancers, coincident with the characteristics of the remote metastasis and recurrence of the cancers. Therefore, we performed a cluster analysis of CSCs, and a total of six sub-clusters were found. After annotating and analyzing all sub-clusters, sub-cluster 0 is ubiquitous in UGICs, including 19 EC stem cells, 356 GC stem cells, and 114 HNSCC stem cells, which proves that sub-cluster 0 preliminarily meets the characteristics of common CSCs (Figures 2C,D). Therefore, we concentrated on sub-cluster 0 in the follow-up analysis.

To verify whether sub-cluster 0 reflects the characteristics of UGIC rather than only GC, we performed a down-sampling process in sub-cluster 0 since there are more than 70% GC stem cells in sub-cluster 0. We randomly selected the same number of GC cells as the HNSCC cells and named new sub-cluster 0. Subsequently, we compared the differential genes between sub-cluster 0 and the new sub-cluster 0 in CSCs. The merge ratio is 77.14% (Figure 2E), which means these two sub-clusters share the same differential gene set. These results indicate that the differentially expressed genes (DEGs) of sub-cluster 0 represent the features of UGICs.

To validate the specificity of UGCSCs, we compared UGCSCs with other tumor cells. We collected 71 samples (215,291 cells) from six types of cancers including glioma (GLM), melanoma (MELA), osteosarcoma (OSTC), breast cancer (BC), ovarian cancer (OVC), and stellate cell cancer (SCC) (Supplementary Table S2). After normalizing the cells and removing the batch effect (see Method), all the cells were gathered into 29 sub-clusters (Figures 3A,C). After annotating all cancer cells in the CellMarker database, we noticed that there are plenty of cell types due to the complexity of tissue types involved. Therefore, we only annotated CSCs by using marker genes. The tumor stem cells were obviously aggregated with CXCR4 markers (Figure 3B; Supplementary Figure S3), which are sub-clusters 4 and 7 and contain 21323 cells, as circled in Figure 3A. We compared CSCs of other cancers with CSCs of UGICs. We re-clustered and obtained 31 sub-clusters in all CSCs, which reveals the differences between CSCs of different tumor types (Figure 3D). But at the same time, the cluster distribution of CSCs from different tissues is uniform, which indicates that there are similarities between different tissues in CSCs (Figures 3E,F). This phenomenon is also coincident with the heterogeneity of tumors. The cluster annotation of cancer types shows that the UGCSC is self-clustering and far away from other tumor CSCs (Figure 3E). Therefore, the UGCSC is the specific cancer stem cell in UGIC while UGCSC does not exist in other cancers.

We comprehensively analyzed the distribution and function of UGCSCs. The cell sources of UGCSC cancers were analyzed and counted (Figure 4A). As shown in Figure 4, UGCSCs are averagely expressed in UGIC patients, including 10 GC patients, four EC patients, and nine HNSCC patients. In summary, the UGCSCs are distributed uniformly, which proves that the UGCSC is common in upper gastrointestinal patients.

We analyzed the expression network of UGCSCs. First, we compared the expression profiles of UGCSCs and all other tumor stem cells and obtained 174 genes with significant differences, including 33 upregulated genes and 141 downregulated genes (Supplementary Data S1). We uncovered that the gene information function reflects the characteristics of UGCSCs as a digestive system and as cancer stem cells by analyzing the function annotation (Ashburner et al., 2000) of DEGs (Figures 4B,D). The upregulated genes are related to antibacterial response, such as “antibacterial humoral response,” “antibacterial humoral immune response-mediated,” and “mucosal immune response,” which are consistent with the function of the digestive tract in the human body. In addition, the functions of downregulated genes mainly focus on reducing the activity of T cells and lymphocytes and downregulating cell killing which could reduce the body’s immune response and enhance the survival rate of tumor cells, and these are also the characteristics of cancer stem cells. The downregulated genes also play a role in the regulation of cell–cell adhesion to facilitate the distant metastasis of tumors, which is in line with the feature of metastasis. Through the analysis of the KEGG pathway (Kanehisa et al., 2021), we observed that the significantly differentially expressed genes are enriched in inflammation-related mucosal infections such as “Staphylococcus aureus infection,” “epithelial cell signaling in Helicobacter pylori infection,” “IL-17 signaling pathway,” and “chemokine signal pathway” (Figure 4C). These results uncovered a potential carcinogenic factor of UGIC, that is, mucosal damage induced activation and mutation in inflammatory pathways. Through gene set enrichment analysis (GSEA), we found that the significantly differentially expressed genes were significantly enriched in ATF2- and ERBB2-related cancer genes (Figures 4E,F). To further confirm the reliability of these genes, we calculated the DEGs between UGCSCs and 30,267 cells of normal tissues in the upper alimentary tract (Cillo et al., 2020; Zhang X. et al., 2021) (Supplementary Figure S7A). The DEGs of UGCSCs and tumor cells are few and share fewer genes with the DEGs of UGCSCs and normal cells (Supplementary Figure S7B). Also, the functional analysis of DEGs of UGCSCs and normal cells indicates that the functional pathways are related to cell development, which is a common feature of tumors (Supplementary Figure S7C). In summary, these results indicate that the DEGs of UGCSCs and tumor cells are oncogenes related to the function of the digestive tract.

In order to study the pathogenic mechanisms that may exist in UGCSCs, we mapped 174 proteins into the human protein–protein interaction network (PPIN). We constructed PPIN and deleted low-quality text mining terms in the String database (Szklarczyk et al., 2021). After mapping 174 DEGs in UGCSCs into PPIN, an interaction network consisting of 144 proteins and 545 edges was obtained (Figure 5A). Through the analysis of the topological properties of the network, we found that the degree of DEGs in PPIN is 362.174, which is significantly higher than 175.418 in the human total network (Figure 5B). This result indicates that the shortest path through different genes is significantly higher than the average value (Figure 5B), which implied that these genes are hub genes in the UGCSC network. Furthermore, another topological property, the clustering coefficient (CC), is significantly higher than the background network, which points out that the 144 genes are closely linked compared with the random gene set in the network. The close interaction means a similar or synergistic function in cells. Through the comprehensive analysis of degree and CC, we inferred that the 144 genes are tightly connected hub genes in PPIN, which means that they play an important function in UGCSCs as a co-operative hub gene set.

We have performed functional annotations on the possible functions of these genes and inferred regulatory pathways with the aim to explore the possible pathogenic mechanisms and potential therapeutic targets in UGIC. We analyzed the regulation pathway of those genes through published articles and proved that the upregulated genes are basically related to cancer (Supplementary Table S3). Here are some exciting discoveries. Some genes are related to inflammatory pathways, such as CXCL8 (Ha et al., 2017), BPIFB1 (Li J. et al., 2020), PIGR (Kakiuchi et al., 2020), CXCL3, and RNASE1 (Wang et al., 2006), and some genes are related to specific functions of the digestive tract, such as GAST (Giraud et al., 2016), REG1A (Sha et al., 2019), and TFF3 (Braga Emidio et al., 2020). These results illustrated that there may have similar pathogenic mechanisms and common regulatory pathways in some UGICs. We speculated that mucosal damage is induced by long-term unhealthy eating habits, which include smoking, drinking, and hot food breed inflammation. Persistent inflammation leads to carcinogenic mutations and early gastrointestinal tumors. These conjectures have been confirmed in the specific regulatory network of UGCSCs. Based on the detected differentially expressed genes and the mining of relevant research literature studies, we speculated the pathogenesis of the disease, as shown in the Figure 6. Inflammation-associated interleukin (CXCL8 and CXCL3) and inflammation defense-related BPIFB1, PIGR, and RNASE1 are activated in UGCSCs. Combined with the epidemiological investigation of gastrointestinal cancer, there is a hypothesis that chronic inflammation is incited by mucosal damage due to long-term bad eating habits. We present that the cancerous chronic inflammation is activated by GAST, REG1A, TFF3, and ZG16B in the Wnt signaling pathway. Upregulated hPG80 and TFF3 induce PI3K/Ras and lead to tumor cell growth and invasion, which may be one reason for the poor prognosis of UGIC. Hence, this resource provides a novel view for the occurrence and development of UGIC and the advancement of gastrointestinal cancer diagnosis and therapy.