96291204 1 Patient-specific cell communication networks associate with disease 1 progression in cancer 2 3 4 David L Gibbs1, Boris Aguilar1, VΓ©steinn Thorsson1, Alexander V Ratushny2, Ilya Shmulevich1 5 6 1Institute for Systems Biology, 401 Terry Avenue North, Seattle, WA 98109, USA; 2Bristol-7 Myers Squibb, 400 Dexter Avenue North, Suite 1200, Seattle, WA 98109, USA 8 9 Correspondence: 10 David L Gibbs 11 david.gibbs@isbscience.org 12 13 Abstract 14 15 The maintenance and function of tissues in health and disease depends on cell-cell communication. This 16 work shows how high-level features, representing cell-cell communication, can be defined and used to 17 associate certain signaling 'axes' with clinical outcomes. Using cell-sorted gene expression data, we 18 generated a scaffold of cell-cell interactions and define a probabilistic method for creating per-patient 19 weighted graphs based on gene expression and cell deconvolution results. With this method, we generated 20 over 9,000 graphs for TCGA patient samples, each representing likely channels of intercellular 21 communication in the tumor microenvironment. It was shown that particular edges were strongly 22 associated with disease severity and progression, in terms of survival time and tumor stage. Within 23 individual tumor types, there are predominant cell types and the collection of associated edges were found 24 to be predictive of clinical phenotypes. Additionally, genes associated with differentially weighted edges 25 were enriched in Gene Ontology terms associated with tissue structure and immune response. Code, data, 26 and notebooks are provided to enable the application of this method to any expression dataset 27 (https://github.com/IlyaLab/Pan-Cancer-Cell-Cell-Comm-Net). 28 Keywords 29 Networks, cell communication, immuno-oncology, computational oncology, bioinformatics, systems 30 biology 31 Introduction 32 The maintenance and function of tissues depends on cell-cell communication (Wilson et al., 2000; Haass 33 and Herlyn, 2005). While cell communication can take place through physically binding cell membrane 34 surface proteins, cells also release ligand molecules that diffuse and bind to receptors on other cells 35 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint https://paperpile.com/c/5TES3g/QCquJ https://paperpile.com/c/5TES3g/QCquJ https://paperpile.com/c/5TES3g/P8KBB https://paperpile.com/c/5TES3g/P8KBB https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ 2 (paracrine or endocrine), or even the same cell (autocrine), triggering a signaling cascade that can 36 potentially activate a gene regulatory program (Cameron and Kelvin, 2013; Heldin et al., 2016; Cohen 37 and Nelson, 2018). More generally, a message is sent and received, transferring some information as part 38 of a large network (Frankenstein et al., 2006). Cells communicate in order to coordinate activity, such as, 39 to correctly (and jointly) respond to environmental changes (Song et al., 2019). 40 Altered cellular communication can cause disease, and conversely diseases can alter 41 communication (Wei et al., 2004). Cancer, once thought of as purely a disease of genetics, is now 42 recognized as being enmeshed in complex cellular interactions within the tumor microenvironment 43 (TME) (Trosko and Ruch, 1998). The cell-cell interactions are important for cell differentiation, tumor 44 growth (West and Newton, 2019), and response to therapeutics (Kumar et al., 2018). 45 Between cells, information transfer is directional in nature, where cells produce molecules that 46 are received by the properly paired, and expressed, receptor. There is often a sender and receiver, which 47 makes the cell-cell networks directionally linked by molecules. The dynamics of the signal is greatly 48 important (Fridman et al., 2012, Behar et al., 2013), but unfortunately is difficult to detect in bulk 49 sequencing experiments. One approach to studying cell interactions is through the use of graphical 50 models of communication networks (Morel et al., 2017). By incorporating experimental data, the 51 graphical models can become quantitative, providing predictions that can be tested and used in 52 discovering novel drug targets and developing optimal intervention strategies. 53 In recent work (Thorsson et al., 2018), we developed a method used to identify cellular 54 communication networks at work in the tumor microenvironment. Given a set of samples with a similar 55 tumor microenvironment, the method identified ligands, receptors and cells meeting certain criteria of 56 abundance and concordance within that set of samples. The method was applied to identify networks 57 playing a role within specific tumor types and molecular subtypes and is available as a workflow and 58 interactive module on the iAtlas portal for immuno-oncology (Eddy et al., 2020). 59 In this work, we have combined multiple sources of data with a new probabilistic method for 60 constructing patient-specific cell-cell communication networks (Figure 1). In total, we built networks for 61 9,234 samples in The Cancer Genome Atlas (TCGA), starting from a network of 64 cell types and 1,894 62 ligand-receptor pairs. This is a rich feature set from which to investigate biological alterations in cell 63 communication within the tumor microenvironment. We identified informative network features that are 64 associated with disease progression. The method can be applied to any cancer type, but in this manuscript 65 we focus on a selection of cancer types with very high mortality rates, including pancreatic 66 adenocarcinoma (PAAD), melanoma (SKCM), lung (LUSC), and cancers of the gastrointestinal tract 67 (ESCA, STAD, COAD, READ) (Cancer Genome Atlas Network, 2015). 68 This represents a new method that provides information on possible modes of intercellular 69 signaling in the TME, something that is currently lacking. While there are many methods on gene set 70 scoring, cellular abundance estimation, differential expression, there are still few ways to investigate cell-71 cell communication diversity in the TME with respect to patient outcomes. Fortunately, new databases of 72 receptor-ligand pairs are becoming available, making work in this area possible (Efremova et al., 2019; 73 Jin et al., 2020; Nath and Leier, 2020; Shao et al., 2020). The methods, code, data, and complete results 74 are available and open to all researchers (https://github.com/IlyaLab/Pan-Cancer-Cell-Cell-Comm-Net). 75 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint https://paperpile.com/c/5TES3g/1F4wh+OtjGl+kvv7o https://paperpile.com/c/5TES3g/1F4wh+OtjGl+kvv7o https://paperpile.com/c/5TES3g/lLbFv https://paperpile.com/c/5TES3g/H1JGR https://paperpile.com/c/5TES3g/vADrs https://paperpile.com/c/5TES3g/RMQGg https://paperpile.com/c/5TES3g/qtDx5 https://paperpile.com/c/5TES3g/nDiJk https://paperpile.com/c/5TES3g/KWOdD https://paperpile.com/c/5TES3g/5wEd8 https://paperpile.com/c/5TES3g/5wEd8 https://paperpile.com/c/5TES3g/CTvjs https://paperpile.com/c/5TES3g/EFysE https://paperpile.com/c/5TES3g/CXcXS https://paperpile.com/c/5TES3g/SO0k5 https://paperpile.com/c/5TES3g/IS5iX+ZgbVl+UyBhS+FfNDP https://paperpile.com/c/5TES3g/IS5iX+ZgbVl+UyBhS+FfNDP https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ 3 Methods 76 Data aggregation and integration 77 Data sources including TCGA and cell-sorted gene expression, bulk tumor expression, cell type scores, 78 cell-ligand and cell-receptor presence estimations were used for network construction and probabilistic 79 weighting on a per-sample basis. 80 81 Each tumor sample is composed of a mixture of cell types including tumor, immune, and stromal cells. 82 Recently, methods have been developed to 'deconvolve' mixed samples into estimated fractions of cell 83 type quantities. For example, xCell, which resembles gene set enrichment, has performed this estimation 84 for 64 cell types across most TCGA samples (Aran et al., 2017). We use these xCell estimates of cellular 85 fractions in this work. 86 87 Ramilowski et al. performed a comprehensive survey of cellular communication, generating a 88 compendium that includes 1,894 ligand-receptor pairs, and a mapping between 144 cell types and 89 expression of ligand or receptor molecules (Ramilowski et al., 2015) The compendium was shared via 5th 90 edition of the FANTOM Project, FANTOM5. These ligand-receptor pairs were adopted for this study. 91 Unfortunately, the FANTOM5 collection of cell types does not overlap well with cell types in xCell. In 92 order to integrate the xCell and FANTOM5 data resources, it was necessary to determine the expressed 93 ligands and receptors for each of the 64 cell types in xCell, using the source gene expression data. 94 95 The xCell project used six public cell sorted bulk gene expression data sets in order to generate gene 96 signatures and score each TCGA sample. Across the data sets, there is some discrepancy in cell type 97 nomenclature, making it necessary to manually curate cell type names to improve alignment across 98 experiments (Supplementary Table 1). Typically, for a given cell type, there are several replicate 99 expression profiles, often across the data sets. 100 Building the cell-cell communication network scaffold 101 102 In the FANTOM5 'draft of cellular communication', an expression threshold of 10 TPMs was used to link 103 a cell type to a ligand or receptor. When considering the distribution of expression in the FANTOM5 104 project, 10 TPMs is close to the median. 105 106 To construct our scaffold, we used a majority voting scheme based on comparing expression levels to 107 median levels. For each cell type, paired with ligands and receptors, if the expression level was greater 108 than the median, it was counted as a vote (i.e., ligand expressed in this cell type). If a ligand or receptor 109 recieved a majority vote across all available data sources, it was accepted, and entered into the cell-cell 110 scaffold. 111 112 With this procedure, a network scaffold is induced, where cells produce ligands that bind to receptors on 113 receiving cells. One edge in the network is composed of components cell - ligand - receptor - cell. This 114 produced a cell-cell communication network with over 1M edges. Each edge represents a possible 115 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint https://paperpile.com/c/5TES3g/YAa5P https://paperpile.com/c/5TES3g/P3pxd https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ 4 interaction in the tumor microenvironment. We subsequently determine the probability that an edge is 116 active in a particular patient sample using a probabilistic method described below. 117 Patient level cell-cell communication network weights 118 With a cell-cell scaffold, expression values and cell type estimations per sample, we can produce a per-119 sample weighted cell-cell communication network (Figure 2). This is done probabilistically, using the 120 following definition: 121 122 𝑃(𝑒𝑖 ) = 𝑃(π‘™π‘Ž , 𝑐𝑙 ) Β· 𝑃(π‘Ÿπ‘ , π‘π‘Ÿ ), (Eq. 1) 123 124 where 𝑒𝑖 is edge i, π‘™π‘Ž is ligand a, π‘Ÿπ‘ is receptor b, and 𝑐𝑙 and π‘π‘Ÿ are cells that can produce ligand a and 125 receptor b respectively. 𝑃(𝑒𝑖 ) represents a probability that edge i is active and is based on the premise that 126 the physical and biochemical link and activation is possible only if all the components are present, and 127 that activity becomes increasingly possible with greater availability of those components. The joint 128 probabilities can be decomposed to: 129 130 𝑃(π‘™π‘Ž , 𝑐𝑙 ) = 𝑃(π‘™π‘Ž | 𝑐𝑙 ) 𝑃(𝑐𝑙 ) and 131 𝑃(π‘Ÿπ‘ , π‘π‘Ÿ ) = 𝑃(π‘Ÿπ‘ | π‘π‘Ÿ ) 𝑃(π‘π‘Ÿ ). (Eq. 2) 132 133 The 𝑃(π‘π‘™π‘˜ ) is short for CDF 𝑃(𝐢𝑙 < π‘π‘™π‘˜ ) which indicates the probability that a randomly sampled value 134 from the empirical 𝐢𝑙 distribution (over all 9K TCGA samples) would be less than the cell estimate for 135 cell type l, in sample k. To do this, for a given cell type, using all samples available, an empirical 136 distribution 𝑃(𝐢𝑙 ) is computed, and for any query, essentially using a value π‘π‘™π‘˜ , the probability can be 137 found by integrating from 0 to π‘π‘™π‘˜. 138 139 To compute 𝑃(π‘™π‘Ž | 𝑐𝑙 ), each 𝐢𝑙 distribution was divided into quartiles, and then (again using the 9K 140 samples) empirical gene expression distributions within each quartile were fit. This expresses the 141 probability that with an observed cell quantity (thus within a quartile), the probability that a randomly 142 selected gene expression value (for gene π‘™π‘Ž) would be lower than what is observed in sample k. 143 144 We refer to "edge weights" to be the probability as shown in Eq. (1). To compute edge weights, each 145 TCGA sample was represented as a column vector of gene expression and a column vector of cell 146 quantities (or enrichments). For each edge in the scaffold (cell-ligand-receptor-cell), data was used to look 147 up probabilities using the defined empirical distributions and taking products for the resulting edge 148 weight probability. This leads to over 9K tumor-specific weighted networks, one for each TCGA 149 participant. 150 151 Probability distributions were precomputed using the R language empirical cumulative distribution 152 function (ecdf). For example, fitting P(CD8 T cells) is done by taking all available estimates across the 153 Pan-Cancer samples and computing the ecdf. Then, for a sample k, we find 𝑃(𝐢𝑙 < π‘π‘™π‘˜ ) using the ecdf. 154 The same technique is used to find the conditional probability functions, where for each gene, the 155 expression values are selected after binning samples using the R function 'quantile', and then used to 156 compute the ecdf. With all distributions precomputed, 9.8 billion joint probability functions were 157 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ 5 computed using an HPC environment, then transferred to a Google BigQuery table where analysis 158 proceeded. This table of network weights was structured so that each row contained one weight from one 159 edge and one tumor sample. Although being a large table of 9.8 billion rows, taking nearly 500GB, 160 BigQuery allows for fast analytical queries that can produce statistics using a selection of standard 161 mathematical functions. 162 Association of network features and survival-based phenotypes 163 The S1 statistic is a robust measure based on the difference of medians (Yahaya et al., 2004; Ahad et al., 164 2016; Babu et al., 1999; Hubert et al. 2012), in this case the median of edge weights for a defined 165 phenotypic group. S1 statistics were computed using the NCI cancer research data commons cloud 166 resource, the ISB-CGC, per tissue type. 167 168 𝑆1 = π‘šπ‘’π‘‘π‘–π‘Žπ‘›(𝑋) βˆ’ π‘šπ‘’π‘‘π‘–π‘Žπ‘›(π‘Œ) √1.4826 𝑀𝐴𝐷(𝑋) + 1.4826 𝑀𝐴𝐷(π‘Œ) 169 170 This statistic allowed for cell-cell interactions to be ranked within a defined context. The results were 171 again saved to BigQuery tables to allow for further cloud-based analysis and integration with underlying 172 data. 173 174 To judge the magnitude of the statistic with respect to a random context (Figure 3), an ensemble of three 175 edge-weight sample-pools were generated, each with 100K weights. Then, for each member of the 176 ensemble, 1 million S1 statistics were generated using sample sizes that match the analyzed data. These 177 random S1 statistic distributions were used to compare to the observed results (i.e., a resampling 178 procedure). 179 180 As an initial examination of the interplay of cell communication and disease, two proxies of disease 181 severity were investigated: progression-free interval (PFI) and tumor stage (Liu et al., 2018). The staging 182 variable used the AJCC pathologic tumor stage. The PFI feature was computed using days until a 183 progression event. The staging variable was binarized by binning stages I-II together (β€œearly stage”), and 184 III-IV together (β€œlate stage”). A binary PFI variable was created by computing the median PFI on non-185 censored samples and then applying the split to all samples. Both clinical features were computed by 186 tissue type (TCGA Study). As Liu et al writes, "The event time is the shortest period from the date of 187 initial diagnosis to the date of an event. The censored time is from the date of initial diagnosis to the date 188 of last contact or the date of death without disease." 189 190 For example, in LUSC, the median time to PFI event was 420 days (14.0 months) and in the censored 191 group, 649 days (21.1 months). After splitting samples at 420 days (14 months), the short PFI group was 192 composed of 67 uncensored samples and 128 censored samples. The long PFI group was composed of 68 193 uncensored samples and 234 censored samples. 194 195 Null distributions, using these same sample sizes (e.g., one group of 68 and another group of 234), were 196 generated by repeatedly drawing from the previously described ensemble of three sample-pools. The 197 distributions, while heavy tailed, were close to Normal (Supplemental Figure 1). The S1 statistics scale 198 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint https://paperpile.com/c/5TES3g/aP9Gd https://paperpile.com/c/5TES3g/aP9Gd https://paperpile.com/c/5TES3g/QrFnS https://paperpile.com/c/5TES3g/QrFnS https://paperpile.com/c/5TES3g/JNgul https://paperpile.com/c/5TES3g/7SC9C https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ 6 with the difference in median values (Supplemental Figure 4). After combining resampled statistics across 199 the ensemble, an edge was selected as a high edge weight if it were in the top 1 millionth percentile when 200 compared to the null. Each tissue and contrast generates a weighted subgraph of the starting scaffold, 201 which is retained for further analysis (e.g., a LUSC-PFI network). 202 203 To identify informative cell-cell edges that relate to disease progression, machine learning models were 204 trained on binarized clinical data as described. With clinical features such as progression free interval 205 (PFI) and tumor stage for each sample, a matrix of patient-specific edge weights was constructed 206 representing each tissue and contrast. Classification of samples was performed with XGBoost classifiers 207 (Chen and Guestrin, 2016) , which are composed of an ensemble of tree classifiers. To avoid overfitting 208 the models, the tree depth was set at maximum of 2 and the early-stopping parameter was set at 2 rounds 209 (training was stopped after the classification error did not improve on a test set for two rounds). XGBoost 210 provides methods for determining the information gain of each feature in the model and was used to rank 211 edges that are most informative for classification. 212 213 Gene ontology (GO) term enrichment was performed using the GONet tool (Pomaznoy et al., 2018). The 214 set of 1,175 genes in the cell-cell scaffold was used as the enrichment background. GONet builds on the 215 "goenrich" software package, which maps genes onto terms and propagates them up the GO graph, 216 performs Fisher's exact tests, and moderates results with FDR. To compare the results, random collections 217 of genes were generated from the cell-cell scaffold and produced no significant results. 218 Results 219 The scaffold network graph is heterogeneous, containing nodes representing cells, ligands (e.g. 220 cytokines), and receptors. Edges are directed, following communication routes from cell to cell. But, to 221 simplify the graph, a cell produces a ligand that binds a receptor found on another cell type, which could 222 make a single edge "LCell-Ligand-Receptor-RCell". In total, there were 1,062,718 cell-cell edges in the 223 network. 224 225 The number of edges for ligand-producing cells varies from 32,910 for osteoblasts to 6,587 for Multi-226 Potent Progenitor (MPPs). For receptor-producing cells, the range spans from 30,225 for platelets to 227 5,763 edges for MPPs. 228 229 Applying the proposed probabilistic framework allowed for the creation of 9,234 weighted networks. The 230 edge weight distributions generally follow approximately exponentially decreasing function 231 (Supplemental Figure 1). There are few edges with strong weights and many with low (near zero) 232 weights. 233 234 We first sought to find communication edges that were most characteristic of an individual tumor type. 235 The S1 statistics comparing one tumor type to all other tissue types was computed, with a high score 236 indicating a substantial difference in edge weights between the two groups. Edges were found that clearly 237 delineated tissues (Figure 4). For example, in SKCM (skin cutaneous melanoma), the top scoring edge is 238 between melanocytes the most cell of origin for cutaneous melanoma (Melanocytes-MIA-CDH19-239 Melanocytes, S1 score 2.5, median edge weight 0.86 higher than in other tumor types). Normal tissue 240 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint https://paperpile.com/c/5TES3g/YrF4J https://paperpile.com/c/5TES3g/x5oMI https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ 7 differences can contribute to differences in edge weights, though in this case the central role of 241 melanocytes in melanomas implies that the high scores are likely due to cancerous cell activity. The study 242 with the most similar edge weights is uveal melanoma (UVA), which arises from melanocytes resident in 243 the uveal tract (Robertson et al., 2017) (Fig. 4A). Additionally, we observed that when a cell type is 244 highly prevalent in a particular tissue, and the scaffold has an autocrine loop, interactions between that 245 type of cell tend to have high weights. If we exclude cell types communicating with self-types, then for 246 SKCM, osteoblasts, natural killer T cells, and mesenchymal stem cells (MSCs) interact with melanocytes 247 in the top 10 scoring edges, consistent with the emerging role of these cell types in melanoma. An 248 important role for osteoblasts is now coming to light for melanoma (Ferguson et al., 2020). Natural killer 249 T cells are being investigated for their applicability in immunotherapy of cancers such as melanoma 250 (Wolf et al., 2018). MSCs appear to interact with melanoma cells, as work by Zhang et al. (Zhang et al., 251 2017) showed the proliferation of A375 cells (a melanoma cell line) was inhibited and the cell cycle of 252 A375 was arrested by MSCs, and cell-cell signaling related to NF-ΞΊB was down-regulated. Overall, the 253 number of high weight edges in each tumor type did not associate with the number of samples, as might 254 be expected (Supplemental Figure 2). 255 256 To identify which elements of cellular communication networks might be associated with clinical 257 progression of particular tumor types, we identified edges associated with disease. 258 Disease progression and severity were examined using dichotomous values of tumor stage and 259 progression free interval (PFI) as described in the methods. Statistical scores were calculated comparing 260 edge weight distributions between the two clinical groups using S1. Results were carried forward if larger 261 than the threshold set by the millionth percentage of resampled statistics (Supplemental Figure 3-5), 262 yielding differentially weighted edges (DWEs). 263 264 Most tumor types showed DWEs for PFI, and fewer for the early to late tumor stage comparison 265 (Supplemental Figure 5). For example, STAD (gastric adenocarcinoma) had several hundred edges in for 266 both comparisons, while PAAD (pancreatic adenocarcinoma) showed fewer DWEs, and only for PFI. 267 Figure 5 shows median edge weights between the two groups for the selected studies. Some tumor types, 268 like SKCM, show much stronger deviations between the medians, compared to the other studies like 269 STAD, ESCA, and LUSC, which may be an indication of a stronger immune response. According to 270 CRI-iAtlas (Eddy et al., 2020), among our example studies, SKCM has the highest estimated level of 271 CD8 T cells and generally has a robust immune response. 272 Tumor stage comparison showed DWEs in 17 of 32 studies and ranged widely from 4 edges for 273 MESO (mesothelioma) to over 63K edges for BLCA (urothelial bladder cancer adenocarcinoma). The 274 PFI comparison showed results in 28 / 32 studies and ranged from 4 edges in READ to over 21K in 275 LIHC. See Table 1 for edge counts from selected studies. The studies with larger numbers of samples had 276 permuted S1 distributions that were narrow compared to studies with few samples (Supp. Fig. 3), but there 277 was not a strong association between DWE counts and sample sizes. The variation thus more likely has to 278 do with clinical factors. 279 280 Within a tumor type and clinical response variable, the set of high scoring edges were usually dominated 281 by a small number of cell-types, ligands, or receptors (Figure 6, Supplemental Figure 6A,6B). For SKCM, 282 in the tumor stage contrast, a majority of ligand-producing cells include GMP cells, Osteoblasts, MSC 283 cells, and Melanocytes, in order of prevalence. The number of edges starting with these four cells 284 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint https://paperpile.com/c/5TES3g/iigUR https://paperpile.com/c/5TES3g/EgnuO https://paperpile.com/c/5TES3g/2INhk https://paperpile.com/c/5TES3g/TD8Vb https://paperpile.com/c/5TES3g/TD8Vb https://paperpile.com/c/5TES3g/CXcXS https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ 8 accounts for 53% of DWEs. Certainly, melanocytes are well known in melanoma, and mesenchymal stem 285 cells are drawn to inflammation, but the role of osteoblasts is less well documented, but still have been 286 associated with melanoma progression (Ferguson et al., 2020). 287 In the PFI contrast with gastroesophageal cancers, megakaryocytes are the most common cell 288 type in STAD DWEs (40 edges out of 142), and the second most common in ESCA (49 edges out of 137, 289 following CD8+ Tcm interactions). The megakaryocyte DWEs include ligands and receptors that 290 represent both interleukins and ECM-associated molecules such as integrins and collagen, but also 291 NOTCH1 and PF4 (platelet factor 4). For STAD, most edge weights are lower with longer PFI. Put 292 another way, the shorter PFI intervals (adverse outcome) were associated with increased megakaryocyte-293 involved edge weights (Supplemental Figure 7). 294 However, the opposite is observed in ESCA, where higher edge weights were generally 295 associated with longer PFI (negative S1 score). In ESCA, edges that show high weights for short PFI 296 include Neutrophils-HMGB1-SDC1-Sebocytes (0.17). Although ESCA has a much lower xCell mean 297 megakaryocyte score than STAD (38% lower), the cell score trends from xCell follow opposite trends 298 with STAD decreasing with longer PFI and ESCA increasing with PFI. STAD is among the tissues with 299 highest megakaryocyte scores (59, 56th rank out of 64 for PFI 1,2 resp.), ESCA is at a respectable rank of 300 49 and 44 out of 64 for short-PFI, respectively. 301 In COAD (colorectal adenocarcinoma), for ligand-producing cells, the DWEs were dominated by 302 astrocytes, MSCs, megakaryocytes, and sebocytes, while receptor-producing cells included astrocytes, 303 chondrocytes, and MSCs in order of counts of DWEs. By summarizing DWEs we can possibly categorize 304 cancer types based on which cells are taking part in potentially active interactions. 305 The above-described edge dominance is related to cells (graph nodes) with high degree. In the 306 language of graphs, the degree is the count of edges connected to a given node or vertex. In STAD the 307 cell types with highest degree are megakaryocytes (degree 50), followed by neutrophils (31), CLP cells 308 (26), and erythrocytes (23)(Supplemental Figure 6A,B). However, if we look at the directionality for the 309 directed graph, we see that while megakaryocytes are split nearly evenly in and out, cells like the Th1 310 have 5 edges in, and only a single edge out, whereas B cells have zero edges in and 4 edges out. The 311 network directionality should be considered in activities such as the modeling of dynamical systems. 312 313 Within the tumor microenvironment, communication between the multitude of cells happens 314 simultaneously through many ligand-receptor axes. By considering a set of differentially weighted edges 315 within a tissue type, we can construct connected networks that potentially represent dynamic 316 communication. DWEs derived by comparing edge weights between clinical groups may indicate which 317 parts of the cell-cell communication network shift together with disease severity. 318 We sought to identify which aspects of intercellular communication could relate to tumor staging 319 or disease severity. The edges making up the differential networks were used to model clinical states for 320 individual tumors. XGBoost models (Chen and Guestrin, 2016) were fit on each clinical feature, using 321 edge weights as predictive variables, to infer which edges carried the most information in classification 322 (Figure 8). 323 The purpose of the modeling was within-data inference rather than classification outside of the 324 TCGA pan-cancer data set. After fitting, it is possible to examine what model features (edges) are most 325 useful in classification. The XGBoost classifiers are regularized models, not all features will be used and 326 often only a small subset of features are retained in the final model. We assess the relative usefulness of a 327 feature by comparing the feature gain -- the improvement in accuracy when a feature is added to a tree. 328 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint https://paperpile.com/c/5TES3g/EgnuO https://paperpile.com/c/5TES3g/YrF4J https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ 9 All classification models had an accuracy between 91% (SKCM, PFI) and 99% (COAD, Stage). As 329 mentioned above, there can be a high degree of correlation between edge values in a data set. While 330 features are selected first based on improving prediction, the machine learning model accounts for 331 correlated features by selecting the one that has best predictive power, leaving out other correlated 332 features. That said, the number of features selected by the model is then related to the correlation 333 structure. In a set of uncorrelated features where all features add to the predictive power, all features will 334 be selected, whereas for correlated features, only a small number will be selected. This is seen in results 335 here in terms of differences in the numbers of features compared with the starting network. 336 In the COAD-PFI case, the number of features was reduced by approximately 20%, keeping 50 337 edges in the model. The STAD-PFI features were reduced by approximately 45%. Other examples are are 338 LUSC-PFI at 60% reduction, ESCA-PFI at 74%, and SKCM-PFI at 95% (12 edges selected) indicating a 339 high degree of internal feature correlation. 340 A similar pattern was observed in the tumor stage contrasts, where SKCM-stage had a 96% 341 reduction in features, STAD-stage 52%, READ-stage 47%. For COAD-stage, feature reduction was 95% 342 reduction, but attributable to the large number of starting edges (1851) compared to the 84 edges selected. 343 A collection of the most predictive edges is given in Table 2. 344 345 The collection of genes from each differential network was used for gene ontology (GO) term enrichment 346 using the GONet tool (Pomaznoy et al., 2018). All tissue-contrast combinations with differentially 347 weighted edges produced enriched GO terms (FDR < 0.05, within tissue contrasts) except the SKCM-348 stage group, which although contained 77 genes in the differential network, produced no enriched terms. 349 Common themes included structural GO terms such as "extracellular structure organization" (for 350 SKCM), cell-substrate adhesion (ESCA, LUSC), cell-cell adhesion (STAD), ECM / extracellular matrix 351 organization (LUSC, COAD, READ, STAD). Cell migration was also a common theme with "cell 352 migration" (STAD), "epithelial cell migration" (SKCM), and "regulation of cell migration" (LUSC, 353 COAD/READ). Among immune related themes, GO terms included "IFNG signaling" and "antigen 354 processing and presentation" (SKCM), "regulation of immune processes" and "IL2" (STAD), and "viral 355 host response" (COAD / READ). See Table 3 for a summary and supplemental table 3 for complete 356 results. 357 Discussion 358 Patient outcome or response to therapy is not necessarily well predicted by tumor stage alone (Kirilovsky 359 et al., 2016). As Fridman et al. wrote, "different types of infiltrating immune cells have different effects 360 on tumour progression, which can vary according to cancer type" (Fridman et al., 2012). This idea has 361 been developed further with the creation of the 'Immunoscore', a prognostic based on the presence and 362 density of particular immune cells in the TME context, expanded to include the peripheral margin as well 363 as tumor core. For example, the Immunoscore in colorectal cancer depends on the density of both CD3+ 364 lymphocytes (any T cell) and specifically, CD8+ cytotoxic T cells in the tumor core and invasive margin 365 (PagΓ¨s et al., 2018). The differences in factors that relate to stage and survival is reflected in the current 366 work in the identification of different cell-cell interactions of importance for each. 367 368 Previous studies have shown that cellular interactions within the tumor microenvironment have an 369 impact on patient survival, drug response, and tumor growth. X. Zhao et al. (Zhou et al., 2017) described 370 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint https://paperpile.com/c/5TES3g/x5oMI https://paperpile.com/c/5TES3g/OVhrZ https://paperpile.com/c/5TES3g/OVhrZ https://paperpile.com/c/5TES3g/KWOdD https://paperpile.com/c/5TES3g/iSTvg https://paperpile.com/c/5TES3g/3xjIz https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ 10 alterations in ligand-receptor pair associations in cancer compared to normal tissue, the cell-cell 371 communication structures thereby becoming a generalized phenotype for malignancy. Using the same 372 foundational database of possible interactions as this work, ligand-receptor pair expression correlation 373 was compared between tumor and normal tissue. Their "aggregate analysis revealed that … tumors of 374 most cancer types generally had reduced (ligand-receptor) correlation compared with the normal tissues." 375 The ligand-receptor pairs that commonly showed such differences across the ten tissue types studied 376 included PLAU-ITGA5, LIPH-LPAR2, SEM14G-PLXNB2, SEMABD-TYROBP, CCL2-CCR5, CCL3-377 CCR5, and CGN-TYROBP. 378 379 Like the Zhao et al. work, we found the collection of associated edges enriched for related biological 380 processes, especially to ECM organization and cell adhesion -- possibly related to the progression towards 381 dysplasia. For example, in Zhao et al., the ligand-receptor pairs COL11A1-ITGA2, COL7A1-ITGA2, 382 MDK-GPC2 and MMP1-ITGA2 were found to be positively correlated in cancer but not in normal tissue. 383 In the current work, integrins and laminins generally have elevated edge weights in late tumor stage. In 384 the PFI contrasts, except for ESCA, such edges have higher weights in shorter PFIs, corresponding to 385 more severe progression. Regarding SEMA7A, found in the PFI STAD results as a predictive feature, 386 previous findings report the collagen gene COL1A1 has been associated with metastasis, and SEMA7A is 387 known to play an important role in integrin-mediated signaling and functions both in regulating cell 388 migration and immune responses. Cancers such as esophageal, gastric, and colorectal all show transitions 389 to metaplasia and dysplasia, a process that breaks down the structural order of a tissue, replacing it with 390 disorder and cell transdifferentiation. 391 392 In our model, a host response is reflected in a change in S1 score, negative if the edge weight is higher 393 with longer PFI times. In the PFI results, Th1 cells appeared in 13 high scoring edges in SKCM, all with 394 negative S1 values. Also, for SKCM and COAD, ligand producing (pro-inflammatory) M1 macrophage 395 edges are present but show both positive and negative S1 scores. Inflammation cytokines IL1B and IL18 396 are both present in the results of ESCA and STAD (Figure 9). In the tumor stage contrasts, we see Th2 397 and NK cells with inflammation cytokines IL1A, IL1B, IL4, TNF in STAD and COAD. So, while certain 398 inflammatory signatures are observed, the absence of well-known canonical edges such as Th1-IL12-399 IL12RB1-M1 macrophages, may be due to essentially no difference, or undetectable differences in the 400 quantity of Th1 cells or IL12A expression between PFI groups (4.9 vs 3.3 TCGA Pan-Cancer RSEM for 401 short vs long PFI). These observations point to possible mechanisms of action for immune cells known to 402 be important for cancer immune response, the CD4+ T helper 1 cells and M1 macrophages, in relation to 403 tumor progression 404 405 In tissues susceptible to dysplasia, such as the tissues explored here, unexpected cell types may be 406 detected. For example, the '...disruption of tissue organization appears to trigger a profound change in 407 cellular commitment, which leads to hepatocyte differentiation in the β€œoval cells” in … the epithelial 408 cells lining the small pancreatic ductules' (Reddy et al., 1991). As another example, pancreatic cancer is 409 known to have desmoplastic stroma, the source of which may include MSCs which are defined by their 410 ability to differentiate into osteoblasts, chondrocytes, and adipocytes (Mathew et al., 2016). In line with 411 that finding, it's been observed that "...stromal cells isolated from the neoplastic pancreas can differentiate 412 into osteoblasts, chondrocytes, and adipocytes" (Mathew et al., 2016). 413 414 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint https://paperpile.com/c/5TES3g/trKJF https://paperpile.com/c/5TES3g/7gKvG https://paperpile.com/c/5TES3g/7gKvG https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ 11 It has been reported that (YÑñez et al., 2017), "granulocyte-monocyte progenitors (GMPs) and monocyte-415 dendritic cell progenitors (MDPs) produce monocytes during homeostasis and in response to increased 416 demand during infection." or as in (Weston et al., 2018), "Granulocyte-monocyte progenitor (GMP) cells 417 play a vital role in the immune system by maturing into a variety of white blood cells, including 418 neutrophils and macrophages, depending on exposure to cytokines such as various types of colony 419 stimulating factors (CSF)." 420 421 In our results for SKCM and COAD, GMPs had negative S1 statistics, meaning the late-stage cases had 422 edges with higher weights. The GMP cells most often interacted with (as receptor bearing cells) MSC, 423 Melanocytes, both M1 and M2 macrophages, and CD8+ Tem (T effector-memory cells). The presence of 424 GMP related edges may be indicative of the commonly observed 'myeloid dysfunction', which "can 425 promote tumor progression through immune suppression, tissue remodeling, angiogenesis or 426 combinations of these mechanisms."(Messmer et al., 2015) Also, "tumors secrete a variety of factors such 427 as G-CSF that act in a systemic way to reduce IRF-8 within progenitor cells, releasing myelopoiesis from 428 IRF-8 control such that the granulocytic lineage (blue cell) undergoes hyperplasia, leading to increased 429 immature suppressive cells to promote tumor growth." This is in line with our observations. 430 431 Megakaryocytes, a multipotent stem cell, are cells that typically reside in the bone marrow and produce 432 platelets. Megakaryocytes are also produced in the liver, kidney, and spleen. Additionally, 433 megakaryocytes have been observed in the lung and circulating blood where they were useful as a 434 biomarker in prostate cancer. Case reports exist showing megakaryocytes in the metaplasia of gastric 435 cancer patients (Chatelain et al., 2004). Megakaryocytes respond to a variety of cytokines such as IL-3, 436 IL-6, IL-11, CXCL5, CXCL7, and CCL5. A majority of interacting cells are leukocytes. In both 437 esophageal and gastric cancers β€œ...thrombocytosis has been reported in general to be associated with 438 adverse clinical outcomes. (Voutsadakis, 2014)" Additionally, there are reports of 'tumor educated 439 platelets' that can be useful as part of a liquid biopsy (Best et al., 2018) (Haemmerle et al., 2018). 440 441 Among the rich literature regarding oncological cytokine networks, there is a strong emphasis on the 442 cancer cell as a central actor. Many of the review articles and research focuses on the cancer cell 443 interactions in the TME. For example, cancer cells producing an overabundance of IL6 or IL10 that has 444 been associated with poor prognosis (Burkholder et al., 2014; Fisher et al., 2014; Lippitz and Harris, 445 2016). 446 447 However, in this work, the focus has been put on the environment and less about the cancer cell itself. 448 This is largely because in performing cell deconvolution on gene expression data to determine the 449 presence and quantity of different cell types in the mixed sample, reliable signatures for cancer cells are 450 not readily available. Because in carcinomas, a cancer cell derives from the epithelium, and in many 451 ways remains similar to epithelial cells. Even in single cell RNA-seq studies, it is often difficult to 452 determine what cells are cancerous and picking this signature out of a mixed expression dataset is difficult 453 and remains an open question. 454 455 This work is based upon gene expression, rather than protein expression, cell-surface expression or 456 secretion measurements. Also, the base expression data is taken from sorted cells, rather than cells in 457 tissue with an assumption that we cannot get β€œnew/non-scaffold edges” in a tissue/cancer context. 458 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint https://paperpile.com/c/5TES3g/X6Njd https://paperpile.com/c/5TES3g/Dtjl8 https://paperpile.com/c/5TES3g/XBz0P https://paperpile.com/c/5TES3g/kSjUZ https://paperpile.com/c/5TES3g/6AFSZ https://paperpile.com/c/5TES3g/OxEbr https://paperpile.com/c/5TES3g/1Nr0r https://paperpile.com/c/5TES3g/QAJ4q+Ew2iX+H4TcE https://paperpile.com/c/5TES3g/QAJ4q+Ew2iX+H4TcE https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ 12 However new data types and methods including scRNA-seq and PIC-seq will provide ways of 459 determining new cell-cell interactions that are context specific (Giladi et al., 2020). 460 461 Importantly, the physical and biochemical process of secretion, binding and activation cannot be 462 identified with the current data and method. By identifying the propensity of edge constituents in 463 particular tumor microenvironments in comparison with others, it becomes more likely that 464 communication with activation can take place, as the presence of those constituents is a prerequisite. 465 466 With the data and results publicly available in a Google BigQuery table (Supplemental Figure 8), this 467 resource is open to researchers to explore and ask questions. It is a low-cost way (with free options) to 468 achieve compute cluster performance for quickly answering such questions. The table is easily joined to 469 clinical and molecular annotations and can be worked with from R and python notebooks. With the 470 addition of resources like GTEx, it should begin to be possible to tease aberrant, cancer specific 471 interactions apart. 472 473 In terms of future work, it could be important to examine communication networks given the immune 474 subtypes of (Thorsson et al., 2018) and communication differences between TCGA tumor molecular 475 subtypes. New data types can be applied to enhance the scaffold with knowledge gained from (for 476 example) single-cell RNA-seq. 477 478 In this work, we have introduced a method and identified lines of communication between cells that may 479 play a role in disease. These lines include both established/recognized cells in the context of cancer, as 480 well as others that should be explored further, with targeted methods. 481 Acknowledgments 482 The authors would like to thank Samuel Danziger, David Reiss, Mark McConnell, Andrew Dervan, 483 Matthew Trotter, Douglas Bassett, Robert Hershberg, the Shmulevich Lab and the Institute for Systems 484 Biology for engaging and informative discussions. This study was supported by Celgene, a wholly owned 485 subsidiary of Bristol-Myers Squibb, in part through a Sponsored Research Award to D.L.G., B.A. and I.S, 486 and by the Cancer Research Institute (D.L.G, V.T, I.S). We thank the ISB-CGC for their ongoing support. 487 ISB-CGC has been funded in whole or in part with Federal funds from the National Cancer Institute, 488 National Institutes of Health, Task Order No. 17X148under Contract No. 75N91019D00024. The content 489 of this publication does not necessarily reflect the views or policies of the Department of Health and 490 Human Services, nor does mention of trade names, commercial products, or organizations imply 491 endorsement by the U.S. Government. 492 Competing Interests. 493 D.L.G., B.A., V.T. and I.S. declare no competing interests. A.V.R.: Bristol-Myers Squibb: Employment, 494 Equity Ownership. 495 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint https://paperpile.com/c/5TES3g/JZSBx https://paperpile.com/c/5TES3g/EFysE https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ 13 Author Contributions 496 D.L.G., B.A., V.T, A.R., I.S. conceived of the idea. D.L.G. developed the method, wrote the code, and 497 performed the computations. D.L.G. wrote the manuscript with contributions from B.A., V.T., A.R., I.S. 498 and A.R. supervised the project. All authors provided critical feedback and helped shape the research, 499 analysis and manuscript. 500 501 Tables 502 Table 1. Counts of differentially weighted edges compared to the number of samples in each study. 503 504 Study N samples PFI short/long PFI DWEs Selected Feat. Model Accuracy GO results? ESCA 170 73/97 137 36 94.7 y STAD 409 155/231 142 78 95.1 y PAAD 178 68/83 8 - - y COAD 281 96/183 63 50 97.1 y READ 91 16/71 4 - - y SKCM 102 27/75 249 12 91.1 y LUSC 494 193/285 304 119 98.7 y Study N samples Stage early/late Stage DWEs Selected Feat. Model Accuracy GO results? ESCA 170 86/63 0 - - - STAD 409 167/198 241 114 99.7 y PAAD 178 142/7 0 - - - COAD 281 151/118 1851 84 99.6 y READ 91 36/44 34 18 97.5 y SKCM 102 68/29 221 8 99 n LUSC 494 390/89 0 - - - Study: tissue type, N samples: number of samples used, PFI short/long: number of samples in each group, 505 PFI DWEs: number of differentially weighted edges, Model Accuracy: accuracy of predicting group, GO 506 results?: if yes, significant GO enrichments. 507 508 509 510 511 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ 14 Table 2. Top 5 most predictive edges from XGBoost models. 512 513 Contrast Study EdgeID LCell Ligand Receptor RCell S1 Median Diff Information Gain PFI COAD 586640 Megakaryocyt es BMP10 ENG Epithelial cells 0.169 0.082 0.109 PFI COAD 50871 astrocytes TNC ITGA5 mv Endothelial cells 0.168 0.061 0.069 PFI COAD 406871 Hepatocytes GDF2 ENG Epithelial cells 0.168 0.082 0.067 PFI COAD 49669 astrocytes EFNB1 EPHB4 Mesangial cells 0.199 0.117 0.066 PFI COAD 632560 MEP TIMP2 ITGB1 MEP 0.167 0.095 0.051 Stage COAD 406579 Hepatocytes CGN TGFBR2 Eosinophils -0.165 -0.077 0.043 Stage COAD 330377 Eosinophils LAMB3 ITGB1 Eosinophils -0.144 -0.048 0.038 Stage COAD 616033 Memory B- cells BMP15 BMPR2 Epithelial cells -0.150 -0.060 0.037 Stage COAD 784400 NK cells TNFSF10 TNFRSF10B CD4+ memory T- cells 0.137 0.043 0.037 Stage COAD 630108 MEP B2M KIR2DL1 iDC 0.138 0.055 0.037 PFI ESCA 457801 Keratinocytes GS ADCY7 CD4+ Tcm 0.167 0.078 0.078 PFI ESCA 182483 CD8+ Tcm RBP3 NOTCH1 pDC -0.184 -0.073 0.071 PFI ESCA 1051114 Th2 cells CALM1 GP6 naive B-cells 0.171 0.085 0.070 PFI ESCA 658080 Mesangial cells SPP1 CD44 Tregs 0.184 0.080 0.064 PFI ESCA 397215 GMP HMGB1 THBD MEP 0.184 0.060 0.059 PFI LUSC 879775 Plasma cells VEGFA ITGB1 GMP 0.120 0.047 0.041 PFI LUSC 451902 iDC VEGFA ITGB1 Plasma cells 0.137 0.067 0.038 PFI LUSC 398971 GMP ADAM17 ITGB1 Plasma cells 0.120 0.059 0.030 PFI LUSC 340857 Epithelial cells COL4A6 ITGB1 CD8+ naive T-cells 0.124 0.054 0.026 PFI LUSC 471558 Keratinocytes THBS1 ITGA6 Plasma cells 0.120 0.068 0.025 Stage READ 632552 MEP TGFB3 TGFBR2 MEP -0.267 -0.134 0.127 Stage READ 795527 NKT GZMB PGRMC1 CD4+ memory T- cells 0.343 0.144 0.115 Stage READ 402754 Hepatocytes CGN TGFBR2 CD4+ Tem -0.274 -0.134 0.108 Stage READ 808308 NKT GZMB IGF2R Plasma cells 0.261 0.101 0.103 Stage READ 800747 NKT IL7 IL2RG GMP 0.264 0.136 0.095 PFI SKCM 1008243 Smooth muscle SEMA7A PLXNC1 pro B-cells 0.438 0.259 0.242 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ 15 PFI SKCM 517677 Macrophages UBA52 NOTCH1 Osteoblast -0.284 -0.145 0.200 PFI SKCM 80934 Basophils VIM CD44 NKT -0.383 -0.254 0.103 PFI SKCM 1007915 Smooth muscle PSAP SORT1 Preadipocytes 0.311 0.175 0.082 PFI SKCM 84049 Basophils CALM1 PTPRA Th1 cells -0.285 -0.151 0.080 Stage SKCM 275306 CLP GI2 CXCR1 Osteoblast 0.353 0.176 0.207 Stage SKCM 399084 GMP TIMP1 CD63 Plasma cells -0.302 -0.147 0.206 Stage SKCM 273727 CLP GI2 F2R MEP 0.290 0.123 0.182 Stage SKCM 182981 CD8+ Tcm GI2 TBXA2R Plasma cells -0.283 -0.095 0.123 Stage SKCM 397545 GMP BST1 CAV1 MSC -0.337 -0.194 0.109 PFI STAD 461765 Keratinocytes CALM3 KCNQ1 Eosinophils -0.136 -0.067 0.062 PFI STAD 644724 Mesangial cells TGFB2 ACVR1 Erythrocytes 0.149 0.061 0.054 PFI STAD 105991 CD4+ T-cells IL1B IL1R2 Megakaryocytes 0.134 0.081 0.047 PFI STAD 269013 CLP ADAM28 ITGA4 CD4+ T-cells 0.145 0.075 0.046 PFI STAD 343620 Epithelial cells VCAN TLR1 CLP 0.134 0.051 0.033 Stage STAD 128412 CD4+ Tem CALM1 KCNQ1 Macrophages 0.140 0.058 0.057 Stage STAD 43832 astrocytes FBN1 ITGB6 Epithelial cells -0.146 -0.058 0.036 Stage STAD 346120 Epithelial cells LAMB1 ITGAV Hepatocytes -0.139 -0.066 0.035 Stage STAD 403540 Hepatocytes SHH PTCH1 CD8+ T-cells -0.138 -0.069 0.034 Stage STAD 648983 Mesangial cells FGB ITGAV Megakaryocytes -0.140 -0.060 0.031 Contrast: the groupwise test performed, Study: tissue type, Edge ID: BigQuery table lookup ID, LCell: 514 cell producing ligands, Ligand: ligand gene symbol , Receptor: receptor gene symbol, R Cell: receptor 515 producing cell, S1: between group S1 statistic, Median Diff: difference in edge weights between groups, 516 Information Gain: Xgboost information gain after adding feature to model. 517 518 Table 3. Enriched GO terms. 519 Tissue Contrast Num GOs ECM Migration Immune Immune2 SKCM PFI 34 extracellular structure organization epithelium cell migration IFNG signaling antigen processing and presentation ESCA PFI 3 cell-substrate adhesion .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ 16 STAD PFI 59 cell-cell adhesion mediated by integrin cell migration regulation of immune system process IL2 LUSC PFI 39 extracellular matrix organization positive regulation of cell migration COAD / READ stage 85 ECM regulation of epithelial cell migration viral host response STAD stage 28 ECM / adhesion cell migration Tissue: TCGA study, Contrast: the groupwise test performed, Num GOs: number of gene ontology terms 520 found significantly enriched, ECM: GO categories involving ECM, Migration: GO terms involving cell 521 migration, Immune: GO terms involving immune response, Immune2: additional GO terms involving 522 immune response. 523 524 References 525 526 Ahad, N. A., Yahaya, S. S. S., and Yin, L. P. (2016). Robustness of S1 statistic with Hodges-Lehmann for 527 skewed distributions. AIP Conf. Proc. 1782, 050002. 528 Aran, D., Hu, Z., and Butte, A. J. (2017). xCell: digitally portraying the tissue cellular heterogeneity 529 landscape. Genome Biol. 18, 220. 530 Babu, G. J., Padmanabhan, A. R., and Puri, M. L. (1999). Robust one-way ANOVA under possibly non-531 regular conditions. Biometrical Journal: Journal of Mathematical Methods in Biosciences 41, 321–532 339. 533 Behar, M., Barken, D., Werner, S. L., and Hoffmann, A. (2013). The dynamics of signaling as a 534 pharmacological target. Cell 155, 448–461. 535 Best, M. G., Wesseling, P., and Wurdinger, T. (2018). Tumor-Educated Platelets as a Noninvasive 536 Biomarker Source for Cancer Detection and Progression Monitoring. Cancer Res. 78, 3407–3412. 537 Burkholder, B., Huang, R.-Y., Burgess, R., Luo, S., Jones, V. S., Zhang, W., et al. (2014). Tumor-induced 538 perturbations of cytokines and immune cell networks. Biochim. Biophys. Acta 1845, 182–201. 539 Cameron, M. J., and Kelvin, D. J. (2013). Cytokines, Chemokines and Their Receptors. Landes 540 Bioscience. 541 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint http://paperpile.com/b/5TES3g/QrFnS http://paperpile.com/b/5TES3g/QrFnS http://paperpile.com/b/5TES3g/QrFnS http://paperpile.com/b/5TES3g/QrFnS http://paperpile.com/b/5TES3g/YAa5P http://paperpile.com/b/5TES3g/YAa5P http://paperpile.com/b/5TES3g/YAa5P http://paperpile.com/b/5TES3g/YAa5P http://paperpile.com/b/5TES3g/JNgul http://paperpile.com/b/5TES3g/JNgul http://paperpile.com/b/5TES3g/JNgul http://paperpile.com/b/5TES3g/JNgul http://paperpile.com/b/5TES3g/JNgul http://paperpile.com/b/5TES3g/5wEd8 http://paperpile.com/b/5TES3g/5wEd8 http://paperpile.com/b/5TES3g/5wEd8 http://paperpile.com/b/5TES3g/5wEd8 http://paperpile.com/b/5TES3g/OxEbr http://paperpile.com/b/5TES3g/OxEbr http://paperpile.com/b/5TES3g/OxEbr http://paperpile.com/b/5TES3g/OxEbr http://paperpile.com/b/5TES3g/H4TcE http://paperpile.com/b/5TES3g/H4TcE http://paperpile.com/b/5TES3g/H4TcE http://paperpile.com/b/5TES3g/H4TcE http://paperpile.com/b/5TES3g/kvv7o http://paperpile.com/b/5TES3g/kvv7o http://paperpile.com/b/5TES3g/kvv7o http://paperpile.com/b/5TES3g/kvv7o https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ 17 Cancer Genome Atlas Network (2015). Genomic Classification of Cutaneous Melanoma. Cell 161, 1681–542 1696. 543 Chatelain, D., Devendeville, A., Rudelli, A., Bruniau, A., Geslin, G., and Sevestre, H. (2004). Gastric 544 myeloid metaplasia: a case report and review of the literature. Arch. Pathol. Lab. Med. 128, 568–545 570. 546 Chen, T., and Guestrin, C. (2016). XGBoost: A Scalable Tree Boosting System. in Proceedings of the 547 22Nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining KDD 548 ’16. (New York, NY, USA: ACM), 785–794. 549 Cohen, D. J., and Nelson, W. J. (2018). Secret handshakes: cell-cell interactions and cellular mimics. 550 Curr. Opin. Cell Biol. 50, 14–19. 551 Eddy, J. A., Thorsson, V., Lamb, A. E., Gibbs, D. L., Heimann, C., Yu, J. X., et al. (2020). CRI iAtlas: an 552 interactive portal for immuno-oncology research. F1000Res. 9, 1028. 553 Efremova, M., Vento-Tormo, M., Teichmann, S. A., and Vento-Tormo, R. (2019). CellPhoneDB v2.0: 554 Inferring cell-cell communication from combined expression of multi-subunit receptor-ligand 555 complexes. bioRxiv. doi:10.1101/680926. 556 Ferguson, J., Wilcock, D. J., McEntegart, S., Badrock, A. P., Levesque, M., Dummer, R., et al. (2020). 557 Osteoblasts contribute to a protective niche that supports melanoma cell proliferation and survival. 558 Pigment Cell Melanoma Res. 33, 74–85. 559 Fisher, D. T., Appenheimer, M. M., and Evans, S. S. (2014). The two faces of IL-6 in the tumor 560 microenvironment. Semin. Immunol. 26, 38–47. 561 Frankenstein, Z., Alon, U., and Cohen, I. R. (2006). The immune-body cytokine network defines a social 562 architecture of cell interactions. Biol. Direct 1, 32. 563 Fridman, W. H., PagΓ¨s, F., SautΓ¨s-Fridman, C., and Galon, J. (2012). The immune contexture in human 564 tumours: impact on clinical outcome. Nat. Rev. Cancer 12, 298–306. 565 Giladi, A., Cohen, M., Medaglia, C., Baran, Y., Li, B., Zada, M., et al. (2020). Dissecting cellular 566 crosstalk by sequencing physically interacting cells. Nat. Biotechnol. 38, 629–637. 567 Haass, N. K., and Herlyn, M. (2005). Normal human melanocyte homeostasis as a paradigm for 568 understanding melanoma. J. Investig. Dermatol. Symp. Proc. 10, 153–163. 569 Haemmerle, M., Stone, R. L., Menter, D. G., Afshar-Kharghan, V., and Sood, A. K. (2018). The Platelet 570 Lifeline to Cancer: Challenges and Opportunities. Cancer Cell 33, 965–983. 571 Heldin, C.-H., Lu, B., Evans, R., and Gutkind, J. S. (2016). Signals and Receptors. Cold Spring Harb. 572 Perspect. Biol. 8, a005900. 573 Hubert M., Pison G., Struyf A., Van Aelst S., editors. Theory and applications of recent robust methods. 574 BirkhΓ€user; 2012 Dec 6. 575 Jin, S., Guerrero-Juarez, C. F., Zhang, L., Chang, I., Myung, P., Plikus, M. V., et al. (2020). Inference and 576 analysis of cell-cell communication using CellChat. Cold Spring Harbor Laboratory, 577 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint http://paperpile.com/b/5TES3g/SO0k5 http://paperpile.com/b/5TES3g/SO0k5 http://paperpile.com/b/5TES3g/SO0k5 http://paperpile.com/b/5TES3g/SO0k5 http://paperpile.com/b/5TES3g/kSjUZ http://paperpile.com/b/5TES3g/kSjUZ http://paperpile.com/b/5TES3g/kSjUZ http://paperpile.com/b/5TES3g/kSjUZ http://paperpile.com/b/5TES3g/kSjUZ http://paperpile.com/b/5TES3g/YrF4J http://paperpile.com/b/5TES3g/YrF4J http://paperpile.com/b/5TES3g/YrF4J http://paperpile.com/b/5TES3g/YrF4J http://paperpile.com/b/5TES3g/YrF4J http://paperpile.com/b/5TES3g/1F4wh http://paperpile.com/b/5TES3g/1F4wh http://paperpile.com/b/5TES3g/1F4wh http://paperpile.com/b/5TES3g/1F4wh http://paperpile.com/b/5TES3g/CXcXS http://paperpile.com/b/5TES3g/CXcXS http://paperpile.com/b/5TES3g/CXcXS http://paperpile.com/b/5TES3g/CXcXS http://paperpile.com/b/5TES3g/IS5iX http://paperpile.com/b/5TES3g/IS5iX http://paperpile.com/b/5TES3g/IS5iX http://paperpile.com/b/5TES3g/IS5iX http://paperpile.com/b/5TES3g/IS5iX http://dx.doi.org/10.1101/680926 http://paperpile.com/b/5TES3g/IS5iX http://paperpile.com/b/5TES3g/EgnuO http://paperpile.com/b/5TES3g/EgnuO http://paperpile.com/b/5TES3g/EgnuO http://paperpile.com/b/5TES3g/EgnuO http://paperpile.com/b/5TES3g/QAJ4q http://paperpile.com/b/5TES3g/QAJ4q http://paperpile.com/b/5TES3g/QAJ4q http://paperpile.com/b/5TES3g/QAJ4q http://paperpile.com/b/5TES3g/lLbFv http://paperpile.com/b/5TES3g/lLbFv http://paperpile.com/b/5TES3g/lLbFv http://paperpile.com/b/5TES3g/lLbFv http://paperpile.com/b/5TES3g/KWOdD http://paperpile.com/b/5TES3g/KWOdD http://paperpile.com/b/5TES3g/KWOdD http://paperpile.com/b/5TES3g/KWOdD http://paperpile.com/b/5TES3g/JZSBx http://paperpile.com/b/5TES3g/JZSBx http://paperpile.com/b/5TES3g/JZSBx http://paperpile.com/b/5TES3g/JZSBx http://paperpile.com/b/5TES3g/P8KBB http://paperpile.com/b/5TES3g/P8KBB http://paperpile.com/b/5TES3g/P8KBB http://paperpile.com/b/5TES3g/P8KBB http://paperpile.com/b/5TES3g/1Nr0r http://paperpile.com/b/5TES3g/1Nr0r http://paperpile.com/b/5TES3g/1Nr0r http://paperpile.com/b/5TES3g/1Nr0r http://paperpile.com/b/5TES3g/OtjGl http://paperpile.com/b/5TES3g/OtjGl http://paperpile.com/b/5TES3g/OtjGl http://paperpile.com/b/5TES3g/OtjGl http://paperpile.com/b/5TES3g/ZgbVl http://paperpile.com/b/5TES3g/ZgbVl http://paperpile.com/b/5TES3g/ZgbVl http://paperpile.com/b/5TES3g/ZgbVl https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ 18 2020.07.21.214387. doi:10.1101/2020.07.21.214387. 578 Kirilovsky, A., Marliot, F., El Sissy, C., Haicheur, N., Galon, J., and PagΓ¨s, F. (2016). Rational bases for 579 the use of the Immunoscore in routine clinical settings as a prognostic and predictive biomarker in 580 cancer patients. Int. Immunol. 28, 373–382. 581 Kumar, M. P., Du, J., Lagoudas, G., Jiao, Y., Sawyer, A., Drummond, D. C., et al. (2018). Analysis of 582 Single-Cell RNA-Seq Identifies Cell-Cell Communication Associated with Tumor Characteristics. 583 Cell Rep. 25, 1458–1468.e4. 584 Lippitz, B. E., and Harris, R. A. (2016). Cytokine patterns in cancer patients: A review of the correlation 585 between interleukin 6 and prognosis. Oncoimmunology 5, e1093722. 586 Liu, J., Lichtenberg, T., Hoadley, K. A., Poisson, L. M., Lazar, A. J., Cherniack, A. D., et al. (2018). An 587 Integrated TCGA Pan-Cancer Clinical Data Resource to Drive High-Quality Survival Outcome 588 Analytics. Cell 173, 400–416.e11. 589 Mathew, E., Brannon, A. L., Del Vecchio, A., Garcia, P. E., Penny, M. K., Kane, K. T., et al. (2016). 590 Mesenchymal Stem Cells Promote Pancreatic Tumor Growth by Inducing Alternative Polarization 591 of Macrophages. Neoplasia 18, 142–151. 592 Messmer, M. N., Netherby, C. S., Banik, D., and Abrams, S. I. (2015). Tumor-induced myeloid 593 dysfunction and its implications for cancer immunotherapy. Cancer Immunol. Immunother. 64, 1–594 13. 595 Morel, P. A., Lee, R. E. C., and Faeder, J. R. (2017). Demystifying the cytokine network: Mathematical 596 models point the way. Cytokine 98, 115–123. 597 Nath, A., and Leier, A. (2020). Improved cytokine-receptor interaction prediction by exploiting the 598 negative sample space. BMC Bioinformatics 21, 493. 599 PagΓ¨s, F., Mlecnik, B., Marliot, F., Bindea, G., Ou, F.-S., Bifulco, C., et al. (2018). International 600 validation of the consensus Immunoscore for the classification of colon cancer: a prognostic and 601 accuracy study. Lancet 391, 2128–2139. 602 Pomaznoy, M., Ha, B., and Peters, B. (2018). GOnet: a tool for interactive Gene Ontology analysis. BMC 603 Bioinformatics 19, 470. 604 Ramilowski, J. A., Goldberg, T., Harshbarger, J., Kloppmann, E., Lizio, M., Satagopam, V. P., et al. 605 (2015). A draft network of ligand–receptor-mediated multicellular signalling in human. Nat. 606 Commun. 6, 7866. 607 Reddy, J. K., Rao, M. S., Yeldandi, A. V., Tan, X. D., and Dwivedi, R. S. (1991). Pancreatic hepatocytes. 608 An in vivo model for cell lineage in pancreas of adult rat. Dig. Dis. Sci. 36, 502–509. 609 Robertson, A. G., Shih, J., Yau, C., Gibb, E. A., Oba, J., Mungall, K. L., et al. (2017). Integrative 610 Analysis Identifies Four Molecular and Clinical Subsets in Uveal Melanoma. Cancer Cell 32, 204–611 220.e15. 612 Shao, X., Liao, J., Li, C., Lu, X., Cheng, J., and Fan, X. (2020). CellTalkDB: a manually curated database 613 of ligand-receptor interactions in humans and mice. Brief. Bioinform. doi:10.1093/bib/bbaa269. 614 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint http://paperpile.com/b/5TES3g/ZgbVl http://dx.doi.org/10.1101/2020.07.21.214387 http://paperpile.com/b/5TES3g/ZgbVl http://paperpile.com/b/5TES3g/OVhrZ http://paperpile.com/b/5TES3g/OVhrZ http://paperpile.com/b/5TES3g/OVhrZ http://paperpile.com/b/5TES3g/OVhrZ http://paperpile.com/b/5TES3g/OVhrZ http://paperpile.com/b/5TES3g/nDiJk http://paperpile.com/b/5TES3g/nDiJk http://paperpile.com/b/5TES3g/nDiJk http://paperpile.com/b/5TES3g/nDiJk http://paperpile.com/b/5TES3g/Ew2iX http://paperpile.com/b/5TES3g/Ew2iX http://paperpile.com/b/5TES3g/Ew2iX http://paperpile.com/b/5TES3g/Ew2iX http://paperpile.com/b/5TES3g/7SC9C http://paperpile.com/b/5TES3g/7SC9C http://paperpile.com/b/5TES3g/7SC9C http://paperpile.com/b/5TES3g/7SC9C http://paperpile.com/b/5TES3g/7SC9C http://paperpile.com/b/5TES3g/7gKvG http://paperpile.com/b/5TES3g/7gKvG http://paperpile.com/b/5TES3g/7gKvG http://paperpile.com/b/5TES3g/7gKvG http://paperpile.com/b/5TES3g/7gKvG http://paperpile.com/b/5TES3g/XBz0P http://paperpile.com/b/5TES3g/XBz0P http://paperpile.com/b/5TES3g/XBz0P http://paperpile.com/b/5TES3g/XBz0P http://paperpile.com/b/5TES3g/XBz0P http://paperpile.com/b/5TES3g/CTvjs http://paperpile.com/b/5TES3g/CTvjs http://paperpile.com/b/5TES3g/CTvjs http://paperpile.com/b/5TES3g/CTvjs http://paperpile.com/b/5TES3g/FfNDP http://paperpile.com/b/5TES3g/FfNDP http://paperpile.com/b/5TES3g/FfNDP http://paperpile.com/b/5TES3g/FfNDP http://paperpile.com/b/5TES3g/iSTvg http://paperpile.com/b/5TES3g/iSTvg http://paperpile.com/b/5TES3g/iSTvg http://paperpile.com/b/5TES3g/iSTvg http://paperpile.com/b/5TES3g/iSTvg http://paperpile.com/b/5TES3g/x5oMI http://paperpile.com/b/5TES3g/x5oMI http://paperpile.com/b/5TES3g/x5oMI http://paperpile.com/b/5TES3g/x5oMI http://paperpile.com/b/5TES3g/P3pxd http://paperpile.com/b/5TES3g/P3pxd http://paperpile.com/b/5TES3g/P3pxd http://paperpile.com/b/5TES3g/P3pxd http://paperpile.com/b/5TES3g/P3pxd http://paperpile.com/b/5TES3g/trKJF http://paperpile.com/b/5TES3g/trKJF http://paperpile.com/b/5TES3g/trKJF http://paperpile.com/b/5TES3g/trKJF http://paperpile.com/b/5TES3g/iigUR http://paperpile.com/b/5TES3g/iigUR http://paperpile.com/b/5TES3g/iigUR http://paperpile.com/b/5TES3g/iigUR http://paperpile.com/b/5TES3g/iigUR http://paperpile.com/b/5TES3g/UyBhS http://paperpile.com/b/5TES3g/UyBhS http://paperpile.com/b/5TES3g/UyBhS http://paperpile.com/b/5TES3g/UyBhS http://dx.doi.org/10.1093/bib/bbaa269 http://paperpile.com/b/5TES3g/UyBhS https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ 19 Song, D., Yang, D., Powell, C. A., and Wang, X. (2019). Cell-cell communication: old mystery and new 615 opportunity. Cell Biol. Toxicol. 35, 89–93. 616 Theory and Applications of Recent Robust Methods | Mia Hubert | Springer Available at: 617 https://www.springer.com/gp/book/9783764370602 [Accessed June 19, 2020]. 618 Thorsson, V., Gibbs, D. L., Brown, S. D., Wolf, D., Bortone, D. S., Ou Yang, T.-H., et al. (2018). The 619 Immune Landscape of Cancer. Immunity 48, 812–830.e14. 620 Trosko, J. E., and Ruch, R. J. (1998). Cell-cell communication in carcinogenesis. Front. Biosci. 3, d208–621 36. 622 Voutsadakis, I. A. (2014). Thrombocytosis as a prognostic marker in gastrointestinal cancers. World J. 623 Gastrointest. Oncol. 6, 34–40. 624 Wei, C.-J., Xu, X., and Lo, C. W. (2004). Connexins and cell signaling in development and disease. 625 Annu. Rev. Cell Dev. Biol. 20, 811–838. 626 West, J., and Newton, P. K. (2019). Cellular interactions constrain tumor growth. Proc. Natl. Acad. Sci. 627 U. S. A. 116, 1918–1923. 628 Weston, B. R., Li, L., and Tyson, J. J. (2018). Mathematical Analysis of Cytokine-Induced Differentiation 629 of Granulocyte-Monocyte Progenitor Cells. Front. Immunol. 9, 2048. 630 Wilson, M. R., Close, T. W., and Trosko, J. E. (2000). Cell Population Dynamics (Apoptosis, Mitosis, 631 and Cell–Cell Communication) during Disruption of Homeostasis. Exp. Cell Res. 254, 257–268. 632 Wolf, B. J., Choi, J. E., and Exley, M. A. (2018). Novel Approaches to Exploiting Invariant NKT Cells in 633 Cancer Immunotherapy. Front. Immunol. 9, 384. 634 Yahaya, S. S. S., Othman, A. R., and Keselman, H. J. (2004). Testing the Equality of Location Parameters 635 for Skewed Distributions Using S1 with High Breakdown Robust Scale Estimators. in Theory and 636 Applications of Recent Robust Methods (BirkhΓ€user Basel), 319–328. 637 YÑñez, A., Coetzee, S. G., Olsson, A., Muench, D. E., Berman, B. P., Hazelett, D. J., et al. (2017). 638 Granulocyte-Monocyte Progenitors and Monocyte-Dendritic Cell Progenitors Independently 639 Produce Functionally Distinct Monocytes. Immunity 47, 890–902.e4. 640 Zhang, J., Hou, L., Zhao, D., Pan, M., Wang, Z., Hu, H., et al. (2017). Inhibitory effect and mechanism of 641 mesenchymal stem cells on melanoma cells. Clin. Transl. Oncol. 19, 1358–1374. 642 Zhou, J. X., Taramelli, R., Pedrini, E., Knijnenburg, T., and Huang, S. (2017). Extracting Intercellular 643 Signaling Network of Cancer Tissues using Ligand-Receptor Expression Patterns from Whole-644 tumor and Single-cell Transcriptomes. Sci. Rep. 7, 8815. 645 646 647 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint http://paperpile.com/b/5TES3g/H1JGR http://paperpile.com/b/5TES3g/H1JGR http://paperpile.com/b/5TES3g/H1JGR http://paperpile.com/b/5TES3g/H1JGR http://paperpile.com/b/5TES3g/Ups9L http://paperpile.com/b/5TES3g/Ups9L https://www.springer.com/gp/book/9783764370602 http://paperpile.com/b/5TES3g/Ups9L http://paperpile.com/b/5TES3g/EFysE http://paperpile.com/b/5TES3g/EFysE http://paperpile.com/b/5TES3g/EFysE http://paperpile.com/b/5TES3g/EFysE http://paperpile.com/b/5TES3g/RMQGg http://paperpile.com/b/5TES3g/RMQGg http://paperpile.com/b/5TES3g/RMQGg http://paperpile.com/b/5TES3g/RMQGg http://paperpile.com/b/5TES3g/6AFSZ http://paperpile.com/b/5TES3g/6AFSZ http://paperpile.com/b/5TES3g/6AFSZ http://paperpile.com/b/5TES3g/6AFSZ http://paperpile.com/b/5TES3g/vADrs http://paperpile.com/b/5TES3g/vADrs http://paperpile.com/b/5TES3g/vADrs http://paperpile.com/b/5TES3g/vADrs http://paperpile.com/b/5TES3g/qtDx5 http://paperpile.com/b/5TES3g/qtDx5 http://paperpile.com/b/5TES3g/qtDx5 http://paperpile.com/b/5TES3g/qtDx5 http://paperpile.com/b/5TES3g/Dtjl8 http://paperpile.com/b/5TES3g/Dtjl8 http://paperpile.com/b/5TES3g/Dtjl8 http://paperpile.com/b/5TES3g/Dtjl8 http://paperpile.com/b/5TES3g/QCquJ http://paperpile.com/b/5TES3g/QCquJ http://paperpile.com/b/5TES3g/QCquJ http://paperpile.com/b/5TES3g/QCquJ http://paperpile.com/b/5TES3g/2INhk http://paperpile.com/b/5TES3g/2INhk http://paperpile.com/b/5TES3g/2INhk http://paperpile.com/b/5TES3g/2INhk http://paperpile.com/b/5TES3g/aP9Gd http://paperpile.com/b/5TES3g/aP9Gd http://paperpile.com/b/5TES3g/aP9Gd http://paperpile.com/b/5TES3g/aP9Gd http://paperpile.com/b/5TES3g/aP9Gd http://paperpile.com/b/5TES3g/X6Njd http://paperpile.com/b/5TES3g/X6Njd http://paperpile.com/b/5TES3g/X6Njd http://paperpile.com/b/5TES3g/X6Njd http://paperpile.com/b/5TES3g/X6Njd http://paperpile.com/b/5TES3g/TD8Vb http://paperpile.com/b/5TES3g/TD8Vb http://paperpile.com/b/5TES3g/TD8Vb http://paperpile.com/b/5TES3g/TD8Vb http://paperpile.com/b/5TES3g/3xjIz http://paperpile.com/b/5TES3g/3xjIz http://paperpile.com/b/5TES3g/3xjIz http://paperpile.com/b/5TES3g/3xjIz http://paperpile.com/b/5TES3g/3xjIz https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ 20 Figure Legends 648 649 Figure 1. Overview of workflow showing the transition from data sources to results. 650 651 Figure 2. Illustration of the probabilistic model and edge weight computations. (A) For a given cell-cell 652 communication edge, (B) per patient values are used to 'look up' probabilities from the distributions 653 learned from all TCGA data. Those probabilities are then used to compute an edge weight. 654 655 Figure 3. Diagram of how differentially weighted edges were determined. Three samples of edge weights 656 were taken from the pool by tissue source. Then matching the sample proportions in the clinical features, 657 permutations were sampled and used for computing randomized S1 statistics. Each sample was used to 658 produce 1 million permuted statistics, and taken together, the millionth percentile was used as a cutoff in 659 determining important edges. 660 661 Figure 4. Top edges (by S1 scores) that can distinguish tissue types. Each point represents a tumor sample 662 and each panel represents one edge. (A) EdgeID 605551, Melanocytes-MIA-CDH19-Melanocyte SKCM 663 red, UVM blue, BRCA purple, PAAD orange. (B) EdgeID 687457, MSC-TFPI-F3-MSC, PAAD red. (C) 664 EdgeID 968128, Sebocytes-WNT5A-FZD6- Sebocytes, LUSC red, LUAD blue, HNSC purple. (D) 665 EdgeID 1049823, Th2 cells-IL4-IL2RG-Megakaryocytes, STAD red, READ blue, COAD purple, ESCA 666 orange. 667 668 Figure 5. (A) Median values for each differentially weighted cell-cell edge (DWE) for the PFI categories 669 (in row, DWE edges in columns). (B) Examples of differentially weighted edges. 670 671 Figure 6. Edge member dominance in DWEs shown by log10 counts of cell types. 672 673 Figure 7. High probability edges (DWEs) from PFI contrasts form predictive connected subnetworks. 674 Color indicates the magnitude and direction of S1 statistics (+ / -). 675 676 Figure 8. Informative edges selected by XGBoost models for prediction within study. Color indicates 677 information gain. 678 679 Figure 9. Cell-cell interaction diagram demonstrating complexity in communication with three cell types 680 that produce the IL1B ligand that have two possible binding partners on the same receptor bearing cell. 681 Edge weight violin plots are shown for two STAD PFI groups, short (left) and long (right) PFI. 682 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/ .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted February 10, 2021. ; https://doi.org/10.1101/2021.02.08.430343doi: bioRxiv preprint https://doi.org/10.1101/2021.02.08.430343 http://creativecommons.org/licenses/by/4.0/