Analysis of cancer-associated fibroblasts in cervical cancer by single-cell RNA sequencing

Selection of research subjects

This research selected surgically resected tumor tissues from three patients with cervical cancer hospitalized at the First Affiliated Hospital of Zhengzhou University for single-cell sequencing. This study was approved by the Ethics Committee of the First Affiliated Hospital of Zhengzhou University (ethics number: 2018-KY-28), and informed consent was obtained from patients. Besides, single-cell sequencing data were attained from GSE168652 in the Gene Expression Omnibus (GEO) database for one case of cervical cancer paracancer tissue (GSM5155197) [8].

Preparation of fresh single-cell suspensions

After surgery, the excised tissues were placed in lyophilized tubes containing tissue protection solution pre-cooled at 4° C and washed with DPBS pre-cooled at 4° C and repeated 2-3 times to remove the residual tissue protection solution. Then, with a small amount of digestion solution added, the tissue was cut into small pieces and placed in centrifuge tubes (50 ml) containing 5 mL of digestion solution and digested by shaking in a constant temperature water bath at 37° C. After digestion, with the digest passed through a pre-wetted 40 μm cell sieve, the filtrate was collected into a new 50 mL centrifuge tube. 2 mL of pre-cooled erythrocyte lysate at 4° C was added to the cell sediment, gently blown and mixed, and incubated at room temperature for 3 min. Further, 10 μL of the blown cell suspension was taken and gently mixed with 10 μL of 0.4% Typan blue. Finally, 10μL of the obtained mixture was used to be quickly added to the Countess® II Cell Counting plate.

Library construction and sequencing

The prepared cell suspensions, together with 10×Barcode gel beads and oil, were added to Chromium Chip B microplates, and the beads and cells with Cell Barcode were encapsulated in droplets by microfluidic techniques. These cell-encapsulated droplets were collected and the cells in the droplets were lysed to form GelBeads-in-emulsion (GEM) through a microfluidic “double-cross” crossover system. GEM was then reversely transcribed using a PCR instrument, and RNA in the cells was reversely transcribed into the first strand of cDNA with the barcode Unique Molecular Identifier (UMI) information using the reverse transcription primers on the gel beads. The first strand of cDNA was further purified utilizing magnetic beads for RT-PCR amplification (12 cycles of amplification). After cDNA amplification was completed, it was fragmented using enzymes and the optimal fragment was screened with magnetic beads to construct a cDNA library by sequencing primers by PCR. The labeled cDNA library mix was sequenced on the Illumina NovaSeq platform. The sequencing was performed by the end-pairing procedure for double-end sequencing, and Illumina provided software to manipulate the sequencing process and to analyze the real-time data generated by sequencing.

Quality control of sequencing data

Sample quality control was conducted using Cell Ranger, the official software of 10x genomics, which integrates STAR software internally [9]. The reads were first compared to the reference genome to obtain quality control results such as high-quality cell count, gene count, and genomic matching rate in the raw data. Subsequently, the quality of each sample was assessed. The Seurat (version 4.3.0) [10] software package was adopted to further QC the experimental data based on the initial QC of Cell Ranger to exclude data that were multicellular, bicellular, or unbound on cells. Cells with a total gene count of less than 7000 and more than 200 and a mitochondrial UMI percentage of less than 10% were retained as high-quality cells. Double cells were removed for downstream analysis with the aid of the DoubletFinder (version 2.0.3) [11] software. Eventually, three cases of cervical cancer sequencing data and one case of single-cell sequencing data from paracancerous tissue were merged, with Harmony (version 0.1.1) [12] employed to remove batch effects during the merging process.

Data downscaling and clustering

The first 50 principal components were selected as statistically significant components for subsequent analysis, and the dimensionality reduction results were visualized in two dimensions using t-distributed random neighborhood embedding (t-SNE) plots [13, 14]. The research utilized the single-cell sequencing cell type annotation software Single R [15] which annotates cell types by calculating the correlation between the single-cell reference expression profile dataset and the cell expression profile to be identified and annotates the cell to be identified as a cell type with the highest correlation with the reference dataset. Further, the cell types were artificially annotated based on Marker genes, in combination with results from published literature (Supplementary File 1).

Detection of copy number variation in tumor cells

InferCNV [16] can be used to identify large-scale chromosomal copy number variations (CNV) in tumor single-cell RNA-Seq data, including amplifications and deletions of whole or large segments of chromosomes. This approach determines the intensity of gene expression by comparing the gene expression of each tumor cell with that of a reference cell (normal cells). Therefore, this study adopted the InferCNV (version 1.10.1) software package to detect CNVs in epithelial cells, identifying real tumor cells. The control group was 1 cluster of epithelial cells derived from normal tissue.

Cellchat analysis of intercellular communication

CellChat refers to an open-source R package (https://github.com/sqjin/CellChat), which can be applied for the analysis and visualization of scRNA-seq data intercellular communication (version1.6.1) [17]. CellChat is capable of characterizing different intercellular communication by a variety of learning methods including pattern recognition methods [18] and can combine the interaction between gene expression and signaling ligands and receptors to construct a probabilistic model of cell-to-cell communication. Such an analytical model can accurately elaborate on specific signals played by each cell population. Consequently, CellChat was used in this paper to analyze the communication between the nine cell types that were identified.

Proposed time series analysis of CAFs subtypes

The clustering of CAFs was downscaled and time series analysis of their subpopulations was proposed to visualize cell movement trajectories and map kinetic expression heat maps.

Transcription factor analysis of CAFs subtypes

PySCENIC [19] is an implementation for SCENIC [20] (Single Cell Regulatory Network Inference and Clustering) on python, which can infer transcription factors, gene regulatory networks, and cell types from single-cell RNA-seq data. Compared to SCENIC (which is based on the R language), pySCENIC offers a dramatic increase in computational speed. Since pySCENIC scales easily to multicore clusters when running on a computer, tens of thousands of cells can be analyzed in a relatively short period. This computing process is mainly achieved through the Dask framework of distributed computing. Hence, pySCENIC was utilized for transcription factor analysis of CAFs subsets in this study. Additionally, to distinguish cell subtypes more efficiently, the Ag20 (Average per 20 cells) [21] method was chosen for single-cell data and repeated three times to eliminate the error caused by random sampling.

Availability of data and materials

The dataset supporting the conclusions of this article is included within the article. The relevant codes can be found in Supplementary File 2.

Patient information

Basic information collected from 3 patients with cervical cancer is as follows (Table 1).

Results of dimensionality reduction clustering for single-cell data

Visualization of the clustering results by t-SNE revealed that all cells from 3 tumor samples and 1 normal sample from the database were clustered into 12 cell populations (Figure 1A). To clarify the type of each cell population, the marker gene was obtained for each cell population. Each cell subpopulation based on the above-identified marker genes was then annotated using Single R. The annotation results were based on marker genes from the Cell Marker database, combined with cell marker genes from the published literature for more detailed identification of each cell subpopulation. The results demonstrated that 12 cell subpopulations were annotated into 9 cell types: T cells, epithelial cells, smooth muscle cells, CAFs, endothelial cells, macrophage cells, B cells, lymphocytes, and plasma cells (Figure 1B). In the t-SNE plots according to tissue type, it was observed that immune cells (T cells, B cells, and plasma cells) occupied a large part of the cells in the tumor tissue (Figure 1C). The proportion of different cells in each sample was statistically analyzed and plotted as a histogram. It was found that the proportion of T cells was significantly higher in tumor samples than that in normal samples; while in normal samples, smooth muscle cells and endothelial cells were more predominant (Figure 1D). Marker genes of various cells were well distinguished for cell types, and each marker gene was abundantly expressed in its corresponding cell type (Figure 1E). Specifically, CD8A was highly expressed in T cells; KRT19 was highly expressed in epithelial cells; ACTA2 was highly expressed in smooth muscle cells; PDGFRA, FAP, and COL1A1 were highly expressed in CAFs; CD74 was highly expressed in macrophage cells; CD37 was highly expressed in B cells; KIT was highly expressed in lymphocytes; JCHAIN was highly expressed in plasma cells (Figure 1F).

Figure 1. Visualization of tSNE for cell descending clustering. (A) tSNE plot of 12 cell clusters; (B) tSNE plot of 9 cell types; (C) tSNE plot of normal and tumor cells; (D) Proportion of each cell type in different samples; (E) Heat map of marker gene expression for 9 types of cells; (F) Bubble map of marker gene expression for 9 types of cells.

GSEA pathway enrichment analysis for each cell type

ALL_OGRAFT_REJECTION, PIK_AKT_mTOR signaling pathway, α interferon response, γ interferon response were enriched in T cells; The IL6_JAK_STAT3 signaling pathway, complement, inflammatory response, and upregulation of KRAS signaling pathway, and TNFa/NF_kb signaling pathway were enriched in macrophages; The pathways of epithelial and mesenchymal transition and coagulation were significantly enriched in CAFs; Wnt/beta-catenin signaling pathway was significantly enriched in endothelial cells (Figure 2).

Figure 2. Enrichment pathways in each cell type. The colours represent the enrichment of the pathways, with the pathway enrichment increasing from blue to red.

CNV profile in tumor

Heterogeneity among tumors is often caused by copy number variation [22]. Cervical squamous carcinoma, the most common type of cervical cancer, is a squamous cell carcinoma of epithelial origin. This study used the gene copy number of normal epithelial cells as control and analyzed the degree of malignancy of epithelial-derived cells by the inferCNV algorithm. Epithelial cell clusters were selected from the cell population and separated according to tumor and normal tissues. As observed from the figure, a significantly higher copy number level occurred in the tumor group than in the normal group (Figure 3).

Figure 3. InferCNV heat map. The colours represent the copy number, which gradually increases from green to purple.

Intercellular communication network

Intercellular receptor ligands and signaling pathways were calculated using CellChat to analyze intercellular interactions. For analysis, the cells were separated into two groups according to normal and tumor tissues and the associations between different cell types in the two subgroups were calculated separately. In normal tissues, CAFs and other cells presented more interactions and stronger effects (Figure 4A); while in tumor tissues, endothelial cells, CAFs, smooth muscle cells, and other cells had more interactions and stronger effects (Figure 4B). This means that there was more than one cell involved in the development of a tumour.

Figure 4. Cellular communication networks. (A) Number of interactions and intensity of action of different cell types in the normal group; (B) Number of interactions and intensity of action of different cell types in the tumor group. The thickness of the line refers to the number of interactions and the strength of the effects.

Changes in intercellular communication in the tumor group relative to the normal group

After differential analysis of cellular communication between the tumor and normal groups, it was found that, relative to the normal group, endothelial cells, smooth muscle cells, CAFs, and macrophages were enhanced in the number of interactions with other cells in the tumor group, while the intensities of the effects of endothelial cells, smooth muscle cells, and other cells were improved, along with largely positive correlations (Figure 5).

Figure 5. Changes in the number of interactions and intensity of the effects of different cell types in the tumor versus normal group, with blue representing weakening and red representing intensification.

Communication between CAFs and other cells

From the obtained ligand-receptor pairs and signaling pathways, the MK (Midkine) pathway was chosen as the object of analysis. In both tumor and normal groups, CAFs could be the main source of MK action on several other cells in the MK pathway. This suggested that CAFs were the main mediators in the MK pathway, acting as gatekeepers for cell-cell communication. (Figure 6A, 6B).

Figure 6. Communication of CAFs with other cells in the MK signaling pathway. (A) Hierarchical diagram of the MK signaling pathway in the communication network of the normal group of cells; (B) Hierarchical diagram of the MK signaling pathway in the communication network of the tumor group of cells. This diagram consists of two parts: the left and right parts highlight the autocrine and paracrine signals of CAFs and several other cells, respectively. The solid and hollow circles represent the source and target, respectively. Different colors of the circles represent different cells, and the width of the connecting lines represents the communication probability. The color of the connecting lines is consistent with the source of the signal.

Role of CAFs in the Midkine signaling pathway

Cellchat takes the centrality metric previously used in graph theory for social network analysis to identify the major sources, targets, key mediators, and key influencers in cellular communication [17]. In contrast to other cells, CAFs acted as important signal sources, targets, mediators, and influencers in the MK signaling pathway in both normal and tumor groups (Figure 7A, 7B).

Figure 7. The role of CAFs in the MK signaling pathway in the normal group (A) and the tumor group (B).

Defining new subtypes of CAFs

To further investigate the biological functions of CAFs in cervical cancer, CAFs were re-clustered to distinguish different subgroups. CAFs were re-clustered into 11 clusters (Figure 8A), with cells identified according to the marker gene of each cluster. Three of these clusters, 4, 5, and 6, with high expression of C7 and CNN1, were named type1 cells; two clusters, 0 and 8, with high expression of SELENOM and RACK, were named type2 cells; and several clusters, 1, 2, 3, 7, 9 and 10, with high expression of CTGF and TNXB, were named type3 cells (Figure 8B).

Figure 8. Distinguishing new subtypes of CAFs. (A) tSNE plot of CAFs with reduced dimensional clustering; (B) tSNE plot of maker genes for 3 cell subpopulations.

Analysis of the cell developmental trajectories of each subpopulation of CAFs

It could be seen from Figure 9A that the three subpopulations were divided into three developmental states. Figure 9B showed the timeline of cell differentiation, with the shades of color representing early and late developmental stages. Figure 9C introduces the developmental stages in which different subpopulations of CAFs were located. As could be seen from Figure 9C, most type2 cells were at the early developmental stage of CAFs, most type3 cells were at the middle developmental stage of CAFs, and most type1 cells were at middle and late developmental stages of CAFs. The expression heat map of genes was drawn according to the developmental trajectory, and these genes were clustered into 5 modules (Figure 9D). From the heat map, it was observed that cluster1, cluster5, and cluster3 represented the early, middle, and late stages of CAFs development, respectively.

Figure 9. The pseudo-time trajectory of CAFs. (A–C) The pseudo-time trajectory of CAFs with gene expression profiles derived from Monocle 2, with each point representing a single cell; (D) Heat map of gene expression changes.

KEGG pathway enrichment analysis of genes in cluster 1, cluster 5, and cluster 3

The proteoglycans in cancer, axon guidance, and lysosome were significantly enriched in CAFs at early developmental stages (cluster 1) (Figure 10A); The pathways in cancer, IL-17 signaling pathway, and cytokine-cytokine receptor interactions were significantly enriched in CAFs at mid-developmental stages (cluster 5) (Figure 10B); The cell cycle, DNA replication and mismatch repair pathways were significantly enriched in CAFs at the late developmental stage (cluster 3) (Figure 10C).

Figure 10. KEGG enrichment analysis of CAFs developing in early (A) to mid (B) to late (C) stages.

Transcription factor prediction of CAFs subtypes

A comprehensive network analysis of pySCENIC results was performed to systematically identify key transcription factors of the cells. The activity of each transcription factor associated with three fibroblast subtypes was evaluated, leading to the definition of a transcription factor specificity score (RSS) based on the Jensen-Shannon scatter [23]. The positively regulated transcription factors with the highest RSS values were then selected and their functional properties were further examined. Therefore, TCF21, TFF3, and FOXN3 were identified as type1 cell-specific transcription factors (Figure 11A) and their specificity was validated using tSNE plots (Figure 11B, 11C); TAF1, ZC3H11A, MEIS2, RREB1, and ZNF471 were identified as type2 cell-specific transcription factors (Figure 11D) and their specificity was validated using tSNE plots (Figure 11E, 11F); FOSL2, MYEF2, GADD45, IES6, ZNF1 were identified as type3 cell-specific transcription factors (Figure 11G), and their specificity was also validated utilizing tSNE plots (Figure 11H, 11I). Furthermore, the accuracy of these transcription factor predictions was systematically evaluated using two methods, namely SEEK [24] and CoCiter [25]. SEEK is a method based on mining public datasets and CoCiter is an approach for mining the literature. Notably, SEEK should be employed first to filter a dataset in which a transcription factor and its target gene present co-expression. Usually, the functions of these co-expressed genes tend to reflect the functions of that class of cells. Consequently, it is effective to search the headings of the dataset to obtain the functions of the cells being studied. Secondly, since descriptions related to functionally relevant genes tend to appear in abstracts in the literature, this research identified the functions of these genes and corresponding cell types in the literature abstracts by CoCiter analysis.

Figure 11. Transcription factor prediction of CAFs subtypes. (A) Ranking of transcription factors of type1 cells based on transcription factor specificity score (RSS); (B) Highlighting type1 in the t-SNE plot (indicated by red dots); (C) Binary conversion of transcription factor activity scores for transcription factors with higher RSS TCF21 (Z-Score normalization for all samples with a 2.5 dividing line to convert values to 0 and 1), highlighting in t-SNE plots (indicated by green dots); (D) Ranking of transcription factors of type2 cells based on transcription factor specificity score; (E) Highlighting type2 in the t-SNE plot (indicated by red dots); (F) Binary conversion of transcription factor activity scores for transcription factors with higher RSS ZC3H11A, highlighting in t-SNE plots (indicated by green dots); (G) Ranking of transcription factors of type3 cells based on transcription factor specificity score; (H) Highlighting type3 in the t-SNE plot (indicated by red dots); (I) Binary conversion of transcription factor activity scores for transcription factors with higher RSS MYEF2, highlighting in t-SNE plots (indicated by green dots); (J) Co-expression results of TCF21 and its target genes in different GEO datasets retrieved using SEEK.; (K) Co-expression results of ZC3H11A and its target genes in different GEO datasets retrieved using SEEK; (L) Co-expression results of MYEF2 and its target genes in different GEO datasets retrieved using SEEK. The x-axis represents different datasets and the y-axis represents the statistical significance of the correlation between target genes in each dataset. Functionally related and significantly correlated (p < 0.01) datasets are indicated by yellow dots.

For type1 cells, SEEK was used to search for GEO datasets where TCF21 and its target genes were significantly co-expressed. Among the more than 3000 datasets detected, those associated with type1 cells were ranked high (Figure 11J, p<0.05). The two terms of “Breast” and “immune” from the descriptions of the datasets were extracted as a description of the function of type1 cells. For type2 cells, in the dataset of more than 3000 ZC3H11A detections showing co-expression with their target genes, most of those associated with type2 cells were ranked high (Figure 11K, p < 0.05). The two terms of “Lung” and “cervical” from the descriptions associated with the dataset were extracted. For type3 cells, SEEK was adopted to search the GEO dataset for significant co-expression of MYEF2 and its target genes. The higher-ranked ones were associated with type3 cells (Figure 11L, p < 0.05). The two terms of “CIN” and “Colon” were extracted from the descriptions of the dataset as the most frequently occurring terms for type3 cells.

Modularization of transcription factors

Transcription factors usually interact with each other to regulate gene expression. To describe this mode of action more systematically, this research compared transcription factor activity scores (RAS) based on the Connection Specificity Index (CSI) [26], allowing transcription factors with similarities to cluster into one category. The results indicated that these transcription factors were mainly distinguished into six modules. The binding sequences of transcription factors were subsequently mapped by selecting representative transcription factors in each module based on their average activity scores. The representative transcription factors of M1 included TCF21, FOSL2, and RREB1, which were three subtypes of cell-specific transcription factors (Figure 12A). When the average activity fraction of each module was expressed using tSNE, it was found that each module occupied a different region, along with the presence of a complementary relationship between all highlighted regions of these modules (Figure 12B). KEGG pathway enrichment analysis of transcription factors in M1 revealed that significantly enriched pathways were connected with cancer, transcriptional misregulation in cancer, viral oncogenic effects, and human papillomavirus infection (Figure 12C).

Figure 12. The modularity of transcription factors and pathway enrichment analysis. (A) Identification of transcription factor modules based on transcription factor ligand specificity indices, with representative transcription factors selected, and binding sequences of transcription factors plotted; (B) tSNE plot of the average transcription factor activity score for each module; (C) KEGG pathway enrichment analysis of genes in module 1.