An Integrative Model of Single Cell Transcriptomic States for Triple-Negative Breast Cancer

An Integrative Model of Single Cell Transcriptomic States for Triple-Negative Breast Cancer | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article An Integrative Model of Single Cell Transcriptomic States for Triple-Negative Breast Cancer Erik Reinhold Samuelsson, Roy Francis, Taobo Hu, Mats Nilsson This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-5974271/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract Triple-negative breast cancer (TNBC), an aggressive and highly heterogeneous subtype of breast cancer, poses significant challenges for treatment due to its molecular diversity and resistance to standard therapies. Accounting for 10–20% of all breast cancer cases, TNBC lacks specific biological markers, making it difficult to classify and treat effectively. Traditional approaches based on bulk RNA sequencing obscure intratumoral heterogeneity and fail to capture distinct cellular states within tumors. In this study, we constructed a comprehensive single-cell transcriptomic map of TNBC by analyzing a cohort of published TNBC patient datasets, identifying nine transcriptomic states, or metaprograms, which capture the core behaviors of TNBC cells, including cancer stem cell properties, epithelial-to-mesenchymal transition (EMT), immune modulation, metabolic adaptation, and vasculogenic mimicry. We observed that these metaprograms are variably expressed within and across patient tumors, underscoring the complexity of TNBC. By integrating TNBC-specific metaprograms with established breast cancer subtypes, we found significant prognostic associations, with specific metaprograms correlating with poor survival outcomes. This study highlights the need for single-cell approaches to uncover TNBC’s molecular heterogeneity and suggests that metaprogram-based classification could facilitate more precise therapeutic interventions. Triple-negative breast cancer Single-cell RNA sequencing Consensus Non-Negative Matrix Factorization Intratumoral heterogeneity Metaprograms Copy Number Variations Tumor microenvironment Prognostic biomarkers. Figures Figure 1 Figure 2 Figure 3 Figure 4 Background Triple-negative breast cancer (TNBC) is a subtype of breast cancer that lies at the forefront of difficulties due to heterogeneity and the need to develop targeted therapies. TNBC is not a biological classification, but an operational term that stemmed from the fact that patients with estrogen receptor-negative (ER-), progesterone receptor-negative (PR-), and human epidermal growth factor receptor 2-negative (HER2-negative) disease are insensitive to endocrine treatment and targeted therapies. TNBC encompasses a collection of distinct epithelial cancers that are vastly different in terms of their therapeutic response, clinical behavior, histologies, and genomic characteristics 1 . It accounts for 10–20% of all breast cancers, and is an aggressive form of breast cancer with a poorer prognosis than other types, with an average survival rate of 10.2 months from diagnosis 2 . According to the WHO 2019 Classification of Breast Cancer, TNBCs can include a range of histological classifications 3 , though the most common include metaplastic carcinoma, medullary carcinoma, invasive lobular carcinoma, and apocrine carcinoma 4 . This heterogeneity increases the difficulty of diagnosis and prognosis 5 , 6 , as well as treatment as there is a high degree of diversity within and between patients 6 , 7 , 8 , 9 , as well as a high degree of treatment resistance 5 . Compounding the issue of heterogeneity is the lack of consensus on stratifying TNBC into clinically actionable subtypes. Several efforts have been made to stratify TNBC into actionable subtypes, including transcriptomic 10 , 11 , 12 , proteomic 13 , epigenetic 14 and metabolic 15 . While there are observed overlaps between classification systems, there is still large variability and inconsistencies between them 16 . The challenges in implementing these systems stem from the diversity of datasets, poor integration between methods, and inconsistencies in classification. A further limitation lies in the bulk nature of such signatures; for these methods a single breast tumor is classified according to the dominant subpopulation into a predefined subtype. In reality, there can exist cell-state diversity both within the tumor and between patients that is masked at the bulk level. The tumor microenvironment's (TME) complexity complicates categorization, as similar transcriptomic and genomic signatures between malignant and nonmalignant cells can confound classification and treatment decisions 5 , 7 . A more precise, single-cell-based classification of TNBC can improve prognosis and better characterize the dynamic process of disease progression driven by intratumoral heterogeneity, which is shaped by diverse malignant breast cancer phenotypes, dynamic phenotype conversion, and clonal evolution. This is crucial, as these phenomena directly contribute to the development of metastatic disease and the emergence of treatment-resistant phenotypes 17 . To address this issue, we have assessed a cohort of TNBC patient single-cell RNA sequencing (scRNA-seq) datasets and identified nine transcriptomic states that underlie range of behaviors and driving forces in this disease; from cancer stem cell differentiation, Epithelial-to-Mesenchymal Transition (EMT), plasticity, metabolic regulation, and stress response, to clonal evolution, selection, immunogenic mimicry, vasculogenic mimicry, and prognosis. Methods Reference Dataset Generation Raw counts and associated metadata from each study were combined into a single dataset. The dataset was analyzed mostly using Scanpy 55 in Python. Cells with less than 200 detected genes and genes with less than 5 cells were discarded. QC metrics such as percentage of mitochondrial expression and ribosomal expression were computed. Data was normalized followed by log transformation and scaling. The effect of the total number of counts was regressed out. Cell cycle state was estimated. Doublets were estimated using Scrublet 56 . Top highly variable genes were identified followed by dimensionality reduction methods: PCA and UMAP. The data was integrated over 'datasets' using Harmony 57 . Cell types were estimated using Celltypist 58 , 59 . Tabula sapiens 11k ( https://cellxgene.cziscience.com/e/2ba40233-8576-4dec-a5f1-2adfa115e2dc.cxg/ ) and Kumar 2023 700k ( https://cellxgene.cziscience.com/collections/4195ab4c-20bd-4cd3-8b3d-65601277e731 ) were used as reference datasets. Datasets were classified into high, low and normal. High, low, and normal reference assignments made based on the site of the sample biopsy according to the study and the enrichment of cells before scRNA-seq; high was from TNBC patients with no enrichment, low was from TNBC patients with non-epithelial enrichment (usually immune) and reference normal were from non-TNBC patients. High and low datasets were used to infer copy number variants (CNV). InferCNV ( https://github.com/broadinstitute/inferCNV ) was run to estimate CNVs using gencode v45 as annotation. High datasets were used as reference. R package SCEVAN 32 was also run to identify CNVs using only the high samples. Metaprogram Generation To identify the heterogeneous transcriptomic programs present across the samples, we have applied an approach similar to the one pursued by Gavish et al 2023 31 , but with the aim of characterizing TNBC tumor-type specific programs rather than pan-tumor programs. We performed consensus Non-negative Matrix factorization (cNMF, https://github.com/dylkot/cNMF ) on the data of each individual patient to identify the programs of each sample. Since application of cNMF requires a “K” parameter that influences the results, we run cNMF using different values (K = 5,6,7,8,9,10), and generate 45 programs for each tumor. Each cNMF program is summarized by the top 100 genes representing that program based on cNMF coefficients. Then, we identified the most robust cNMF programs across the patient cohort as those that recur within the tumor (gene lists have at least 70% overlap), recur across tumor (have at least 20% overlap with any other cNMF program in other patients analyzed), and are then non-redundant within the tumor (rank programs by similarity with programs from other tumors, remove programs that have at least a 20% overlap with other programs within the same patient). From the robust cNMF programs, we clustered them together based on their Jaccard similarity, and identified the most consistent genes present across the grouped cNMF programs. This is how we end up with our meta-programs. The clustering process resulted in 30 metaprograms (clusters), each containing samples with high internal similarity. Assessment of Clustered cNMF's Comprising Each Metaprogram To validate the resulting metaprograms, we assessed the overall similarity via Jaccard indices of each cNMF assigned to a cluster and their tendency to be more similar to cNMFs within said cluster than outside of it. The Jaccard index, defined as the size of the intersection divided by the size of the union of two gene sets, was computed using the formula: $$\:Jaccard\:Index\:=\:\frac{\left|A\:\cap\:\:B\right|}{\left|A\:\cup\:\:B\right|}$$ , where A and B are two gene sets. First, we calculated the Jaccard similarity indices for all pairs of cNMFs within a dataset, resulting in a similarity matrix. This matrix represents the proportion of shared features between pairs of samples. Next, we grouped the cNMFs together based on the metaprogram clustering which grouped the samples into distinct clusters based on their similarity indices. The similarity matrix could then be visualized as a heatmap to provide a clear visual representation of the internal similarity of the cNMFs within clusters and the separation between different clusters. Cumulative Expression of Metaprograms Malignant cells from the TNBCMap were identified based on the iCNV attribute, specifically selecting those labeled as 'aneuploid'. Metaprograms were defined by the gene sets from each of the 9 metaprograms (Supplementary Table 3 ). The gene expression data for the malignant cells was normalized using StandardScaler and log-transformed to emphasize relative differences in expression levels. For each metaprogram, cumulative expression was calculated by summing the expression levels of its constituent genes. To quantify the relative expression score for each cell, the cumulative expression values were centered around zero, with positive values indicating higher expression and negative values indicating lower expression relative to the mean. Hierarchical clustering was performed on the cells using the average linkage method, and the resulting order was applied to the cumulative expression data. Scoring of Metaprograms in Cells To analyze the distribution of cells across different metaprograms, we used scanpy's tl.score_genes function to calculate gene scores for each metaprogram, using the list of genes present within the data.Each cell was then assigned to the MP with the highest score. Correlation of Metaprograms Across Cells: The generation of the correlation matrix involved the dataset utilized for this analysis comprising metaprogram (MP) scores for various cells. The correlation coefficients, specifically Pearson correlation coefficients, were calculated using the .corr() method in pandas. Comparison of TNBC Metaprograms to Reference Metaprograms The analysis involved three datasets: the nonmalignant TNBC cells, normal reference breast cells (Reed 2024), and Gavish 2023 Hallmarks of Cancer metaprograms. To determine the similarities between these gene sets, the Jaccard similarity index was calculated for each pair of columns (representing gene sets) from the metaprograms dataset and each of the nonmalignant TNBC and Reed2024 datasets. To extract the overlapping genes between the metaprograms, we then developed a custom function to calculate the overlap between genes in each pair of MPs across different datasets. This function iteratively compares each MP from one dataset with every MP from another, identifying common genes and recording the number and identities of overlapping genes. To ensure clarity, self-comparisons within the same dataset were excluded, and only unique comparisons were considered to avoid redundancy. For each comparison, we generated two DataFrames: one to capture the count of overlapping genes and another to store the actual overlapping gene names. These DataFrames were then used to filter overlaps with five or more genes, which were considered significant for further analysis. The results were saved into CSV files for documentation and further review. Finally, the significant overlaps were formatted and printed to provide a comprehensive view of the gene overlaps across different metaprogram datasets. Enrichment Assessment of TNBC Metaprograms to Lehmann and Burstein Subtypes To assess the enrichment of the metaprograms within defined breast cancer subtypes, we conducted an enrichment analysis using a custom Python script. We imported gene expression datasets for TNBC metaprograms, Lehmann and Burstein subtypes, nonmalignant TNBC clusters, and Reed 2024 clusters and developed a function to calculate gene set enrichment by identifying shared genes between each metaprogram and subtype. This function performed set intersections of non-null values in each column pair, quantifying both the count and specific list of enriched genes. To focus on meaningful enrichment results, we filtered pairs with five or more shared genes using a threshold-based function. We further evaluated the statistical significance of the enrichment for each gene set using Fisher's exact test, constructing contingency tables that included the size of each metaprogram (assumed to contain 50 genes), the size of the subtype (predefined per subtype), and an estimated total of 20,000 protein-coding genes in the human genome. These tables were processed with the fisher_exact function from the SciPy library, yielding p-values that quantify the likelihood of observing such enrichment by chance. CNV Assessment in TNBC Patient Data To assess copy number variations (CNVs) across patient samples, we utilized output data from SCEVAN. We began by loading the CNV data from a pre-processed pickle file and defined a list of patient IDs alongside corresponding sample keys. These mappings enabled accurate cross-referencing between patient data in the CNV file and additional clinical annotations. For each patient, we extracted CNV data based on their specific sample key, which allowed for the identification and categorization of subclonal CNV patterns. The CNV heatmaps were generated with a custom function that first filtered the CNV DataFrame to retain only relevant values, excluding columns with chromosomal coordinates and gene annotations. To provide additional clarity, chromosome labels were added at the midpoint of each boundary. Assessment of Metaprogram Residuals Across Patient and Subclones To investigate the association between chromosomal copy number variations (CNVs) and metaprogram (MP) assignments in individual patients, we calculated residuals using a chi-squared test on a patient-by-patient basis. We began by identifying each unique patient in our dataset based on the donor_id field from single-cell RNA sequencing (scRNA-seq) data annotations. For each patient, we grouped their data by subclone and assigned MP, constructing a contingency table that captured the distribution of MP assignments within each subclone. To prevent issues arising from zero values in the table, a small constant (0.5) was added to each cell. A chi-squared test of independence was performed on each patient's contingency table to determine if a significant association existed between subclone and MP assignment. This test provided the observed chi-squared statistic, degrees of freedom, and p-value, indicating if the observed MP distribution within subclones significantly differed from the expected distribution. Following this, we calculated residuals by subtracting the expected counts from the observed counts in each cell. These residuals reflect the extent to which each MP was over- or under-represented within subclones, with positive values indicating over-representation and negative values indicating under-representation. TCGA BRCA Analysis This study utilized publicly available gene expression and clinical data from The Cancer Genome Atlas (TCGA) for breast invasive carcinoma (BRCA) to perform a survival analysis focused on triple-negative breast cancer (TNBC) samples. The methodology involved several steps, including data acquisition, processing, and statistical analysis. The analysis involved retrieving, processing, and analyzing RNA-Seq gene expression data, combined with clinical data, to generate survival curves for metaprograms identified in TNBC. Cox Proportional Hazard Regression Analysis Gene expression data for the entire TCGA-BRCA cohort was obtained using the TCGAbiolinks package. The query targeted the "Transcriptome Profiling" data category, specifically RNA-Seq data processed using the STAR workflow, focusing on "Gene Expression Quantification." Data was downloaded and prepared using the GDCprepare function to obtain the FPKM-UQ (Fragments Per Kilobase of transcript per Million mapped reads - Upper Quartile normalized) gene expression matrix. Clinical data corresponding to the TCGA-BRCA cohort was retrieved using the GDCquery_clinic function, and relevant survival information (patient barcode, days to last follow-up, and vital status) was selected for further analysis. TNBC samples were identified using an external reclassification dataset provided in an Excel file, which contained barcodes of TNBC patients. These barcodes were extracted and cleaned to match the format of the barcodes in the gene expression matrix. The cleaned barcodes were then used to subset both the gene expression data and the clinical data, ensuring that only TNBC samples with corresponding expression data were included in the subsequent analysis. Univariate Cox proportional hazards regression was performed for each gene in the matched expression data using the coxph function from the survival package. The analysis evaluated the association between gene expression levels and patient survival (time to last follow-up and vital status). The hazard ratio (HR), 95% confidence intervals, and p-values for each gene were calculated and stored in a results table. The significant results ( p < 0.05) were filtered and merged with metaprogram gene lists derived from previous studies, categorizing genes into specific metaprograms associated with malignant, nonmalignant, and reference normal breast tissue samples. The final significant results were saved for further analysis and interpretation. The significant genes identified through Cox regression were used to filter and refine metaprogram gene lists. These filtered gene lists were restructured to exclude non-significant genes and saved in CSV format for further exploration. Kaplan-Meier Survival Curve Analysis Gene lists representing various metaprograms were loaded from a pre-filtered CSV file that contained genes previously associated with specific biological processes in TNBC. Clinical data for the TCGA-BRCA cohort was retrieved using the GDCquery_clinic function, providing vital status and survival information for each patient. Reclassification data from an external Excel file was used to identify and extract TNBC-specific patient barcodes, which were then matched with the clinical data. Necessary columns from the clinical data were selected and processed, including converting time-related variables to numeric format and categorizing patient vital status as "Alive" or "Dead." RNA-Seq gene expression data for the TCGA-BRCA cohort was queried and downloaded using the TCGAbiolinks package 60 . This data was processed to generate a gene expression matrix, focusing on "Primary Tumor" samples. The matrix was normalized using variance-stabilizing transformation (VST) from the DESeq2 package, ensuring that downstream analyses were conducted on homogeneously scaled data. For each metaprogram, the expression of genes within the metaprogram was extracted from the normalized gene expression matrix. The expression data for these genes was combined and averaged to compute a "module expression" score for each patient. This score was then stratified into "LOW" or "HIGH" groups based on median expression levels. The clinical and gene expression data were merged to facilitate survival analysis. The merged dataset included the calculated module expression score and clinical variables such as overall survival time and vital status. A Cox proportional hazards model was fitted to this dataset to evaluate the impact of each metaprogram's expression on patient survival. Survival curves were generated using the Kaplan-Meier method, with significance evaluated via log-rank tests. The survfit function from the survival package was employed to fit the survival curves, and the ggsurvplot function from the survminer package was used for visualization. The final survival plots were saved to provide a graphical summary of the association between metaprogram expression and patient survival outcomes. Comparison of Scored Metaprograms to TCGA PAM50 Subtype Classification We began by acquiring clinical and RNA-Seq gene expression data from the TCGA-BRCA cohort using the TCGAbiolinks package. Clinical data were retrieved through the GDCquery_clinic function, while RNA-Seq data were accessed via the GDCquery function, specifically focusing on "Transcriptome Profiling" with the "STAR - Counts" workflow. Additionally, we extracted patient barcode data from an external Excel file that contained reclassified sample information, specifically identifying TNBC cases. These barcodes were carefully matched with the clinical dataset to ensure accurate identification and alignment of TNBC samples. Gene expression data was downloaded and prepared using the GDCprepare function, followed by normalization through variance-stabilizing transformation (VST) using the DESeq2 package. To enhance data robustness, genes with low expression (fewer than 10 reads across all samples) were filtered out. Metaprogram-specific gene lists were then loaded from a predefined CSV file, and the normalized gene expression matrix was subset accordingly. For each metaprogram, cumulative expression scores were calculated by averaging the expression levels of all constituent genes. We further integrated the gene expression data with clinical and subtype information. Subtype data for the TCGA-BRCA cohort was retrieved using the TCGAquery_subtype function and merged with the clinical data and metaprogram scores. Relevant clinical variables, such as vital status and survival times, were selected and combined with the metaprogram expression scores and subtype information, resulting in a comprehensive dataset for survival analysis. The final dataset, which included metaprogram scores, clinical data, and subtype information, was saved as a CSV file. This dataset served as the foundation for subsequent survival analysis, providing key insights into the potential biological roles of metaprograms in TNBC and their association with patient outcomes. Hierarchical Clustering of TCGA Patients by Metaprogram Scores From the preceding analysis, clinical data and RNA-Seq gene expression data were obtained from the TCGA-BRCA project using the TCGAbiolinks package. Clinical data, including patient survival information and subtype classification, were merged with gene expression data, which had been variance-stabilizing transformed using the DESeq2 package. The combined dataset included MP scores for each TNBC patient, which were used for further analysis. To investigate the potential stratification of patients based on their MP scores, hierarchical clustering was performed on TNBC patients classified as having the Basal subtype using the Ward's method. The MP scores were standardized before clustering, and a specific cutoff distance was applied to define clusters. The resulting clusters were then merged back into the main dataset to associate cluster labels with patient survival data. Kaplan-Meier survival curves were generated for each identified cluster, and log-rank tests were conducted to assess the statistical significance of survival differences between clusters. The Kruskal-Wallis test was used to evaluate the differences in MP scores across clusters. Significant cluster comparisons were identified and further analyzed to calculate the mean and median survival times for each cluster, and these results were saved for further interpretation. Results Construction of the TNBC Map Given the high heterogeneity of TNBC, it is necessary to combine the datasets from a representative population of TNBC patients in order to establish common transcriptomic definitions of the malignant phenotypes across a sufficient number of patients’ scRNA-seq samples. We identified 9 publicly available TNBC scRNA-seq datasets that include data derived from TNBC patient biopsies or reference healthy breast tissue and included more than 500 cells per donor (Fig. 1 , Supplementary Table 1, Supplementary Table 2 ). After integrating the datasets, the TNBC Map consisted of 302,509 cells with 15,097 genes from 178 donors (Fig. 2 a). To identify malignant cells in the TNBC Map, we performed inferred copy number variation (CNV) assignment on the confirmed TNBC patient biopsy donors via inferCNV ( https://github.com/broadinstitute/inferCNV ). Aneuploidy is common in human cancers, with their resultant extensive genome-wide CNVs in chromosomal arms being a universal feature of human cancer 18 , 19 , 20 . We identified a total of 57,183 malignant cells from 46 donors. After excluding donors with fewer than 10 identified malignant cells, we were left with 57,104 malignant cells from 21 patients. As seen in Fig. 2 a, the aneuploid and diploid cells form largely distinct clusters along with their distinct CNV profiles. Some overlap between clusters suggests regions where gene expression profiles might share common features or transitional states (Fig. 2 b ) . Genomic and Transcriptomic Heterogeneity within TNBC To start, we assessed the CNV patterns present in each patient to assess aneuploidy. Not only can different cancers adopt different CNV patterns 21 , 22 , 23 , but CNV can be used to assess patient prognosis and therapeutic resistance 24 , 25 . In general, CNVs that have been consistently observed in TNBC are gains of 1q, 8q, and 10q, with losses of 5q and 8p 1 , 26 , 27 , 28 , 29 . Significantly, TNBC in and of itself does not have any CNVs used to classify it clinically 3 . Examination of the CNV patterns within patients reveals considerable heterogeneity in the gains and losses within donors ( Supplementary Fig. 1 ). The variability in CNV patterns observed across patients underscores the genomic complexity and heterogeneity within TNBC, and the fact that TNBC itself is a collection of different types of tumor. To better structure and facilitate interpretation of transcriptional ITH, we conducted an analysis of the metaprograms using the identified malignant TNBC cells. A metaprogram refers to a higher-level collection of gene expression patterns that capture more complex, often broader biological states or processes. Metaprograms represent overarching themes or states of cellular activity that can encompass multiple pathways and broader regulatory networks. Each metaprogram represents a group of genes that are co-expressed and potentially co-regulated, reflecting specific biological processes or cellular states. The identification of metaprograms has been used effectively to characterize transcriptomic states in tumors before, either in the case of Neftel et al 2019 which identified 6 metaprograms in glioblastoma 30 , or Gavish et al 2023 which identified hallmark metaprograms present across many types of cancer 31 . The heatmap in Fig. 2 c reveals several distinct clusters of Non-negative Matrix factorization (NMF) programs, each represented by blocks along the diagonal of the matrix. Given that these metaprograms were identified specifically from malignant, aneuploid cells in scRNA-seq data, they reflect the intrinsic properties of the cancer cells themselves (see Supplementary Table 3 for the genes in each metaprogram). To assess the functions indicated by each metaprogram, we first assigned each malignant cell into a metaprogram based on the highest scoring metaprogram per cell. We then extracted the top genes by log fold expression in the population of cells assigned to each metaprogram, and used these genes to perform gene set enrichment analysis (GSEA). We extracted the most significant ( p -value < 0.05) upregulated pathways and functions to assess the functions most associated with each malignant metaprogram (Fig. 3 , Supplementary Fig. 2 ). TNBC_MP1, which we term Immune-like 1 (IM-like 1), is characterized by immune-related functions, primarily involving macrophage and monocyte activation, inflammatory signaling, cytokine production, and interactions within the tumor microenvironment, with a notable role in immune checkpoint pathways such as PD-1, suggesting immune suppression and potential responsiveness to immunotherapy. TNBC_MP2, which we term Proliferating, is functionally associated with high cell proliferation, DNA repair, and aggressive cancer traits, characterized by pathways involved in cell cycle regulation, mitotic control, and poor prognosis, suggesting a metaprogram that drives rapid tumor growth and potential metastatic behavior in TNBC. TNBC_MP3 is associated with metabolic activity and sees high expression across most malignant cells in our cohort. This suggests that TNBC_MP3 is associated with high metabolic activity, which is often necessary for supporting the rapid growth and energy demands of cancer cells. TNBC_MP4, which we term mesenchymal-like (MES-like), is enriched for genes involved in ECM remodeling, mesenchymal and stromal characteristics, consistent with roles in EMT and stromal interactions, particularly resembling signatures seen in certain breast cancer subtypes (like metaplastic carcinoma) and developmental stromal programs. TNBC_MP5, which we term IM-like 2 (IM-like 2) likely represents an immune-related metaprogram characterized by T-cell and NK-cell activation, cytotoxic functions, and checkpoint pathways, suggesting an immune-active TME focused on adaptive immunity and inflammation. TNBC_MP6 likely represents a hypoxia-driven, basal-like phenotype characterized by EMT, metabolic reprogramming, and pro-inflammatory signaling associated with aggressive and invasive tumor behavior (Fig. 3 c ) . The TNBC_MP7 metaprogram demonstrates luminal-like characteristics, with significant enrichment in luminal-associated pathways and a lack of alignment with basal or mesenchymal traits (Fig. 3 d ) . TNBC_MP8, which we term immune-like 3 (IM-like 3), is primarily associated with immune and interferon responses, including antiviral defense, inflammatory signaling, and potential roles in EMT and immune cell recruitment, which may support immune modulation and cancer cell invasion. TNBC_MP9, which we term vascular mimicry-like (VM-like) is primarily associated with endothelial and vascular functions, including angiogenesis, vascular integrity, extracellular matrix organization, and potentially facilitating tumor-associated vascular mimicry to support tumor growth and metastasis. The absence of significant GSEA results may suggest that TNBC_MP9 could reflect a vascular phenotype in TNBC that does not closely align with common reference datasets, indicating its specific relevance to tumor vasculature and possibly niche-specific endothelial roles. These results indicate that while malignant cells in TNBC exhibit significant genomic and transcriptomic heterogeneity, it is possible to identify consistent states present across this cohort. The distinct metaprograms—ranging from immune-modulating to proliferative, metabolic, and vascular phenotypes—highlight the diverse functional landscapes within TNBC tumors, each potentially contributing differently to tumor progression, treatment response, and immune evasion. Characterizing ITH in TNBC: The Relationship between Malignant cells, Subclones and Metaprograms Our next step, after identifying a range of metaprograms present in our TNBC cohort, was to assess the distribution of scores across all malignant cells in our cohort (Fig. 3 ). We did this via scoring all individual malignant cells for each metaprogram and assessing their distribution across various conditions. To start, we scored each cell for eight metaprograms (MPs), excluding TNBC_MP3 (metabolic activity), which showed high expression across all cells and were thus not informative for distinguishing cellular states. We then assigned each cell to the MP for which it scored the highest. This approach is specifically to determine which metaprograms are most prominent within malignant cells and provide a clear view of the dominant biological processes in play. The majority of cells were assigned to TNBC_MP6 (Basal-like) and TNBC_MP7 (Luminal-like), with 31,621 and 13,894 cells, respectively (Fig. 3 a). This suggests that these two metaprograms are key states driving the malignant phenotype in TNBC. In contrast, TNBC_MP8 (IM-like 3, interferon signaling) and TNBC_MP5 (IM-like 2, cytotoxic immune activity) were assigned to fewer cells, suggesting that immune-related processes are less prevalent within the bulk malignant population. Similarly, TNBC_MP4 (MES-like) exhibited moderate representation, underscoring its role in tumor-stroma interactions, albeit in a less dominant manner across the malignant cells. Notably, TNBC_MP1 (IM-like 1, immune mimicry) and TNBC_MP9 (vasculogenic mimicry) were found in the fewest cells, indicating that these processes might be niche-specific rather than widespread drivers of malignant activity. These patterns suggest that while immune evasion, immune mimicry, and vasculogenic mimicry are present, they may be more localized to specific niches, such as perivascular regions or tumor margins, rather than predominating within the tumor core. Overall, these results underscore the central role of epithelial lineage development-like programs in TNBC progression, while also reflecting the complexity and heterogeneity of the tumor microenvironment. To assess whether the metaprograms could reflect genetic subclones in the malignant cells, we used SCEVAN 32 to identify distinct subclones within patient donor malignant cells based on their CNV profile. From the 19 patients that passed quality control for SCEVAN subclone analysis, we identified a total of 102 subclones, with 1–7 distinct subclones in each tumor. Next, we assigned each malignant cell into a metaprogram based on the highest scoring metaprogram per cell. Notably, each of the subclones contained cells in multiple metaprograms (Fig. 2 d). Out of the 19 donors, 8 donors contain cells in all assigned metaprograms. To further investigate whether the MPs are associated with specific genetic subclones, we assessed the residuals in MP-assigned cells across these subclones. Residuals, which represent the difference between the observed and expected counts of cells assigned to each metaprogram within each subclone, allow us to quantify patterns of over- or under-representation of transcriptional states in the context of genetic subclones. Expected counts are determined by assuming a random distribution of metaprogram assignments across all subclones within each patient, based on the overall metaprogram proportions of that patient ( Supplementary Fig. 3 ). There are some trends in the representation of metaprograms within specific subclones that underscore the substantial degree of intratumoral heterogeneity across patients, with metaprograms activity varying across subclones. We observed that both TNBC_MP6 (basal-like) and TNBC_MP7 (luminal-like) states are commonly present within the same patient. Moreover, the basal-like and luminal-like states frequently display an inverse relationship comparing sub-clones within the same patient, with one being overrepresented where the other is underrepresented, hinting at distinct functional roles or selective pressures acting on these metaprograms, causing switches between these two programs. This detailed residual analysis provides insights into patient-specific subclonal evolution, highlighting key metaprograms that could influence tumor behavior and inform targeted treatment responses on a subclone-specific basis. Breast cancer, including TNBC, is known for its phenotypic plasticity 33 , 34 , allowing malignant cells to shift states as the tumor develops and encounters different conditions 35 , 36 , 37 , 38 . While our work identifies discrete transcriptomic states, TNBC operates as a continuum, with plasticity enabling multiple metaprogram expressions within individual cells. To gauge the expression of identified metaprograms, Fig. 3 b illustrates the cumulative gene expression of TNBC_MP1 to TNBC_MP9 in malignant cells. The heatmap reveals substantial heterogeneity, with distinct clusters of cells expressing specific metaprograms. For instance, TNBC_MP2, linked to cell cycle progression, is strongly expressed in a proliferative cell cluster, while TNBC_MP6 (basal-like, EMT/hypoxia) and TNBC_MP9 (VM-like, angiogenesis) are elevated in separate clusters, indicating specialized roles. Some metaprograms, like TNBC_MP3 (metabolic) and TNBC_MP4 (MES-like, ECM remodeling), show broader expression across cells, suggesting general importance in tumor maintenance. Co-expression patterns, such as TNBC_MP2 (cell cycle) and TNBC_MP7 (luminal-like), indicate that certain subpopulations, like luminal-like cells, are also engaged in robust proliferation. Additionally, TNBC_MP8 (IM-like 3, interferon signaling) appears functionally segregated, marking a subset potentially involved in immune response. Overall, the distinct patterns across metaprograms highlight the diversity of malignant cell states, reflecting common functional responses or adaptations to environmental cues. To identify metaprograms with similar correlation patterns, we performed hierarchical clustering on a correlation matrix of the metaprogram scores in the malignant cells of the TNBC Map (Fig. 4 a), revealing how different biological processes are coordinated or segregated within the tumor. TNBC_MP2, associated with cell cycle progression, shows minimal correlation with other metaprograms, suggesting that proliferative cells operate independently from immune or metabolic functions. In contrast, TNBC_MP9 (VM-like), linked to angiogenesis, has a strong positive correlation with TNBC_MP4 (MES-like, ECM remodeling), indicating that cells involved in vascular development are also engaged in modifying the extracellular matrix, supporting structural changes that facilitate tumor invasion. TNBC_MP1 (IM-like 1, immune mimicry) is strongly correlated with TNBC_MP5 (IM-like 2, cytotoxic activity) and TNBC_MP8 (interferon signaling), suggesting that immune-mimicking cells participate in cytotoxic and inflammatory responses, potentially enhancing immune modulation within the tumor. Notably, TNBC_MP7, characterized by luminal-like gene expression and associated with cell adhesion, shows negative correlations with nearly all other metaprograms, suggesting it represents a distinct subset focused on structural roles. Overall, these correlation patterns highlight the functional specialization among different malignant cell subsets, with proliferative, immune-modulatory, and structural roles operating in tandem to support the tumor’s complexity and adaptability. We also assessed the correlation patterns patient-by-patient ( Supplementary Results 1: Sample Specific MP correlations ) to identify cohort and patient-specific patterns. The sample-by-sample comparison to the overall correlation matrix reveals key trends in the relationships between MPs across the dataset. Some MP pairs, such as TNBC_MP1 and TNBC_MP5 (IM-like 1 and 2) exhibit stable positive correlations across both individual samples and the overall dataset. These consistent relationships suggest that the co-activation of immune-related processes are fundamental aspects of tumor biology across the cohort. Similarly, the stable negative correlation between TNBC_MP7 (Luminal-like) and TNBC_MP3 (metabolic) implies distinct cellular states that may rarely overlap. In contrast, other MP pairs show more variability in their correlations between individual samples and the overall dataset. For instance, the correlation between TNBC_MP8 (IM-like 3) with TNBC_MP6 (basal-like) and TNBC_MP9 (VM-like) varies across samples, indicating that their co-activation is more context-dependent, potentially influenced by specific tumor characteristics or environmental factors. Additionally, while the overall correlation matrix tends to smooth out variations, closer inspection of individual samples reveals important differences that highlight tumor heterogeneity and patient-specific differences. These insights underscore the importance of considering both overall trends and sample-specific variations to fully understand the complexity and diversity within TNBC tumors, which may have significant implications for therapeutic strategies. Characterization of TNBC Metaprograms via Comparison to Other Data Sets As there are parallels between normal breast and malignant development in breast cancer, we conducted a comprehensive analysis using gene expression data from various studies to characterize the range of shared features between our derived TNBC MPs and MPs derived from other relevant data sets. Specifically, we examined the overlap among TNBC_MP1-9 derived from malignant cells in the TNBCMap, MPs from nonmalignant cells from TNBC patients in the TNBCMap, MPs from the Reed et al. 2024 atlas of breast cells 39 , and MPs from the Gavish et al. 2023 study on biological Hallmarks across many varieties of cancer 31 ( Supplementary Table 3 ). The heatmap ( Supplementary Fig. 7 ) of the Jaccard similarity indices between gene sets of various metaprograms and clusters from the nonmalignant, Reed2024 and Hallmark datasets reveals shared features between normal breast programs, general cancer programs and our TNBC metaprograms. Several key observations can be made from the heatmap. TNBC_MP1, TNBC_MP5, TNBC_MP7, and TNBC_MP8 exhibit low similarity but some overlap with specific clusters, indicating connections to immune response (TNBC_MP1), immune signaling (TNBC_MP5), luminal cells (TNBC_MP7), and inflammation (TNBC_MP8). TNBC_MP2, TNBC_MP3, and TNBC_MP6 display moderate overlap with clusters related to cell cycle, metabolism, and basal cells, reflecting shared processes across malignant and general cancer contexts. TNBC_MP9 shows high similarity with Reed2024_MP4 and NC9, highlighting the similarity of vascular processes in malignant cells with nonmalignant counterparts. These observations are further strengthened when assessing the overlapping genes. See Supplementary Results 2 for a full list of overlapping genes with the TNBC Metaprograms. Next, we validated the single-cell metaprograms against the bulk RNA subtypes found in the Burstein 40 and Lehmann 10 , 11 . Given that bulk RNA expression profiles represent mixtures of diverse malignant and nonmalignant cell types and states, the expression of each metaprogram in bulk samples provides an estimate of the abundance of the corresponding cellular state within the TME. The analysis comparing TNBC metaprograms derived from single-cell RNA sequencing with TNBC subtypes defined by bulk RNA sequencing data reveals key insights into the relationships between cellular states and the tumor microenvironment (TME) across different TNBC subtypes. Given that bulk RNA expression profiles represent mixtures of diverse malignant and nonmalignant cell types and states, the expression of each metaprogram in bulk samples provides an estimate of the abundance of the corresponding cellular state within the TME (see Supplementary Table 5 for extracted Lehmann and Burstein subtype genes). In Fig. 4 b, the heatmap of the Jaccard similarity indices between gene sets of various metaprograms and bulk RNA subtypes reveals shared features between our single-cell and the established bulk RNA signatures. Several key observations can be made from the heatmap. It is likely that the bulk subtypes masked the malignant MPs that are most similar to immune and stroma components. TNBC_MP3 and TNBC_MP7 show minimal overlap with TNBC subtypes, suggesting they capture rare or specialized cell states underrepresented in bulk RNA data. TNBC_MP2 overlaps with Lehmann subtypes related to cell cycle, while TNBC_MP1, 5 and 8 overlap with immune-modulatory subtypes. TNBC_MP4 and 6 align with basal-like and mesenchymal features in Lehmann subtypes. TNBC_MP8 overlaps with Burstein subtypes, reflecting immune and mesenchymal cell states. These findings highlight the importance of integrating single-cell and bulk RNA data to better understand the tumor microenvironment and its impact on TNBC subtype biology, and the enrichment analysis reveals key overlaps between TNBC MPs and cancer-related processes, emphasizing the diverse cellular states contributing to TNBC heterogeneity. MPs like TNBC_MP1 and TNBC_MP4 capture immune and mesenchymal features critical to TNBC subtype classification, while others, like TNBC_MP6 and TNBC_MP9, highlight distinct biological processes such as EMT and angiogenesis essential for tumor progression. Malignant and Nonmalignant Metaprograms Influence TNBC Patient Survival To systematically examine the association between transcriptomic states and patient outcome, we next performed an analysis of 192 identified TNBC bulk specimens from Cancer Genome Atlas Breast Invasive Carcinoma (TCGA-BRCA), which we identified via reclassification of the TCGA-BRCA samples performed in another study 41 . Bulk RNA expression profiles reflect mixtures of diverse malignant and nonmalignant cell types and states present in the TME, and therefore, the expression of each metaprogram defines a rough estimate for the abundance of the corresponding cellular state in bulk samples. Before generating the survival curves for the extracted TCGA samples based on low or high cumulative expression of the genes in the metaprograms, we perform univariate Cox proportional hazard regression analysis on each gene individually. The point of this is to evaluate the association between the expression of each gene and survival, then select the genes with the strongest association (e.g., based on the p-value) for inclusion in the metaprogram survival analysis ( see Supplementary Table 6 .) We then scored each bulk sample for expression of the survival-associated genes in each metaprogram from the 9 malignant TNBC, 14 nonmalignant TNBC and 32 reference normal breast metaprograms and examined the association between high or low expression and survival (Fig. 4 C, Supplementary Figs. 4, 5 ). The survival analyses for various metaprograms, including those derived from malignant TNBC cells, nonmalignant cells, and normal breast tissue (Reed2024 metaprograms), reveal interesting insights into their prognostic relevance in TNBC. For TNBC metaprograms ( Supplementary Fig. 4C ), TNBC_MP3 shows that higher expression correlates with significantly worse survival outcomes ( p = 0.018). TNBC_MP4 also associates higher gene expression with poorer survival ( p = 0.015). In contrast, TNBC_MP8 indicates that lower expression of its genes is linked to poorer survival ( p = 0.019). These findings underscore the distinct biological impacts of malignant metaprograms associated with enhanced metabolism, a mesenchymal transcriptomic state and immune interferon response can have on patient outcome. The nonmalignant cell (NC) metaprograms further highlight the critical role of nonmalignant cells in the tumor microenvironment and their influence on disease outcomes ( Supplementary Fig. 4 ). These NC metaprograms demonstrate the significant impact of nonmalignant TME on TNBC outcomes. The Reed2024 metaprograms, derived from normal breast tissue, provide additional insights into how these normal tissue signatures are associated with TNBC survival ( Supplementary Fig. 5 ). Although these Reed2024 metaprograms originate from normal breast tissue, their dysregulation in TNBC is associated with processes that drive tumor aggressiveness and reduced survival outcomes. Overall, these analyses underscore the diverse biological impacts of the malignant, nonmalignant, and normal breast metaprograms on TNBC prognosis. Next, we sought to validate the biological relevance of the metaprograms within the broader context of breast cancer subtypes, assess how well these TNBC-specific metaprograms align with or diverge from established molecular classifications, and potentially uncover novel insights into the heterogeneity and underlying biology of TNBC. This analysis aimed to explore whether these metaprograms could refine our understanding of TNBC beyond the traditional PAM50 framework, with implications for more personalized treatment strategies. The PAM50 panel is a list of 50 genes used to classify breast cancers into five intrinsic subtypes in FFPE tissue sections via real time polymerase chain reaction (RT-PCR) 42 . Their utility lies in that PAM50 subtypes have been found independently prognostic for long-term breast cancer survival. This classification system was furthermore developed with the major types of breast cancer (HER2+, ER, and PR positive) in mind, not for TNBC. To compare the malignant metaprograms to the PAM50 classification, we scored each metaprogram in the identified TNBC TCGA patient samples and hierarchically clustered them to assess patient stratification (Supplementary Fig. 6). After that, we compared them to their assigned PAM50 classification ( Supplementary Table 7 ). In our cohort, 87% of the samples were Basal, 8% Her2, 3% Normal, and 2% LumA. This is in line with what is known about PAM50 classification of TNBC 43 . After scoring each patient in the TCGA TNBC cohort for the malignant MP mean expression, we performed hierarchical clustering to identify five groups of patients in the Basal, and for the significantly different clusters, observed a median survival time of 1398 days in Cluster 2, consisting of 33 patients, compared to 1813 days for Cluster 5, consisting of 43 patients ( Supplementary Table 8 ) differences in survival time of defined patient groups based on MP. Overall, these results underscore the critical role of metaprogram expression levels in predicting survival outcomes for TNBC patients, highlighting the potential for these metaprograms to serve as valuable biomarkers for clinical prognosis and personalized treatment strategies. Discussion & Conclusions In this study, we used single-cell RNA sequencing data to construct a comprehensive molecular map of TNBC and to address the significant challenge of intratumoral heterogeneity. TNBC’s aggressive nature and poor prognosis is driven by its heterogeneity, which complicates diagnosis, treatment selection, and patient outcomes. Using our TNBCMap, we identified distinct transcriptomic signatures, revealing critical insights into the molecular diversity of TNBC. However, as our analysis is based on publicly available datasets, the full spectrum of TNBC heterogeneity may not be entirely captured. While the publically available datasets included 53 TNBC patients, only 21 patients had sufficient cells for further analysis of the malignant component after filtering, which may impact the generalizability of our findings. Future studies incorporating a larger and more diverse cohort, will be necessary to further validate these transcriptomic states. Our findings reinforce the idea that TNBC is not a single disease but a collection of tumors with varying molecular and clinicopathological characteristics, based on the differing CNV profiles between the patients in our cohort. However, the identification of nine key transcriptomic metaprograms highlights the common biological processes that malignant TNBC cells co-opt to enhance tumor progression, including immune evasion and mimicry, EMT, metabolic reprogramming, cell cycle regulation, and ECM remodeling, vasculogenic mimicry, and finally, luminal- and basal-like transcriptomic states. Our findings underscore the pivotal roles of vasculogenic and immunogenic mimicry in TNBC progression, highlighting their contributions to intratumoral heterogeneity and poor patient outcomes. The vasculogenic mimicry metaprogram (TNBC_MP9) shows that TNBC cells can adopt vascular-like properties, facilitating alternative vascular networks that correlate with more aggressive phenotypes and poorer survival, as supported by the literature 44 , 45 , 46 . Similarly, the metaprograms associated with immune mimicry (TNBC_MP1 and 5) reveals how TNBC cells co-opt immune signatures to evade immune surveillance, consistent with reports that immunogenic mimicry aids cancer in escaping immune detection 47 , 48 . By associating these mimicry-driven metaprograms with patient outcomes, we offer a unique single-cell perspective that validates and extends previous findings, suggesting that targeting these mimicry processes could be a promising therapeutic strategy in TNBC. The expression of an inflammatory metaprogram (TNBC_MP8) highlights the multifaceted role of inflammatory pathways within the tumor microenvironment, particularly in driving immune-modulatory behaviors. The strong enrichment of interferon-stimulated genes within this metaprogram suggests that TNBC cells actively utilize inflammatory signaling to shape immune responses. Studies have shown that interferon-gamma (IFNγ) can modulate antigen processing and presentation in TNBC cells, enhancing immune visibility 49 . At the same time, chronic inflammation, as seen in inflamed TNBC subpopulations, has been linked to chemotherapy resistance and genomic instability 50 . Although some research has associated inflammation with improved outcomes in TNBC 51 , our analysis did not reveal a direct survival association for TNBC_MP8, underscoring the complex and context-dependent role of inflammation. These findings highlight the dualistic nature of inflammatory signaling, both activating immune responses and promoting tumor persistence, and suggest that targeting inflammatory pathways in TNBC could reduce pro-tumorigenic inflammation and potentially improve therapeutic outcomes. The mesenchymal-like (TNBC_MP8, MES-like) phenotype, characterized by EMT, has significant implications for the progression and treatment of TNBC. EMT is linked to enhanced cellular motility, invasiveness, and resistance to apoptosis, making MES-like TNBC particularly aggressive and chemoresistant. MES-like states are often associated with poor patient outcomes due to their ability to evade standard therapies that target epithelial-like cells, so targeting the EMT process, either by reversing it through mesenchymal-epithelial transition (MET) or inhibiting key EMT drivers, can sensitize MES-like cells to treatment. These findings highlight the need for therapeutic strategies that disrupt the plasticity between EMT and MET states and cut off an avenue of evasion in TNBC, offering the potential to improve treatment outcomes by preventing metastasis and overcoming chemoresistance. However, EMT is a highly plastic and dynamic process, and static transcriptomic snapshots may not fully capture the temporal transitions between epithelial and mesenchymal states, which could influence therapeutic strategies. In terms of metabolism, the fact that we observe a metaprogram associated with upregulated metabolism (TNBC_MP3) and enhanced glycolysis is supported by literature, as TNBC generally exhibits a higher rate of glycolysis in general and compared to other types of breast cancer. One study identified a metabolic signature of enhanced glycolysis and lactate secretion. The identification of a proliferative metaprogram (TNBC_MP2) emphasizes the role enhanced proliferation has in TNBC, corroborating our results by earlier work that identified proliferation highly associated with luminal transcriptomic signature. The preferential association of certain subclones with specific metaprograms in TNBC suggests that these subpopulations occupy distinct niches within the tumor, driven by selective pressures that promote diverse adaptive traits. However, our dataset includes a limited number of patients, and larger studies incorporating longitudinal samples could provide further insights into the evolution and stability of these transcriptomic states over time. Karaayvaz et al. 2018 demonstrated that these subclonal populations not only coexist but also display gene expression profiles linked to poor patient outcomes, indicating a predisposition to resist therapy and drive metastasis 52 . Kim et al. 2021 found that TNBC cells under chemotherapeutic stress often harbor resistant subclones with metabolic adaptations like enhanced oxidative phosphorylation, allowing them to survive treatment and contribute to recurrence and metastasis 53 . Mavrommati et al . 2021 further emphasized that the persistence of distinct, resilient subclones complicates treatment, as single-target therapies may fail to eliminate the full spectrum of tumor heterogeneity, leaving resistant clones to propagate 54 . Consequently, effective TNBC treatment strategies must adopt multi-targeted approaches to address the range of subclonal states and prevent relapse by inhibiting multiple pathways concurrently. Notably, certain metaprograms, such as those related to cell proliferation, basal-like characteristics and vasculogenic mimicry, were strongly associated with poorer survival outcomes, underscoring their role in TNBC malignancy. However, as our survival analysis is based on bulk RNA-seq data from TCGA, it does not fully capture the single-cell-level heterogeneity within tumors, highlighting the need for validation using single-cell or spatial transcriptomics approaches. This finding is further supported by recent research demonstrating that targeting the basal-like to luminal-like state transition offers a potential therapeutic strategy. A recent study by Schade et al. revealed that combining AKT and EZH2 inhibitors drives basal-like TNBC cells into a more differentiated, luminal-like state, resulting in substantial tumor regression 38 . This work underscores the importance of these two cell states in TNBC progression and suggests that therapies targeting this transition could improve patient outcomes by disrupting aggressive basal-like phenotypes and promoting differentiation. Additionally, our analysis extends beyond malignant cells to the tumor microenvironment, revealing that non-malignant cells, including immune and stromal components, play a significant role in tumor development. These findings support the notion that TNBC progression is shaped by both intrinsic tumor characteristics, extrinsic interactions within the microenvironment and evolutionary selection. This study underscores the necessity of embracing TNBC’s heterogeneity to develop more personalized therapeutic strategies based on the varying proportions of common characteristics in this disease. By identifying molecular subpopulations within TNBC, our work provides a foundation for the future integration of multi-omic approaches, potentially guiding more precise clinical interventions. Ultimately, understanding the complexity of TNBC at the single-cell level will be pivotal in overcoming current treatment challenges and improving patient outcomes. Abbreviations TNBC: Triple-Negative Breast Cancer scRNA-seq: Single-Cell RNA Sequencing ITH: Intratumoral Heterogeneity TME: Tumor Microenvironment EMT: Epithelial-to-Mesenchymal Transition VM-like: Vasculogenic Mimicry-like NMF: Non-Negative Matrix Factorization CNV: Copy Number Variation MP: Metaprogram GSEA: Gene Set Enrichment Analysis PAM50: Prediction Analysis of Microarray 50 (breast cancer classification system) TCGA-BRCA: The Cancer Genome Atlas - Breast Invasive Carcinoma IM-like: Immune-like MES-like: Mesenchymal-like cNMF: Consensus Non-Negative Matrix Factorization SCEVAN: Single-Cell Evolutionary Analysis (CNV-based subclone analysis) VST: Variance-Stabilizing Transformation RT-PCR: Real-Time Polymerase Chain Reaction FFPE: Formalin-Fixed Paraffin-Embedded QC: Quality Control UMAP: Uniform Manifold Approximation and Projection PCA: Principal Component Analysis HUGO: Human Genome Organization FPKM-UQ: Fragments Per Kilobase of transcript per Million mapped reads - Upper Quartile NC: Nonmalignant Cluster GDC: Genomic Data Commons DESeq2: Differential Expression Analysis for Sequence Count Data FDR: False Discovery Rate cBioPortal: Cancer Genomics Data Portal CoxPH: Cox Proportional Hazards Model MPs: Metaprograms Declarations Ethics approval and consent to participate Not applicable. Consent for publication Not applicable. Availability of data and materials The processed data is stored in the Figshare (https://doi.org/10.6084/m9.figshare.27996128). The code used in this study is available in the GitHub repositories at https://github.com/eriksamuelsson1/TNBC. Competing interests The authors declare that they have no competing interests. Funding This work is supported by funds from by Chan Zuckerberg Initiative; an advised fund of the Silicon Valley Community Foundation; the Erling-Persson Family Foundation (Erling-Perssons Stiftelse; the Human Developmental Cell Atlas); the Knut and Alice Wallenberg Foundation (Knut och Alice Wallenbergs Stiftelse; KAW 2018.0172); the Swedish Research Council (Vetenskapsrådet; 2019-01238); and the Swedish Cancer Society (Cancerfonden; CAN 2021/1726). Author contributions Conceptualization: E.S., M.N.; Methodology: E.S.; Software:E.S., R.F.; Validation: E.S., R.F., T.H; Formal analysis: , E.S., R.F.; Investigation: E.S.; Resources: E.S., R.F.; Data curation: E.S. R.F. ; Writing - original draft: E.S.; Writing - review & editing: E.S, T.H., M.N. ; Visualization: E.S, T.H.; Supervision: M.N., E.S.; Project administration: M.N., E.S.; Funding acquisition: M.N., E.S.. Acknowledgements This work was supported by the National Bioinformatics Infrastructure Sweden (NBIS) at SciLifeLab. The computations and data storage was enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS), partially funded by the Swedish Research Council through grant agreement no. 2022-06725. References Derakhshan F, Reis-Filho JS. Pathogenesis of Triple-Negative Breast Cancer. Annu Rev Pathol 17 , 181-204 (2022). Almansour NM. Triple-Negative Breast Cancer: A Brief Review About Epidemiology, Risk Factors, Signaling Pathways, Treatment and Role of Artificial Intelligence. Front Mol Biosci 9 , 836417 (2022). Tan PH , et al. The 2019 World Health Organization classification of tumours of the breast. Histopathology 77 , 181-185 (2020). Zhao S, Zuo WJ, Shao ZM, Jiang YZ. Molecular subtypes and precision treatment of triple-negative breast cancer. Ann Transl Med 8 , 499 (2020). Guo L , et al. Breast cancer heterogeneity and its implication in personalized precision therapy. Exp Hematol Oncol 12 , 3 (2023). Luond F, Tiede S, Christofori G. Breast cancer as an example of tumour heterogeneity and tumour cell plasticity during malignant progression. Br J Cancer 125 , 164-175 (2021). Polyak K. Heterogeneity in breast cancer. J Clin Invest 121 , 3786-3788 (2011). Rivenbark AG, O'Connor SM, Coleman WB. Molecular and cellular heterogeneity in breast cancer: challenges for personalized medicine. Am J Pathol 183 , 1113-1124 (2013). Turashvili G, Brogi E. Tumor Heterogeneity in Breast Cancer. Front Med (Lausanne) 4 , 227 (2017). Lehmann BD , et al. Identification of human triple-negative breast cancer subtypes and preclinical models for selection of targeted therapies. J Clin Invest 121 , 2750-2767 (2011). Lehmann BD , et al. Refinement of Triple-Negative Breast Cancer Molecular Subtypes: Implications for Neoadjuvant Chemotherapy Selection. PLoS One 11 , e0157368 (2016). Wang DY, Jiang Z, Ben-David Y, Woodgett JR, Zacksenhaus E. Molecular stratification within triple-negative breast cancer subtypes. Sci Rep 9 , 19107 (2019). Masuda H , et al. Reverse phase protein array identification of triple-negative breast cancer subtypes and comparison with mRNA molecular subtypes. Oncotarget 8 , 70481-70495 (2017). DiNome ML , et al. Clinicopathological Features of Triple-Negative Breast Cancer Epigenetic Subtypes. Ann Surg Oncol 26 , 3344-3353 (2019). Gong Y , et al. Metabolic-Pathway-Based Subtyping of Triple-Negative Breast Cancer Reveals Potential Therapeutic Targets. Cell Metab 33 , 51-64 e59 (2021). Ensenyat-Mendez M , et al. Current Triple-Negative Breast Cancer Subtypes: Dissecting the Most Aggressive Form of Breast Cancer. Front Oncol 11 , 681476 (2021). Turner KM, Yeo SK, Holm TM, Shaughnessy E, Guan JL. Heterogeneity within molecular subtypes of breast cancer. Am J Physiol Cell Physiol 321 , C343-C354 (2021). Beroukhim R , et al. The landscape of somatic copy-number alteration across human cancers. Nature 463 , 899-905 (2010). Hanahan D. Hallmarks of Cancer: New Dimensions. Cancer Discov 12 , 31-46 (2022). Stratton MR, Campbell PJ, Futreal PA. The cancer genome. Nature 458 , 719-724 (2009). Cai H, Kumar N, Ai N, Gupta S, Rath P, Baudis M. Progenetix: 12 years of oncogenomic data curation. Nucleic Acids Res 42 , D1055-1062 (2014). Steele CD , et al. Signatures of copy number alterations in human cancer. Nature 606 , 984-991 (2022). Cai H, Gupta S, Rath P, Ai N, Baudis M. arrayMap 2014: an updated cancer genome resource. Nucleic Acids Res 43 , D825-830 (2015). Kim TM, Xi R, Luquette LJ, Park RW, Johnson MD, Park PJ. Functional genomic analysis of chromosomal aberrations in a compendium of 8000 cancer genomes. Genome Res 23 , 217-227 (2013). Zack TI , et al. Pan-cancer patterns of somatic copy number alteration. Nat Genet 45 , 1134-1140 (2013). Curtis C , et al. The genomic and transcriptomic architecture of 2,000 breast tumours reveals novel subgroups. Nature 486 , 346-352 (2012). Shiu KK, Natrajan R, Geyer FC, Ashworth A, Reis-Filho JS. DNA amplifications in breast cancer: genotypic-phenotypic correlations. Future Oncol 6 , 967-984 (2010). Turner N , et al. Integrative molecular profiling of triple negative breast cancers identifies amplicon drivers and potential therapeutic targets. Oncogene 29 , 2013-2023 (2010). Berger AC , et al. A Comprehensive Pan-Cancer Molecular Study of Gynecologic and Breast Cancers. Cancer Cell 33 , 690-705 e699 (2018). Neftel C , et al. An Integrative Model of Cellular States, Plasticity, and Genetics for Glioblastoma. Cell 178 , 835-849 e821 (2019). Gavish A , et al. Hallmarks of transcriptional intratumour heterogeneity across a thousand tumours. Nature 618 , 598-606 (2023). De Falco A, Caruso F, Su XD, Iavarone A, Ceccarelli M. A variational algorithm to detect the clonal copy number substructure of tumors from scRNA-seq data. Nat Commun 14 , 1074 (2023). Kong D, Hughes CJ, Ford HL. Cellular Plasticity in Breast Cancer Progression and Therapy. Front Mol Biosci 7 , 72 (2020). Kvokackova B, Remsik J, Jolly MK, Soucek K. Phenotypic Heterogeneity of Triple-Negative Breast Cancer Mediated by Epithelial-Mesenchymal Plasticity. Cancers (Basel) 13 , (2021). Pujals M , et al. RAGE/SNAIL1 signaling drives epithelial-mesenchymal plasticity in metastatic triple-negative breast cancer. Oncogene 42 , 2610-2628 (2023). Heilala M , et al. Fibrin Stiffness Regulates Phenotypic Plasticity of Metastatic Breast Cancer Cells. Adv Healthc Mater 12 , e2301137 (2023). Guo Z, Han S. Targeting cancer stem cell plasticity in triple-negative breast cancer. Explor Target Antitumor Ther 4 , 1165-1181 (2023). Schade AE , et al. AKT and EZH2 inhibitors kill TNBCs by hijacking mechanisms of involution. Nature , (2024). Reed AD , et al. A single-cell atlas enables mapping of homeostatic cellular shifts in the adult human breast. Nat Genet 56 , 652-662 (2024). Burstein MD , et al. Comprehensive genomic analysis identifies novel subtypes and targets of triple-negative breast cancer. Clin Cancer Res 21 , 1688-1698 (2015). Lehmann BD , et al. Multi-omics analysis identifies therapeutic vulnerabilities in triple-negative breast cancer subtypes. Nat Commun 12 , 6276 (2021). Perou CM , et al. Molecular portraits of human breast tumours. Nature 406 , 747-752 (2000). Lehmann BD, Pietenpol JA. Identification and use of biomarkers in treatment strategies for triple-negative breast cancer subtypes. J Pathol 232 , 142-150 (2014). Zheng S , et al. Vasculogenic mimicry regulates immune infiltration and mutational status of the tumor microenvironment in breast cancer to influence tumor prognosis. Environ Toxicol 39 , 2948-2960 (2024). Liang X , et al. Identification of new subtypes of breast cancer based on vasculogenic mimicry related genes and a new model for predicting the prognosis of breast cancer. Heliyon 10 , e36565 (2024). Andonegui-Elguera MA, Alfaro-Mora Y, Caceres-Gutierrez R, Caro-Sanchez CHS, Herrera LA, Diaz-Chavez J. An Overview of Vasculogenic Mimicry in Breast Cancer. Front Oncol 10 , 220 (2020). Gao R , et al. Cancer cell immune mimicry delineates onco-immunologic modulation. iScience 24 , 103133 (2021). Timar J, Honn KV, Hendrix MJC, Marko-Varga G, Jalkanen S. Newly identified form of phenotypic plasticity of cancer: immunogenic mimicry. Cancer Metastasis Rev 42 , 323-334 (2023). Goncalves G , et al. IFNgamma Modulates the Immunopeptidome of Triple Negative Breast Cancer Cells by Enhancing and Diversifying Antigen Processing and Presentation. Front Immunol 12 , 645770 (2021). Jacobo Jacobo M, Donnella HJ, Sobti S, Kaushik S, Goga A, Bandyopadhyay S. An inflamed tumor cell subpopulation promotes chemotherapy resistance in triple negative breast cancer. Sci Rep 14 , 3694 (2024). Oshi M , et al. Inflammation Is Associated with Worse Outcome in the Whole Cohort but with Better Outcome in Triple-Negative Subtype of Breast Cancer Patients. J Immunol Res 2020 , 5618786 (2020). Karaayvaz M , et al. Unravelling subclonal heterogeneity and aggressive disease states in TNBC through single-cell RNA-seq. Nat Commun 9 , 3588 (2018). Kim C , et al. Chemoresistance Evolution in Triple-Negative Breast Cancer Delineated by Single-Cell Sequencing. Cell 173 , 879-893 e813 (2018). Mavrommati I, Johnson F, Echeverria GV, Natrajan R. Subclonal heterogeneity and evolution in breast cancer. NPJ Breast Cancer 7 , 155 (2021). Wolf FA, Angerer P, Theis FJ. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol 19 , 15 (2018). Wolock SL, Lopez R, Klein AM. Scrublet: Computational Identification of Cell Doublets in Single-Cell Transcriptomic Data. Cell Syst 8 , 281-291 e289 (2019). Korsunsky I , et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat Methods 16 , 1289-1296 (2019). Xu C , et al. Automatic cell-type harmonization and integration across Human Cell Atlas datasets. Cell 186 , 5876-5891 e5820 (2023). Dominguez Conde C , et al. Cross-tissue immune cell analysis reveals tissue-specific features in humans. Science 376 , eabl5197 (2022). Colaprico A , et al. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res 44 , e71 (2016). Tables Tables 1-8 are available in the Supplementary Files section. Additional Declarations No competing interests reported. Supplementary Files SupplementaryResults1.docx SupplementaryFigures.docx TablesFormattedForBCR.xlsx Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-5974271","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":416093383,"identity":"4b38152f-44db-46ec-bcbf-a899127d8b7c","order_by":0,"name":"Erik Reinhold Samuelsson","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAwUlEQVRIiWNgGAWjYFACxmYQKQciDjwgRYsxWEsCkdYwg4jEBhBJlBZ+6cPNxhV/6tLnhx1+CLTFTk63gYAWyb7E5sSzbYdzN95OMwBqSTY2O0BAi8EZxuaDjQ0HcjfOTgBpOZC4jZAWe5CWBqDDDGenfyBOiwEPY3NiAxtzgrx0DpG2SABtMWxsO2y4QTqn4ECCARF+4e9hfywJdJi8/Oz0zR8+VNjJEdSCcCFYpQGxykFAvoEU1aNgFIyCUTCiAAArOkWIVb1jkAAAAABJRU5ErkJggg==","orcid":"","institution":"Stockholm University","correspondingAuthor":true,"prefix":"","firstName":"Erik","middleName":"Reinhold","lastName":"Samuelsson","suffix":""},{"id":416093384,"identity":"37eba024-d9ed-49b7-8f2a-85fc64e051dc","order_by":1,"name":"Roy Francis","email":"","orcid":"","institution":"Uppsala University","correspondingAuthor":false,"prefix":"","firstName":"Roy","middleName":"","lastName":"Francis","suffix":""},{"id":416093385,"identity":"e20ae499-35a2-4eaf-a126-3be0b15515a0","order_by":2,"name":"Taobo Hu","email":"","orcid":"","institution":"Stockholm University","correspondingAuthor":false,"prefix":"","firstName":"Taobo","middleName":"","lastName":"Hu","suffix":""},{"id":416093386,"identity":"13084fa4-8562-4b23-a8f7-5387fd115855","order_by":3,"name":"Mats Nilsson","email":"","orcid":"","institution":"Stockholm University","correspondingAuthor":false,"prefix":"","firstName":"Mats","middleName":"","lastName":"Nilsson","suffix":""}],"badges":[],"createdAt":"2025-02-06 14:08:35","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-5974271/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-5974271/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":76663177,"identity":"b7a9b42a-0b5b-487c-84f6-54d539d3022d","added_by":"auto","created_at":"2025-02-19 12:14:54","extension":"jpg","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":401127,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eStudy design and computational analysis summary. \u003c/strong\u003eRepresentation of the study design and computational analysis summary. This included collection of the publically available datasets, merging and annotation, metaprogram extraction and computational analysis.\u003c/p\u003e","description":"","filename":"Figure1.tiff.jpg","url":"https://assets-eu.researchsquare.com/files/rs-5974271/v1/3be711885314c04ab7256cd6.jpg"},{"id":76663603,"identity":"bd877d1d-58b6-4486-bcf6-59bfe7baef68","added_by":"auto","created_at":"2025-02-19 12:22:54","extension":"jpg","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":925092,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eGenomic and Transcriptomic Heterogeneity in TNBC Assessed via scRNA-seq\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ea) UMAP of cells of the TNBC-MAP colored by aneuploid or diploid status.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eb) UMAP of malignant cells of the TNBC-MAP colored by donor ID.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ec) Metaprograms of ITH in Malignant TNBC Cells: Heatmap of Jaccard Similarity Indices with Cluster Annotations. \u003c/strong\u003eThis heatmap represents the Jaccard similarity indices between various NMF programs, visualized to highlight the clustering of similar programs. Each cell indicates the Jaccard similarity between the corresponding row and column programs, with darker shades representing higher similarity. Clusters identified through hierarchical clustering are annotated with distinct colors along the rows and columns. The legend indicates the cluster assignments, and the colorbar provides a scale for interpreting the similarity values.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ed) Bar plot of MPs per Subclone. \u003c/strong\u003eX axis, proportion of cells assigned to MPs. Y axis, donor subclones.\u003c/p\u003e","description":"","filename":"Figure2.tiff.jpg","url":"https://assets-eu.researchsquare.com/files/rs-5974271/v1/9807894c9663c6251ff9192b.jpg"},{"id":76661792,"identity":"65f860c0-0c40-4ee7-bd10-3b949c8a0be0","added_by":"auto","created_at":"2025-02-19 12:06:54","extension":"jpg","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":534249,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eAssessment of Metaprogram Scoring in Malignant Cells\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ea)\u003c/strong\u003e \u003cstrong\u003eProportions of Malignant Cells with Highest scoring Metaprogram. \u003c/strong\u003eThe chart illustrates the distribution of cells assigned to various MPs, highlighting that MP 6 and MP 7 represent the largest portions at 55.3% and 24.3%, respectively. The accompanying table lists the exact number of cells assigned to each MP, providing a detailed view of the data distribution.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eb) Upregulation of MP Signatures in Malignant Cells: Cumulative Gene Expression Heatmap of Metaprograms in Malignant Cells.\u003c/strong\u003e This heatmap visualizes the cumulative gene expression of metaprograms across malignant cells. Each row represents a metaprogram, while each column represents a single cell. The colors indicate the relative expression levels (log-ratio), with red shades showing higher expression and blue shades showing lower expression. Cells have been hierarchically clustered to reveal patterns in metaprogram activity across the cell population.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ec) GSEA results for TNBC metaprograms TNBC_MP6. \u003c/strong\u003eThe plot displays the top enriched pathways with positive GeneRatios and lowest adjusted \u003cem\u003ep\u003c/em\u003e-values. The pathways are plotted on the x-axis by their GeneRatio, indicating the extent of enrichment, while the size of the points reflects pathway significance.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ed) GSEA results for TNBC metaprograms TNBC_MP7. \u003c/strong\u003eThe plot displays the top enriched pathways with positive GeneRatios and lowest adjusted \u003cem\u003ep\u003c/em\u003e-values. The pathways are plotted on the x-axis by their GeneRatio, indicating the extent of enrichment, while the size of the points reflects pathway significance.\u003c/p\u003e","description":"","filename":"Figure3.tiff.jpg","url":"https://assets-eu.researchsquare.com/files/rs-5974271/v1/5eeb4b56f1dd044f2305df6d.jpg"},{"id":76661793,"identity":"4c224f0c-1348-4796-975f-770897252550","added_by":"auto","created_at":"2025-02-19 12:06:54","extension":"jpg","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":453359,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eComparison of TNBC Metaprograms Across Themselves and Bulk RNA TNBC Subtypes\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ea) The heatmap illustrates the correlation between MPs across all malignant cells in the TNBCMap. \u003c/strong\u003eThe correlation matrix in Figure 2E displays the relationships between the expression patterns of the nine metaprograms (TNBC_MP1-9) across malignant cells. The correlation values range from -1 to 1, where positive values (closer to 1) indicate a positive correlation, and negative values (closer to -1) indicate a negative correlation.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eb) Comparison of Metaprograms to Lehmann \u0026amp; Burstein Subtypes. \u003c/strong\u003eHeatmap representation of Jaccard similarity indices comparing TNBC metaprograms (TNBC_MP1–9) with molecular subtypes defined by Lehmann and Burstein. The indices highlight the degree of overlap between gene sets associated with TNBC metaprograms and subtype classifications, illustrating potential relationships between metaprograms and established subtype features. Warmer colors indicate higher similarity.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ec) Kaplan-Meier survival plots for TNBC patients grouped by high and low expression levels of malignant metaprograms. \u003c/strong\u003eThe survival analysis reveals significant differences between high- and low- expression groups, with p-values below 0.05. Low expression levels are represented by the cyan line, while high expression levels are shown in red. Underneath is a risk table displaying the number of patients per condition at each time point.\u003c/p\u003e","description":"","filename":"Figure4.jpg","url":"https://assets-eu.researchsquare.com/files/rs-5974271/v1/a522a5693a37182cce863bb0.jpg"},{"id":91944766,"identity":"164a28cf-6b5c-4354-9d58-7cfbb9102874","added_by":"auto","created_at":"2025-09-23 05:01:47","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":4156450,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-5974271/v1/81e3c306-b391-4b12-b368-20a502488bbc.pdf"},{"id":76661807,"identity":"5a0e9be9-7bb9-4391-979e-1ae424452b75","added_by":"auto","created_at":"2025-02-19 12:06:56","extension":"docx","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":470741,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryResults1.docx","url":"https://assets-eu.researchsquare.com/files/rs-5974271/v1/447402446ae5813759b19a1f.docx"},{"id":76663604,"identity":"8f177e40-1025-4e4d-a6b3-9c9c678c2040","added_by":"auto","created_at":"2025-02-19 12:22:54","extension":"docx","order_by":2,"title":"","display":"","copyAsset":false,"role":"supplement","size":4091629,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryFigures.docx","url":"https://assets-eu.researchsquare.com/files/rs-5974271/v1/8795020c0bc65a61effdc890.docx"},{"id":76661805,"identity":"228ced86-311b-4e62-af53-5115ae48a2bc","added_by":"auto","created_at":"2025-02-19 12:06:56","extension":"xlsx","order_by":3,"title":"","display":"","copyAsset":false,"role":"supplement","size":369531,"visible":true,"origin":"","legend":"","description":"","filename":"TablesFormattedForBCR.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-5974271/v1/16ab780eff83ea40ad18ccdd.xlsx"}],"financialInterests":"No competing interests reported.","formattedTitle":"An Integrative Model of Single Cell Transcriptomic States for Triple-Negative Breast Cancer","fulltext":[{"header":"Background","content":"\u003cp\u003eTriple-negative breast cancer (TNBC) is a subtype of breast cancer that lies at the forefront of difficulties due to heterogeneity and the need to develop targeted therapies. TNBC is not a biological classification, but an operational term that stemmed from the fact that patients with estrogen receptor-negative (ER-), progesterone receptor-negative (PR-), and human epidermal growth factor receptor 2-negative (HER2-negative) disease are insensitive to endocrine treatment and targeted therapies. TNBC encompasses a collection of distinct epithelial cancers that are vastly different in terms of their therapeutic response, clinical behavior, histologies, and genomic characteristics\u003csup\u003e\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e\u003c/sup\u003e. It accounts for 10\u0026ndash;20% of all breast cancers, and is an aggressive form of breast cancer with a poorer prognosis than other types, with an average survival rate of 10.2 months from diagnosis\u003csup\u003e\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e\u003c/sup\u003e. According to the WHO 2019 Classification of Breast Cancer, TNBCs can include a range of histological classifications\u003csup\u003e\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e\u003c/sup\u003e, though the most common include metaplastic carcinoma, medullary carcinoma, invasive lobular carcinoma, and apocrine carcinoma\u003csup\u003e\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e\u003c/sup\u003e. This heterogeneity increases the difficulty of diagnosis and prognosis\u003csup\u003e\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e,\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e\u003c/sup\u003e, as well as treatment as there is a high degree of diversity within and between patients\u003csup\u003e\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e,\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e,\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e,\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e\u003c/sup\u003e, as well as a high degree of treatment resistance\u003csup\u003e\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e \u003cp\u003eCompounding the issue of heterogeneity is the lack of consensus on stratifying TNBC into clinically actionable subtypes. Several efforts have been made to stratify TNBC into actionable subtypes, including transcriptomic\u003csup\u003e\u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e,\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e,\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e\u003c/sup\u003e, proteomic\u003csup\u003e\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e\u003c/sup\u003e, epigenetic\u003csup\u003e\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e\u003c/sup\u003e and metabolic\u003csup\u003e\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e\u003c/sup\u003e. While there are observed overlaps between classification systems, there is still large variability and inconsistencies between them\u003csup\u003e\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e\u003c/sup\u003e. The challenges in implementing these systems stem from the diversity of datasets, poor integration between methods, and inconsistencies in classification. A further limitation lies in the bulk nature of such signatures; for these methods a single breast tumor is classified according to the dominant subpopulation into a predefined subtype. In reality, there can exist cell-state diversity both within the tumor and between patients that is masked at the bulk level. The tumor microenvironment's (TME) complexity complicates categorization, as similar transcriptomic and genomic signatures between malignant and nonmalignant cells can confound classification and treatment decisions\u003csup\u003e\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e,\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e\u003c/sup\u003e. A more precise, single-cell-based classification of TNBC can improve prognosis and better characterize the dynamic process of disease progression driven by intratumoral heterogeneity, which is shaped by diverse malignant breast cancer phenotypes, dynamic phenotype conversion, and clonal evolution. This is crucial, as these phenomena directly contribute to the development of metastatic disease and the emergence of treatment-resistant phenotypes\u003csup\u003e\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e \u003cp\u003eTo address this issue, we have assessed a cohort of TNBC patient single-cell RNA sequencing (scRNA-seq) datasets and identified nine transcriptomic states that underlie range of behaviors and driving forces in this disease; from cancer stem cell differentiation, Epithelial-to-Mesenchymal Transition (EMT), plasticity, metabolic regulation, and stress response, to clonal evolution, selection, immunogenic mimicry, vasculogenic mimicry, and prognosis.\u003c/p\u003e"},{"header":"Methods","content":"\u003cp\u003e \u003cb\u003eReference Dataset Generation\u003c/b\u003e \u003c/p\u003e \u003cp\u003eRaw counts and associated metadata from each study were combined into a single dataset. The dataset was analyzed mostly using Scanpy\u003csup\u003e\u003cspan citationid=\"CR55\" class=\"CitationRef\"\u003e55\u003c/span\u003e\u003c/sup\u003e in Python. Cells with less than 200 detected genes and genes with less than 5 cells were discarded. QC metrics such as percentage of mitochondrial expression and ribosomal expression were computed. Data was normalized followed by log transformation and scaling. The effect of the total number of counts was regressed out. Cell cycle state was estimated. Doublets were estimated using Scrublet\u003csup\u003e\u003cspan citationid=\"CR56\" class=\"CitationRef\"\u003e56\u003c/span\u003e\u003c/sup\u003e. Top highly variable genes were identified followed by dimensionality reduction methods: PCA and UMAP. The data was integrated over 'datasets' using Harmony\u003csup\u003e\u003cspan citationid=\"CR57\" class=\"CitationRef\"\u003e57\u003c/span\u003e\u003c/sup\u003e. Cell types were estimated using Celltypist\u003csup\u003e\u003cspan citationid=\"CR58\" class=\"CitationRef\"\u003e58\u003c/span\u003e,\u003cspan citationid=\"CR59\" class=\"CitationRef\"\u003e59\u003c/span\u003e\u003c/sup\u003e. Tabula sapiens 11k (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://cellxgene.cziscience.com/e/2ba40233-8576-4dec-a5f1-2adfa115e2dc.cxg/\u003c/span\u003e\u003cspan address=\"https://cellxgene.cziscience.com/e/2ba40233-8576-4dec-a5f1-2adfa115e2dc.cxg/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003e)\u003c/span\u003e and Kumar 2023 700k (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://cellxgene.cziscience.com/collections/4195ab4c-20bd-4cd3-8b3d-65601277e731\u003c/span\u003e\u003cspan address=\"https://cellxgene.cziscience.com/collections/4195ab4c-20bd-4cd3-8b3d-65601277e731\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003e)\u003c/span\u003e were used as reference datasets. Datasets were classified into high, low and normal. High, low, and normal reference assignments made based on the site of the sample biopsy according to the study and the enrichment of cells before scRNA-seq; high was from TNBC patients with no enrichment, low was from TNBC patients with non-epithelial enrichment (usually immune) and reference normal were from non-TNBC patients. High and low datasets were used to infer copy number variants (CNV). InferCNV (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://github.com/broadinstitute/inferCNV\u003c/span\u003e\u003cspan address=\"https://github.com/broadinstitute/inferCNV\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) was run to estimate CNVs using gencode v45 as annotation. High datasets were used as reference. R package SCEVAN\u003csup\u003e\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e\u003c/sup\u003e was also run to identify CNVs using only the high samples.\u003c/p\u003e \u003cdiv id=\"Sec3\" class=\"Section2\"\u003e \u003ch2\u003eMetaprogram Generation\u003c/h2\u003e \u003cp\u003eTo identify the heterogeneous transcriptomic programs present across the samples, we have applied an approach similar to the one pursued by Gavish et al 2023\u003csup\u003e31\u003c/sup\u003e, but with the aim of characterizing TNBC tumor-type specific programs rather than pan-tumor programs. We performed consensus Non-negative Matrix factorization (cNMF, \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://github.com/dylkot/cNMF\u003c/span\u003e\u003cspan address=\"https://github.com/dylkot/cNMF\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) on the data of each individual patient to identify the programs of each sample. Since application of cNMF requires a \u0026ldquo;K\u0026rdquo; parameter that influences the results, we run cNMF using different values (K\u0026thinsp;=\u0026thinsp;5,6,7,8,9,10), and generate 45 programs for each tumor. Each cNMF program is summarized by the top 100 genes representing that program based on cNMF coefficients. Then, we identified the most robust cNMF programs across the patient cohort as those that recur within the tumor (gene lists have at least 70% overlap), recur across tumor (have at least 20% overlap with any other cNMF program in other patients analyzed), and are then non-redundant within the tumor (rank programs by similarity with programs from other tumors, remove programs that have at least a 20% overlap with other programs within the same patient). From the robust cNMF programs, we clustered them together based on their Jaccard similarity, and identified the most consistent genes present across the grouped cNMF programs. This is how we end up with our meta-programs. The clustering process resulted in 30 metaprograms (clusters), each containing samples with high internal similarity.\u003c/p\u003e \u003c/div\u003e\n\u003ch3\u003eAssessment of Clustered cNMF's Comprising Each Metaprogram\u003c/h3\u003e\n\u003cp\u003eTo validate the resulting metaprograms, we assessed the overall similarity via Jaccard indices of each cNMF assigned to a cluster and their tendency to be more similar to cNMFs within said cluster than outside of it. The Jaccard index, defined as the size of the intersection divided by the size of the union of two gene sets, was computed using the formula:\u003cdiv id=\"Equa\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equa\" name=\"EquationSource\"\u003e\n$$\\:Jaccard\\:Index\\:=\\:\\frac{\\left|A\\:\\cap\\:\\:B\\right|}{\\left|A\\:\\cup\\:\\:B\\right|}$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003cp\u003e, where A and B are two gene sets. First, we calculated the Jaccard similarity indices for all pairs of cNMFs within a dataset, resulting in a similarity matrix. This matrix represents the proportion of shared features between pairs of samples. Next, we grouped the cNMFs together based on the metaprogram clustering which grouped the samples into distinct clusters based on their similarity indices. The similarity matrix could then be visualized as a heatmap to provide a clear visual representation of the internal similarity of the cNMFs within clusters and the separation between different clusters.\u003c/p\u003e\n\u003ch3\u003eCumulative Expression of Metaprograms\u003c/h3\u003e\n\u003cp\u003eMalignant cells from the TNBCMap were identified based on the iCNV attribute, specifically selecting those labeled as 'aneuploid'. Metaprograms were defined by the gene sets from each of the 9 metaprograms \u003cb\u003e(Supplementary Table\u0026nbsp;3\u003c/b\u003e). The gene expression data for the malignant cells was normalized using StandardScaler and log-transformed to emphasize relative differences in expression levels. For each metaprogram, cumulative expression was calculated by summing the expression levels of its constituent genes. To quantify the relative expression score for each cell, the cumulative expression values were centered around zero, with positive values indicating higher expression and negative values indicating lower expression relative to the mean. Hierarchical clustering was performed on the cells using the average linkage method, and the resulting order was applied to the cumulative expression data.\u003c/p\u003e\n\u003ch3\u003eScoring of Metaprograms in Cells\u003c/h3\u003e\n\u003cp\u003eTo analyze the distribution of cells across different metaprograms, we used scanpy's tl.score_genes function to calculate gene scores for each metaprogram, using the list of genes present within the data.Each cell was then assigned to the MP with the highest score.\u003c/p\u003e\n\u003ch3\u003eCorrelation of Metaprograms Across Cells:\u003c/h3\u003e\n\u003cp\u003eThe generation of the correlation matrix involved the dataset utilized for this analysis comprising metaprogram (MP) scores for various cells. The correlation coefficients, specifically Pearson correlation coefficients, were calculated using the .corr() method in pandas.\u003c/p\u003e \u003cdiv id=\"Sec8\" class=\"Section2\"\u003e \u003ch2\u003eComparison of TNBC Metaprograms to Reference Metaprograms\u003c/h2\u003e \u003cp\u003eThe analysis involved three datasets: the nonmalignant TNBC cells, normal reference breast cells (Reed 2024), and Gavish 2023 Hallmarks of Cancer metaprograms. To determine the similarities between these gene sets, the Jaccard similarity index was calculated for each pair of columns (representing gene sets) from the metaprograms dataset and each of the nonmalignant TNBC and Reed2024 datasets.\u003c/p\u003e \u003cp\u003eTo extract the overlapping genes between the metaprograms, we then developed a custom function to calculate the overlap between genes in each pair of MPs across different datasets. This function iteratively compares each MP from one dataset with every MP from another, identifying common genes and recording the number and identities of overlapping genes. To ensure clarity, self-comparisons within the same dataset were excluded, and only unique comparisons were considered to avoid redundancy.\u003c/p\u003e \u003cp\u003eFor each comparison, we generated two DataFrames: one to capture the count of overlapping genes and another to store the actual overlapping gene names. These DataFrames were then used to filter overlaps with five or more genes, which were considered significant for further analysis. The results were saved into CSV files for documentation and further review. Finally, the significant overlaps were formatted and printed to provide a comprehensive view of the gene overlaps across different metaprogram datasets.\u003c/p\u003e \u003c/div\u003e\n\u003ch3\u003eEnrichment Assessment of TNBC Metaprograms to Lehmann and Burstein Subtypes\u003c/h3\u003e\n\u003cp\u003eTo assess the enrichment of the metaprograms within defined breast cancer subtypes, we conducted an enrichment analysis using a custom Python script. We imported gene expression datasets for TNBC metaprograms, Lehmann and Burstein subtypes, nonmalignant TNBC clusters, and Reed 2024 clusters and developed a function to calculate gene set enrichment by identifying shared genes between each metaprogram and subtype. This function performed set intersections of non-null values in each column pair, quantifying both the count and specific list of enriched genes. To focus on meaningful enrichment results, we filtered pairs with five or more shared genes using a threshold-based function. We further evaluated the statistical significance of the enrichment for each gene set using Fisher's exact test, constructing contingency tables that included the size of each metaprogram (assumed to contain 50 genes), the size of the subtype (predefined per subtype), and an estimated total of 20,000 protein-coding genes in the human genome. These tables were processed with the fisher_exact function from the SciPy library, yielding p-values that quantify the likelihood of observing such enrichment by chance.\u003c/p\u003e\n\u003ch3\u003eCNV Assessment in TNBC Patient Data\u003c/h3\u003e\n\u003cp\u003eTo assess copy number variations (CNVs) across patient samples, we utilized output data from SCEVAN. We began by loading the CNV data from a pre-processed pickle file and defined a list of patient IDs alongside corresponding sample keys. These mappings enabled accurate cross-referencing between patient data in the CNV file and additional clinical annotations. For each patient, we extracted CNV data based on their specific sample key, which allowed for the identification and categorization of subclonal CNV patterns. The CNV heatmaps were generated with a custom function that first filtered the CNV DataFrame to retain only relevant values, excluding columns with chromosomal coordinates and gene annotations. To provide additional clarity, chromosome labels were added at the midpoint of each boundary.\u003c/p\u003e \u003cdiv id=\"Sec11\" class=\"Section2\"\u003e \u003ch2\u003eAssessment of Metaprogram Residuals Across Patient and Subclones\u003c/h2\u003e \u003cp\u003eTo investigate the association between chromosomal copy number variations (CNVs) and metaprogram (MP) assignments in individual patients, we calculated residuals using a chi-squared test on a patient-by-patient basis. We began by identifying each unique patient in our dataset based on the donor_id field from single-cell RNA sequencing (scRNA-seq) data annotations. For each patient, we grouped their data by subclone and assigned MP, constructing a contingency table that captured the distribution of MP assignments within each subclone. To prevent issues arising from zero values in the table, a small constant (0.5) was added to each cell. A chi-squared test of independence was performed on each patient's contingency table to determine if a significant association existed between subclone and MP assignment. This test provided the observed chi-squared statistic, degrees of freedom, and p-value, indicating if the observed MP distribution within subclones significantly differed from the expected distribution. Following this, we calculated residuals by subtracting the expected counts from the observed counts in each cell. These residuals reflect the extent to which each MP was over- or under-represented within subclones, with positive values indicating over-representation and negative values indicating under-representation.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec12\" class=\"Section2\"\u003e \u003ch2\u003eTCGA BRCA Analysis\u003c/h2\u003e \u003cp\u003eThis study utilized publicly available gene expression and clinical data from The Cancer Genome Atlas (TCGA) for breast invasive carcinoma (BRCA) to perform a survival analysis focused on triple-negative breast cancer (TNBC) samples. The methodology involved several steps, including data acquisition, processing, and statistical analysis. The analysis involved retrieving, processing, and analyzing RNA-Seq gene expression data, combined with clinical data, to generate survival curves for metaprograms identified in TNBC.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec13\" class=\"Section2\"\u003e \u003ch2\u003eCox Proportional Hazard Regression Analysis\u003c/h2\u003e \u003cp\u003eGene expression data for the entire TCGA-BRCA cohort was obtained using the TCGAbiolinks package. The query targeted the \"Transcriptome Profiling\" data category, specifically RNA-Seq data processed using the STAR workflow, focusing on \"Gene Expression Quantification.\" Data was downloaded and prepared using the GDCprepare function to obtain the FPKM-UQ (Fragments Per Kilobase of transcript per Million mapped reads - Upper Quartile normalized) gene expression matrix. Clinical data corresponding to the TCGA-BRCA cohort was retrieved using the GDCquery_clinic function, and relevant survival information (patient barcode, days to last follow-up, and vital status) was selected for further analysis. TNBC samples were identified using an external reclassification dataset provided in an Excel file, which contained barcodes of TNBC patients. These barcodes were extracted and cleaned to match the format of the barcodes in the gene expression matrix. The cleaned barcodes were then used to subset both the gene expression data and the clinical data, ensuring that only TNBC samples with corresponding expression data were included in the subsequent analysis. Univariate Cox proportional hazards regression was performed for each gene in the matched expression data using the coxph function from the survival package. The analysis evaluated the association between gene expression levels and patient survival (time to last follow-up and vital status). The hazard ratio (HR), 95% confidence intervals, and p-values for each gene were calculated and stored in a results table. The significant results (\u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.05) were filtered and merged with metaprogram gene lists derived from previous studies, categorizing genes into specific metaprograms associated with malignant, nonmalignant, and reference normal breast tissue samples. The final significant results were saved for further analysis and interpretation. The significant genes identified through Cox regression were used to filter and refine metaprogram gene lists. These filtered gene lists were restructured to exclude non-significant genes and saved in CSV format for further exploration.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec14\" class=\"Section2\"\u003e \u003ch2\u003eKaplan-Meier Survival Curve Analysis\u003c/h2\u003e \u003cp\u003eGene lists representing various metaprograms were loaded from a pre-filtered CSV file that contained genes previously associated with specific biological processes in TNBC. Clinical data for the TCGA-BRCA cohort was retrieved using the GDCquery_clinic function, providing vital status and survival information for each patient. Reclassification data from an external Excel file was used to identify and extract TNBC-specific patient barcodes, which were then matched with the clinical data. Necessary columns from the clinical data were selected and processed, including converting time-related variables to numeric format and categorizing patient vital status as \"Alive\" or \"Dead.\" RNA-Seq gene expression data for the TCGA-BRCA cohort was queried and downloaded using the TCGAbiolinks package\u003csup\u003e\u003cspan citationid=\"CR60\" class=\"CitationRef\"\u003e60\u003c/span\u003e\u003c/sup\u003e. This data was processed to generate a gene expression matrix, focusing on \"Primary Tumor\" samples. The matrix was normalized using variance-stabilizing transformation (VST) from the DESeq2 package, ensuring that downstream analyses were conducted on homogeneously scaled data. For each metaprogram, the expression of genes within the metaprogram was extracted from the normalized gene expression matrix. The expression data for these genes was combined and averaged to compute a \"module expression\" score for each patient. This score was then stratified into \"LOW\" or \"HIGH\" groups based on median expression levels. The clinical and gene expression data were merged to facilitate survival analysis. The merged dataset included the calculated module expression score and clinical variables such as overall survival time and vital status. A Cox proportional hazards model was fitted to this dataset to evaluate the impact of each metaprogram's expression on patient survival. Survival curves were generated using the Kaplan-Meier method, with significance evaluated via log-rank tests. The survfit function from the survival package was employed to fit the survival curves, and the ggsurvplot function from the survminer package was used for visualization. The final survival plots were saved to provide a graphical summary of the association between metaprogram expression and patient survival outcomes.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec15\" class=\"Section2\"\u003e \u003ch2\u003eComparison of Scored Metaprograms to TCGA PAM50 Subtype Classification\u003c/h2\u003e \u003cp\u003eWe began by acquiring clinical and RNA-Seq gene expression data from the TCGA-BRCA cohort using the TCGAbiolinks package. Clinical data were retrieved through the GDCquery_clinic function, while RNA-Seq data were accessed via the GDCquery function, specifically focusing on \"Transcriptome Profiling\" with the \"STAR - Counts\" workflow. Additionally, we extracted patient barcode data from an external Excel file that contained reclassified sample information, specifically identifying TNBC cases. These barcodes were carefully matched with the clinical dataset to ensure accurate identification and alignment of TNBC samples. Gene expression data was downloaded and prepared using the GDCprepare function, followed by normalization through variance-stabilizing transformation (VST) using the DESeq2 package. To enhance data robustness, genes with low expression (fewer than 10 reads across all samples) were filtered out. Metaprogram-specific gene lists were then loaded from a predefined CSV file, and the normalized gene expression matrix was subset accordingly. For each metaprogram, cumulative expression scores were calculated by averaging the expression levels of all constituent genes. We further integrated the gene expression data with clinical and subtype information. Subtype data for the TCGA-BRCA cohort was retrieved using the TCGAquery_subtype function and merged with the clinical data and metaprogram scores. Relevant clinical variables, such as vital status and survival times, were selected and combined with the metaprogram expression scores and subtype information, resulting in a comprehensive dataset for survival analysis. The final dataset, which included metaprogram scores, clinical data, and subtype information, was saved as a CSV file. This dataset served as the foundation for subsequent survival analysis, providing key insights into the potential biological roles of metaprograms in TNBC and their association with patient outcomes.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec16\" class=\"Section2\"\u003e \u003ch2\u003eHierarchical Clustering of TCGA Patients by Metaprogram Scores\u003c/h2\u003e \u003cp\u003eFrom the preceding analysis, clinical data and RNA-Seq gene expression data were obtained from the TCGA-BRCA project using the TCGAbiolinks package. Clinical data, including patient survival information and subtype classification, were merged with gene expression data, which had been variance-stabilizing transformed using the DESeq2 package. The combined dataset included MP scores for each TNBC patient, which were used for further analysis. To investigate the potential stratification of patients based on their MP scores, hierarchical clustering was performed on TNBC patients classified as having the Basal subtype using the Ward's method. The MP scores were standardized before clustering, and a specific cutoff distance was applied to define clusters. The resulting clusters were then merged back into the main dataset to associate cluster labels with patient survival data. Kaplan-Meier survival curves were generated for each identified cluster, and log-rank tests were conducted to assess the statistical significance of survival differences between clusters. The Kruskal-Wallis test was used to evaluate the differences in MP scores across clusters. Significant cluster comparisons were identified and further analyzed to calculate the mean and median survival times for each cluster, and these results were saved for further interpretation.\u003c/p\u003e \u003c/div\u003e"},{"header":"Results","content":"\u003cdiv id=\"Sec18\" class=\"Section2\"\u003e \u003ch2\u003eConstruction of the TNBC Map\u003c/h2\u003e \u003cp\u003eGiven the high heterogeneity of TNBC, it is necessary to combine the datasets from a representative population of TNBC patients in order to establish common transcriptomic definitions of the malignant phenotypes across a sufficient number of patients’ scRNA-seq samples. We identified 9 publicly available TNBC scRNA-seq datasets that include data derived from TNBC patient biopsies or reference healthy breast tissue and included more than 500 cells per donor (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e, \u003cb\u003eSupplementary Table\u0026nbsp;1, Supplementary Table\u0026nbsp;2\u003c/b\u003e). After integrating the datasets, the TNBC Map consisted of 302,509 cells with 15,097 genes from 178 donors (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ea).\u003c/p\u003e \u003cp\u003eTo identify malignant cells in the TNBC Map, we performed inferred copy number variation (CNV) assignment on the confirmed TNBC patient biopsy donors via inferCNV (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://github.com/broadinstitute/inferCNV\u003c/span\u003e\u003cspan address=\"https://github.com/broadinstitute/inferCNV\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e). Aneuploidy is common in human cancers, with their resultant extensive genome-wide CNVs in chromosomal arms being a universal feature of human cancer\u003csup\u003e\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e,\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e,\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e\u003c/sup\u003e. We identified a total of 57,183 malignant cells from 46 donors. After excluding donors with fewer than 10 identified malignant cells, we were left with 57,104 malignant cells from 21 patients. As seen in Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ea, the aneuploid and diploid cells form largely distinct clusters along with their distinct CNV profiles. Some overlap between clusters suggests regions where gene expression profiles might share common features or transitional states (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eb\u003cb\u003e)\u003c/b\u003e.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec19\" class=\"Section2\"\u003e \u003ch2\u003eGenomic and Transcriptomic Heterogeneity within TNBC\u003c/h2\u003e \u003cp\u003eTo start, we assessed the CNV patterns present in each patient to assess aneuploidy. Not only can different cancers adopt different CNV patterns\u003csup\u003e\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e,\u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e,\u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e\u003c/sup\u003e, but CNV can be used to assess patient prognosis and therapeutic resistance\u003csup\u003e\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e,\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e\u003c/sup\u003e. In general, CNVs that have been consistently observed in TNBC are gains of 1q, 8q, and 10q, with losses of 5q and 8p\u003csup\u003e\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e,\u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e,\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e,\u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e,\u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e\u003c/sup\u003e. Significantly, TNBC in and of itself does not have any CNVs used to classify it clinically\u003csup\u003e\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e\u003c/sup\u003e. Examination of the CNV patterns within patients reveals considerable heterogeneity in the gains and losses within donors (\u003cb\u003eSupplementary Fig.\u0026nbsp;1\u003c/b\u003e). The variability in CNV patterns observed across patients underscores the genomic complexity and heterogeneity within TNBC, and the fact that TNBC itself is a collection of different types of tumor.\u003c/p\u003e \u003cp\u003eTo better structure and facilitate interpretation of transcriptional ITH, we conducted an analysis of the metaprograms using the identified malignant TNBC cells. A metaprogram refers to a higher-level collection of gene expression patterns that capture more complex, often broader biological states or processes. Metaprograms represent overarching themes or states of cellular activity that can encompass multiple pathways and broader regulatory networks. Each metaprogram represents a group of genes that are co-expressed and potentially co-regulated, reflecting specific biological processes or cellular states. The identification of metaprograms has been used effectively to characterize transcriptomic states in tumors before, either in the case of Neftel \u003cem\u003eet al\u003c/em\u003e 2019 which identified 6 metaprograms in glioblastoma\u003csup\u003e\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e\u003c/sup\u003e, or Gavish \u003cem\u003eet al\u003c/em\u003e 2023 which identified hallmark metaprograms present across many types of cancer\u003csup\u003e\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e\u003c/sup\u003e. The heatmap in Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ec reveals several distinct clusters of Non-negative Matrix factorization (NMF) programs, each represented by blocks along the diagonal of the matrix. Given that these metaprograms were identified specifically from malignant, aneuploid cells in scRNA-seq data, they reflect the intrinsic properties of the cancer cells themselves (see \u003cb\u003eSupplementary Table\u0026nbsp;3\u003c/b\u003e for the genes in each metaprogram).\u003c/p\u003e \u003cp\u003eTo assess the functions indicated by each metaprogram, we first assigned each malignant cell into a metaprogram based on the highest scoring metaprogram per cell. We then extracted the top genes by log fold expression in the population of cells assigned to each metaprogram, and used these genes to perform gene set enrichment analysis (GSEA). We extracted the most significant (\u003cem\u003ep\u003c/em\u003e-value \u0026lt; 0.05) upregulated pathways and functions to assess the functions most associated with each malignant metaprogram (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e, \u003cb\u003eSupplementary Fig.\u0026nbsp;2\u003c/b\u003e). TNBC_MP1, which we term Immune-like 1 (IM-like 1), is characterized by immune-related functions, primarily involving macrophage and monocyte activation, inflammatory signaling, cytokine production, and interactions within the tumor microenvironment, with a notable role in immune checkpoint pathways such as PD-1, suggesting immune suppression and potential responsiveness to immunotherapy. TNBC_MP2, which we term Proliferating, is functionally associated with high cell proliferation, DNA repair, and aggressive cancer traits, characterized by pathways involved in cell cycle regulation, mitotic control, and poor prognosis, suggesting a metaprogram that drives rapid tumor growth and potential metastatic behavior in TNBC. TNBC_MP3 is associated with metabolic activity and sees high expression across most malignant cells in our cohort. This suggests that TNBC_MP3 is associated with high metabolic activity, which is often necessary for supporting the rapid growth and energy demands of cancer cells. TNBC_MP4, which we term mesenchymal-like (MES-like), is enriched for genes involved in ECM remodeling, mesenchymal and stromal characteristics, consistent with roles in EMT and stromal interactions, particularly resembling signatures seen in certain breast cancer subtypes (like metaplastic carcinoma) and developmental stromal programs. TNBC_MP5, which we term IM-like 2 (IM-like 2) likely represents an immune-related metaprogram characterized by T-cell and NK-cell activation, cytotoxic functions, and checkpoint pathways, suggesting an immune-active TME focused on adaptive immunity and inflammation. TNBC_MP6 likely represents a hypoxia-driven, basal-like phenotype characterized by EMT, metabolic reprogramming, and pro-inflammatory signaling associated with aggressive and invasive tumor behavior (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ec\u003cb\u003e)\u003c/b\u003e. The TNBC_MP7 metaprogram demonstrates luminal-like characteristics, with significant enrichment in luminal-associated pathways and a lack of alignment with basal or mesenchymal traits (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ed\u003cb\u003e)\u003c/b\u003e. TNBC_MP8, which we term immune-like 3 (IM-like 3), is primarily associated with immune and interferon responses, including antiviral defense, inflammatory signaling, and potential roles in EMT and immune cell recruitment, which may support immune modulation and cancer cell invasion. TNBC_MP9, which we term vascular mimicry-like (VM-like) is primarily associated with endothelial and vascular functions, including angiogenesis, vascular integrity, extracellular matrix organization, and potentially facilitating tumor-associated vascular mimicry to support tumor growth and metastasis. The absence of significant GSEA results may suggest that TNBC_MP9 could reflect a vascular phenotype in TNBC that does not closely align with common reference datasets, indicating its specific relevance to tumor vasculature and possibly niche-specific endothelial roles. These results indicate that while malignant cells in TNBC exhibit significant genomic and transcriptomic heterogeneity, it is possible to identify consistent states present across this cohort. The distinct metaprograms—ranging from immune-modulating to proliferative, metabolic, and vascular phenotypes—highlight the diverse functional landscapes within TNBC tumors, each potentially contributing differently to tumor progression, treatment response, and immune evasion.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec20\" class=\"Section2\"\u003e \u003ch2\u003eCharacterizing ITH in TNBC: The Relationship between Malignant cells, Subclones and Metaprograms\u003c/h2\u003e \u003cp\u003eOur next step, after identifying a range of metaprograms present in our TNBC cohort, was to assess the distribution of scores across all malignant cells in our cohort (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e). We did this via scoring all individual malignant cells for each metaprogram and assessing their distribution across various conditions. To start, we scored each cell for eight metaprograms (MPs), excluding TNBC_MP3 (metabolic activity), which showed high expression across all cells and were thus not informative for distinguishing cellular states. We then assigned each cell to the MP for which it scored the highest. This approach is specifically to determine which metaprograms are most prominent within malignant cells and provide a clear view of the dominant biological processes in play. The majority of cells were assigned to TNBC_MP6 (Basal-like) and TNBC_MP7 (Luminal-like), with 31,621 and 13,894 cells, respectively (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ea). This suggests that these two metaprograms are key states driving the malignant phenotype in TNBC. In contrast, TNBC_MP8 (IM-like 3, interferon signaling) and TNBC_MP5 (IM-like 2, cytotoxic immune activity) were assigned to fewer cells, suggesting that immune-related processes are less prevalent within the bulk malignant population. Similarly, TNBC_MP4 (MES-like) exhibited moderate representation, underscoring its role in tumor-stroma interactions, albeit in a less dominant manner across the malignant cells. Notably, TNBC_MP1 (IM-like 1, immune mimicry) and TNBC_MP9 (vasculogenic mimicry) were found in the fewest cells, indicating that these processes might be niche-specific rather than widespread drivers of malignant activity. These patterns suggest that while immune evasion, immune mimicry, and vasculogenic mimicry are present, they may be more localized to specific niches, such as perivascular regions or tumor margins, rather than predominating within the tumor core. Overall, these results underscore the central role of epithelial lineage development-like programs in TNBC progression, while also reflecting the complexity and heterogeneity of the tumor microenvironment.\u003c/p\u003e \u003cp\u003eTo assess whether the metaprograms could reflect genetic subclones in the malignant cells, we used SCEVAN\u003csup\u003e\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e\u003c/sup\u003e to identify distinct subclones within patient donor malignant cells based on their CNV profile. From the 19 patients that passed quality control for SCEVAN subclone analysis, we identified a total of 102 subclones, with 1–7 distinct subclones in each tumor. Next, we assigned each malignant cell into a metaprogram based on the highest scoring metaprogram per cell. Notably, each of the subclones contained cells in multiple metaprograms (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ed). Out of the 19 donors, 8 donors contain cells in all assigned metaprograms.\u003c/p\u003e \u003cp\u003eTo further investigate whether the MPs are associated with specific genetic subclones, we assessed the residuals in MP-assigned cells across these subclones. Residuals, which represent the difference between the observed and expected counts of cells assigned to each metaprogram within each subclone, allow us to quantify patterns of over- or under-representation of transcriptional states in the context of genetic subclones. Expected counts are determined by assuming a random distribution of metaprogram assignments across all subclones within each patient, based on the overall metaprogram proportions of that patient (\u003cb\u003eSupplementary Fig.\u0026nbsp;3\u003c/b\u003e). There are some trends in the representation of metaprograms within specific subclones that underscore the substantial degree of intratumoral heterogeneity across patients, with metaprograms activity varying across subclones. We observed that both TNBC_MP6 (basal-like) and TNBC_MP7 (luminal-like) states are commonly present within the same patient. Moreover, the basal-like and luminal-like states frequently display an inverse relationship comparing sub-clones within the same patient, with one being overrepresented where the other is underrepresented, hinting at distinct functional roles or selective pressures acting on these metaprograms, causing switches between these two programs. This detailed residual analysis provides insights into patient-specific subclonal evolution, highlighting key metaprograms that could influence tumor behavior and inform targeted treatment responses on a subclone-specific basis.\u003c/p\u003e \u003cp\u003eBreast cancer, including TNBC, is known for its phenotypic plasticity\u003csup\u003e\u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e33\u003c/span\u003e,\u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e\u003c/sup\u003e, allowing malignant cells to shift states as the tumor develops and encounters different conditions\u003csup\u003e\u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e35\u003c/span\u003e,\u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e36\u003c/span\u003e,\u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e37\u003c/span\u003e,\u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e38\u003c/span\u003e\u003c/sup\u003e. While our work identifies discrete transcriptomic states, TNBC operates as a continuum, with plasticity enabling multiple metaprogram expressions within individual cells. To gauge the expression of identified metaprograms, Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eb illustrates the cumulative gene expression of TNBC_MP1 to TNBC_MP9 in malignant cells. The heatmap reveals substantial heterogeneity, with distinct clusters of cells expressing specific metaprograms. For instance, TNBC_MP2, linked to cell cycle progression, is strongly expressed in a proliferative cell cluster, while TNBC_MP6 (basal-like, EMT/hypoxia) and TNBC_MP9 (VM-like, angiogenesis) are elevated in separate clusters, indicating specialized roles. Some metaprograms, like TNBC_MP3 (metabolic) and TNBC_MP4 (MES-like, ECM remodeling), show broader expression across cells, suggesting general importance in tumor maintenance. Co-expression patterns, such as TNBC_MP2 (cell cycle) and TNBC_MP7 (luminal-like), indicate that certain subpopulations, like luminal-like cells, are also engaged in robust proliferation. Additionally, TNBC_MP8 (IM-like 3, interferon signaling) appears functionally segregated, marking a subset potentially involved in immune response. Overall, the distinct patterns across metaprograms highlight the diversity of malignant cell states, reflecting common functional responses or adaptations to environmental cues.\u003c/p\u003e \u003cp\u003eTo identify metaprograms with similar correlation patterns, we performed hierarchical clustering on a correlation matrix of the metaprogram scores in the malignant cells of the TNBC Map (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ea), revealing how different biological processes are coordinated or segregated within the tumor. TNBC_MP2, associated with cell cycle progression, shows minimal correlation with other metaprograms, suggesting that proliferative cells operate independently from immune or metabolic functions. In contrast, TNBC_MP9 (VM-like), linked to angiogenesis, has a strong positive correlation with TNBC_MP4 (MES-like, ECM remodeling), indicating that cells involved in vascular development are also engaged in modifying the extracellular matrix, supporting structural changes that facilitate tumor invasion. TNBC_MP1 (IM-like 1, immune mimicry) is strongly correlated with TNBC_MP5 (IM-like 2, cytotoxic activity) and TNBC_MP8 (interferon signaling), suggesting that immune-mimicking cells participate in cytotoxic and inflammatory responses, potentially enhancing immune modulation within the tumor. Notably, TNBC_MP7, characterized by luminal-like gene expression and associated with cell adhesion, shows negative correlations with nearly all other metaprograms, suggesting it represents a distinct subset focused on structural roles. Overall, these correlation patterns highlight the functional specialization among different malignant cell subsets, with proliferative, immune-modulatory, and structural roles operating in tandem to support the tumor’s complexity and adaptability.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eWe also assessed the correlation patterns patient-by-patient (\u003cb\u003eSupplementary Results 1: Sample Specific MP correlations\u003c/b\u003e) to identify cohort and patient-specific patterns. The sample-by-sample comparison to the overall correlation matrix reveals key trends in the relationships between MPs across the dataset. Some MP pairs, such as TNBC_MP1 and TNBC_MP5 (IM-like 1 and 2) exhibit stable positive correlations across both individual samples and the overall dataset. These consistent relationships suggest that the co-activation of immune-related processes are fundamental aspects of tumor biology across the cohort. Similarly, the stable negative correlation between TNBC_MP7 (Luminal-like) and TNBC_MP3 (metabolic) implies distinct cellular states that may rarely overlap. In contrast, other MP pairs show more variability in their correlations between individual samples and the overall dataset. For instance, the correlation between TNBC_MP8 (IM-like 3) with TNBC_MP6 (basal-like) and TNBC_MP9 (VM-like) varies across samples, indicating that their co-activation is more context-dependent, potentially influenced by specific tumor characteristics or environmental factors. Additionally, while the overall correlation matrix tends to smooth out variations, closer inspection of individual samples reveals important differences that highlight tumor heterogeneity and patient-specific differences. These insights underscore the importance of considering both overall trends and sample-specific variations to fully understand the complexity and diversity within TNBC tumors, which may have significant implications for therapeutic strategies.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec21\" class=\"Section2\"\u003e \u003ch2\u003eCharacterization of TNBC Metaprograms via Comparison to Other Data Sets\u003c/h2\u003e \u003cp\u003eAs there are parallels between normal breast and malignant development in breast cancer, we conducted a comprehensive analysis using gene expression data from various studies to characterize the range of shared features between our derived TNBC MPs and MPs derived from other relevant data sets. Specifically, we examined the overlap among TNBC_MP1-9 derived from malignant cells in the TNBCMap, MPs from nonmalignant cells from TNBC patients in the TNBCMap, MPs from the Reed et al. 2024 atlas of breast cells\u003csup\u003e\u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e39\u003c/span\u003e\u003c/sup\u003e, and MPs from the Gavish et al. 2023 study on biological Hallmarks across many varieties of cancer\u003csup\u003e\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e\u003c/sup\u003e (\u003cb\u003eSupplementary Table\u0026nbsp;3\u003c/b\u003e). The heatmap (\u003cb\u003eSupplementary Fig.\u0026nbsp;7\u003c/b\u003e) of the Jaccard similarity indices between gene sets of various metaprograms and clusters from the nonmalignant, Reed2024 and Hallmark datasets reveals shared features between normal breast programs, general cancer programs and our TNBC metaprograms. Several key observations can be made from the heatmap. TNBC_MP1, TNBC_MP5, TNBC_MP7, and TNBC_MP8 exhibit low similarity but some overlap with specific clusters, indicating connections to immune response (TNBC_MP1), immune signaling (TNBC_MP5), luminal cells (TNBC_MP7), and inflammation (TNBC_MP8). TNBC_MP2, TNBC_MP3, and TNBC_MP6 display moderate overlap with clusters related to cell cycle, metabolism, and basal cells, reflecting shared processes across malignant and general cancer contexts. TNBC_MP9 shows high similarity with Reed2024_MP4 and NC9, highlighting the similarity of vascular processes in malignant cells with nonmalignant counterparts. These observations are further strengthened when assessing the overlapping genes. See \u003cb\u003eSupplementary Results 2\u003c/b\u003e for a full list of overlapping genes with the TNBC Metaprograms.\u003c/p\u003e \u003cp\u003eNext, we validated the single-cell metaprograms against the bulk RNA subtypes found in the Burstein\u003csup\u003e\u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e40\u003c/span\u003e\u003c/sup\u003e and Lehmann\u003csup\u003e\u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e,\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e\u003c/sup\u003e. Given that bulk RNA expression profiles represent mixtures of diverse malignant and nonmalignant cell types and states, the expression of each metaprogram in bulk samples provides an estimate of the abundance of the corresponding cellular state within the TME. The analysis comparing TNBC metaprograms derived from single-cell RNA sequencing with TNBC subtypes defined by bulk RNA sequencing data reveals key insights into the relationships between cellular states and the tumor microenvironment (TME) across different TNBC subtypes. Given that bulk RNA expression profiles represent mixtures of diverse malignant and nonmalignant cell types and states, the expression of each metaprogram in bulk samples provides an estimate of the abundance of the corresponding cellular state within the TME (see \u003cb\u003eSupplementary Table\u0026nbsp;5\u003c/b\u003e for extracted Lehmann and Burstein subtype genes). In Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eb, the heatmap of the Jaccard similarity indices between gene sets of various metaprograms and bulk RNA subtypes reveals shared features between our single-cell and the established bulk RNA signatures. Several key observations can be made from the heatmap. It is likely that the bulk subtypes masked the malignant MPs that are most similar to immune and stroma components. TNBC_MP3 and TNBC_MP7 show minimal overlap with TNBC subtypes, suggesting they capture rare or specialized cell states underrepresented in bulk RNA data. TNBC_MP2 overlaps with Lehmann subtypes related to cell cycle, while TNBC_MP1, 5 and 8 overlap with immune-modulatory subtypes. TNBC_MP4 and 6 align with basal-like and mesenchymal features in Lehmann subtypes. TNBC_MP8 overlaps with Burstein subtypes, reflecting immune and mesenchymal cell states.\u003c/p\u003e \u003cp\u003eThese findings highlight the importance of integrating single-cell and bulk RNA data to better understand the tumor microenvironment and its impact on TNBC subtype biology, and the enrichment analysis reveals key overlaps between TNBC MPs and cancer-related processes, emphasizing the diverse cellular states contributing to TNBC heterogeneity. MPs like TNBC_MP1 and TNBC_MP4 capture immune and mesenchymal features critical to TNBC subtype classification, while others, like TNBC_MP6 and TNBC_MP9, highlight distinct biological processes such as EMT and angiogenesis essential for tumor progression.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec22\" class=\"Section2\"\u003e \u003ch2\u003eMalignant and Nonmalignant Metaprograms Influence TNBC Patient Survival\u003c/h2\u003e \u003cp\u003eTo systematically examine the association between transcriptomic states and patient outcome, we next performed an analysis of 192 identified TNBC bulk specimens from Cancer Genome Atlas Breast Invasive Carcinoma (TCGA-BRCA), which we identified via reclassification of the TCGA-BRCA samples performed in another study\u003csup\u003e\u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e41\u003c/span\u003e\u003c/sup\u003e. Bulk RNA expression profiles reflect mixtures of diverse malignant and nonmalignant cell types and states present in the TME, and therefore, the expression of each metaprogram defines a rough estimate for the abundance of the corresponding cellular state in bulk samples. Before generating the survival curves for the extracted TCGA samples based on low or high cumulative expression of the genes in the metaprograms, we perform univariate Cox proportional hazard regression analysis on each gene individually. The point of this is to evaluate the association between the expression of each gene and survival, then select the genes with the strongest association (e.g., based on the p-value) for inclusion in the metaprogram survival analysis (\u003cb\u003esee Supplementary Table\u0026nbsp;6\u003c/b\u003e.) We then scored each bulk sample for expression of the survival-associated genes in each metaprogram from the 9 malignant TNBC, 14 nonmalignant TNBC and 32 reference normal breast metaprograms and examined the association between high or low expression and survival (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eC, \u003cb\u003eSupplementary Figs.\u0026nbsp;4, 5\u003c/b\u003e). The survival analyses for various metaprograms, including those derived from malignant TNBC cells, nonmalignant cells, and normal breast tissue (Reed2024 metaprograms), reveal interesting insights into their prognostic relevance in TNBC.\u003c/p\u003e \u003cp\u003eFor TNBC metaprograms (\u003cb\u003eSupplementary Fig.\u0026nbsp;4C\u003c/b\u003e), TNBC_MP3 shows that higher expression correlates with significantly worse survival outcomes (\u003cem\u003ep\u003c/em\u003e = 0.018). TNBC_MP4 also associates higher gene expression with poorer survival (\u003cem\u003ep\u003c/em\u003e = 0.015). In contrast, TNBC_MP8 indicates that lower expression of its genes is linked to poorer survival (\u003cem\u003ep\u003c/em\u003e = 0.019). These findings underscore the distinct biological impacts of malignant metaprograms associated with enhanced metabolism, a mesenchymal transcriptomic state and immune interferon response can have on patient outcome. The nonmalignant cell (NC) metaprograms further highlight the critical role of nonmalignant cells in the tumor microenvironment and their influence on disease outcomes (\u003cb\u003eSupplementary Fig.\u0026nbsp;4\u003c/b\u003e). These NC metaprograms demonstrate the significant impact of nonmalignant TME on TNBC outcomes. The Reed2024 metaprograms, derived from normal breast tissue, provide additional insights into how these normal tissue signatures are associated with TNBC survival (\u003cb\u003eSupplementary Fig.\u0026nbsp;5\u003c/b\u003e). Although these Reed2024 metaprograms originate from normal breast tissue, their dysregulation in TNBC is associated with processes that drive tumor aggressiveness and reduced survival outcomes. Overall, these analyses underscore the diverse biological impacts of the malignant, nonmalignant, and normal breast metaprograms on TNBC prognosis.\u003c/p\u003e \u003cp\u003eNext, we sought to validate the biological relevance of the metaprograms within the broader context of breast cancer subtypes, assess how well these TNBC-specific metaprograms align with or diverge from established molecular classifications, and potentially uncover novel insights into the heterogeneity and underlying biology of TNBC. This analysis aimed to explore whether these metaprograms could refine our understanding of TNBC beyond the traditional PAM50 framework, with implications for more personalized treatment strategies. The PAM50 panel is a list of 50 genes used to classify breast cancers into five intrinsic subtypes in FFPE tissue sections via real time polymerase chain reaction (RT-PCR)\u003csup\u003e\u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e42\u003c/span\u003e\u003c/sup\u003e. Their utility lies in that PAM50 subtypes have been found independently prognostic for long-term breast cancer survival. This classification system was furthermore developed with the major types of breast cancer (HER2+, ER, and PR positive) in mind, not for TNBC.\u003c/p\u003e \u003cp\u003eTo compare the malignant metaprograms to the PAM50 classification, we scored each metaprogram in the identified TNBC TCGA patient samples and hierarchically clustered them to assess patient stratification \u003cb\u003e(Supplementary Fig.\u0026nbsp;6).\u003c/b\u003e After that, we compared them to their assigned PAM50 classification (\u003cb\u003eSupplementary Table\u0026nbsp;7\u003c/b\u003e). In our cohort, 87% of the samples were Basal, 8% Her2, 3% Normal, and 2% LumA. This is in line with what is known about PAM50 classification of TNBC\u003csup\u003e\u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e43\u003c/span\u003e\u003c/sup\u003e. After scoring each patient in the TCGA TNBC cohort for the malignant MP mean expression, we performed hierarchical clustering to identify five groups of patients in the Basal, and for the significantly different clusters, observed a median survival time of 1398 days in Cluster 2, consisting of 33 patients, compared to 1813 days for Cluster 5, consisting of 43 patients (\u003cb\u003eSupplementary Table\u0026nbsp;8\u003c/b\u003e) differences in survival time of defined patient groups based on MP. Overall, these results underscore the critical role of metaprogram expression levels in predicting survival outcomes for TNBC patients, highlighting the potential for these metaprograms to serve as valuable biomarkers for clinical prognosis and personalized treatment strategies.\u003c/p\u003e \u003cdiv id=\"Sec23\" class=\"Section3\"\u003e \u003c/div\u003e \u003c/div\u003e"},{"header":"Discussion \u0026 Conclusions","content":"\u003cp\u003eIn this study, we used single-cell RNA sequencing data to construct a comprehensive molecular map of TNBC and to address the significant challenge of intratumoral heterogeneity. TNBC’s aggressive nature and poor prognosis is driven by its heterogeneity, which complicates diagnosis, treatment selection, and patient outcomes. Using our TNBCMap, we identified distinct transcriptomic signatures, revealing critical insights into the molecular diversity of TNBC. However, as our analysis is based on publicly available datasets, the full spectrum of TNBC heterogeneity may not be entirely captured. While the publically available datasets included 53 TNBC patients, only 21 patients had sufficient cells for further analysis of the malignant component after filtering, which may impact the generalizability of our findings. Future studies incorporating a larger and more diverse cohort, will be necessary to further validate these transcriptomic states. Our findings reinforce the idea that TNBC is not a single disease but a collection of tumors with varying molecular and clinicopathological characteristics, based on the differing CNV profiles between the patients in our cohort. However, the identification of nine key transcriptomic metaprograms highlights the common biological processes that malignant TNBC cells co-opt to enhance tumor progression, including immune evasion and mimicry, EMT, metabolic reprogramming, cell cycle regulation, and ECM remodeling, vasculogenic mimicry, and finally, luminal- and basal-like transcriptomic states.\u003c/p\u003e\u003cp\u003eOur findings underscore the pivotal roles of vasculogenic and immunogenic mimicry in TNBC progression, highlighting their contributions to intratumoral heterogeneity and poor patient outcomes. The vasculogenic mimicry metaprogram (TNBC_MP9) shows that TNBC cells can adopt vascular-like properties, facilitating alternative vascular networks that correlate with more aggressive phenotypes and poorer survival, as supported by the literature\u003csup\u003e\u003cspan citationid=\"CR44\" class=\"CitationRef\"\u003e44\u003c/span\u003e,\u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e45\u003c/span\u003e,\u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e46\u003c/span\u003e\u003c/sup\u003e. Similarly, the metaprograms associated with immune mimicry (TNBC_MP1 and 5) reveals how TNBC cells co-opt immune signatures to evade immune surveillance, consistent with reports that immunogenic mimicry aids cancer in escaping immune detection\u003csup\u003e\u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e47\u003c/span\u003e,\u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e48\u003c/span\u003e\u003c/sup\u003e. By associating these mimicry-driven metaprograms with patient outcomes, we offer a unique single-cell perspective that validates and extends previous findings, suggesting that targeting these mimicry processes could be a promising therapeutic strategy in TNBC.\u003c/p\u003e\u003cp\u003eThe expression of an inflammatory metaprogram (TNBC_MP8) highlights the multifaceted role of inflammatory pathways within the tumor microenvironment, particularly in driving immune-modulatory behaviors. The strong enrichment of interferon-stimulated genes within this metaprogram suggests that TNBC cells actively utilize inflammatory signaling to shape immune responses. Studies have shown that interferon-gamma (IFNγ) can modulate antigen processing and presentation in TNBC cells, enhancing immune visibility\u003csup\u003e\u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e49\u003c/span\u003e\u003c/sup\u003e. At the same time, chronic inflammation, as seen in inflamed TNBC subpopulations, has been linked to chemotherapy resistance and genomic instability\u003csup\u003e\u003cspan citationid=\"CR50\" class=\"CitationRef\"\u003e50\u003c/span\u003e\u003c/sup\u003e. Although some research has associated inflammation with improved outcomes in TNBC\u003csup\u003e\u003cspan citationid=\"CR51\" class=\"CitationRef\"\u003e51\u003c/span\u003e\u003c/sup\u003e, our analysis did not reveal a direct survival association for TNBC_MP8, underscoring the complex and context-dependent role of inflammation. These findings highlight the dualistic nature of inflammatory signaling, both activating immune responses and promoting tumor persistence, and suggest that targeting inflammatory pathways in TNBC could reduce pro-tumorigenic inflammation and potentially improve therapeutic outcomes.\u003c/p\u003e\u003cp\u003eThe mesenchymal-like (TNBC_MP8, MES-like) phenotype, characterized by EMT, has significant implications for the progression and treatment of TNBC. EMT is linked to enhanced cellular motility, invasiveness, and resistance to apoptosis, making MES-like TNBC particularly aggressive and chemoresistant. MES-like states are often associated with poor patient outcomes due to their ability to evade standard therapies that target epithelial-like cells, so targeting the EMT process, either by reversing it through mesenchymal-epithelial transition (MET) or inhibiting key EMT drivers, can sensitize MES-like cells to treatment. These findings highlight the need for therapeutic strategies that disrupt the plasticity between EMT and MET states and cut off an avenue of evasion in TNBC, offering the potential to improve treatment outcomes by preventing metastasis and overcoming chemoresistance. However, EMT is a highly plastic and dynamic process, and static transcriptomic snapshots may not fully capture the temporal transitions between epithelial and mesenchymal states, which could influence therapeutic strategies.\u003c/p\u003e\u003cp\u003eIn terms of metabolism, the fact that we observe a metaprogram associated with upregulated metabolism (TNBC_MP3) and enhanced glycolysis is supported by literature, as TNBC generally exhibits a higher rate of glycolysis in general and compared to other types of breast cancer. One study identified a metabolic signature of enhanced glycolysis and lactate secretion. The identification of a proliferative metaprogram (TNBC_MP2) emphasizes the role enhanced proliferation has in TNBC, corroborating our results by earlier work that identified proliferation highly associated with luminal transcriptomic signature.\u003c/p\u003e\u003cp\u003eThe preferential association of certain subclones with specific metaprograms in TNBC suggests that these subpopulations occupy distinct niches within the tumor, driven by selective pressures that promote diverse adaptive traits. However, our dataset includes a limited number of patients, and larger studies incorporating longitudinal samples could provide further insights into the evolution and stability of these transcriptomic states over time. Karaayvaz et al. 2018 demonstrated that these subclonal populations not only coexist but also display gene expression profiles linked to poor patient outcomes, indicating a predisposition to resist therapy and drive metastasis\u003csup\u003e\u003cspan citationid=\"CR52\" class=\"CitationRef\"\u003e52\u003c/span\u003e\u003c/sup\u003e. Kim \u003cem\u003eet al.\u003c/em\u003e 2021 found that TNBC cells under chemotherapeutic stress often harbor resistant subclones with metabolic adaptations like enhanced oxidative phosphorylation, allowing them to survive treatment and contribute to recurrence and metastasis\u003csup\u003e\u003cspan citationid=\"CR53\" class=\"CitationRef\"\u003e53\u003c/span\u003e\u003c/sup\u003e. Mavrommati \u003cem\u003eet al\u003c/em\u003e. 2021 further emphasized that the persistence of distinct, resilient subclones complicates treatment, as single-target therapies may fail to eliminate the full spectrum of tumor heterogeneity, leaving resistant clones to propagate\u003csup\u003e\u003cspan citationid=\"CR54\" class=\"CitationRef\"\u003e54\u003c/span\u003e\u003c/sup\u003e. Consequently, effective TNBC treatment strategies must adopt multi-targeted approaches to address the range of subclonal states and prevent relapse by inhibiting multiple pathways concurrently.\u003c/p\u003e\u003cp\u003eNotably, certain metaprograms, such as those related to cell proliferation, basal-like characteristics and vasculogenic mimicry, were strongly associated with poorer survival outcomes, underscoring their role in TNBC malignancy. However, as our survival analysis is based on bulk RNA-seq data from TCGA, it does not fully capture the single-cell-level heterogeneity within tumors, highlighting the need for validation using single-cell or spatial transcriptomics approaches. This finding is further supported by recent research demonstrating that targeting the basal-like to luminal-like state transition offers a potential therapeutic strategy. A recent study by Schade et al. revealed that combining AKT and EZH2 inhibitors drives basal-like TNBC cells into a more differentiated, luminal-like state, resulting in substantial tumor regression\u003csup\u003e\u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e38\u003c/span\u003e\u003c/sup\u003e. This work underscores the importance of these two cell states in TNBC progression and suggests that therapies targeting this transition could improve patient outcomes by disrupting aggressive basal-like phenotypes and promoting differentiation.\u003c/p\u003e\u003cp\u003eAdditionally, our analysis extends beyond malignant cells to the tumor microenvironment, revealing that non-malignant cells, including immune and stromal components, play a significant role in tumor development. These findings support the notion that TNBC progression is shaped by both intrinsic tumor characteristics, extrinsic interactions within the microenvironment and evolutionary selection. This study underscores the necessity of embracing TNBC’s heterogeneity to develop more personalized therapeutic strategies based on the varying proportions of common characteristics in this disease. By identifying molecular subpopulations within TNBC, our work provides a foundation for the future integration of multi-omic approaches, potentially guiding more precise clinical interventions. Ultimately, understanding the complexity of TNBC at the single-cell level will be pivotal in overcoming current treatment challenges and improving patient outcomes.\u003c/p\u003e"},{"header":"Abbreviations","content":"\u003cp\u003eTNBC: Triple-Negative Breast Cancer\u003cbr\u003e\u0026nbsp;scRNA-seq: Single-Cell RNA Sequencing\u003cbr\u003e\u0026nbsp;ITH: Intratumoral Heterogeneity\u003cbr\u003e\u0026nbsp;TME: Tumor Microenvironment\u003cbr\u003e\u0026nbsp;EMT: Epithelial-to-Mesenchymal Transition\u003cbr\u003e\u0026nbsp;VM-like: Vasculogenic Mimicry-like\u003cbr\u003e\u0026nbsp;NMF: Non-Negative Matrix Factorization\u003cbr\u003e\u0026nbsp;CNV: Copy Number Variation\u003cbr\u003e\u0026nbsp;MP: Metaprogram\u003cbr\u003e\u0026nbsp;GSEA: Gene Set Enrichment Analysis\u003cbr\u003e\u0026nbsp;PAM50: Prediction Analysis of Microarray 50 (breast cancer classification system)\u003cbr\u003e\u0026nbsp;TCGA-BRCA: The Cancer Genome Atlas - Breast Invasive Carcinoma\u003cbr\u003e\u0026nbsp;IM-like: Immune-like\u003cbr\u003e\u0026nbsp;MES-like: Mesenchymal-like\u003cbr\u003e\u0026nbsp;cNMF: Consensus Non-Negative Matrix Factorization\u003cbr\u003e\u0026nbsp;SCEVAN: Single-Cell Evolutionary Analysis (CNV-based subclone analysis)\u003cbr\u003e\u0026nbsp;VST: Variance-Stabilizing Transformation\u003cbr\u003e\u0026nbsp;RT-PCR: Real-Time Polymerase Chain Reaction\u003cbr\u003e\u0026nbsp;FFPE: Formalin-Fixed Paraffin-Embedded\u003cbr\u003e\u0026nbsp;QC: Quality Control\u003cbr\u003e\u0026nbsp;UMAP: Uniform Manifold Approximation and Projection\u003cbr\u003e\u0026nbsp;PCA: Principal Component Analysis\u003cbr\u003e\u0026nbsp;HUGO: Human Genome Organization\u003cbr\u003e\u0026nbsp;FPKM-UQ: Fragments Per Kilobase of transcript per Million mapped reads - Upper Quartile\u003cbr\u003e\u0026nbsp;NC: Nonmalignant Cluster\u003cbr\u003e\u0026nbsp;GDC: Genomic Data Commons\u003cbr\u003e\u0026nbsp;DESeq2: Differential Expression Analysis for Sequence Count Data\u003cbr\u003e\u0026nbsp;FDR: False Discovery Rate\u003cbr\u003e\u0026nbsp;cBioPortal: Cancer Genomics Data Portal\u003cbr\u003e\u0026nbsp;CoxPH: Cox Proportional Hazards Model\u003cbr\u003e\u0026nbsp;MPs: Metaprograms\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eEthics approval and consent to participate\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eNot applicable.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eConsent for publication\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eNot applicable.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAvailability of data and materials\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe processed data is stored in the Figshare (https://doi.org/10.6084/m9.figshare.27996128). The code used in this study is available in the GitHub repositories at https://github.com/eriksamuelsson1/TNBC.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCompeting interests\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe authors declare that they have no competing interests.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eFunding\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThis work is supported by funds from by Chan Zuckerberg Initiative; an advised fund of the Silicon Valley Community Foundation; the Erling-Persson Family Foundation (Erling-Perssons Stiftelse; the Human Developmental Cell Atlas); the Knut and Alice Wallenberg Foundation (Knut och Alice Wallenbergs Stiftelse; KAW 2018.0172); the Swedish Research Council (Vetenskapsr\u0026aring;det; 2019-01238); and the Swedish Cancer Society (Cancerfonden; CAN 2021/1726).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAuthor contributions\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eConceptualization: E.S., M.N.; Methodology: E.S.; Software:E.S., R.F.; Validation: E.S., R.F., T.H; Formal analysis: , E.S., R.F.; Investigation: E.S.; Resources: E.S., R.F.; Data curation: E.S. R.F. ; Writing - original draft: E.S.; Writing - review \u0026amp; editing: E.S, T.H., M.N. ; Visualization: E.S, T.H.; Supervision: M.N., E.S.; Project administration: M.N., E.S.; Funding acquisition: M.N., E.S..\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAcknowledgements\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThis work was supported by the National Bioinformatics Infrastructure Sweden (NBIS) at SciLifeLab. The computations and data storage was enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS), partially funded by the Swedish Research Council through grant agreement no. 2022-06725.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eDerakhshan F, Reis-Filho JS. Pathogenesis of Triple-Negative Breast Cancer. \u003cem\u003eAnnu Rev Pathol\u003c/em\u003e \u003cstrong\u003e17\u003c/strong\u003e, 181-204 (2022).\u003c/li\u003e\n\u003cli\u003eAlmansour NM. Triple-Negative Breast Cancer: A Brief Review About Epidemiology, Risk Factors, Signaling Pathways, Treatment and Role of Artificial Intelligence. \u003cem\u003eFront Mol Biosci\u003c/em\u003e \u003cstrong\u003e9\u003c/strong\u003e, 836417 (2022).\u003c/li\u003e\n\u003cli\u003eTan PH\u003cem\u003e, et al.\u003c/em\u003e The 2019 World Health Organization classification of tumours of the breast. \u003cem\u003eHistopathology\u003c/em\u003e \u003cstrong\u003e77\u003c/strong\u003e, 181-185 (2020).\u003c/li\u003e\n\u003cli\u003eZhao S, Zuo WJ, Shao ZM, Jiang YZ. Molecular subtypes and precision treatment of triple-negative breast cancer. \u003cem\u003eAnn Transl Med\u003c/em\u003e \u003cstrong\u003e8\u003c/strong\u003e, 499 (2020).\u003c/li\u003e\n\u003cli\u003eGuo L\u003cem\u003e, et al.\u003c/em\u003e Breast cancer heterogeneity and its implication in personalized precision therapy. \u003cem\u003eExp Hematol Oncol\u003c/em\u003e \u003cstrong\u003e12\u003c/strong\u003e, 3 (2023).\u003c/li\u003e\n\u003cli\u003eLuond F, Tiede S, Christofori G. Breast cancer as an example of tumour heterogeneity and tumour cell plasticity during malignant progression. \u003cem\u003eBr J Cancer\u003c/em\u003e \u003cstrong\u003e125\u003c/strong\u003e, 164-175 (2021).\u003c/li\u003e\n\u003cli\u003ePolyak K. Heterogeneity in breast cancer. \u003cem\u003eJ Clin Invest\u003c/em\u003e \u003cstrong\u003e121\u003c/strong\u003e, 3786-3788 (2011).\u003c/li\u003e\n\u003cli\u003eRivenbark AG, O\u0026apos;Connor SM, Coleman WB. Molecular and cellular heterogeneity in breast cancer: challenges for personalized medicine. \u003cem\u003eAm J Pathol\u003c/em\u003e \u003cstrong\u003e183\u003c/strong\u003e, 1113-1124 (2013).\u003c/li\u003e\n\u003cli\u003eTurashvili G, Brogi E. Tumor Heterogeneity in Breast Cancer. \u003cem\u003eFront Med (Lausanne)\u003c/em\u003e \u003cstrong\u003e4\u003c/strong\u003e, 227 (2017).\u003c/li\u003e\n\u003cli\u003eLehmann BD\u003cem\u003e, et al.\u003c/em\u003e Identification of human triple-negative breast cancer subtypes and preclinical models for selection of targeted therapies. \u003cem\u003eJ Clin Invest\u003c/em\u003e \u003cstrong\u003e121\u003c/strong\u003e, 2750-2767 (2011).\u003c/li\u003e\n\u003cli\u003eLehmann BD\u003cem\u003e, et al.\u003c/em\u003e Refinement of Triple-Negative Breast Cancer Molecular Subtypes: Implications for Neoadjuvant Chemotherapy Selection. \u003cem\u003ePLoS One\u003c/em\u003e \u003cstrong\u003e11\u003c/strong\u003e, e0157368 (2016).\u003c/li\u003e\n\u003cli\u003eWang DY, Jiang Z, Ben-David Y, Woodgett JR, Zacksenhaus E. Molecular stratification within triple-negative breast cancer subtypes. \u003cem\u003eSci Rep\u003c/em\u003e \u003cstrong\u003e9\u003c/strong\u003e, 19107 (2019).\u003c/li\u003e\n\u003cli\u003eMasuda H\u003cem\u003e, et al.\u003c/em\u003e Reverse phase protein array identification of triple-negative breast cancer subtypes and comparison with mRNA molecular subtypes. \u003cem\u003eOncotarget\u003c/em\u003e \u003cstrong\u003e8\u003c/strong\u003e, 70481-70495 (2017).\u003c/li\u003e\n\u003cli\u003eDiNome ML\u003cem\u003e, et al.\u003c/em\u003e Clinicopathological Features of Triple-Negative Breast Cancer Epigenetic Subtypes. \u003cem\u003eAnn Surg Oncol\u003c/em\u003e \u003cstrong\u003e26\u003c/strong\u003e, 3344-3353 (2019).\u003c/li\u003e\n\u003cli\u003eGong Y\u003cem\u003e, et al.\u003c/em\u003e Metabolic-Pathway-Based Subtyping of Triple-Negative Breast Cancer Reveals Potential Therapeutic Targets. \u003cem\u003eCell Metab\u003c/em\u003e \u003cstrong\u003e33\u003c/strong\u003e, 51-64 e59 (2021).\u003c/li\u003e\n\u003cli\u003eEnsenyat-Mendez M\u003cem\u003e, et al.\u003c/em\u003e Current Triple-Negative Breast Cancer Subtypes: Dissecting the Most Aggressive Form of Breast Cancer. \u003cem\u003eFront Oncol\u003c/em\u003e \u003cstrong\u003e11\u003c/strong\u003e, 681476 (2021).\u003c/li\u003e\n\u003cli\u003eTurner KM, Yeo SK, Holm TM, Shaughnessy E, Guan JL. Heterogeneity within molecular subtypes of breast cancer. \u003cem\u003eAm J Physiol Cell Physiol\u003c/em\u003e \u003cstrong\u003e321\u003c/strong\u003e, C343-C354 (2021).\u003c/li\u003e\n\u003cli\u003eBeroukhim R\u003cem\u003e, et al.\u003c/em\u003e The landscape of somatic copy-number alteration across human cancers. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e463\u003c/strong\u003e, 899-905 (2010).\u003c/li\u003e\n\u003cli\u003eHanahan D. Hallmarks of Cancer: New Dimensions. \u003cem\u003eCancer Discov\u003c/em\u003e \u003cstrong\u003e12\u003c/strong\u003e, 31-46 (2022).\u003c/li\u003e\n\u003cli\u003eStratton MR, Campbell PJ, Futreal PA. The cancer genome. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e458\u003c/strong\u003e, 719-724 (2009).\u003c/li\u003e\n\u003cli\u003eCai H, Kumar N, Ai N, Gupta S, Rath P, Baudis M. Progenetix: 12 years of oncogenomic data curation. \u003cem\u003eNucleic Acids Res\u003c/em\u003e \u003cstrong\u003e42\u003c/strong\u003e, D1055-1062 (2014).\u003c/li\u003e\n\u003cli\u003eSteele CD\u003cem\u003e, et al.\u003c/em\u003e Signatures of copy number alterations in human cancer. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e606\u003c/strong\u003e, 984-991 (2022).\u003c/li\u003e\n\u003cli\u003eCai H, Gupta S, Rath P, Ai N, Baudis M. arrayMap 2014: an updated cancer genome resource. \u003cem\u003eNucleic Acids Res\u003c/em\u003e \u003cstrong\u003e43\u003c/strong\u003e, D825-830 (2015).\u003c/li\u003e\n\u003cli\u003eKim TM, Xi R, Luquette LJ, Park RW, Johnson MD, Park PJ. Functional genomic analysis of chromosomal aberrations in a compendium of 8000 cancer genomes. \u003cem\u003eGenome Res\u003c/em\u003e \u003cstrong\u003e23\u003c/strong\u003e, 217-227 (2013).\u003c/li\u003e\n\u003cli\u003eZack TI\u003cem\u003e, et al.\u003c/em\u003e Pan-cancer patterns of somatic copy number alteration. \u003cem\u003eNat Genet\u003c/em\u003e \u003cstrong\u003e45\u003c/strong\u003e, 1134-1140 (2013).\u003c/li\u003e\n\u003cli\u003eCurtis C\u003cem\u003e, et al.\u003c/em\u003e The genomic and transcriptomic architecture of 2,000 breast tumours reveals novel subgroups. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e486\u003c/strong\u003e, 346-352 (2012).\u003c/li\u003e\n\u003cli\u003eShiu KK, Natrajan R, Geyer FC, Ashworth A, Reis-Filho JS. DNA amplifications in breast cancer: genotypic-phenotypic correlations. \u003cem\u003eFuture Oncol\u003c/em\u003e \u003cstrong\u003e6\u003c/strong\u003e, 967-984 (2010).\u003c/li\u003e\n\u003cli\u003eTurner N\u003cem\u003e, et al.\u003c/em\u003e Integrative molecular profiling of triple negative breast cancers identifies amplicon drivers and potential therapeutic targets. \u003cem\u003eOncogene\u003c/em\u003e \u003cstrong\u003e29\u003c/strong\u003e, 2013-2023 (2010).\u003c/li\u003e\n\u003cli\u003eBerger AC\u003cem\u003e, et al.\u003c/em\u003e A Comprehensive Pan-Cancer Molecular Study of Gynecologic and Breast Cancers. \u003cem\u003eCancer Cell\u003c/em\u003e \u003cstrong\u003e33\u003c/strong\u003e, 690-705 e699 (2018).\u003c/li\u003e\n\u003cli\u003eNeftel C\u003cem\u003e, et al.\u003c/em\u003e An Integrative Model of Cellular States, Plasticity, and Genetics for Glioblastoma. \u003cem\u003eCell\u003c/em\u003e \u003cstrong\u003e178\u003c/strong\u003e, 835-849 e821 (2019).\u003c/li\u003e\n\u003cli\u003eGavish A\u003cem\u003e, et al.\u003c/em\u003e Hallmarks of transcriptional intratumour heterogeneity across a thousand tumours. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e618\u003c/strong\u003e, 598-606 (2023).\u003c/li\u003e\n\u003cli\u003eDe Falco A, Caruso F, Su XD, Iavarone A, Ceccarelli M. A variational algorithm to detect the clonal copy number substructure of tumors from scRNA-seq data. \u003cem\u003eNat Commun\u003c/em\u003e \u003cstrong\u003e14\u003c/strong\u003e, 1074 (2023).\u003c/li\u003e\n\u003cli\u003eKong D, Hughes CJ, Ford HL. Cellular Plasticity in Breast Cancer Progression and Therapy. \u003cem\u003eFront Mol Biosci\u003c/em\u003e \u003cstrong\u003e7\u003c/strong\u003e, 72 (2020).\u003c/li\u003e\n\u003cli\u003eKvokackova B, Remsik J, Jolly MK, Soucek K. Phenotypic Heterogeneity of Triple-Negative Breast Cancer Mediated by Epithelial-Mesenchymal Plasticity. \u003cem\u003eCancers (Basel)\u003c/em\u003e \u003cstrong\u003e13\u003c/strong\u003e, (2021).\u003c/li\u003e\n\u003cli\u003ePujals M\u003cem\u003e, et al.\u003c/em\u003e RAGE/SNAIL1 signaling drives epithelial-mesenchymal plasticity in metastatic triple-negative breast cancer. \u003cem\u003eOncogene\u003c/em\u003e \u003cstrong\u003e42\u003c/strong\u003e, 2610-2628 (2023).\u003c/li\u003e\n\u003cli\u003eHeilala M\u003cem\u003e, et al.\u003c/em\u003e Fibrin Stiffness Regulates Phenotypic Plasticity of Metastatic Breast Cancer Cells. \u003cem\u003eAdv Healthc Mater\u003c/em\u003e \u003cstrong\u003e12\u003c/strong\u003e, e2301137 (2023).\u003c/li\u003e\n\u003cli\u003eGuo Z, Han S. Targeting cancer stem cell plasticity in triple-negative breast cancer. \u003cem\u003eExplor Target Antitumor Ther\u003c/em\u003e \u003cstrong\u003e4\u003c/strong\u003e, 1165-1181 (2023).\u003c/li\u003e\n\u003cli\u003eSchade AE\u003cem\u003e, et al.\u003c/em\u003e AKT and EZH2 inhibitors kill TNBCs by hijacking mechanisms of involution. \u003cem\u003eNature\u003c/em\u003e, (2024).\u003c/li\u003e\n\u003cli\u003eReed AD\u003cem\u003e, et al.\u003c/em\u003e A single-cell atlas enables mapping of homeostatic cellular shifts in the adult human breast. \u003cem\u003eNat Genet\u003c/em\u003e \u003cstrong\u003e56\u003c/strong\u003e, 652-662 (2024).\u003c/li\u003e\n\u003cli\u003eBurstein MD\u003cem\u003e, et al.\u003c/em\u003e Comprehensive genomic analysis identifies novel subtypes and targets of triple-negative breast cancer. \u003cem\u003eClin Cancer Res\u003c/em\u003e \u003cstrong\u003e21\u003c/strong\u003e, 1688-1698 (2015).\u003c/li\u003e\n\u003cli\u003eLehmann BD\u003cem\u003e, et al.\u003c/em\u003e Multi-omics analysis identifies therapeutic vulnerabilities in triple-negative breast cancer subtypes. \u003cem\u003eNat Commun\u003c/em\u003e \u003cstrong\u003e12\u003c/strong\u003e, 6276 (2021).\u003c/li\u003e\n\u003cli\u003ePerou CM\u003cem\u003e, et al.\u003c/em\u003e Molecular portraits of human breast tumours. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e406\u003c/strong\u003e, 747-752 (2000).\u003c/li\u003e\n\u003cli\u003eLehmann BD, Pietenpol JA. Identification and use of biomarkers in treatment strategies for triple-negative breast cancer subtypes. \u003cem\u003eJ Pathol\u003c/em\u003e \u003cstrong\u003e232\u003c/strong\u003e, 142-150 (2014).\u003c/li\u003e\n\u003cli\u003eZheng S\u003cem\u003e, et al.\u003c/em\u003e Vasculogenic mimicry regulates immune infiltration and mutational status of the tumor microenvironment in breast cancer to influence tumor prognosis. \u003cem\u003eEnviron Toxicol\u003c/em\u003e \u003cstrong\u003e39\u003c/strong\u003e, 2948-2960 (2024).\u003c/li\u003e\n\u003cli\u003eLiang X\u003cem\u003e, et al.\u003c/em\u003e Identification of new subtypes of breast cancer based on vasculogenic mimicry related genes and a new model for predicting the prognosis of breast cancer. \u003cem\u003eHeliyon\u003c/em\u003e \u003cstrong\u003e10\u003c/strong\u003e, e36565 (2024).\u003c/li\u003e\n\u003cli\u003eAndonegui-Elguera MA, Alfaro-Mora Y, Caceres-Gutierrez R, Caro-Sanchez CHS, Herrera LA, Diaz-Chavez J. An Overview of Vasculogenic Mimicry in Breast Cancer. \u003cem\u003eFront Oncol\u003c/em\u003e \u003cstrong\u003e10\u003c/strong\u003e, 220 (2020).\u003c/li\u003e\n\u003cli\u003eGao R\u003cem\u003e, et al.\u003c/em\u003e Cancer cell immune mimicry delineates onco-immunologic modulation. \u003cem\u003eiScience\u003c/em\u003e \u003cstrong\u003e24\u003c/strong\u003e, 103133 (2021).\u003c/li\u003e\n\u003cli\u003eTimar J, Honn KV, Hendrix MJC, Marko-Varga G, Jalkanen S. Newly identified form of phenotypic plasticity of cancer: immunogenic mimicry. \u003cem\u003eCancer Metastasis Rev\u003c/em\u003e \u003cstrong\u003e42\u003c/strong\u003e, 323-334 (2023).\u003c/li\u003e\n\u003cli\u003eGoncalves G\u003cem\u003e, et al.\u003c/em\u003e IFNgamma Modulates the Immunopeptidome of Triple Negative Breast Cancer Cells by Enhancing and Diversifying Antigen Processing and Presentation. \u003cem\u003eFront Immunol\u003c/em\u003e \u003cstrong\u003e12\u003c/strong\u003e, 645770 (2021).\u003c/li\u003e\n\u003cli\u003eJacobo Jacobo M, Donnella HJ, Sobti S, Kaushik S, Goga A, Bandyopadhyay S. An inflamed tumor cell subpopulation promotes chemotherapy resistance in triple negative breast cancer. \u003cem\u003eSci Rep\u003c/em\u003e \u003cstrong\u003e14\u003c/strong\u003e, 3694 (2024).\u003c/li\u003e\n\u003cli\u003eOshi M\u003cem\u003e, et al.\u003c/em\u003e Inflammation Is Associated with Worse Outcome in the Whole Cohort but with Better Outcome in Triple-Negative Subtype of Breast Cancer Patients. \u003cem\u003eJ Immunol Res\u003c/em\u003e \u003cstrong\u003e2020\u003c/strong\u003e, 5618786 (2020).\u003c/li\u003e\n\u003cli\u003eKaraayvaz M\u003cem\u003e, et al.\u003c/em\u003e Unravelling subclonal heterogeneity and aggressive disease states in TNBC through single-cell RNA-seq. \u003cem\u003eNat Commun\u003c/em\u003e \u003cstrong\u003e9\u003c/strong\u003e, 3588 (2018).\u003c/li\u003e\n\u003cli\u003eKim C\u003cem\u003e, et al.\u003c/em\u003e Chemoresistance Evolution in Triple-Negative Breast Cancer Delineated by Single-Cell Sequencing. \u003cem\u003eCell\u003c/em\u003e \u003cstrong\u003e173\u003c/strong\u003e, 879-893 e813 (2018).\u003c/li\u003e\n\u003cli\u003eMavrommati I, Johnson F, Echeverria GV, Natrajan R. Subclonal heterogeneity and evolution in breast cancer. \u003cem\u003eNPJ Breast Cancer\u003c/em\u003e \u003cstrong\u003e7\u003c/strong\u003e, 155 (2021).\u003c/li\u003e\n\u003cli\u003eWolf FA, Angerer P, Theis FJ. SCANPY: large-scale single-cell gene expression data analysis. \u003cem\u003eGenome Biol\u003c/em\u003e \u003cstrong\u003e19\u003c/strong\u003e, 15 (2018).\u003c/li\u003e\n\u003cli\u003eWolock SL, Lopez R, Klein AM. Scrublet: Computational Identification of Cell Doublets in Single-Cell Transcriptomic Data. \u003cem\u003eCell Syst\u003c/em\u003e \u003cstrong\u003e8\u003c/strong\u003e, 281-291 e289 (2019).\u003c/li\u003e\n\u003cli\u003eKorsunsky I\u003cem\u003e, et al.\u003c/em\u003e Fast, sensitive and accurate integration of single-cell data with Harmony. \u003cem\u003eNat Methods\u003c/em\u003e \u003cstrong\u003e16\u003c/strong\u003e, 1289-1296 (2019).\u003c/li\u003e\n\u003cli\u003eXu C\u003cem\u003e, et al.\u003c/em\u003e Automatic cell-type harmonization and integration across Human Cell Atlas datasets. \u003cem\u003eCell\u003c/em\u003e \u003cstrong\u003e186\u003c/strong\u003e, 5876-5891 e5820 (2023).\u003c/li\u003e\n\u003cli\u003eDominguez Conde C\u003cem\u003e, et al.\u003c/em\u003e Cross-tissue immune cell analysis reveals tissue-specific features in humans. \u003cem\u003eScience\u003c/em\u003e \u003cstrong\u003e376\u003c/strong\u003e, eabl5197 (2022).\u003c/li\u003e\n\u003cli\u003eColaprico A\u003cem\u003e, et al.\u003c/em\u003e TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. \u003cem\u003eNucleic Acids Res\u003c/em\u003e \u003cstrong\u003e44\u003c/strong\u003e, e71 (2016).\u003c/li\u003e\n\u003c/ol\u003e"},{"header":"Tables","content":"\u003cp\u003eTables 1-8 are available in the Supplementary Files section.\u003c/p\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"Triple-negative breast cancer, Single-cell RNA sequencing, Consensus Non-Negative Matrix Factorization, Intratumoral heterogeneity, Metaprograms, Copy Number Variations, Tumor microenvironment, Prognostic biomarkers.","lastPublishedDoi":"10.21203/rs.3.rs-5974271/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-5974271/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eTriple-negative breast cancer (TNBC), an aggressive and highly heterogeneous subtype of breast cancer, poses significant challenges for treatment due to its molecular diversity and resistance to standard therapies. Accounting for 10\u0026ndash;20% of all breast cancer cases, TNBC lacks specific biological markers, making it difficult to classify and treat effectively. Traditional approaches based on bulk RNA sequencing obscure intratumoral heterogeneity and fail to capture distinct cellular states within tumors. In this study, we constructed a comprehensive single-cell transcriptomic map of TNBC by analyzing a cohort of published TNBC patient datasets, identifying nine transcriptomic states, or metaprograms, which capture the core behaviors of TNBC cells, including cancer stem cell properties, epithelial-to-mesenchymal transition (EMT), immune modulation, metabolic adaptation, and vasculogenic mimicry. We observed that these metaprograms are variably expressed within and across patient tumors, underscoring the complexity of TNBC. By integrating TNBC-specific metaprograms with established breast cancer subtypes, we found significant prognostic associations, with specific metaprograms correlating with poor survival outcomes. This study highlights the need for single-cell approaches to uncover TNBC\u0026rsquo;s molecular heterogeneity and suggests that metaprogram-based classification could facilitate more precise therapeutic interventions.\u003c/p\u003e","manuscriptTitle":"An Integrative Model of Single Cell Transcriptomic States for Triple-Negative Breast Cancer","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-02-19 12:06:49","doi":"10.21203/rs.3.rs-5974271/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"2d3b788e-2f85-4b98-a93b-c9dbeb240b77","owner":[],"postedDate":"February 19th, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[],"tags":[],"updatedAt":"2025-09-23T04:53:28+00:00","versionOfRecord":[],"versionCreatedAt":"2025-02-19 12:06:49","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-5974271","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-5974271","identity":"rs-5974271","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}