MatriCom, a single-cell RNA-sequencing data mining tool to infer cell–extracellular matrix interactions Open Access

The extracellular matrix (ECM) is a complex and dynamic meshwork of proteins that forms the framework of all multicellular organisms (Karamanos et al., 2021; Naba, 2024). We have previously defined the ‘matrisome’ as the compendium of genes encoding the structural components of the ECM (‘core matrisome’, e.g. collagens, fibrillar glycoproteins, proteoglycans) and associated modulatory components (e.g. proteases or cross-linking enzymes that modulate ECM architecture, growth factors or morphogens that bind core ECM components) (Naba et al., 2012). This compendium comprises ∼1000 genes in the mammalian genomes (Gebauer and Naba, 2020; Hynes and Naba, 2012).

Importantly, protein–protein interactions are critical for ECM protein functions. Indeed, matrisome–matrisome protein interactions that can occur intracellularly in the biosynthetic and secretory pathway or extracellularly are essential for the proper assembly of the ECM meshwork (Naba, 2024). In addition, ECM proteins exert signaling functions through interactions with transmembrane receptors. These interactions trigger molecular cascades orchestrating most cellular phenotypes, including cell proliferation and survival, stemness, differentiation and migration (Hastings et al., 2019; Hynes, 2009; Pally and Naba, 2024). As a result, alterations of ECM protein interactions compromise the structural integrity and signaling functions of the ECM and have significant consequences (Naba, 2024; Theocharis et al., 2019), as illustrated by the plethora of diseases arising from ECM gene variants (Lamandé and Bateman, 2020) or observations made in the context of cancer (Cox, 2021) and fibrosis (Walraven and Hinz, 2018). Yet, as of now, we lack the tools to probe these interactions with high throughput in the context of a tissue microenvironment and ask fundamental questions pertaining to ECM biology. For example, although we and others have previously shown that matrisome genes present cell-type-specific and organ-specific expression patterns (Nieuwenhuis et al., 2021; Tellman et al., 2022), we have yet to determine which cell populations within a tissue or organ contribute to the building of the ECM meshwork and whether similar populations do so across organs. Moreover, recent proteomic studies have revealed that the ECM of any given tissue is composed of ∼100–200 distinct proteins and ∼25% of these proteins are ECM regulators, i.e. enzymes involved in ECM remodeling (Naba, 2023; Shao et al., 2023), yet we do not know which cell populations express these regulators and contribute to the remodeling of the assembled ECM meshwork. Last, we have yet to build a comprehensive atlas of pairs of sender and receiver populations of ECM signals to map communications in a tissue microenvironment precisely.

The past decade has seen an acceleration of atlasing efforts aimed at mapping biomolecules with high throughput and at the single-cell resolution in health and disease conditions (Ando et al., 2020). This acceleration was made possible thanks to the development of accessible technologies such as single-cell RNA sequencing (scRNA-Seq) and also by the establishment of large consortia, including those supported by the National Institutes of Health (NIH) in the United States, such as the Human Biomolecular Atlas Program (HuBMAP), focused on the profiling of biomolecules in adult healthy tissues (Jain et al., 2023); the Cellular Senescence Network (SenNet), focused on the profiling of senescent cells across human organs (SenNet Consortium, 2022); the Human Tumor Atlas Network (HTAN), focused on cancers (Rozenblatt-Rosen et al., 2020); or organ-specific efforts such as the Kidney Precision Medicine Project (de Boer et al., 2021) or LungMap (Gaddis et al., 2024). In conjunction with these atlasing efforts, tools have been developed to leverage the wealth of data generated to infer protein–protein interactions and signaling events from scRNA-Seq (i.e. gene expression) datasets. Such tools include CellChat (Jin et al., 2021, 2024), NicheNet (Browaeys et al., 2020), CellTalkDB (Shao et al., 2021) and SingleCellSignalR (Cabello-Aguilar et al., 2020). These tools broadly focus on ligand–receptor interactions and rely on custom interaction databases. Interrogation of the interaction databases supporting these tools has revealed a significant under-representation of matrisome interactions, especially with regard to the core matrisome. For example, CellChatDB lists 4870 interactions involving human proteins, of which 2344 involve at least one matrisome component, but only 111 (∼2.2%) are between two matrisome components, and none are between two core matrisome components. Although ∼40% of the 12,659 interactions reported in NicheNet involve at least one matrisome component, only 3% involve one core matrisome component. Similarly, CellTalkDB compiles 3398 interactions, of which 1877 include a matrisome gene either as ligand or receptor, but only 458 (∼13%) involve at least one core matrisome component, 98 (∼2.9%) are between two matrisome components, and only six interactions are between two core matrisome components.

Importantly, matrisome–matrisome protein interactions obey certain rules. For example, functional collagens, the most abundant core ECM proteins, are triple helical proteins. Although 44 genes encode α chains of collagens in the human genome, the assembly of these chains into functional trimers is highly selective and results in 28 functional collagen proteins. Collagen protomers, formed of three collagen α chains, must assemble intracellularly before being secreted extracellularly (Naba, 2024; Ricard-Blum, 2011). Similar rules govern the assembly of trimers of laminins that are structural components of the basement membrane ECM and any other multimeric core matrisome proteins (Naba, 2024). These examples demonstrate that stipulations need to be imposed when inferring protein interactions from scRNA-Seq datasets: (1) that not all collagen chains can assemble with each other and (2) that genes encoding the individual components of functional protomers should be expressed by the same cell and not by different cells. Of note, selective chain assembly does not only apply to collagens but also to trimeric laminins (12 genes in the human genome), key structural components of basement membrane ECMs (Hohenester, 2019) or to the largest class of ECM receptor, the integrin heterodimers (Hynes, 2002).

To account for the specificities governing interactions involving matrisome components, we developed MatriCom, a web application (https://matrinet.shinyapps.io/matricom) and a companion R package (https://github.com/Izzilab/MatriCom), devised to mine scRNA-Seq dataset and infer ECM–ECM and cell–ECM communication systems. MatriCom relies on a unique database, MatriComDB, that includes over 25,000 curated interactions involving matrisome components sourced from seven interaction databases, resulting in data on 80% of the matrisome. In addition, MatriCom offers a model-maximization functionality to account for rules governing ECM interactions (e.g. constraints on homocellular versus heterocellular co-expressions or selective partner assemblies) and a ranking of the reliability of the interactions. We then illustrate the usability of MatriCom through two examples using publicly available scRNA-Seq datasets: (1) the building of the kidney matrisome communication system and (2) the integration of scRNA-Seq datasets from 46 datasets leading to the identification of ubiquitous and tissue-specific ECM communication patterns. These two examples aim to highlight how MatriCom can be used as a hypothesis-generating resource to gain novel insights into the roles of different cell populations in the building of the ECM and cell–ECM interactions and their dysregulations in the context of diseases such as cancer or fibrosis.

To build MatriComDB, we sourced interactions from seven databases, including two ECM-focused databases – MatrixDB and MatrixDB-curated interactions of the IMEX database, and basement membraneBASE – and more generalist interaction databases – Kyoto Encyclopedia of Genes and Genomes (KEGG), BioGRID, STRING and OmniPath (Fig. 1A, Table 1; Fig. S1). In aggregate, MatriComDB comprises 26,571 unique interactions involving at least one matrisome component. Only 8.2% of these interactions are listed by more than one source database (Fig. 1A, right; Table S1A), highlighting the importance of having sourced interactions from multiple databases. Of the 26,571 interactions listed in MatriComDB, over 20,000 are between a matrisome protein and a non-matrisome protein and 6373 (∼24%) are between two matrisome proteins (Fig. 1B; Table S1A). The interactions between a matrisome protein and a non-matrisome protein include interactions with intracellular proteins (9390 interactions), proteins found at the cell surface (6518 interactions) and proteins found in the extracellular space but not the ECM (4290 interactions) (Fig. 1C; Table S1A). Of the 6373 interactions engaging two matrisome proteins, 1442 involve two core matrisome proteins, 1309 involve a core matrisome and a matrisome-associated protein (examples of these interactions include interactions between a collagen substrate and a collagen-crosslinking enzyme such as lysyl oxidase), and 3622 are between two matrisome-associated proteins (Fig. 1D; Table S1A). The contribution of the seven source databases to these different types of interactions included in MatriComDB is provided (Fig. S2).

With the matrisome in focus, MatriComDB includes interactions for the proteins encoded by nearly 6000 genes, of which 995 are matrisome genes (Fig. 1E, left; Table S1B), covering 81% of the in-silico-predicted matrisome and all categories of matrisome components, including structural core matrisome proteins, such as collagens, proteoglycans and ECM glycoproteins, as well as matrisome-associated proteins, such as ECM-affiliated proteins, ECM regulators and ECM-bound secreted factors (Fig. 1E, right; Table S1B).

MatriCom is made available to users as an online web application (available at https://matrinet.shinyapps.io/matricom/) and an offline version served through the MatriCom package (https://github.com/Izzilab/MatriCom). The only differences between the two versions are the availability of preprocessed open-access (OA) sample datasets, which are only available on the online version, and the size of permitted uploads (maximum 1 GB online, unrestricted offline). We will use the online version to illustrate MatriCom functionalities throughout this paper.

Users can upload their own datasets from Seurat (Hao et al., 2021) or Bioconductor's Single Cell Experiment (SCE) pipelines in .rds or .qs format, or from ScanPy/Loom in .h5ad format via the ‘Data Input section’. In the online version, users can also select a dataset from a list of OA sample datasets available from the Tabula Sapiens OA collection (24 healthy human organ datasets; Table S2A), The Human Protein Atlas (THPA) (31 healthy human tissue datasets; Table S2B) and other studies (four organs; Table S2C).

Users can select the metadata column to be used to label cell identities from their dataset using the first dropdown menu of the Query Parameters tab. For OA sample datasets, both original and standardized cell identities are available (see Materials and Methods).

Once cell identity labels are set, users should choose the proportion of cells (per label) that express any gene involved in matrisome communication and the minimum mean level of gene expression to be met to consider a gene as ‘expressed’. Of note, lower values encourage more granular results at the expense of computational time, while higher values select for the most abundant and more stable matrisome features of a dataset but mask more subtle, and potentially important, communications systems in less abundant cell types.

MatriCom implements a set of filters to account for specific features of ECM biology that might affect the precise inference of communication systems. In particular, we have devised options to (1) maximize the output of a MatriCom model by only listing the most reliable communicating pairs in case of duplicates identified in multiple source databases of different reliability; (2) exclude ‘impossible’ pairs between different cells, such as collagen protomers which strictly assemble within the same cells (Naba, 2024); (3) remove pairs where both elements are the same gene (e.g. the same collagen chains), as the likeliness of these would inflect on the recovery of heterotypical ones; and (4) filter the results based on the reliability score of each pair, the type of communication (between the same or different cell type, or both), and the subcellular location of the gene product (matrisome, cell surface, extracellular and intracellular space). Of note, these filters are activated by default but can be deactivated upon relaunching a MatriCom run.

At completion of the analytical operations and initial filtering, a rich output is returned to users in both graphical and tabular formats in a few seconds to minutes, depending on the dataset size. To facilitate navigation, results are grouped into tabs related to a specific functional part of the output:

This tab provides users with a ‘Global Communication Cluster Map’ (Fig. 2B) derived from the dataset queried. The map represented as a bubble plot reports the number of matrisome communication pairs for each cell population pair. The map is interactive: users can click on any bubble to restrict the results and downstream analyses to a particular population pair.

In the same tab, users will have access to a bar graph summarizing the number of matrisome interactions (‘matrisome–matrisome’ or ‘matrisome–non-matrisome’; Fig. 2C) and a bubble plot reporting the number of matrisome interactions involving the different categories of matrisome components (Fig. 2D). Last, in this tab, users have access to the full output in a tabular format. The table lists, for each communication pair, the partners involved (gene symbols), the type of communication (‘matrisome–matrisome’ or ‘matrisome–non-matrisome’), the identities of the cell populations involved, the proportion of each cell population expressing the genes of the communicating pair, and the mean gene expression levels of the two partners (Fig. 2E).

In this tab, MatriCom results are processed into a non-directed, acyclic graph, and the relative importance of any graph node (any gene) in directing ‘traffic’ towards its neighboring nodes is calculated, thus enabling users to identify the most noticeable (and potentially more actionable or interesting) genes within their results (Fig. S3A).

In this tab, MatriCom provides users with ad hoc enrichment analysis against a set of manually curated matrisome-specific signatures (Fig. S3B).

To illustrate the usability of MatriCom, we reanalyzed a previously published scRNA-Seq dataset of the healthy adult human kidney (https://www.ebi.ac.uk/biostudies/studies/S-SUBS7) using the default query parameters, settings and post-run filters. MatriCom analysis returns a total of 12,528 matrisome communications established by 793 distinct pairs established between the 33 cell types represented in the original sample dataset (Fig. 3A; Table S3A). The full MatriCom output for this dataset with original cell identity annotation, including network influencer and enrichment analyses, are provided in Table S3A–C. To demonstrate the cell annotation functionality embedded in MatriCom, we also provide users with the same dataset and query parameters but generated by selecting the ‘Census’ reannotation option (Table S3D–F).

Communications between genes expressed by the same population – i.e. homocellular pairs – account for only 6.5% of the full network, while most communications, 93.5%, are established by heterocellular pairs (Table S3A). Analysis of the full communication network reveals that non-matrisome–matrisome pairs comprise the majority (78.4%) of all communications pairs (Fig. 3A). Of the pairs involving two matrisome components, only 333 (2.7%) involve two core matrisome components (Fig. 3A). Yet we know that core matrisome–core matrisome component interactions are foundational to ECM assembly. This observation reflects, to a certain extent, the content of MatriComDB, which is biased toward non-matrisome and matrisome-associated proteins and highlights the need to develop tools to identify interactions between core matrisome components (see Discussion).

To determine the degree to which individual populations contribute to the kidney matrisome communication network, we determined a number of expressions per population (Fig. S4A), which we defined as the number of times a given population appears in the MatriCom communication network table as either Population1 or Population2. Together, seven populations make up over 50% of all expressions: distinct proximal tubule 1, connecting tubule, epithelial progenitor cell, pelvic epithelium, principal cell and fibroblast (Fig. 3B; Table S4A). Because cell count in the original dataset varies between populations (Fig. S4B), we corrected for population size and found that, despite comprising less than 0.25% of the original sample dataset, fibroblasts are the largest contributors to the kidney matrisome communication network (Fig. 3B; Table S4A). Interestingly, the relative contribution of expression by fibroblasts is nearly twice that of distinct proximal tubule 1 cells, the next largest contributor (Fig. 3B; Table S4A). Because fibroblasts are the primary cell type responsible for ECM deposition in most tissues (Naba, 2024), we next sought to focus on matrisome communications established by fibroblasts.

We found that 1663 communication pairs (corresponding to ∼13% of the kidney matrisome communication network) involved fibroblasts (Fig. 3C; Table S4B). These include 618 (37.2%) matrisome–matrisome pairs and 1045 (62.8%) non-matrisome–matrisome pairs (Fig. 3C). Heterocellular pairs account for over 96% of the fibroblast network, such that less than 4% of communications are between fibroblasts (Table S4B). Only 143 (8.6%) communication pairs involve two core matrisome components (Fig. 3C; Table S4C–E, Fig. S4C), yet the proportion of communications between two core matrisome genes in the fibroblast network is over threefold higher than that of the full kidney network (Fig. 3A,C).

To evaluate the contribution of the different cell populations to the fibroblast-specific matrisome communication network, we computed all population pairs and determined the combined (e.g. fibroblast–partner and partner–fibroblast) communication frequency per partner (Table S4C). We found that the top five contributing populations – distinct proximal tubule 1, principal cell, connecting tubule, pelvic epithelium and epithelial progenitor cell – establish more non-matrisome–matrisome communications with fibroblasts than matrisome–matrisome communications (Fig. 3D; Table S4C, Fig. S4D).

Overall, nearly half (46.6%) of all fibroblast communications involve at least one core matrisome gene, including 41.9% of matrisome-matrisome pairs and 81.8% of matrisome–surfaceome pairs (Table S4B–D). Therefore, we further sought to explore the subnetwork of communication pairs between ECM and ECM receptors involving fibroblasts. To do so, we extracted from the global fibroblast network the communication pairs represented in the KEGG ECM–receptor interaction pathway (hsa04512) for which fibroblasts are the ‘sender’ population – i.e. express the genes encoding ECM ligands (Table S5A). The fibroblast ECM-receptor subnetwork comprises 408 communication pairs. Of these, 324 (79.4%) are non-matrisome–matrisome pairs and 84 (20.6%) are matrisome–matrisome pairs (Table S5).

Inspection of the subnetwork revealed that ten unique ligand genes and 16 unique receptor genes account for all pairs of the fibroblast ECM–receptor communication subnetwork (Fig. 4A,B; Table S5B). All ten ligand genes belong to the core matrisome, including five ECM glycoprotein genes (FN1, THBS1, SPP1, LAMB2 and LAMC3) and five collagen genes (COL1A1, COL1A2, COL6A1, COL6A2, COL6A3) (Fig. 4A; Table S5B). Of the 16 receptor genes, two are the ECM-affiliated transmembrane proteoglycans syndecans 1 (SDC1) and 4 (SDC4) (Fig. 4B; Table S5B). Communications involving genes encoding integrins, the main class of ECM receptors (Campbell and Humphries, 2011), comprise 62.0% of the fibroblast ECM–receptor network, while those involving other cell surface receptor genes make up the other 38.0% (Table S5B). ITGB1, which encodes the β1 integrin subunit, participates in 130 (31.9%) fibroblast ECM–receptor communications. ITGAV, which encodes the αV integrin subunit, is involved in 42 of the 408 pairs, constituting a little over 10% of all fibroblast ECM–receptor communications (Fig. 4B). Additionally, CD44, SDC1 and SDC4 account for 56 (13.7%), 30 (7.4%) and 54 (13.3%) pairs, respectively (Fig. 4B).

We next evaluated the contribution of fibroblast–receiver population pairs to the fibroblast ECM–receptor communication network (Table S5C; Fig. 4C) and found that distinct proximal tubule 1 cells are the largest contributors, followed by principal cells, epithelial progenitor cells, pelvic epithelium and fibroblasts (Fig. 4C). Homocellular pairs contribute only collagen–non-matrisome and ECM glycoprotein–non-matrisome communications, such that no expression of SDC1 or SDC4 receptor genes by fibroblasts within the ECM–receptor network was detected by MatriCom analysis at the default thresholds (Fig. 4C).

Notably, fibroblasts express all three genes – COL6A1, COL6A2 and COL6A3 – of the [α1(VI)α2(VI)α3(VI)] alpha-chain trimer, the functional protomer that assembles extracellularly to form collagen VI beaded filaments at the interface of the basement membrane ECM and the interstitial ECM (Cescon et al., 2015; Ricard-Blum, 2011). Together, COL6A1, COL6A2 and COL6A3 participate in over half, 52.%, of all fibroblast ECM–receptor communications returned by MatriCom (Fig. 4A). Of the 216 communications involving one of these three collagen VI genes, integrins make up 135 (62.5%) of receptor genes, while non-integrin cell surface receptor genes participate in the other 81 (37.5%) communication pairs (Table S5D). Of the kidney tubule cell populations, distinct proximal tubule 1 cells express three distinct sets of integrin genes encoding heterodimers that interact with collagen VI trimers – ITGA1–ITGB1, ITGA3–ITGB1 and ITGAV–ITGB8 – while connecting tubule and principal cells are the only populations with detectable expression of genes encoding the ITGA2–ITGB1 heterodimer, in addition to ITGAV–ITGB8 (Table S5F). Nine of the 12 immune cell populations express the CD44 receptor gene and establish communications with all three collagen VI ligand genes (Table S5E). Altogether, MatriCom results can be leveraged to infer the interaction network of any given matrisome component (here, the collagen VI molecule) with its receptors in the context of any organ's cell population make-up (Fig. 4D).

The ECM is foundational to tissue organization, development, growth, maintenance and repair (Naba, 2024). We thus hypothesized that parts of its communication systems may be under conservative pressure and reused by different organs facing similar biochemical demands. To test this, we investigated the expression of stringent matrisome communication pairs or ‘patterns’ in the entire OA dataset collection made available in MatriCom.

The MatriCom analysis was run with default parameters (mean average gene expression of 1, mean percentage positive population of 30 and reliability of 3), combining and comparing results from ∼680,000 cells across the 24 tissues of the Tabula Sapiens collection and 430,000 cells across the 22 tissues of the Census-reannotated THPA collection (Fig. 5A; Table S6). To capture strong signals of pattern conservation, we then selected a minimal set of pairs that were expressed in at least 50% of all organs and tissues from Tabula Sapiens and that were also present in THPA (Fig. 5B). The 113 resulting patterns involved approximately the same number of genes in the core and associated divisions of the matrisome and non-matrisome (typically, cell receptors) (Fig. 5C), and covered all combinations of matrisome-to-matrisome categories (Fig. 5D).

In keeping with the fundamental contribution of structural cells (fibroblasts, epithelial and endothelial cells) to the total ECM output of tissues, we found that these conserved matrisome patterns largely connected stromal cells to epithelial and endothelial compartments, although all combinations are represented, with a significant set of pairs mediating connections between structural and immune cells and some even impacting on rare cells such as stem cells in Tabula Sapiens (Fig. 5E,F; Table S7A–C).

Previous reports have demonstrated that pairs and even sets of genes are often reused by multiple, different cell types to regulate diverse physiological processes in both health and disease, showcasing the versatility of gene regulatory networks, where specific combinations of genes can be co-opted for distinct functions depending on the cellular context (Fazilaty and Basler, 2023). In line with these reports, we found that the matrisome patterns show high pleiotropy and reuse, as we found them often and in multiple cell types within and across organs (Fig. S6). In particular, we noticed trends that connect the mosaic of pattern diffusion across cells and tissues with high-order functional roles of the matrisome (Fig. 6A,B), with differences in the balance of reuse versus specialization that depend largely on the functions of the given pair. For example, the interaction between collagen VI (COL6A1) and the CD44 receptor (CD44) is fundamental for cell adhesion, migration and survival, and influences processes such as tissue repair, immune response and maintenance of cell homeostasis (Cescon et al., 2015). Thus, it is not surprising to find this pair being reused by all cellular compartments but the stem cells (Fig. S6C). In contrast, the interaction between selectin-P (SELP) and collagen XVIII (COL18A1) (Izzi et al., 2020) is restricted to the endothelial and the immune compartment (Fig. S6C), as expected for the highly specialized roles of both SELP and COL18A1 in vascular biology, angiogenesis and leukocyte adhesion to blood vessel walls (Dupas et al., 2024; Izzi et al., 2020).

Last, we sought to identify transcription factors (TFs) potentially regulating the co-expression of genes involved in matrisome communication pairs and patterns identified above. To do so, we interrogated two independent databases compiling transcription factor targets, TF2DNA and TFTG (Zhou et al., 2024), and identified 56 TFs as potential regulators of matrisome communication patterns, i.e. as capable of regulating both genes of a communication pair and common to both databases (Fig. S7A, Table S8). We observed that multiple TFs already reported to govern ECM transcription in health and disease are among the top hits [e.g. EGR1, the androgen receptor (AR) and FLI1; Fig. S7B]. Further grouping of TFs into TF families (Lambert et al., 2018) revealed enrichment of C2H2 zinc finger TFs (Lambert et al., 2018), basic helix-loop-helix (bHLH) TFs (Jones, 2004) and nuclear receptors (Lambert et al., 2018) (Fig. S7C; Table S8). As observed at the matrisome gene expression level (see above), we also identified a trend towards the functional specialization of TFs with respect to the cell compartment expressing the pairs, in line with the fundamental role of TFs in directing and maintaining cell-of-origin identity. Specifically, we identified three TF clusters (cluster 1, comprising SP3, ETV4 and NR4A1; cluster 2, comprising SMAD4, SOX2 and FOXM1; and cluster 3, comprising ARNT, GATA1, TFAP2C, NFKB1, FOXP3 and HOXB7) driving the expression of genes of different matrisome communications pairs and patterns (Fig. S7D). Interestingly, most of these TFs were already reported to govern ECM expression in healthy tissues and cancers (Chiquet et al., 2009; Izzi et al., 2019; Trojanowska, 2000), suggesting that ECM regulation at the transcriptional level is strongly cell-type-dependent and not greatly affected by network rewiring.

Altogether, our data demonstrate how imputing interactions from expression data can shed light on the mechanisms of conservation and diversification of ECM assemblies and cell–ECM interactions that result in imparting unique functions to each organ.

MatriCom is a one-of-a-kind tool with a framework designed to capture the rules governing the interactions involving ECM and ECM-related proteins with each other, and with cell surface receptors, and infer these interactions from scRNA-Seq datasets. The datasets used to illustrate the usability of MatriCom were generated from healthy tissue samples, but we envision that, applied to datasets generated on disease-relevant samples, MatriCom will help uncover altered matrisome and cell–matrisome communication systems and patterns and help advance our understanding of the roles of the ECM in disease.

MatriCom output relies on the robustness of MatriComDB, a newly developed database of matrisome protein interactions. Although MatriComDB has been meticulously curated from diverse sources, it is worth noting that not all offer stringent experimental validation or satisfy the technical criteria necessary for identifying matrisome interactions. To address the potential for erroneous results, we introduced the concept of source reliability. By default, MatriCom is configured to launch at the highest reliability level, set at 3, ensuring that all identified communication pathways are supported by robust experimental evidence. MatriCom users are therefore encouraged to carefully consider the trade-off between reliability and yield. The assignment of reliability scores to each source reflects the current state of the art in matrisome-focused and broader biological interaction databases and will be subject to future adjustments as sources are updated or as interaction data from additional databases are included. Consequently, interactions deemed to have low reliability today may be reassessed and upgraded in the future. Of note, the current content of MatriComDB, which is biased toward non-matrisome and matrisome-associated proteins, also highlights the need to develop tools to identify interactions between core matrisome proteins, a persisting challenge owing to the high insolubility of core matrisome components (Bains and Naba, 2024; Naba, 2023).

Imputing matrisome communications from scRNA-Seq data presents a significant challenge, one that current computational tools can only partially address. As demonstrated, matrisome gene counts are among the lowest in typical scRNA-Seq datasets (Fig. S5), leading to a higher prevalence of zero counts compared to genes encoding components of the intracellular proteome or of the surfaceome.

Matrisome genes are, therefore, less likely to be comprehensively captured by algorithms that rely on genetic co-expression patterns to represent communications. In other words, matrisome genes will tend to be under-represented in communication analyses. For example, prior knowledge has shown that fibroblasts express a larger number of collagen genes, including the fibrillar collagens III and V (encoded by COL3A1 and COL5A1–3, respectively) not found in our reanalysis of the OA kidney. The dissemination of higher-resolution techniques and velocity assessment should help capture under-represented genes, but, at the moment, the overall low expression of matrisome genes remains a significant challenge preventing a more comprehensive definition of ECM communication networks.

Last, the spatial dimension of physical matrisome interactions remains largely unmodeled and, to date, unmodelable. Unlike soluble ligands, many ECM components are localized to specific regions of the pericellular or intercellular space, stabilized by interactions with cell surface receptors and other matrisome constituents. In addition, and as discussed, some interactions between different chains of ECM protomers can exclusively take place within a given cell and between defined partners. Although the precise localization of these elements is critical for ECM organization (Naba, 2024) and multicellularity (Bich et al., 2019; Hynes, 2012; Rokas, 2008), the actual physical distances between communicating matrisome elements are currently unknown. We anticipate that the broad adoption of methods such as spatially resolved transcriptomics (Marx, 2021) and proteomics (Method of the Year, 2024: spatial proteomics, 2024) will help fill this gap and add a new degree of precision to our understanding of the extracellular space and the ECM.

MatriComDB was assembled by manual curation and compilation of the interactions between matrisome components (‘matrisome–matrisome’) and between matrisome components and cellular components (‘matrisome–non-matrisome’) across seven databases.

To guide curation, we used the latest list of matrisome genes available for download on the Matrisome Project website (https://matrisome.org) across all our tools, such as MatrixDB (Clerc et al., 2019) and Matrisome AnalyzeR (Petrov et al., 2023) to ensure standardization.

The following databases were interrogated in February 2024 (see also Table 1 and Table S1):

• MatrixDB (http://matrixdb.univ-lyon1.fr/) core collection [‘Core MatrixDB dataset’, PSI-MI TAB 2.7 file (Clerc et al., 2019)] augmented with two published interactomes: ‘Interaction network of the pro-peptide of lysyl oxidase’ (Vallet et al., 2020) and ‘Interaction network of the four syndecans’ (Gondelaud and Ricard-Blum, 2019). Interactions were subset to human and mouse (taxid:9606 and taxid:10090, respectively) and to UniProt and gene symbol identifiers only, for a total of 308 binary interactions from the core collection, 329 for the syndecan interactome and 27 for the lysyl oxidase interactome. Finally, we integrated the interactome of membrane collagens (Professor Sylvie Ricard-Blum, University Lyon I, France, personal communication), reporting 111 interactions for a total of 758 unique interactions. These were finally filtered against the matrisome gene list (see above) to ensure that at least one interactor of each pair is a matrisome component, bringing the augmented core MatrixDB collection to a total of 746 interactions.

• MatrixDB – IMEx extended collection [‘IMEx extended MatrixDB dataset’, PSI-MI TAB 2.7 file (Clerc et al., 2019)]. Interactions were subset to human and mouse (taxid:9606 and taxid:10090, respectively) and to UniProt and gene symbol identifiers only and filtered against the matrisome gene list to ensure that at least one interactor of each pair is a matrisome component, for a total of 7367 interactions.

• Basement membraneBASE (https://bmbase.manchester.ac.uk/) (Jayadev et al., 2022). Entries were manually curated to harmonize nomenclature across the dataset (e.g. turning ‘native collagen’ and ‘procollagen’ to ‘collagen’ or turning laminin trimers into their individual constituents) and checked against the matrisome list for a final total of 3534 interactions.

• The KEGG (https://www.genome.jp/kegg/pathway.html) ‘focal adhesion’ (hsa04510) and ‘ECM–receptor interactions’ (hsa04512) pathways (Kanehisa et al., 2023). Gene identifiers were converted to gene symbols using the org.Hs.eg.db package in R and the interactions were run against the matrisome list, which retrieved a combined total of 1179 interactions.

• The OmniPath database was accessed programmatically via the OmnipathR package in R (v3.10.1; https://r.omnipathdb.org/) (Türei et al., 2021), filtered against the key ‘is_stimulation==1’, with multi-subunit complexes simplified by separating each complex into its binary interactions, before the results were filtered against the matrisome gene list for a total of 1650 interactions.

• STRING (https://string-db.org/), physical subnetwork of the whole protein network only. Data were downloaded for human and mouse (taxid:9606 and taxid:10090, respectively), gene identifiers were converted to gene symbols, and the unique results were filtered against the matrisome gene list for a total of 12,191 interactions.

• BioGRID (https://thebiogrid.org/; version 4.4.230), ‘multi-validated physical interaction’ subset only. Data were filtered for human and mouse (taxid:9606 and taxid:10090, respectively) and filtered against the matrisome gene list for a total of 2491 interactions.

The individual datasets were merged into MatriComDB, an omnibus database of 26,571 unique interactions involving at least one matrisome partner. We assigned a reliability score to each resource based on the level of experimental validation of each protein interaction (entire database versus only part of it) and validation of the proposed interactions by independent groups in the ECM community. This led us to assign the maximum reliability score of 3 to interactions retrieved from MatrixDB and KEGG, a score of 2 to basement membraneBASE and MatrixDB-IMEx extended, and 1 to all other sources.

The online and offline versions of the MatriCom Shiny application were built with the R Project for Statistical Computing and Shiny language (https://shiny.rstudio.com/) and share a common set of functions and ‘logic’. The online version is freely accessible at https://matrinet.shinyapps.io/matricom/, while the offline version can be installed from https://github.com/Izzilab/MatriCom. The difference between the online and local (offline) versions of the application is in their intended use: the online version is suitable to handle the analysis of small- to medium-size datasets (maximum upload size 1 GB) and for the retrieval and analysis of large, curated OA sample datasets (see below). The offline version does not offer access to the OA sample datasets but allows the analysis of larger datasets (the size limit is dictated by the system memory of the local machine running the application). Both versions source the same functions for the upload and preprocessing of data, the actual analysis, the filtering of the results according to user actions, and the graphical and tabular reporting of the results.

The MatriCom function is at the core of the application. This function, given an input scRNA-Seq dataset and a set of options (see below), interfaces with MatriComDB (see above) and a custom-built exclusion database and runs the following analysis:

1. all genes expressed at a certain minimum level (set by the user, defaults to 1 in the app) by each cell identity are fetched;

2. the percentage of each cell identity expressing each of the genes meeting the previous criterion is calculated, and only genes for which expression levels are equal or greater than a threshold (set by the user, defaults to 30 in the app) are kept;

3. using MatriComDB, all possible interactions between the same and different cell identities (homocellular and heterocellular co-expressions, respectively) are identified.

The intended file format to be used with MatriCom is the Seurat standard for scRNA-Seq data (Satija et al., 2015). Both the online and local versions of the application can internally convert the ‘AnnData’ format (.h5ad) to a suitable Seurat version. Upon launching, the MatriCom function will represent the Seurat object as version 3 (eventually converting from the current version 5) and search for available assays into the ‘data’ (normalized counts) slot first, then into the ‘counts’ slot, then eventually into the slots meant for integration of multiple datasets (typically, ‘SCT’ and ‘integrated’). MatriCom can also process datasets generated through the SCE pipeline, which will be of use to researchers processing their data via the Bioconductor package.

As customary, MatriCom users will be expected to have preprocessed their data for quality assessment and cell type identification. We also recommend that users provide objects with log- or otherwise normalized data to reduce inter-cellular bias.

Importantly, MatriCom implements a strict session-specific data policy: data uploaded by users are neither stored in our server nor can they leak through sessions. User data are purged upon user disconnection or at session timeout.

Large OA scRNA-Seq datasets have been processed through MatriCom to offer users a streamlined gateway to matrisome communication atlases. In the current release, we offer an interface to Tabula Sapiens (The Tabula Sapiens Consortium et al., 2022), THPA (Karlsson et al., 2021) and additional OA datasets retrieved from Azimuth (Hao et al., 2021) (Table S2). For Tabula Sapiens, we downloaded the full dataset (‘Tabula Sapiens – All Cells’) available through the CZ Biohub, CZI Single-Cell Biology portal (https://cellxgene.cziscience.com/collections). For THPA, the full dataset (‘rna_single_cell_read_count.zip’) was downloaded from the THPA portal (https://www.proteinatlas.org/about/download). Additional OA datasets were downloaded from the ‘References for scRNA-seq Queries’ section of Azimuth (https://azimuth.hubmapconsortium.org/), focusing on human tissues and organs that could easily be integrated with the larger atlases above, namely bone marrow, kidney, lung (V1) and peripheral blood mononuclear cells. In all cases, data were imported with original cell and tissue/organ identities provided in the respective publications. At data upload, metadata columns become available for the users to choose cell identity labels.

To facilitate the comparison of MatriCom outputs across datasets, we are offering users the possibility of having cell populations within their datasets reannotated using the Census package (Ghaddar and De, 2023). We have also used this package to reannotate cell identities in THPA datasets and other OA datasets based on those present in Tabula Sapiens. This resulted in the reannotation of 22 out of 30 THPA datasets (with brain, bronchus, esophagus, fallopian tube, ovary, placenta, spleen and stomach having no matching samples in Tabula Sapiens) and the four other OA datasets available in MatriCom.

Upon completion of the analysis, MatriCom returns a ‘Global Communication Cluster Map’ in the form of a bubble plot reporting the number of communication pairs for each population pair and the full dataset in tabular format (see Results).

Additional outputs from MatriCom include compartmental position and topological and functional enrichments. For compartmental positions, all non-matrisome components are annotated as part of the ‘surfaceome’ (Bausch-Fluck et al., 2018) or as ‘extracellular’ or ‘intracellular’ according to Gene Ontology annotations (The Gene Ontology Consortium, 2019).

At the topological level, we identify influential nodes by converting the results into an acyclic undirected graph and calculating the degree of each node, then subsetting to the nodes with the highest degree (maximum 100, or all the available nodes if the network size is less than 100) and fetching all their first-degree neighbors (i.e. the targets). Then, the contribution (i.e. the influence) of each influencer to the target is calculated as the degree of the influencer divided by the sum of the degrees of all the influencers impinging on the same target and reported as the relative fraction.

At the functional level, we calculate enrichment for 29 matrisome signatures derived from experimental proteomic studies of ECM-enriched samples from a plethora of preclinical and clinical samples previously reported by the Naba laboratory. These signatures are available via the Matrisome Project website (https://matrisome.org) and the Molecular Signature Database (MsigDB; https://www.gsea-msigdb.org/gsea/msigdb) (Liberzon et al., 2011). To calculate the enrichment values, all the matrisome genes found per cell pair in the results are collated and tested against the total number of genes in the matrisome list using a hypergeometric distribution test.

To facilitate the use of MatriCom, we have developed an interactive guided tour accessible from the MatriCom home page. The user interface (UI) also contains multiple help buttons that trigger pop-up windows to assist users in choosing the right options for their analysis. Users in need of support or encountering technical difficulties, can reach us using the ‘Contact Us’ button located on the MatriCom home page.

We will append the content of MatriComDB periodically with new validated interactions involving at least one matrisome component. We will also strive to periodically expand the list of sample datasets as they become available from large sequencing consortia. In particular, future content expansions will aim to include datasets representative of different disease states.

We imported the ‘Kidney’ sample dataset (Stewart et al., 2019) to MatriCom with original annotations, ran analysis using default query parameters, and exported tabular data with all default settings and post-run filters active (Fig. 2A; Table S3).

The code written to perform this analysis is available at https://github.com/Izzilab/MatriCom-analyses/tree/main/CS.

We identified all gene pairs in the fibroblast-specific MatriCom communication network that are represented in the KEGG ECM–receptor interactions list (hsa04512), extracted the subset of ECM–receptor pairs for which fibroblasts are the ‘sender’ population (i.e. express the ligand gene), and determined the proportion of communications involving each unique ligand or receptor gene. Relative contribution of communications to the fibroblast ECM–receptor network per fibroblast–receiver population pair was determined, as described above.

From the subset of fibroblast ECM–receptor communications involving one of three collagen VI ligand genes – COL6A1, COL6A2 or COL6A3 – we extrapolated unique interactions between collagen VI [α1(VI)α2(VI)α3(VI)] trimers and non-integrin cell surface receptors by removing pairs including an integrin receptor gene, concatenating the three collagen VI gene symbols into a shared trimer label and filtering duplicate pairs. For integrin receptors, we separated each collagen VI trimer–integrin heterodimer complex in the KEGG ECM–receptor reference list into its six binary interactions. For each set of integrin genes that encodes a functional heterodimer (e.g. ITGA1–ITGB1), we screened the fibroblast collagen VI–receptor network for gene pairs established between either the α or β subunit gene and one of the three collagen VI α-chain genes. Multi-subunit interactions between [α1(VI)α2(VI)α3(VI)] trimers and integrin heterodimers are only reported if at least one of the six binary interactions was detected in the fibroblast collagen VI–receptor network.

The code written to perform this analysis is available at https://github.com/Izzilab/MatriCom-analyses/tree/main/OA. MatriCom was used to reanalyze open-access scRNA-Seq datasets from Tabula Sapiens (original annotations) and THPA (Census annotation; see above). Tabula Sapiens, having greater depth, was chosen as the primary discovery cohort, with THPA datasets being used to validate findings. Pairs discovered in at least 50% of the available tissues from Tabula Sapiens and present in THPA datasets were termed ‘matrisome communication patterns’. A pair-by-tissue similarity matrix was calculated using Pearson correlation, and exemplary clusters were identified by visual inspection and greedy modularity optimization of the correlation matrix after thresholding for strong correlation (Pearson R>0.7) only.

Enrichment analysis was performed on matrisome co-expression patterns using gprofiler2. Significant results (P<0.05) were further filtered to focus on results from the Reactome database (Milacic et al., 2024) and reported in bubble plots in which the diameter of each bubble is proportional to the antilogarithm base 10 (−log10) of the enrichment P-value.

We thank the members of the Izzi and Naba laboratories for their feedback on MatriCom and the manuscript. We also thank Professor Sylvie Ricard-Blum for feedback on ECM interactions and the construction of MatriComDB.

The peer review history is available online at https://journals.biologists.com/jcs/lookup/doi/10.1242/jcs.263927.reviewer-comments.pdf