Single-cell analysis reveals the pan-cancer invasiveness-associated transition of adipose-derived stromal cells into COL11A1-expressing cancer-associated fibroblasts

During the last ten years, many research results have been referring to a particular type of cancer-associated fibroblasts associated with poor prognosis, invasiveness, metastasis and resistance to therapy in multiple cancer types, characterized by a gene expression signature with prominent presence of genes COL11A1, THBS2 and INHBA. Identifying the underlying biological mechanisms responsible for their creation may facilitate the discovery of targets for potential pan-cancer therapeutics. Using a novel computational approach for single-cell gene expression data analysis identifying the dominant cell populations in a sequence of samples from patients at various stages, we conclude that these fibroblasts are produced by a pan-cancer cellular transition originating from a particular type of adipose-derived stromal cells naturally present in the stromal vascular fraction of normal adipose tissue, having a characteristic gene expression signature. Focusing on a rich pancreatic cancer dataset, we provide a detailed description of the continuous modification of the gene expression profiles of cells as they transition from APOD-expressing adipose-derived stromal cells to COL11A1-expressing cancer-associated fibroblasts, identifying the key genes that participate in this transition. These results also provide an explanation to the well-known fact that the adipose microenvironment contributes to cancer progression.


Introduction
This work investigates, using computational analysis of rich single-cell datasets from many patients, the nature and origin of a particular type of cancer-associated fibroblasts (CAFs) that has been found to be strongly associated with invasiveness, metastasis, resistance to therapy, and poor prognosis, in multiple cancer types. These fibroblasts can be identified by their characteristic signature with prominent presence of collagen COL11A1 and several other coexpressed genes such as THBS2 and INHBA. There are indications that the generation of those CAFs is part of a universal biological process in cancer that plays essential roles in cancer progression. Therefore, the driving vision for this research has been that it may provide testable hypotheses for the development of pan-cancer therapeutics targeting the biological mechanisms responsible for the creation of those CAFs. As described below, to achieve this task we used both established techniques for studying the dynamic changes in gene expression of cells associated with lineages, such as trajectory inference, as well as complementary computational approaches with novel application in single-cell data analysis. These techniques allowed the precise identification of the expression profile of the origin of the underlying cellular transition as a particular cell type of adipose derived stromal/stem cells (ASCs). We also independently validated the presence of those ASCs as naturally occurring, by applying the same computational methods in other available datasets of normal adipose tissue. In the remaining part of this section we provide introductory information about the COL11A1-expressing CAFs, explain the motivation for our choice of computational methods, and provide evidence for their advantages and unique capabilities analyzing the particular data sets that we used.
These CAFs were first identified in 2010 [1] by their cancer stage-associated signature. Specifically, a COL11A1/INHBA/THBS2-expressing gene signature was found to be present only after a particular staging threshold, different in each cancer type, was reached. For example, it only appeared in ovarian cancer of at least stage III; in colon cancer of at least stage II; and in breast cancer of at least invasive stage I (but not in carcinoma in situ). We had observed the striking consistency of that signature across cancer types, which was obvious at that time from bulk microarray data. For example, Table 1 shows the top 15 genes ranked in terms of fold change for three different cancer types (breast [2], ovarian [3], pancreatic [4]) using data provided in papers published independently. The breast cancer data compare invasive ductal carcinoma with ductal carcinoma in situ (supplementary data 3, "up in IDC" of the paper [2]); the ovarian cancer data compare metastatic tissue in the omentum with primary tumor (Table 2 of the paper [3]); and the pancreatic data compare whole tumor tissue with normal pancreatic tissue (Table 1 of the paper [4]). The four genes COL11A1, INHBA, THBS2, COL5A2 appear among the top 15 in all three sets (P = 6×10 −23 by multi-set intersection test [5]). The actual P value is much lower than that, because, in addition to the above overlap, ten additional genes (COL10A1, COL1A1, COL5A1, FAP, FBN1, FN1, LOX, MFAP5, POSTN, SULF1) appear among the top 15 in at least two of the three sets (and are highly ranked in all three sets anyway). This similarity demonstrates that the signature is well-defined and associated with a universal biological mechanism in cancer.
We had also found that gene COL11A1 serves as a proxy of the full signature, in the sense that it is the only gene from which all other genes of the signature are consistently top-ranked in terms of the correlation of their expression with that of COL11A1. Accordingly, we had identified a COL11A1-correlated pan-cancer gene signature, listed in table 4 of [1], which we deposited in the Molecular Signatures Database (MSigDB). We had referred to those CAFs as MAFs ("metastasis-associated fibroblasts"), because their presence suggests that metastasis is imminent. To avoid any inaccurate interpretation of the term as implying that such fibroblasts are markers of metastasis that has occurred already, here we refer to them as "COL11A1expressing CAFs." Since then, many research results were published connecting one of the genes COL11A1, INHBA, THBS2 with poor prognosis, invasiveness, metastasis, or resistance to therapy, in various cancer types [6][7][8][9][10][11][12][13][14][15].
Furthermore, several designated tumor subtypes were identified in individual cancer types as a result of the presence of those pan-cancer CAFs. For example, the top 15 genes distinguishing the ovarian "mesenchymal subtype" according to [16] Table 1. Top 15 ranked genes in terms of fold change (FC) for three different cancer types revealing the signature of the COL11A1-expressing cancer-associated fibroblasts. Shown are the rankings, reported by the authors, for breast, ovarian and pancreatic cancers. We eliminated multiple entries of the same gene (keeping the one that appears first) and dashes. Genes shared in all three cancer types are highlighted in green, while genes appearing twice are highlighted in yellow.

Breast
Ovarian Pancreatic are clearly due to the presence of the COL11A1/INHBA/THBS2-expressing CAFs and therefore these are not cancer-type specific subtype signatures.
To computationally investigate the origin of those CAFs, we reasoned that analysis of rich datasets from single-cell RNA sequencing (scRNA-seq) provides unique opportunities for tracking the trajectories of cell differentiation lineages. There are several single-cell trajectory inference methods [18] performing "trajectory inference analysis," ordering cells along a trajectory based on similarities in expression patterns.
In particular, we identified one exceptionally rich dataset [19] from pancreatic ductal adenocarcinoma, containing gene expression profiles from 24 tumor samples and 11 normal control samples. We found that several among the 24 tumor samples contained populations of cells strongly co-expressing COL11A1, THBS2 and INHBA, while none of the normal samples contained such cells. We also observed that the prominence of this co-expression signature varied significantly among the tumor samples, having only hints of their presence in some of them, suggesting that the corresponding patients were at various stages of the generation of COL11A1-expressing CAFs. This provides an opportunity to perform additional complementary computational analysis by comparing the prevalent fibroblastic cell populations across the tumor samples, and comparing them with those in the normal samples.
Therefore, in this paper we also used attractor analysis (Materials and Methods) in a novel manner for the analysis of rich scRNA-seq data. The unsupervised attractor algorithm [20] iteratively finds co-expression signatures converging to "attractor metagenes" pointing to the core ("heart") of co-expression. Each attractor metagene is defined by a ranked set of genes along with scores determining their corresponding strengths within the signature, so the topranked genes are the most representative of the signature. The attractor algorithm has previously been used successfully for identifying features useful for breast cancer prognosis [21,22]. When applied on single cell data from a sample, it identifies the gene expression profiles of the dominant cell populations in the sample, and the algorithm is designed to ensure that all the top-ranked genes are co-expressed in the same cells. The purpose of the attractor algorithm is not to classify cells into mutually exclusive subsets. Instead, it identifies the genes at the core of co-expression signatures representing cellular populations from single-cell data, and it provides information that cannot be deduced with traditional clustering methods (see Discussion).
When we applied the attractor algorithm separately in each of the normal samples, we identified a set of nearly identical attractor signatures, corresponding to a type of adipose-derived stromal/stem cells (ASCs) naturally present in the stromal vascular fraction (SVF) of normal adipose tissue, expressing a unique characteristic signature containing fibroblastic markers such as LUM and DCN as well as adipose-related genes, such as APOD, CFD and MGP.
When we applied the algorithm in each of the tumor samples, we found a set of signatures that were changing in a remarkably continuous manner across the samples, some of them being very similar to those of the normal samples, while others are similar to the COL11A1based signature. This suggests that the signatures undergo a gradual change as the transition proceeds, starting from the state of the normal ASCs and passing through a continuum of intermediate states. These results were consistent with those found by applying trajectory inference analysis, but they provided additional significant information based on their unique capabilities. Accordingly, this method demonstrated that there is a continuous "ASC to COL11A1-expressing CAF transition." This finding explains the stage association of the COL11A1-expressing signature as resulting from the interaction of tumor cells with the adipose microenvironment: Indeed, adipose tissue is encountered when ovarian cancer cells reach the omentum (stage III); after colon cancer has grown outside the colon (stage II); and in breast cancer from the beginning of the spread (stage I, but not in situ stage 0).
Finally, we validated our results in other cancer types (head and neck, ovarian, lung, breast), suggesting the pan-cancer nature of the ASC to COL11A1-expressing CAF transition.

ASC to COL11A1-expressing CAF transition identified in pancreatic ductal adenocarcinoma (PDAC)
The PDAC dataset [19] consists of 57,530 scRNA-seq profiles from 24 PDAC tumor samples (T1-T24) and 11 normal samples (N1-N11). To find the expression profile of the dominant fibroblastic population in each sample, we applied the attractor algorithm on the set of identified mesenchymal cells (Materials and Methods). All samples (11 normal and 23 tumor samples, excluding sample T20 as it did not contain identified fibroblasts) yielded strong coexpression signatures involving many genes with big overlap among them. Genes LUM, DCN, FBLN1, MMP2, SFRP2 and COL1A2 appear in the top 100 genes in at least 33 out of the 34 samples (S1 Table), revealing a strong similarity shared by all those fibroblastic expression profiles. This strong overlap is consistent with the continuous transition process, as described below.
Dominant fibroblastic population in the normal pancreatic samples is adiposederived. There is a striking similarity among the attractor profiles (Materials and Methods) of the eleven normal pancreatic samples, indicating that they represent a stable and normally occurring cell population. Specifically, there are 12 genes commonly shared among the top 30 genes in the attractors of at least ten of the eleven normal samples (Table 2), of which four genes are shared among all the samples (P = 3×10 −113 by multi-set intersection test [5]). In addition to fibroblastic markers, there are several strongly expressed adipose-related or stemness-related genes in the list, such as APOD, CXCL12, and DPT, revealing that they are ASCs. Consistently, Gene Set Enrichment Analysis (GSEA) of these 12 commonly shared genes identified the most significant enrichment (FDR q value = 2.16 ×10 −19 ) in the "BOQUEST_-STEM_CELL_UP" dataset of genes upregulated in stromal stem cells from adipose tissue versus the non-stem counterparts [23].
To investigate the nature of this ASC population, we referred to recent results from singlecell analysis of general human adipose tissue [24]. We applied the attractor algorithm on the dataset with the single-cell expression profiles of all 26,350 cells taken from the SVF of normal adipose tissue from 25 samples, and compared the identified attractor with the "consensus attractor" (Materials and Methods) of the 11 normal pancreatic samples, which represented the main state of the normal fibroblastic population (Table 3). There are 14 overlapping genes between the top 30 gene lists (P = 10 −33 by hypergeometric test), and most of the nonhighlighted genes in each column are still ranked highly in the other column. This extreme similarity of the two gene expression profiles indicates that they correspond to the same naturally occurring cell population. Furthermore, excluding the general fibroblastic markers LUM and DCN, we found that gene APOD (Apolipoprotein D) has the highest average ranking in Table 3, and is top-ranked in the independently found SVF fibroblastic population of cluster VP4 (supplementary file 20) of [24]. Therefore, we selected APOD as the representative marker for the ASC population.
Establishing the presence of COL11A1-expressing CAFs in PDAC tumor samples. Because COL11A1 serves as proxy of the full signature [1], a reliable test for determining if a sample contains the COL11A1-expressing CAFs is to rank all genes in terms of their association, measured by mutual information (Materials and Methods), with COL11A1 and see if INHBA and THBS2 are top ranked. Indeed, this happens in several tumor samples, as shown in Table 4 for some of them (T23, T11, T6, T15, T18). For each sample, the shown genes are co-expressed in the same cells, because of the high correlations in a single-cell dataset.
Dominant fibroblastic populations in the tumor PDAC samples exhibits a continuous transition from ASCs to COL11A1-expressing CAFs. Based on the selection of APOD as a representative marker for the ASC population as described previously, we rearranged the attractors of the PDAC tumor samples in terms of descending order of the rank of APOD (Table 5) from left to right. There is a remarkable continuity in the shown expression profiles. The samples at the right side of the table include COL11A1 at increasingly high ranks. The intermediate tumor samples shown in the middle have cells expressing genes that are topranked in both the lists on the left as well as on the right. In other words, these cells are in a genuine intermediate state, rather than being a mixture of distinct subtypes (see detailed discussion in Materials and Methods).
Further demonstration of the continuity of the transition. As an additional confirmation of the continuity of the transition (as opposed to the presence of a mixture of distinct Table 2. Top 30 genes of the identified attractors for each pancreatic normal sample (N1-N11). 12 commonly shared genes in at least ten of the eleven normal samples are highlighted.

Table 3. Comparison of the attractors (top 30 genes) identified in the SVF of normal adipose tissue (Dataset 1) and in the normal pancreatic samples (Dataset 2).
Common genes are highlighted in yellow.

PLOS COMPUTATIONAL BIOLOGY
The role of adipose-derived stromal cells in cancer invasiveness To further investigate the continuous transition, we partitioned the 34 pancreatic samples into three groups. Group 1 includes the eleven normal samples (N1 to N11). For tumor samples, we divided the rearranged samples in Table 5 into two groups (Group 2 and Group 3). Group 2 contains all samples to the left of and including T22, so that APOD is ranked before COL11A1 in the attractors of that Group, representing a relatively earlier stage of this transition. We then applied the consensus version of the attractor finding algorithm (Materials and Methods) and identified the signatures representing the main state of the fibroblasts for each of the above three sample groups (Table 6). Although there are many shared genes, the groups have distinct gene rankings. Group 1 (normal samples) contains many adipose-related genes, consistent with Table 2. Group 3 contains, in addition to COL11A1, many among the other CAF genes, such as THBS2, INHBA, AEBP1, MFAP5 and COL10A1. Group 2 displays an intermediate state, including markers of both ASCs as well as CAFs.
To find potential critical genes at the initiation phase of the cellular transition, we focused on the first tumor samples (with highest APOD ranking) in Table 5, so we can compare them with those of the normal ASCs.
We observe that gene SFRP4 stands out, as it appears for the first time remarkably among the top genes in all the first samples T2, T13, T14, T19, ranked 4th, 6th, 4th 8th, respectively. This suggests that the Wnt pathway is involved in the initiation of the cellular transition, because SFRP4 is a Wnt pathway regulator whose expression has been found associated with various cancer types [25,26]. Interestingly, SFRP4 disappears from the list of the attractors, indicating that it is downregulated in the final stage of the transition.
It is also known that gene RARRES1 (aka TIG1) plays an important role in regulating the proliferation and differentiation of ASCs [27]. Consistently, Table 6 reveals that RARRES1 appears for the first time in the attractors of the initial tumor samples. Just like SFRP4, RARRES1 is downregulated in the final stage, related to the fact that it has been suggested as a tumor suppressor [28,29].
We also performed differential expression (DE) analysis comparing the normal samples with the first samples (T2, T13, T14, T19) of Table 6 (Materials and Methods; S2 Table). The results of such DE analysis represent the full population of fibroblasts and not necessarily reflect the expression changes in the particular cells undergoing the ASC to COL11A1 expressing CAF transition. Gene CFD was found to be most downregulated, consistent with the expected downregulation of adipose-related genes as they differentiate into fibroblasts. Genes SFRP4 and RARRES1 are upregulated consistent with their appearance in the attractors.
On the other hand, the top upregulated gene is phospholipase A2 group IIA (PLA2G2A), which is not among the top genes of any attractors we identified, indicating that it is not expressed by cells undergoing the ASC to COL11A1-expressing CAF transition. It probably still plays, however, an important related parallel role and many previous studies referred to its effects on prognosis of multiple cancer types [30][31][32]. The PLA2G2A protein is a member of a family of enzymes catalyzing the hydrolysis of phospholipids into free fatty acid. We hypothesize that this process leads to fatty acid oxidation, which may facilitate metastatic progression. Indeed, it has been recognized that fatty acid oxidation is associated with the final COL11A1expressing stage of the transition [33]. These results suggest that lipid metabolic reprogramming plays an important role in the metastasis-associated biological mechanism [34], by potentially providing energy for the metastasizing tumor cells.
Validation with trajectory inference. We independently applied trajectory inference (TI) analysis on the PDAC fibroblasts by using the Slingshot [35] method in an unsupervised manner. We first performed unsupervised clustering on the identified fibroblasts (Materials and Methods), resulting in four subgroups X1, X2, X3, X4 (S1A Fig) with the top differentially expressed genes shown in S1B Fig. One of these clusters (X4) was discarded from further TI analysis, because it mainly expressed the IL1 CAF marker HAS1 (Hyaluronan Synthase 1), which is not expressed by either ASCs or COL11A1-expressing CAFs (and does not appear at

PLOS COMPUTATIONAL BIOLOGY
The role of adipose-derived stromal cells in cancer invasiveness all in S1 Table), and contained only 3% of fibroblasts resulting almost exclusively from patient T11 (S1C Fig). As seen from the list of top differentially expressed genes of each cluster, X1 contains CAF genes top ranked (including MMP11, COL11A1, THBS2, INHBA), X2 has RARRES1 at the top, and X3 has ASC genes top ranked, including DPT, C7, CXCL12 and CFD. Consistently, S2A and S2B Table shows the top 100 genes with zero P value, ranked by their variances, resulting from pseudotime-based differential gene expression analysis (Materials and Methods). We can clearly identify as top-ranked several ASC genes, as well as CAF genes, while some general fibroblastic markers, such as DCN, are missing, consistent with the continuity of the ASC to COL11A1-expressing CAF transition. We then used a generalized additive model (GAM) fit to pseudotime-ordered expression data to visualize the trend of gene expressions (Fig 2A).
There was a prominent difference between adipose-related genes and COL11A1-associated genes. The expression of the adipose-related genes steadily fell across the process (Fig 2B), while the expression of COL11A1-associated genes gradually increased (Fig 2C). There is a significant negative correlation between these two groups of genes, e.g., COL11A1 (the last among those genes to increase its expression) was exclusively overexpressed in the mature CAFs, which did not express C7. Of particular interest, genes SFRP4 and RARRES1 (Fig 2D) increased consistently at the beginning and then decreased after reaching a peak, suggesting that they may play important roles in the differentiation path.

Validation in other cancer types
Next, we validated the ASC to COL11A1-expressing CAF transition in other solid cancer types. Although we could not find currently available datasets as rich as the PDAC dataset, we selected those containing a large (at least 100) number of fibroblasts and separately analyzed each of them, obtaining consistent results. Specifically, we used four scRNA-seq datasets from head and neck cancer (HNSCC) [36], ovarian cancer [37], lung cancer [38] and breast cancer [39].
The COL11A1-expressing CAF signature has been confirmed to be a pan-cancer signature [40][41][42]. Therefore, the most important validation task would be to confirm the existence of the APOD/CFD/CXCL12/MGP/PTGDS-expressing ASCs as the starting point of the transition, and to also confirm that some samples are at an intermediate stage, expressing genes such as SFRP4, RARRES1 and THBS2, in addition to the core ASC genes, demonstrating that they are at an intermediate stage of the transition.
Head and neck squamous cell carcinoma. For the HNSCC dataset, the authors of the paper presenting the data [36] reported that the cancer-associated fibroblasts in the dataset can be partitioned into two subsets, which they name CAF1 and CAF2. In S5 Table of that paper, the top three differentially expressed genes of the CAF2 group are CFD, APOD and CXCL12, while the full gene list for CAF2 presented in the same S5 Table also includes genes MGP, C3, C7, DPT, PTGDS. This strongly suggests that the partitioning used in the paper was influenced by the presence of an ASC cell subpopulation, identical, or at least very similar to, those discovered in the PDAC. Similarly, the list of differentially expressed genes for CAF1 in S5 Table includes genes INHBA, THBS2, CTHRC1, POSTN, MMP11, COL5A2, COL12A1, suggesting that the identified CAF1 subpopulation was influenced by the presence of differentiated CAFs, which would eventually express COL11A1. Finally, gene RARRES1 also appears among the list of CAF2 genes, suggesting that it was captured among cells which had started the process of ASC to COL11A1-expressing CAF transition.
Among the individual patients, we found that the most prominent case is sample HNSCC28, which contains a rich set of cells undergoing differentiation. Applying the attractor finding algorithm on the fibroblasts of that sample (S5 Table) resulted in genes LUM, APOD, COL6A3, PDGFRA, DCN, and CFD being among the top-ranked, revealing that it represents an ASC population. Furthermore, the presence of genes THBS2, MFAP5 and VCAN in the same attractor reveals that these cells have already started undergoing the transition.
Ovarian cancer. For the ovarian dataset, the clustering results showed two clusters (X6 and X9) expressing COL11A1-associated genes and ASC-associated genes, respectively  Table; Materials and Methods). Among the individual patients, we found that the ones validating our hypotheses most are HG2F and LG2, both of whose datasets, consistently, contain cells from the fatty omental tissue. S5 Table includes the corresponding two attractors identified in the cells of each patient. Among the top ranked genes for HG2F are DCN, LUM, C1S, C7, and C3, but also RARRES1, suggesting that they represent fibroblasts undergoing the transition, while the LG2-based attractor contains highly ranked all three genes COL11A1, INHBA, THBS2.

Potential therapeutic targets inhibiting the invasiveness-associated transition
This work provides opportunities for identifying therapeutic targets inhibiting the cellular transition. For example, targeting of gene MFAP5 was recently found to enhance chemosensitivity in ovarian and pancreatic cancers [43]. Specifically, the author states that "MFAP5 blockade suppresses fibrosis through downregulating of fibrosis-related genes such as COL11A1." Consistently, we found MFAP5 as one of the most highly associated genes with COL11A1 ( Table 4).
As mentioned earlier, genes SFRP4 and RARRES1 are transiently expressed in Group 2 of Table 6, suggesting that they can be investigated for inhibiting the cellular transition.
Of particular interest as potential drivers are noncoding RNAs due to their typical regulatory role. Because the expression of these genes is not accurately captured by scRNA-seq technology, we did a thorough analysis of the full set of The Cancer Genome Atlas (TCGA) pancancer data. For the RNA sequencing and miRNA sequencing dataset of each cancer type, we removed the genes in which more than 50% of the samples have zero counts. Then we performed quantile normalization using the limma package [44] (v3.40.6) on log2 transformed counts. In each of the 33 cancer types, we ranked all protein-coding genes in terms of the association (using the metric of mutual information) of their expression with that of gene COL11A1. We excluded the 11 cancer types (LGG, SKCM, SARC, LAML, PCPG, GBM, TGCT, THYM, ACC, UVM, UCS) in which neither THBS2 nor INHBA was among the 50 top-ranked genes, because of the absence of significant amounts of COL11A1-expressing CAFs in those samples (1st sheet in S9 Table). In each of the remaining 22 cancer types, we then ranked all long noncoding RNAs (lncRNAs) and microRNAs (miRNAs) in terms of their association with COL11A1 (2nd and 3rd sheets in S9 Table). Finally, we did pan-cancer sorting of all lncRNAs and miRNAs in terms of the median rank of all lncRNAs and miRNAs (4th sheet in S9 Table).
We found that LINC01614 represents a particularly promising therapeutic target. It had a perfect score of 1 in the pan-cancer sorting list, being strikingly the top-ranked gene in 14 (BRCA, UCEC, KIRC, HNSC, LUAD, LUAD, LUSC, OV, STAD, ESCA, PAAD, MESO, DLBC, CHOL) out of the 22 cancer types (2nd sheet in S9 Table). In fact, the association of LINC01614 was even higher than that of marker protein-coding gene INHBA. The pan-cancer consensus ranking of protein-coding genes in terms of LINC01614 (5th sheet in S9 Table) corresponds precisely to the COL11A1-expressing CAF signature. These rankings, in which marker genes unique to the original and intermediate stages are missing, indicate that LINC01614 is involved in the very final stage of the creation of the COL11A1-expressing CAFs. Therefore, we hypothesize that therapeutics targeting LINC01614 specifically in patients' CAFs may inhibit the final metastasis-facilitating stage of the transition.

Discussion
Our results indicate that the cancer invasiveness-associated COL11A1-expressing CAFs are produced as a result of the interaction of tumor cells with the adipose microenvironment. Therefore, one contribution of our work is that it provides a potential explanation to the wellknown fact that adipose tissue contributes to the development and progression of cancer [45][46][47].
Another contribution is that it precisely identifies the ASC population, as evidenced by the consistent presence of its marker genes among the top-ranked attractor genes in each of the eleven columns of Table 2. The identification of those particular marker genes (APOD prominent among them) cannot be due to chance, because these were eleven totally independent unbiased experiments, and also because the attractor algorithm applied on the SVF of normal adipose tissue in another independent dataset identified precisely the same genes. This finding could not have been achieved with traditional methods.
There is consensus agreement that CAFs are a promising potential target for optimizing therapeutic strategies against cancer, but such developments are restricted by our current limitations in our understanding of the origin of CAFs and heterogeneity in CAF function [48]. Therefore, there is an urgent need to enhance our understanding of those matters. Our results provide clarity on one important particular component (out of several) of the heterogeneous fibroblast tumor microenvironment. To avoid potential erroneous conclusions after applying bioinformatics algorithms, single-cell data analysis provides an unprecedented capability to validate results, including those resulting from the attractor algorithm, by "seeing" individual cells in color-coded scatter plots, such as the one shown in Fig 1, observing and confirming the presence or absence of distinct populations characterized by the combined presence of particular marker genes.
In particular, there are several published papers relying on the application of clustering algorithms following dimensionality reduction on the particular datasets they use, and concluding that there exist a number of distinct and mutually exclusive CAF subpopulations. These reported fibroblastic subpopulations occasionally have gene expression profiles that are conflicting with each other in significant ways among these publications. Examples include the hC1 and hC0 clusters in [49], the C9 and C10 clusters in [42], the CAF2 and CAF1 clusters in [36], the iCAF and myCAF clusters in [50,51] and the iCAF an mCAF clusters in [52]. A review of such results in pancreatic cancer appears in [53].
As an example of conflicting results, the "iCAFs" identified in [52] have significant differences from those identified in other papers and are, in fact, identical to the normal ASCs ( Fig  3B of [52]) identified in this paper, as evidenced by the list of its differentially expressed genes (PTGDS, LUM, CFD, FBLN1, APOD, DCN, CXCL14, SFRP2, MMP2, all of which appear in Table 3, further validating the ASC signature. Therefore, this identified cluster contains mainly normal cells at the origin of the transition, which should not even be called CAFs. Similarly, a recent single-cell data analysis [54] identified two clusters "touching" each other in a UMAP plot (Fig 2A of [54]), C0 and C3, which are precisely the two endpoints of the ASC to COL11A1-expressing CAF transition. Indeed, as identified in Table S6-1 of [54], C0 cluster has the marker genes APOD, PTGDS, C7, C3, MGP. . . which the attractor algorithm had identified and validated in this paper. On the other hand, the marker genes of cluster C3 are precisely those of the COL11A1-expressing CAFs, in which all three genes COL11A1, INHBA and THBS2 are top-ranked (because the metastatic process was already underway). Importantly, as shown in Fig 2B of [54], the ASC marker genes APOD and PTGDS (top ranked in C0 and unrelated to CAFs) are significantly expressed even in the COL11A1-expressing cluster C3 of the paper, providing further evidence of the presence of intermediate states consistent with the transition-and the separating line between C0 and C3 in the diagram is not generated by any biologically reliable manner, consistent with the continuity.
On the other hand, our work is consistent with, and complementary to the results of [49] focusing on the immunotherapy response, in which the presence of the "TGF-beta CAFs" was inferred by an 11-gene signature consisting of MMP11, COL11A1, C1QTNF3, CTHRC1, COL12A1, COL10A1, COL5A2, THBS2, AEBP1, LRRC15, ITGA11. This population apparently represents the COL11A1-expressing CAF endpoint of the transition, and gene LRRC15 was selected as the representative gene based on the fact that it was found to be the most differentially expressed gene between CAFs and normal tissue fibroblasts in mouse models. Indeed, LRRC15 is a key member of the COL11A1-expressing CAF signature (Table 4 of [1]) and we also found that COL11A1 is the highest associated gene to LRRC15 in the Group 3 PDAC patients.
In our work we used a detailed gene association-based scrutiny of all our results, including numerous color-coded scatter plots, rather than blindly accepting clustering results. We believe that this nontraditional computational methodology, when used on rich single-cell data, represents a paradigm shift in which systems biology alone can be trusted, by itself, for producing reliable results. We hope that our results will give rise to testable hypotheses that could eventually lead to the development of pan-cancer therapeutics inhibiting the ASC to COL11A1-expressing CAF transition.

Data processing and cell identification
We selected the Seurat R toolkit [55] for data processing and cell identification. Seurat implements the entire clustering workflow and has an advantage in speed and scalability to analyze large datasets [56]. We applied the Seurat (v3.1.4) to process the gene expression matrix and characterize the cell type identity for each scRNA-seq dataset. The count matrix was normalized and log transformed by using the NormalizeData function. We selected the 2,000 most variable features and then performed principal component analysis (PCA) followed by applying an unsupervised graph-based clustering approach. We used default parameter settings in all the above steps except that the resolution parameter in the FindCluster function is set to 1.0 to increase the granularity of downstream clustering. To identify differentially expressed genes for each cluster, we used the FindMarkers function. To characterize the identity of mesenchymal cells in each dataset, we made use of the expression of known markers: LUM, DCN, COL1A1 for fibroblasts, and RGS5, ACTA2, PDGFRB and ADIRF for pericytes.
For the smaller-size datasets (ovarian, breast), we performed clustering once on all cells for mesenchymal cell identification. For datasets of larger size (PDAC, HNSCC, lung), we applied 'two-step clustering' to ensure accuracy: The first step was initial clustering within individual samples. Then we excluded samples with very few (< 20) detected fibroblasts and pooled the mesenchymal cells of the remaining samples together for a second clustering, which resulted in the final set of mesenchymal cells for the dataset. For the PDAC dataset, we included an additional step to remove low-quality cells, by retaining cells for which at least one of the corresponding markers had expression levels � 3.

Mutual information
Mutual information (MI) is a general measure of the association between two random variables [57]. We used a spline based estimator [58] to estimate MI values and normalized so the maximum possible value is 1. The MI value is clipped to zero if the Pearson correlation between the two variables is negative. The details of the estimation method are described in the paper introducing the attractor algorithm [20]. We used the getMI or getAllMIWz function implemented in the cafr R package with parameter negateMI = TRUE.

Attractor-based analysis
The attractor algorithm was first proposed for identifying co-expression signatures from bulk expression values in samples [20]. In this study, we use the attractor algorithm for the first time for the purpose of scrutinizing cell populations in single-cell data. Compared to conventional single-cell methods, the attractor algorithm features the unique capability of discovering precise profiles of cell populations, which other methods cannot achieve (see Discussion).
Briefly, the algorithm iteratively finds mutually associated genes from an expression matrix, converging to the core of the co-expression mechanism. The association measure used is the normalized mutual information (as described above), which captures the general relationships (including nonlinear effects) between variables. Using the expression vector corresponding to a seed gene as input, the algorithm converges to an "attractor" in the form of a list of ranked genes, together with scores (ranging from 0 to 1) for each of these genes measuring the strength of the membership of that gene in the signature. It has a characteristic property that using different "attractee" genes belonging to a co-expression signature as seeds leads to the identical attractor.
The attractor algorithm had previously been used to find co-expression signatures in bulk gene expression data, in which case a converged attractor could represent a mixture of contributions from distinct cell subpopulations. When using single-cell data, however, the characteristic genes of particular distinct subpopulations will have high expression values only in the cells from those subpopulations and low values in other cells. These genes will have pairwise positive and large correlations, and therefore they will be highly ranked in attractor signatures representing such individual subpopulations. On the other hand, two characteristic marker genes belonging to two different distinct subpopulations will have reverse-associated expression values across those cells, which will contribute negatively to the overall correlation between these two genes. Only if two genes are co-expressed across individual cells will they appear highly ranked in the same attractor.
For single dataset, we applied the attractor finding algorithm using the findAttractor function implemented in the cafr (v0.312) R package [20] with the general fibroblastic marker gene LUM as seed. Identical results in all samples will be found, with very rare exceptions, if other fibroblastic markers, such as DCN, are used. The exponent (a) was set to different values for scRNA-seq datasets profiled from different protocols. For the analysis of UMI based (e.g. 10x) and full-length-based (e.g. Smart-seq2) datasets, we used a = 3 and a = 5, respectively. To find the consensus attractor for multiple datasets, we applied the consensus version of the attractor finding algorithm as described in [59]. In the consensus version, the association measures between genes are evaluated as the weighted median of the corresponding measures taken from the individual datasets. The weights are proportional to the number of samples included in each individual dataset in log scale.

Trajectory inference (TI) analysis
We selected the Slingshot [35] method for TI analysis, based on its robustness and suggestions made by the benchmarking pipeline dynverse [18]. We used the raw counts as input and followed the Slingshot lineage analysis workflow (v1.4.0). To begin this process, Slingshot chose robustly expressed genes if it has at least 10 cells that have at least 1 read for each. After gene filtering, we proceeded to full quantile normalization. Following diffusion map dimensionality reduction, Gaussian mixture modelling was performed to classify cells, where the number of clusters in the Mclust function was set to 3 based on the fact that there were three clusters in our previous Seurat clustering results. The final step of lineage inference analysis used the slingshot wrapper function in an unsupervised manner. A cluster based minimum spanning tree was subjected to describe the lineage. After analyzing the global lineage structure, we fitted a generalized additive model (GAM) for pseudotime and computed P values. Genes were ranked by P values and variances. After running Slingshot, we identified genes whose expression values significantly vary over the derived pseudotime by using a GAM, allowing us to detect non-linear patterns in gene expression.

Statistical analysis
P value evaluation for overlapping genes from different sets. We applied the hypergeometric test for evaluating the significance of genes shared by different sets. If there are two sets to compare, we used the phyper R function. If there are more than two sets to compare, we employed the multi-set intersection test [5] by applying the cpsets function implemented in the SuperExactTest R package. Regarding the background universe size of genes, we used the total number of genes analyzed in the specific expression matrix. In the case of comparing sets coming from different studies, we used 20,000 as the universe size.
Differential expression analysis. We used a Wilcoxon Rank Sum test by applying the FindMarkers function in Seurat to identify the differentially expressed (DE) genes between fibroblasts of different groups. DE genes with |log fold change| > 0. 25 Table. Ranked genes lists in terms of their association (mutual information) with gene COL11A1 by using the full set of pan-cancer TCGA bulk RNA-seq data. (XLSX)