dataesr/scientific-paragraphs-categorization huggingface dataset:
BulkSignalR approach and design Most single-cell tools infer LRIs by relying only on the ligand L and the receptor R abundance because they have access to separate data for each cell population. With bulk data, the observed abundance of transcripts (or proteins) is the net contribution of different cell types, each of which represents an unknown proportion of the analyzed tissue. Therefore, we hypothesized that LRI inference using bulk data should rely on additional modeling steps. Specifically, we modeled the consequence of a putative LRI by considering the participation of R in biological pathways and the regulation of the genes targeted by these pathways. Furthermore, as we did not assume any knowledge about the different samples, such as clusters harboring similar profiles, we modeled relationships between ligands, receptors, and downstream pathways based on Spearman correlations across the whole dataset.
Consequently, BulkSignalR searches for triples (L, R, pw) where pw is a pathway downstream of R (Figure  1A ). L, R and pw must display correlated activities to be deemed significant by our statistical model. Potential LRIs are taken from the LRdb database  (7) , while pathways originate from Reactome  (18)  or GO biological processes  (19). To identify genes reporting on a pathway activity, which we called target genes, we exploited the topology of a reference network combining Reactome and KEGG  (17)  interactions. BulkSignalR uses a statistical model to assess the significance of all possible triples (L, R, pw) based on the null distributions of L-R and also R-T correlations, pathway sizes, and total number of target genes. Importantly, BulkSignalR combines all correlation values independently from the sample number.
This avoids the issue of very small, but highly significant correlations in large datasets and allows the analysis of small cohorts. BulkSignalR statistical analysis output is summarized in a table that contains putative LRIs, pathways including the receptors, and statistical parameters (Figures  1B  and S2 ). Due to redundancies in pathway definitions, the same interaction can appear in multiple pathways. Moreover, some ligands can have more than one receptor and vice versa. Therefore, we introduced different reduction operations to limit redundancy, or to emphasize the ligand, the receptor, or the pathway (Figures  1C  and S3 ). Such operations can be chained. Reduction to the pathway that yields the best P-value for each LRI results in a table with unique interactions. Using this reduction, we estimated the TP and FP rates in order to obtain pseudo-ROC curves (Figure  1D ).
Since our statistical model allows varying the depth at which target gene correlations are assessed, we explored different choices. Although considering all the target genes might be an option, large pathways may comprise multiple branches that are not all activated. As a compromise, we did not go deeper than the top 45% of the correlation (55th percentile in Figure  1D ). This revealed that it was advantageous to exploit the pathway thoroughly. The L-R correlation on its own did not provide a useful inference mechanism. Indeed, the performance of sole L-R correlation was either limited or did not detect any significant LRI. As explained in Introduction, single-cell scores have been proposed to infer LRIs from bulk data, although they require separate datasets for the ligand-secreting cells and the receptor-expressing cells.
We investigated whether such scores could be adapted when no upfront cell separation is available. To this aim, we selected a salivary duct carcinoma (SDC) dataset (  6 ) that included well-separated clusters of samples with limited intra-cluster variability. The scoring of common immune cell gene signatures (Table  S1 ) revealed such clusters (Figure  1E ). We then compared the clusters to identify LRIs that were enriched in one cluster as a selection mechanism. Due to the necessary presence of both ligand and receptor in the dataset, and the default requirements of BulkSignalR analysis (L-R correlation > 0. 25, minimum four target genes found in pathways of sizes between 5 and 200), we could evaluate 777 LRIs with BulkSignalR. We submitted those 777 LRIs to adapted single-cell scores. We used a score based on ICELLNET (a product of average values) with specific data normalization  (12) , and another one based on CellPhoneDB (P-values of the ligand and receptor average expression by sample random permutations across clusters).
By default, the average of the ligand and the receptor z-scores accounts for one half of the score, while the average of the target gene z-scores included in the statistical model (correlations above chosen percentile) account for the other half with equal individual weights. Gene signatures and scores remain compatible with any combination of reduction operations due to BulkSignalR software design. To facilitate the analysis of non-human datasets, we integrated a generic ortholog mapping mechanism that allows users to virtually work with any species (Figure  1G ). We designed BulkSignalR with the aim of proposing a user-friendly tool. Therefore, only few lines of code, accessible to basic R users, are sufficient to perform a complete analysis and generate informative outputs (Figure  S4 ) including graphical representations (Fig-  ures 1HI  and S5 ).
Potential biases in bulk-identified LRIs The number of LRIs found by BulkSignalR was lower compared to what was obtained when comparable single-cell datasets were available. As expected, the LRIs found in bulk data only were enriched in low-abundance molecules (Figure  S6 ). To better characterize a potential bias introduced by the bulk approach, we generated synthetic bulk data  (27)  from a cohort of 29 primary, untreated breast tumors analyzed in single-cell transcriptomics (28) (Figure  S7A ). LRIs were inferred from single-cell data with SingleCellSignalR and synthetic bulk with BulkSignalR (Figure  S7B ). Results were qualitatively similar to what was reported in Figure  S6. By relaxing BulkSignalR filters we could recover most of the single-cell-inferred LRIs (Figure  S7C ). This indicates that what was missed in bulk was mostly due to signal loss because of averaging in the bulk readout.
Autocrine communications in lung cancer cells For the first application of BulkSignalR, we analyzed bulk transcriptomic data from 206 lung cancer cell lines from DepMap  (23). We reduced BulkSignalR output to the best pathway for each LRI to obtain only unique interactions, and we imposed a FDR <0. 1% (full output in Table  S1 ). We obtained and scored gene signatures (Figure  2A ). Cell lines originating from the two main histologic subgroups, smallcell lung cancer (SCLC) and NSCLC, harbored distinct autocrine communication patterns. Mesothelioma-derived cell lines were close to those of NSCLC-derived cell lines, although they constitute a different entity. This similarity may simply reflect the absence of microenvironment in cultured cell lines. The presence of only two lung carcinoid cell lines in the dataset prevented drawing any conclusion on this subgroup.
Summarizing at the pathway level In the analysis of lung cancer bulk data, we found LRIs that implicated a single receptor or related to a single pathway. BulkSignalR capacity to reduce data can be exploited to investigate a dataset at the pathway level. We illustrate this procedure using transcriptomic data from human brain and heart during development and adult life  (43). We started by reducing the output at the pathway level. To this aim, we pooled together all receptors implicated in a given pathway, followed by pooling of all their ligands. In this way, that pathway was associated with a meta-ligand and a meta-receptor instead of many single interactors (Figure  S3 ): a standard triple (L, R, pw) be- comes ({L 1 , • • • , L N }, {R 1 , • • • , R M }, pw). If two different pathways were associated with the same meta-ligand and -receptor, we performed a reduction to maintain only the best pathway chaining BulkSignalR chained reduction operations (Figure  S3 ).
Relating ligand-receptor interactions to cell types In the single-cell paradigm, LRI inference comes with the knowledge of which cell types express the receptor and the ligand. As such information is not directly accessible in bulk datasets, we implemented an algorithm to predict cell type-LRI associations. This requires scoring a set of cell type gene signatures in all samples. It can be achieved with a simple z-score average-based function provided by BulkSig-nalR, or with more advanced cellular deconvolution tools the output of which can then be imported in BulkSignalR. By comparing cell type and LRI gene signature scores (Figure  3A ), we build a sparse linear model in which LRI activity is linked to cell type abundances by the following equation: LRI l activity I l ∼ = j ∈J α j CT j for a small set J of cell types CT j.
See Materials and Methods for details. Due to the intrinsic limitations of bulk data and the absence of complete LRI reference, BulkSignalR cannot predict whether a given cell type expresses the ligand, the receptor, or both. Using the SDC dataset, we scored cell type signatures for common tumor-infiltrating cell types: immune cells, endothelial cells, and fibroblasts (Table  S1 ). We summarized cell communications by summing the weights of the cell types in all LRI activity models as above, i. e. by summing the α j values (Figure  3B ). We noted stronger communications within the stromal (endothelial cells and fibroblasts) and immune components of the tumor, as expected. By selecting all LRIs that were exclusively associated with stromal cells or with immune cells, we found that they were involved in very relevant pathways (Figure  3B ).
2 and 0. 4. We computed the TP rates (TPR) and true negative rates (TNR) using the small SDC and the large TCGA BRCA cohorts (Figure  3D ). We also generated randomized data that should not be assigned to any cell type by picking LRI signature scores and shuffling 10, 25 or 50% of the values (Figure  3D ). The obtained TNR and TPR showed that the algorithm avoided wrong assignments, but its sensitivity to detect complex relationships with three cell types decreased with increasing noise, especially in the smaller dataset.
Other functionalities and applications Our model of LRI (Figure  1A ) implies significant positive correlations between the receptor and a downstream pathway as well as correlations between ligand and receptor. By default, BulkSignalR does not consider L-R correlations <0. 25. However, it is possible to impose a different minimum, particularly -1, and to find a number of (L, R, pw) triples with strong P-values, but negative L-R correlations. This case is properly handled by our statistical model and suggests that such ligands may have an inhibitory action. We reanalyzed the SDC dataset by imposing FDR <0. 1% and L-R correlation >0. 25 in absolute value. We identified 424 unique LRIs, 361 positive and 63 negative (Table  S1 ). By focusing on the Notch pathway, which is deregulated in many tumors, we found common activators (DLL1, DLL4, DLK1, DLK2, JAG1 and JAG2) with strong P-values and positive correlations with one or several of the four Notch receptors.
We also found several interactors that were negatively correlated with Notch receptors, for instance MFNG (correlation -0. 52 with NOTCH1 and NOTCH2), a glycosyltransferase that modulates Notch activity by modifying O-fucose residues at specific EGF-like domains  (44). DLL3, which can Nucleic Acids Research, 2023 11 inhibit Notch  (45) , also displayed negative correlations (-0. 26 with NOTCH3 and, just below our 0. 25 threshold, -0. 245 with NOTCH4) as well as PSEN1 (-0. 50 with NOTCH1, -0. 27 with NOTCH2 and -0. 41 with NOTCH3) and UBA52 (-0. 51 with NOTCH1). It is worth noting that searching for negatively correlated L-R pairs is much more prone to FP (Figure  S12 ) because L-R databases are dominated by activating LRIs. Therefore, this procedure should be considered as exploratory and the output needs to be experimentally validated.
NicheNet  (46) , a single-cell tool, exploits an integrated molecular interaction network that includes LRIs to relate user-chosen gene sets to ligands that might drive their expression according to the network. The authors tested this functionality using the 100-gene signature proposed by Puram (  5 ) for a partial EMT transcriptional program taking place at the invasive front of head and neck squamous cell carcinoma (HNSC). As NicheNet integrated reference network and our method to link ligands to receptor-pathway target genes are conceptually similar, we use BulkSignalR inferences to provide a similar functionality, but for bulk data. Figure  S13  describes its application to investigate the partial EMT program using the TCGA HNSC bulk transcriptome dataset (n = 500). In Supplementary Information, we describe the analysis of DepMap transcriptomic and proteomic data of breast cancer cell lines (Figures  S14-S17 ), and we use a mouse dataset to illustrate the built-in ortholog mapping functionality (Figures  S18-S20 ).
Investigating spatial transcriptomics data with BulkSignalR We then applied BulkSignalR to several multicellular resolution datasets (Visium). BulkSignalR functions were used unchanged, we only relaxed some parameters to adapt to much reduced data dynamics and the large number of spots (minimum correlation, number of target genes, expression level, and proportion of samples with expression, see Figure  S21 ). ST spots were hence regarded as independent samples and we relied on the correlation analysis over the whole tissue to ensure meaningful LRIs would be identified. Attempts to better exploit spatial proximity are detailed in the comparison section below. We started with a triple-negative breast cancer (TNBC) dataset (47) the overall structure of which is presented in Figure  4A. We identified 224 LRIs (Table  S1  and Figure  S22 ), some of which we briefly present here to show that BulkSignalR identified relevant interactions.
Most LRIs were associated with the stroma or with the invasive cancer tissue. In the stroma, we found LRIs that implicated adhesion molecules, e. g. CAM1, VCAM1, and several integrins, i. e. ITGB1, ITGB2, ITGB7, ITGAL, IT-GAX, as well as immune-related L-R complexes, e. g. B2M-CD3D, B2M-(HLA-F), (HLA-A)-CD3D and IL16-CD4. The presence of activated interactions related to antigen presentation (B2M; HLA-A/B) was consistent with findings by the authors of this dataset  (47). Figure  4B  shows the spatial distribution of the LRI between the metalloproteinase MMP9 and the integrin ITGB2. MMPs are important TME regulators that promote EMT, apoptosis, resis-tance, angiogenesis, and tissue remodeling  (48). Moreover, MMP9 has been associated with aggressive and metastatic breast cancer  (49). In the invasive cancer tissue, we observed the activation of the Notch signaling pathway, a TNBC hallmark (50), through interaction with several ligands, for instance JAG2.
Comparison with existing spatial transcriptomic data analysis packages We identified three recent or widely used tools that offer analysis at the multicellular resolution: stLearn (  14 ), Cell-PhoneDB  (11)  and SpaTalk  (15). We compared these tools with BulkSignalR using three human datasets: the previously used TNBC dataset  (47) , a HER2+ breast cancer dataset  (52)  and a dorsolateral prefrontal cortex dataset  (53). To obtain comparable results, we employed LRdb (our LRI database) with stLearn and SpaTalk. Despite our efforts, we could not replace the CellPhoneDB LRI database. Therefore, we used its native database that combines individual molecules and also complexes. When one interaction involved a complex (or two), we generated all possible pairwise LRIs. Then, we intersected this LRI list with LRdb to be as close as possible to the other three tools.
Nucleic Acids Research, 2023 13 PhoneDB requires the definition of the tissue regions (like in Figure  4A ). They exploit this information to reconstitute cell type-specific transcriptomic profiles, restricted to each region for CellPhoneDB, or overall for SpaTalk. As SpaTalk reports LRIs by specifying the cell types, we selected all LRIs based on their P-values and eliminated redundant LRIs if they were significant in several cell type pairs. CellPhoneDB gives a P-value for each pair of cell types occurring in each region. We took the minimum P-value for each LRI and imposed a minimum Pvalue <0. 1%. The TNBC, cortex, and HER2 + breast tumor datasets included 4895, 3636 and 306 spots, respectively. BulkSig-nalR, stLearn, and CellPhoneDB computing times scaled with the dataset size, while SpaTalk times were very long and difficult to explain (Figure  4D ).
CellPhoneDB was the fastest tool. For the TNBC dataset, the four tools reported heterogeneous numbers of LRIs with limited overlap (Figures  4E  and S25A ). For CellPhoneDB and SpaTalk, we defined the dominant cell type at each spot according to the dataset authors who used single-cell data and a deconvolution software to determine them. For CellPhoneDB, we used the tumor regions defined by the authors and shown in Figure  4A. We obtained the largest number of LRIs with SpaTalk (twice as many as with the other tools), and the smallest number with stLearn (first selection). To relate the numbers of identified LRIs to their quality in the absence of a complete reference, we defined two objective quality criteria. First, we considered that the product of the ligand and the receptor, indicating co-presence at a spot, should display a statistical association with the tissue regions.
Next, we compared the results obtained with the HER2+ breast tumor dataset. Tissue regions and dominant cell types were defined by the dataset authors. The number of significant LRIs identified by each tool was even more heterogeneous than for the TNBC dataset (Figures  4H  and S25B ). SpaTalk and CellPhoneDB found the largest number of interactions. Associations of the L × R product with tissue regions showed that BulkSignalR outperformed significantly SpaTalk and CellPhoneDB and also stLearn (first selection), but not stLearn (second selection). We obtained similar results when we allowed negative L-R correlations in BulkSignalR (Figure  4I ). In terms of L-R correlations (Figure  4J ), BulkSignalR outperformed all the other tools. CellPhoneDB and SpaTalk returned a very large number of negatively correlated LRIs, which might lead to substantial FP rates.
For the last comparison we used a dorsolateral prefrontal cortex dataset for which the authors defined regions, but no dominant cell types at each spot. As CellPhoneDB and SpaTalk compare remote spots anyway, we decided to apply them by using the regions as cell type definitions. stLearn (first selection) and CellPhoneDB returned the smallest lists of LRIs, while the other tools gave comparable numbers of LRIs (Figure  4K ). Inter-tool heterogeneity remained substantial (Figures  4K  and S25C ). Region association was better with BulkSignalR than CellPhoneDB and stLearn (first selection), but not compared with stLearn (second selection) and SpaTalk (Figure  4L ). L-R correlations indicated that BulkSignalR performed better than the other tools. We made an attempt to exploit spatial proximity through local averaging or ligand signal maximization, but it did not result in quality improvement.
Application to colorectal cancer liver metastasis spatial data To experimentally confirm some LRIs identified by BulkSignalR by means of IF, we generated a new ST dataset for four CRC-LM samples. Analysis of this dataset by BulkSignalR gave 173, 251, and 241 unique LRIs for the first three CRC-LM samples, respectively (FDR < 1%, L-R correlation > 0. 02 in absolute value), and only 84 LRIs for the fourth sample (Table  S1 ). As the obtained number of reads per spot was also smaller for the last sample (∼30% less than the mean of the other three samples; data not shown), we excluded it. Application of clustering in Seurat combined with sample analysis by two pathologists allowed defining different regions in each sample (Figure  5A ). A complete analysis of the CRC-LM data with biological interpretation would obviously be out of the scope of this paper.
Spatial configurations of all the LRIs appearing in Figure  5B  are depicted in Figure  S27. We then experimentally validated CDH1-EGFR and CDH1-cMET colocalization in CRCLM1 and CRCLM2 to illustrate BulkSignalR capacity to identify interactions that are not necessarily sample region-specific. CDH1-EGFR and CDH1-cMET interactions have been described in many cell types; however, their relationship remains complex and context-dependent  (57). The CDH1-MET LRI was not identified by BulkSignalR because it was not included in the LRdb database. We assessed EGFR, cMET and CDH1 abundance, and ERK1/2 nuclear localization in the two CRC-LM samples by IF (Figures  6AB  and S28AB ). We used ERK1/2 as downstream reporter of EGFR and cMET signaling since upon stimulation of these receptors, activation of ERK1/2 has been well-documented  (58, (59).
DISCUSSION Invaluable transcriptomic and proteomic datasets are available and continue to be generated with a bulk methodology. Here, we showed that the R library BulkSignalR offers researchers a solution to exploit these datasets to unravel cellular networks. BulkSignalR includes a rich set of functionalities, comparable to the best libraries for single-cell data analysis. The infrastructure that supports BulkSignalR computations to link LRIs and downstream pathways allows data analysis at both the pathway and individual LRI levels. This infrastructure also allows network visualization for relating LRIs to target genes. Although bulk data do not directly convey information about the specific transcriptomes of individual cell populations, we propose a simple machine learning model that can infer what populations are likely to participate in each LRI.
By comparing LRIs inferred by BulkSignalR from several bulk datasets and corresponding single-cell datasets, we found that the single-cell data analysis typically identified 2-3 times more LRIs mostly because of signal suppression in bulk. We also found that LRIs identified in bulk only tended to involve low-abundance ligands and receptors, which reflects the higher sensitivity of bulk technologies. We found no functional category bias in bulk-identified LRIs, differences were driven by sample cohorts. Clearly, single-cell analysis is the approach of choice to map cellular networks when few representative samples are sufficient. In the other cases, bulk approaches can be used for cellular network inference, but with reduced details. The two formats, bulk and single-cell, could be eventually combined as well as data analysis methods in large studies.
By reviewing the literature, we discovered that most existing ST software tools have been developed for single-cell or subcellular resolution data. Some nevertheless claim to be compatible with Visium-type data, for instance CellPhoneDB Nucleic Acids Research, 2023 17  (11) , SpaTalk  (15) , and stLearn  (14). In stLearn, an initial CellPhoneDB application is combined with a local posterior analysis to compare a spot with its neighbors. Conversely, CellPhoneDB alone and SpaTalk consider LRIs between spots that may be far from each other. Although some LRIs have a longer range, many occur at a short range or only at cell contacts. We think that this renders the output of these tools difficult to interpret. Nevertheless, we compared CellPhoneDB, SpaTalk, and stLearn with BulkSig-nalR using three ST datasets.