Globally, there were approximately 2.3 million new cases of and 685,000 deaths due to breast cancer (BC) in 2020 (Lei et al., 2021). BC is also the leading cause of cancer-related deaths among women (Beiki et al., 2012; World Health Organization, 2018) and the second leading cause of cancer deaths globally, worldwide (Beiki et al., 2012). The tumorigenesis involves uncontrolled growth of cells in breast tissues which can be either benign or malignant (Liu et al., 2013). Several studies on breast cancer patients have revealed the different anti- and pro-tumorigenic roles of the CAFs involved (Chang et al., 2012; Brechbuhl et al., 2017; Su et al., 2018).

Extensive studies over the past few decades have uncovered a variety of cell populations in tumors, thus leading to the active research area of intratumor heterogeneity (ITH) (Marusyk et al., 2020). In 2010, Hanahan and Weinberg noted that tumors exhibit an additional dimension of complexity through their “tumor microenvironment” that contributes to the acquisition of the so-called hallmark traits of cancer. ITH is attributed to genetic, epigenetic, and microenvironmental factors (McGranahan and Swanton, 2017; Marusyk et al., 2020) and associated with poor prognosis, therapeutic resistance and treatment failure leading to poor overall survival in cancer patients (Landau et al., 2013; Patel et al., 2014; Zhang et al., 2014; Jamal-Hanjani et al., 2015; Jamal-Hanjani et al., 2017). Indeed, the persistence of some of the drug-tolerant intratumor cell populations could be attributed to their high phenotypic plasticity (Flavahan et al., 2017).

Interestingly, hierarchies of differentiation also exist among normal cells in healthy tissues, but the populations of tumor cells display far greater cell-to-cell variability and the resulting phenotypic instability (Landau et al., 2014; Jenkinson et al., 2017). Such ITH could be attributed to genetic causes ranging from aneuploidy to other factors such as complex contextual signals in the highly aberrant tumor microenvironments, or even global alterations in cancer cell epigenomes (Senovilla et al., 2012). ITH also involves immune cell infiltration, which is important to immunotherapies. Tumor antigen diversity could be determined by the T Cell clonality in the different regions of the same tumor (Senovilla et al., 2012). Studies have shown spatially complex interactions between tumor microenvironments and the patient’s immune system (Hanahan and Weinberg, 2011; Vogelstein and Kinzler, 2015).

While heterogeneous cell types are prevalent within the tumor microenvironment, some of which may account for cancer development and progression, it also contains different non-malignant components, including the cancer-associated fibroblasts (CAFs) (Pietras and Östman, 2010; Cortez et al., 2014; Kalluri, 2016). Although the origin and activation mechanism of CAFs remains an area of active research (Anderberg and Pietras, 2009; Shiga et al., 2015; LeBleu and Kalluri, 2018; Chen and Song, 2019), studies have attributed the processes of formation and derivation of CAFs to various precursor cells (Anderberg and Pietras, 2009; Shiga et al., 2015; LeBleu and Kalluri, 2018; Chen and Song, 2019), which may be the source of the well-known heterogeneity among the CAFs (Du and Che, 2017; Öhlund et al., 2017; Costa et al., 2018; Raz et al., 2018; Lee et al., 2020). Indeed, in certain tumors, such as in the breast, in which the prevalence of CAFs could be as high as 80%, they can play both anti-as well as pro-tumorigenic roles (Chang et al., 2012; Brechbuhl et al., 2017; Su et al., 2018). Importantly, CAFs can facilitate drug resistance dynamically by altering the cell-matrix interactions that control the outer layer of cells’ sensitivity to apoptosis, producing proteins that control cell survival and proliferation, assisting with cell-cell communications, and activating epigenetic plasticity in neighboring cells (Cuiffo and Karnoub, 2012; Junttila and De Sauvage, 2013). CAF-targeted treatments can have dual effects depending on the target and the tissue under consideration (Özdemir et al., 2014; Koliaraki et al., 2015; Wagner, 2016). For instance, spatial proximity to CAFs has been shown to impact molecular features and therapeutic sensitivity of breast cancer cells influencing clinical outcomes (Marusyk et al., 2016).

In recent years, higher resolution, tissue-specific gene expression analysis is made possible by using new platforms such as single-cell RNA sequencing (scRNA-seq), which has rapidly evolved as a powerful and popular tool (Kalisky et al., 2018; Sun et al., 2021). Unlike previous transcriptomic studies that assayed a “bulk” sample, scRNA-seq data can provide a detailed characterization of each tumor. Indeed, the Human Tumor Atlas Network [https://humantumoratlas.org] is increasingly enriched with data on human cancers based on scRNA-seq assays. The high-resolution transcriptomic platform has led to several scRNA-seq studies of the composition of CAFs in different stages of cancer (Bernardo and Fibbe, 2013; Li et al., 2017; Puram et al., 2017; Lambrechts et al., 2018; Elyada et al., 2019; Hosein et al., 2019; Davidson et al., 2020; Dominguez et al., 2020; Friedman et al., 2020). For focused understanding of the heterogeneous expressions of genes, different sites of the same tumor were analyzed with multiregional RNA sequencing for different cancers (Gerlinger et al., 2012; Zhang et al., 2014; Yates et al., 2015; Thrane et al., 2018).

Despite the advancements and efficacy of scRNA-seq, the lack of spatial information in scRNA-seq analysis is a significant shortcoming for typical scRNA-seq methods to capture cellular heterogeneity. For a tumor sample, the presence of spatial contexts might play a major role which could be combined with scRNA-seq data with the explicit aim to capture microenvironmental heterogeneity. Spatial cell-to-cell communication in a given tissue image can be recovered from a spatial scRNA sequencing data via computational spatial re-mapping (Teves and Won, 2020). Alternatively, integration of high-resolution gene expression data with spatial coordinates can resolve such experimental shortcomings (Eng et al., 2019). While imaging the transcriptome in situ with high accuracy has been a major challenge in single-cell biology, development of high-throughput platforms for sequential fluorescence in situ hybridization such as RNA seqFISH+ and algorithms such as CELESTA can identify cell populations and their spatial organization in intact tissues (Zhang et al., 2014; Eng et al., 2019). Towards this, many recent efforts have developed methods to analyze spatial information in single-cell studies (Lee et al., 2014; McKenna et al., 2016; Shah et al., 2016; Frieda et al., 2017; Alemany et al., 2018; Codeluppi et al., 2018; Raj et al., 2018; Spanjaard et al., 2018; Wang et al., 2018).

High-throughput transcriptomic data are useful not only for identifying genes that are differentially expressed, but also to test for co-regulation of multiple genes, i.e., a gene-set, based on existing empirical knowledge of biological pathways and gene signatures, e.g., the well-known hallmarks of cancer. In this direction, several methods for gene-set analysis (GSA) were introduced by (Goeman et al., 2004; Mansmann and Meister, 2005; Subramanian et al., 2005; Kong et al., 2006; Dinu et al., 2007; Efron and Tibshirani, 2007). Since the genes within such gene-sets share a common biological function, considering the correlations within each set is a key aspect of a useful GSA method. However, it was shown by (Tsai and Chen, 2009) that the above GSA methods were affected by large type II errors.

An important limitation of many GSA methods is that they can only accommodate binary outcomes, such as disease versus control. Our method, Linear Combination Test (LCT) is a GSA method that was designed to address these limitations by taking into account correlations across genes and outcomes, and dealing with binary, univariate or multivariate continuous outcomes, measured either at a single point in time or at multiple time points, and therefore, allow us to analyze a wider range of studies involving complex study designs (Wang et al., 2014). Studies have shown that LCT can overcome difficulties such as small sample size, large gene-sets, and can accommodate correlations across gene-sets, time points, and multiple correlated continuous phenotypes (Dinu et al., 2013). Thus, while a specific gene may not show consistent expression across individual cells, LCT is more likely than traditional approaches to detect the regulation of a functional process or biological pathway associated with the intercellular diversity of outcomes in a single cell level experiment.

Recently, we have extended LCT beyond any other “bulk” GSA method for application to single cell experiments (Dinu et al., 2021). However, GSA is considerably more complicated in the presence of spatial information since the analyzed gene-sets need not have a uniform impact over the entire area of a spatially continuous phenotype. In fact, the significance of association between a selected gene-set and a particular phenotypic context at various microenvironmental neighborhoods could be different. Yet, variable as they may be, since spatial effects are generally continuous in nature, proximity may determine more correlated associations than those across distant locations within the same tumor space. Notably, this alludes to Tobler’s First Law of Geography, which states that “everything is related to everything else, but near things are more related than distant things.” Traditional testing of such relationships involves global or “aspatial” regression, with the implicit assumption that the impact of the genes in a gene-set (covariates) on the phenotype (spatial outcome) is constant across the tumor space (study area). In the presence of ITH, such stationarity assumption is unlikely to be valid. Geographically weighted regression (GWR) is a well-known method (Brunsdon et al., 1996) that avoids this problem by performing the regression within local windows and each observation is weighted according to its proximity to the center of the window. Adaptive kernel bandwidths allow for heterogeneity among densities of gene expression over the windows in different parts of the study region. Local regression coefficients and associated statistics are mapped to visualize how the explanatory power of a gene-set on the associated phenotypes changes spatially.

In the present study, we combined gene-set analysis of LCT with the local spatial modeling of GWR with the aim to develop geographically weighted LCT (GWLCT) as a statistical test. We demonstrated it on spatial scRNA-seq data from a real breast tumor sample and obtained key insights into its molecular heterogeneity across different spatially continuous phenotypic contexts defined by five well-known markers of CAFs. We note that GWLCT has several distinct advantages. While the popular GSA methods are aspatial and use only bulk gene expression data, GWLCT is developed for spatial single cell gene expression data. The geographical weighting scheme allows nearby neighborhoods to contribute more to each local model, and the regions with significant association of a selected gene-set and a corresponding phenotype are detected using scan statistics on the local test scores and illustrated as maps. At each location, the combined significance of such associations for the selected gene-sets is computed and visualized with a spatial heatmap. We also present new 3D interactive tools for insightful visualization of the tumor space. In the next section, we describe the data and methods, followed by the results of real tumor data analysis and simulations of different association scenarios using GWLCT, and end with discussion, including future work.

The three aspatial methods, LCT, RF-GSEA, and GSEA, were used to identify the “hallmark” cancer gene-sets that are significantly associated with five spatially continuous CAF phenotypes represented by their known markers C3, COL11A1, CXCL12, FBLN1, and S100A4 in the single cell breast cancer spatial transcriptomic data. The three aspatial methods, the phenotype, the size of each gene-set, p-value of the test, and the corresponding q-value are shown in Table 1. We note that such global p-values cannot be obtained from our proposed GWLCT, as this method is spatial in nature and assesses significance at every coordinate point rather than an overall significance measure. Using GSEA, the only significant gene-set was P53 (p-value = 0.009, q-value = 0.010). The results of LCT and RF-GSEA revealed that all the gene-sets are strongly associated with the five continuous phenotypes with p-value and q-value less than 0.001. This was expected since the candidate gene-sets represent commonly known “hallmarks” of cancer (Liberzon et al., 2011). The breast cancer data analysis interpretation resulting from the three aspatial methods considered is limited to the global significance values. This is an important limitation of aspatial methods. For the remainder of this section, we emphasize the advantages of GWLCT results interpretation in the context of spatial data analysis. Since the main advantage of GWLCT consists of an analysis beyond global significance, we study specifically at locations across the tumor space. Scan statistics provide a well-established computational method for detecting spatial clusters based on point count data. We computed scan statistics to detect the putative clusters of spatial regulation based on the number of cells with significantly enriched gene-sets occurring at locations where such counts exceeded what may be expected from an underlying Poisson distribution defined over the tumor space. The clusters are demarcated as white regions in Figure 1. Notably, the 3D plots, which are instrumental in showcasing the advantages of GWLCT over existing aspatial GSA methods, are presented here as well as at the following links on GitHub: https://mortezahaji.github.io/GWLCT-Project/. By clicking on any point in the plot, one can find the coordinate of the cell, corresponding number and name of significant gene-sets at three different levels: Low, Moderate, and High expressions of the selected CAF marker gene.

In addition, for local GWLCT, three different CAF categories were identified as Low (CAF gene expression less than 0.5), Moderate (CAF gene expression between 0.5 and 1), and High (CAF gene expression exceeding 1). Figure 2 demonstrates a snapshot of the 3D plots for the 5 phenotypes at three CAF levels. A snapshot of the 3D plot for the phenotype COL11A1 at high CAF level is also demonstrated in Figure 3. One is able to detect the frequency as well as the names of significant gene-sets at each location (based on the 8 gene-sets) by clicking on each dot, and rotating and zooming into the 3D interactive plot available at the above-mentioned website. The 3D plot in Figure 3 is divided into eight 2-dimension plots in Figure 4 so that one can evaluate the distribution of one to eight significant gene-sets across the regions with COL11A1 expressed at high CAF level.

Finally, Figure 5 shows the combined significance ( C S ) heatmap of the 8 hallmarks for COL11A1 phenotype at a high CAF level. Higher values of C S in Figure 5 represent locations with a combined significant gene-set.

Observations of heterogeneity of cell subpopulations in a tumor and the complex interplay of functions involved in the diverse morphological and phenotypic profiles of cancers have a long history. Even in the 19th century, pleomorphism of cancer cells within tumors was observed by the “father of modern pathology”, Rudolph Virchow. More recently, in the 1970s, G.H. Heppner, I.J. Fidler and others showed the existence of distinct subpopulations of cancer cells in tumors, which differed in terms of their tumorigenicity, resistance to treatment, and potential to metastasize. Heppner reviewed the concept of tumor heterogeneity in 1984, and recognized cancers as being composed of multiple subpopulations (Heppner and Miller, 1983), which leads to heterogeneity of cellular morphology, gene expression, metabolism, motility, proliferation, etc (Marusyk et al., 2020). Importantly, ITH has been shown to be associated with poor outcome and decreased response to cancer treatment multiple human cancer types implying a general role in therapeutic resistance (Landau et al., 2013; Patel et al., 2014; Zhang et al., 2014).

The past decade has revealed the immense potential of immunotherapy in cancer. Therapies that promote anti-tumor immune responses have resulted in marked and durable responses in subsets of patients in several cancers (Egen et al., 2020). For instance, abundance of tumor-infiltrating lymphocytes (TILs) and absence of lymphovascular invasion were found to be useful prognostic factors for disease-free survival in patients with HR-/HER2+ breast cancer who were treated using adjuvant trastuzumab (Lee et al., 2015). Spatial transcriptomic approach (Ståhl et al., 2016) was used to identify a type I interferon response overlapping with regions of T Cell and macrophage subset co-localization in HER2+ breast tumors (Andersson et al., 2021). To address the complex interplay between different molecular backgrounds that can characterize ITH with spatial precision, we used GWLCT at a given location in the tissue space to test for the association between a phenotype of interest and different selected gene-sets. The C S score to summarize the overall significance of such associations is computed and visualized with a spatial heatmap.

Gene-set analysis (GSA) is a well-established methodological approach in bioinformatics to test for significant regulation of a selected collection of genes across given samples that represent distinct outcomes. At the level of single cells, GSA could be extended to samples that are individual cells which admit to different phenotypes of interest (Dinu et al., 2021). Furthermore, for spatial single cell analysis, such phenotypes would ideally have a spatially correlated and continuous representation. The gene-sets used in GSA are typically curated based on existing experimentally obtained knowledge of genes and their involvements in molecular pathways. In the present study, for illustrative purposes, we selected a collection of 8 gene-sets that represent certain distinctive hallmarks of cancer (Hanahan, 2022). To test for their enrichment in relevant intratumor contexts, we selected 5 different CAF phenotypes of interest since CAFs are well-known for their contribution to heterogeneity and plasticity in the tumor microenvironment (Ping et al., 2021).

The usual methods for GSA involve one of the two major approaches: (a) competitive, which examines if the correlation of a gene-set with the phenotype is the same as the other gene-sets, and (b) self-contained hypothesis, which investigates if the expression of a gene-set changes by the experimental condition. Our LCT method belongs to the former approach which is more likely than traditional methods to detect the regulation of a functional process or biological pathway that is significantly associated with the gene expression results of a given SCA experiment (Dinu et al., 2021). Interestingly, LCT also extends to longitudinal (Khodayari Moez et al., 2019), multivariate and continuous outcomes (Wang et al., 2014), which are capabilities that we built upon here for providing more accurate representation of single cell level stochasticity of the transcriptomic behavior than that of the univariate and discrete class labels typically used in traditional bulk sample studies.

Simulations, along with real omic data analysis, have served as a powerful and effective tool for establishing the performance of new GSA methods. Past studies have thus used simulation for comparative analysis of different criteria of performance of LCT and other major GSA methods. It was found that LCT has type I error and power that are comparable to MANOVA-GSA (Wang et al., 2014) and superior to SAM-GS (Dinu et al., 2007), particularly at higher magnitudes of correlation values across gene-sets (as is commonly noted during GSA). In terms of computational efficiency, LCT outperformed both methods. In another simulation study, LCT also outperformed GSEA (Moez et al., 2018). Along this direction, therefore, in the present study, we conducted a large number of simulations to compare the performance of GWLCT against multiple known GSA techniques based on a variety of well-defined criteria under different experimental assumptions (or scenarios). Interestingly, statistical power did not change with a variation in set size, number of coordination points or bandwidth, or total number of genes in the simulation dataset, for any of the methods considered, which represent desirable properties for sound GSA methods. Larger spatial correlation across gene expression measurements and between genes and phenotype are key aspects of improved statistical power across our simulation experiments.

In the present study, we introduced GWLCT as a new computational platform that presents a fusion of ideas from spatial data analysis (GWR) and bioinformatics (GSA). We understand that the dual modes—both spatial and single-cell—in which GWLCT provides a joint extension to other GSA approaches places it in a unique category thus making it difficult to compare with the existing methods. Yet, we conducted extensive simulation studies which revealed better performance of GWLCT based on several criteria as compared to many known GSA methods that work either on bulk transcriptomics for different scenarios or aspatial version of single-cell transcriptomics. In particular, the use of multiple different kernels and flexible (adaptive) choice of corresponding bandwidths for geographical weighting allows the linear combination test to test for local associations between selected gene-sets and phenotypic contexts within a tumor sample. Thus, GWLCT provides a novel spatial version of gene-set analysis using high-resolution spatial scRNA-seq data. It does have some limitations that will be addressed in our future work. For instance, as it is difficult to determine a priori the precise spatial scale at which a gene-set or pathway may be regulated in a given phenotypic context, new multi-scale geographical weighting techniques (Fotheringham et al., 2017) may prove to be useful. We will also extend GWLCT to other omic data as we have previously demonstrated with LCT (Khodayari Moez et al., 2019). As spatial single cell omic platforms become increasingly popular, GWLCT will enrich the ongoing efforts in this rapidly emerging area of research (Zhang et al., 2014; Eng et al., 2019; Hajihosseini et al., 2022). Clinical verification of such new analytical methods will require follow-up studies that must be systematically designed for that specific purpose.