Single-cell RNA-seq dissects the intratumoral heterogeneity of triple-negative breast cancer based on gene regulatory networks

Abstract

Triple-negative breast cancer (TNBC) is a subtype of breast cancer with high intratumoral heterogeneity. Recent studies revealed that TNBC patients might comprise cells with distinct molecular subtypes. In addition, gene regulatory networks (GRNs) constructed based on single-cell RNA sequencing (scRNA-seq) data have demonstrated the significance for decoding the key regulators. We performed a comprehensive analysis of the GRNs for the intrinsic subtypes of TNBC patients using scRNA-seq. The copy number variations (CNVs) were inferred from scRNA-seq data and identified 545 malignant cells. The subtypes of the malignant cells were assigned based on the PAM50 model. The cell-cell communication analysis revealed that the macrophage plays a dominant role in the tumor microenvironment. Next, the GRN for each subtype was constructed through integrating gene co-expression and enrichment of transcription-binding motifs. Then, we identified the critical genes based on the centrality metrics of genes. Importantly, the critical gene ETV6 was ubiquitously upregulated in all subtypes, but it exerted diverse roles in each subtype through regulating different target genes. In conclusion, the construction of GRNs based on scRNA-seq data could help us to dissect the intratumoral heterogeneity and identify the critical genes of TNBC.

Keywords: cell-cell communication; gene regulatory network; intratumoral heterogeneity; key regulator; single-cell RNA sequencing; triple-negative breast cancer; tumor microenvironment.

Graphical abstract

Figure 1

Cell composition of TNBC patients (A) Heatmap of the inferred CNV across 868 epithelial cells, in which genes were sorted by genomic location. (B) Standard deviation of the CNV of the cells in each cluster. Cluster 1 and cluster 3 showed a relatively high variable of CNV. (C) tSNE plot of the epithelial cells using the inferred CNV. (D) Heatmap of the gene-expression profiles of 545 malignant epithelial cells. (E) Breast cancer subtype composition of each TNBC patient based on the PAM50 model. The point in the petal represents a single cell. The number beside the patient identification (ID) is the tumor size.

Figure 2

Cell-cell interaction and hallmark gene set activity analysis (A) Cell-cell interaction network of different cell types. The node size represents the number of interactions. The width of edge represents the number of significant ligand-receptor interactions in two cell types. (B) Differences of the enrichment of the hallmark gene sets across the five molecular subtypes. The colors are encoded by the mean values of the GSVA enrichment scores in the molecular subtypes (one-way ANOVA, ∗∗∗p < 0.001; ∗∗p < 0.01; ∗p < 0.05; NS, p > 0.05). (C) Violin plots of GSVA enrichment scores of the G2/M checkpoint for five molecular subtypes.

Figure 3

GRNs of normal epithelial cells and five molecular subtypes Colored nodes imply the critical genes. (A−F) The GRNs of normal epithelial (A), normal-like (B), basal-like (C), Her2+ (D), LumA (E), and LumB (F).

Figure 4

Functional diversity of ETV6 (A) ETV6 regulates diverse genes in each subtype. The line colors indicate the ETV6 regulations in different subtypes. (B) GO annotation of ETV6-regulated genes in each subtype. The colored lattice showed the ETV6-regulated genes annotated in the corresponding GO term. (C) Expression of ETV6 in normal epithelial cells and five molecular subtypes. (D) Expression of ETV6 in TCGA normal and TNBC samples. (E) Survival analysis of TNBC patients in the METABRIC dataset based on ETV6 expression. The lower expression of ETV6 showed better clinical outcome.

Figure 5

Schematic diagram of this study First, the scRNA-seq data were used to infer CNV. The malignant cells were identified based on the inferred CNV. Then the subtype of each cell was assigned by using the PAM50 model, and the GRN was constructed for each subtype. Second, five centrality metrics were calculated to measure the importance of nodes. Then the Q statistic was used to integrate these centrality metrics. Finally, the diverse roles of the common critical genes were assessed. Differential expression analysis and survival analysis were performed for the critical gene.