Delineation of colorectal cancer ligand-receptor interactions and their roles in the tumor microenvironment and prognosis

Immunotherapies targeting ligand-receptor interactions (LRIs) are advancing rapidly in the treatment of colorectal cancer (CRC), and LRIs also affect many aspects of CRC development. However, the pattern of LRIs in CRC and their effect on tumor microenvironment and clinical value are still unclear. We delineated the pattern of LRIs in 55,539 single-cell RNA sequencing (scRNA-seq) samples from 29 patients with CRC and three bulk RNA-seq datasets containing data from 1411 CRC patients. Then the influence of tumor microenvironment, immunotherapy and prognosis of CRC patients were comprehensively investigated. We calculated the strength of 1893 ligand-receptor pairs between 25 cell types to reconstruct the spatial structure of CRC. We identified tumor subtypes based on LRIs, revealed the relationship between the subtypes and immunotherapy efficacy and explored the ligand-receptor pairs and specific targets affecting the abundance of tumor-infiltrating lymphocytes. Finally, a prognostic model based on ligand-receptor pairs was constructed and validated. Overall, through the comprehensive and in-depth investigation of the existing ligand-receptor pairs, this study provides new ideas for CRC subtype classification, a new risk screening tool for CRC patients, and potential ligand-receptor pair targets and pathways for CRC therapy.


Introduction
According to the latest estimates of the global cancer burden released by the International Agency for Research on Cancer (IARC) in 2020, colorectal cancer (CRC) is the third most common cancer globally, with the second highest case fatality rate. In recent years, research on advanced CRC has advanced rapidly, and preoperative neoadjuvant therapy and total neoadjuvant therapy have been widely studied in clinical trials. After the publication of the KEYNOTE 177 study in 2020, PD-1 blockade therapy based on the ligand-receptor interaction (LRI) between tumor cells and T cells became the first-line therapy for MSI-H metastatic CRC [1]. Increasing attention has been given to research on intercellular LRIs in CRC treatment.
As an essential component of cell-cell communication, LRIs play a vital role in the development and treatment of cancer. LRIs involve various cells in the tumor microenvironment, and interactions in different cells produce different effects. Current research is mainly limited to the impact of a single ligand-receptor pair between two cell types, but the comprehensive effect of multiple ligandreceptor pairs between different cells is still unclear. With the widespread use of single-cell RNA sequencing (scRNA-seq), several studies have developed algorithms and tools to investigate the effects of LRIs via scRNA-seq data [2,3], which makes it possible to explore LRIs in a holistic way. However, the entire landscape of LRIs in the colorectal cancer ecological environment has not been fully clarified.
We delineated the existing LRIs between several cell types in the tumor microenvironment by analyzing scRNA-seq data from human CRC samples. After clustering and annotating scRNA-seq data, we calculated the strength of the LRIs between different cells, described the general situation of the LRIs, and reconstructed the spatial structure of CRC by LRI strength, which will allow us to explore the pattern of cell-cell communication. Furthermore, based on the strength of the LRIs, we identified malignant epithelial cell subtypes and CRC subtypes and revealed differences in the tumor microenvironment and changes in immunotherapy-related factors. In addition, a survival model and ligand-receptor network were established to explore the impact of ligand-receptor interactions on the survival of CRC patients.

RNA-seq data acquisition
The scRNA-seq data of 29 human CRC samples containing 55,543 cells was obtained from the Gene Expression Omnibus database (http:// www. ncbi. nlm. nih. gov/ geo/) from two studies, and the accession numbers were GSE132465 and GSE144735 [3]. Three bulk RNA-seq profiles of CRC samples were accessed from the TCGA database (https:// portal. gdc. cancer. gov/) and the Gene Expression Omnibus database (Accession Numbers GSE39582 and GSE17538). The RNA-seq data of CRC cell lines was downloaded from CCLE (https:// sites. broad insti tute. org/ ccle).

scRNA-seq data processing
The 'Seurat' package was used to analyze the matrix of unique molecular identifier (UMI) counts per gene. Cell selection was based on the following criteria: cells with > 1000 UMI counts, > 200 genes, < 6000 genes and < 5% mitochondrial gene expression in UMI counts. The 'scTransform' package and 'harmony' package were used for batch correction and scRNA-seq data integration. The t-distributed stochastic neighbor embedding (tSNE) algorithm was used for dimensionality reduction and cluster classification analysis. The number of principal components (NPCs) was selected based on the elbow plot, and the resolution was determined by the number of cells. Adjusted P-value < 0.05 and log2(fold change) > 0.5 were considered the cutoff thresholds for identifying marker genes. Clusters were annotated by the 'singleR' package, marker genes and annotation file of the original research. The 'Monocle' package was utilized to construct singlecell pseudotime trajectories of dendritic cells (DCs).

Ligand-receptor pair interaction strength calculation and analysis
The 1893 currently confirmed ligand-receptor pairs were downloaded from the published database and research [4]. The strength of the LRIs was calculated by referencing the research of Lei Zhang et al. [5]. The LRI strength calculation is briefly described as follows. For singlecell samples, the average ligand expression value of cell type A is multiplied by the average receptor expression value of cell type B. For bulk samples, the ligand expression value of the same sample is multiplied by the receptor expression value. The statistical significance of the ligand-receptor pairs between different pairwise cells was achieved by 1000 permutation tests. The reconstruction of intercellular spatial structure in the scRNA-seq data was performed by Cellular Spatial Organization mapper (CSOmap) based on the LRIs [6]. Cells that were located away from the LRI were defined as those that are 40 units away from the origin. Differentially expressed pairs (DEPs) were screened by the 'edgeR' package, and significantly DEPs were defined as having adjusted P-value < 0.05 and log2(fold change) > 0.5. Database for Annotation, Visualization and Integrated Discovery (DAVID, https:// david. ncifc rf. gov/ home. jsp) v6.8 was used to perform functional enrichment analyses of the differentially expressed genes, and significantly enriched terms were defined as having Benjamini-corrected P-value < 0.05. Whole genome expression profile data were analyzed by gene set enrichment analysis (GSEA 4.1.0).

Ligand-receptor subtypes classification of CRC
An unsupervised k-means subclustering algorithm was applied to classify LRI subtypes in the three bulk RNAseq datasets through ligand-receptor expression patterns. The 'ConsensusClusterPlus' package was used to perform unsupervised consensus clustering. The clustering procedure was iterated 1000 times, with 80% sampling in each iteration. The optimal number of clusters was determined by CDF curves and a consensus heatmap. The 50 most significantly upregulated ligand-receptor pairs were selected by the 'edgeR' package as the characteristic ligand-receptor pairs of each subtype. A silhouette value was used to evaluate the clustering effect, and the closer the silhouette value was to 1, the better the clustering effect was. The unsupervised subclass mapping method (SubMap; https:// cloud. genep attern. org/ gp/) was applied to evaluate the similarity between the different subtypes in the three bulk RNA-seq datasets [7]. A Bonferroni-corrected P-value < 0.05 was considered statistically significant. Then, Kaplan-Meier (KM) survival curves were used to evaluate the overall survival (OS) of patients with different LRI subtypes. The log-rank test was used to analyze the significance of survival differences between groups. The characteristic ligand-receptor pairs of each subtype were used to identify the LRI subtype of the single-cell samples. The abundance of tumorinfiltrating lymphocytes (TILs) and efficacy of immune checkpoint blockade therapy in the bulk RNA-seq  samples were estimated by ImmuCellAI (http:// bioin fo. life. hust. edu. cn/ ImmuC ellAI# !/), a tool used to estimate the abundance of immune cells and predict the response to immunotherapy from a gene expression dataset [8].
The correlation between the LRI intensity and the abundance of TILs was analyzed using Pearson correlation. Submap was utilized to analyze the consistency of the expression profiles of immunotherapy response and the LRI subtypes. The "maftools" R package was utilized to analyze the single nucleotide polymorphism (SNPs) and copy number variations (CNVs) of ligand-receptor interaction subtypes in TCGA-COREAD cohort. The mutation information of CRC cell lines was downloaded from CCLE. The GISTIC 2.0 was applied to define the amplified and deleted regions of each subtype [9]. The G-score is calculated by weighing regions of aberration against the likelihood for random occurrence through permutation test [10].

Co-culture and reactivity assay
Peripheral blood was collected from healthy volunteer from the First Affiliated Hospital of Fujian Medical University. Peripheral blood mononuclear cells (PBMC) were isolated with Ficoll-Hypaque by density gradient centrifugation. Cell culture dish were coated with 10 μg/ mL anti-CD3 (novoprotein, GMP-A018) and anti-CD28 (GMP-A063) and kept overnight at 4 °C. PBMC were culture in RPMI 1640 with 10% fetal bovine serum, 1:100 Penicillin-Streptomycin and 500 U/ml IL-2 (novoprotein, c013). Half of the medium, including IL-2, was refreshed every other day. To evaluate the activation of CD8 + T-cells, cells were washed with PBS and stained with human anti-CD3-APC(Biolegend), anti-CD4-FITC(Biolegend), anti-CD8-PE(Biolegend) for 30 min at 4 °C. Cells were washed twice with PBS and recording at flow cytometer (BD FACSCanto II).
Activated CD8 + T-cells were seeded at a 10:1 effector:target ratio in 96-well plates. Co-culture was performed in T cell medium for 24 h. Lactate dehydrogenase (LDH) kit (Beyotime) was used to measure the cytotoxicity of CD8 + T cells to CRC cells.

Generation and validation of the LRI prognostic risk score model
Univariate analysis was performed to determine the survival-related ligand-receptor pairs in the GSE39582 dataset. The survival-related pairs with P-values < 0.05 in univariate analysis were further filtered by least absolute shrinkage and selection operator (LASSO) regression analysis. Multivariate Cox regression was used to calculate the regression coefficients and the prognostic risk score. The CRC patients in the bulk RNA-seq datasets were classified into either a low-risk (low score) group or a high-risk (high score) group according to the median value. Thirty percent of the GSE39582 dataset was used as the internal validation group, and the TCGA colorectal adenocarcinoma and rectum adenocarcinoma (TCGA-COREAD) and GSE17538 datasets were used as the external validation group to jointly verify the survival prediction ability of the ligand-receptor pair risk score. The risk score combined with clinicopathological factors was used to construct an LRI prognostic model. The predictive accuracy of the prognostic model was assessed by time-dependent receiver operating characteristic (ROC) curve analysis within 3 years and 5 years. A nomogram was drawn using the 'rms' package. Calibration plots for 3-year and 5-year OS predictions were constructed to assess the calibration and discrimination of the survival model.

Discovery of LRI patterns in human CRC using scRNA-seq data
A schematic diagram of our research process is presented in Fig. 1. After batch calibration and quality control of two CRC scRNA-seq datasets, 53,537 CRC cells were further studied. The tSNE algorithm was applied to classify all cells into 33 clusters, and six main types and 25 kinds of cells were annotated by 'singleR' and the original annotation files ( Fig. 2A, Additional file 4: Fig. S1). Then, we calculated the strength of the LRIs from these cells. The heatmap and density plot of CRC tissue cells illustrated that the strength of ligand-receptor pairs was high between myeloid cells and stromal cells (Fig. 2B).
Through the circle plot of LRIs in CRC, we found that myeloid cells and stromal cells are at the center of the interaction pattern, especially endothelial cells, fibroblasts, and stromal cells (Fig. 2C). Then, we used CSOmap to reconstruct the spatial relationships in CRC tissue. We found that the interaction between myeloid and stromal cells and malignant epithelial cells constitutes the center of the tumor, and T cells and B cells do not interact closely with other cells (Fig. 2D-F). In addition, malignant epithelial cells have substantial heterogeneity. Some of the malignant epithelial cells were clustered in the spatial center of the tumor, while some were separated from the interaction center. Analysis of the DEGs in these two types of malignant epithelial cells revealed that epithelial cells far from the LRI center showed downregulation in certain terms, including cell adhesion, negative regulation of apoptotic process, antigen processing and presentation via MHC class I and positive regulation of neutrophil chemotaxis (Fig. 2G, H, Additional file 1).

Identification of LRI subtypes in colorectal cancer
To distinguish the LRI subtypes of patients with CRC, we calculated the LRI strength in the three bulk RNAseq datasets (GSE39582, GSE144735, and TCGA-COREAD). An unsupervised analysis was performed on CRC patients from the GSE39582 dataset using consensus clustering. According to the consensus heatmap, the optimal number of clusters was determined to be three (Fig. 3A). The top 50 upregulated marker pairs of each subtype in GSE39582 were identified, and subtype2 could be well distinguished in the three datasets ( Fig. 3B-G, Additional file 2). Each subtype in the different datasets was mapped through a submap, and the results showed that the LRI patterns of subtype2 in the three datasets were significantly matched (Fig. 3I). Functional enrichment analysis of each subtype's marker pairs was performed. We found that the GO-BP terms mainly enriched in subtype1 included cell adhesion, integrin receptor and G protein-coupled receptor pathway. Subtype2 was significantly enriched in immune-, inflammation-and cell adhesion-related terms. Subtype3 was enriched in terms related to the fibroblast growth factor receptor pathway, MAPK pathway and ERK cascade (Fig. 3H). In addition, we found that the prognosis of patients with subtype2 was worse than that of patients with the other two subtypes (Fig. 3J). We found that the expression of immunotherapy-related genes (PDL1, PDL2, CTLA4, CD80, and CD86) were significantly higher in subtype2 than in the other two subtypes (Additional file 5: Fig. S2A-C). Furthermore, the expression pattern of subtype2 patients matched that of patients with a poor response to immune checkpoint blockade therapy (Additional file 5: Fig. S2D). We further analyzed the differences of genomic alterations between the ligand-receptor subtypes in the TCGA-COREAD dataset. We separately showed the significantly mutated genes among the three subtypes (Additional file 6: Fig. S3A). In addition, we determined the mutations associated with targeted therapy for CRC among subtypes (Additional file 6: Fig. S3B). The mutation of oncogenic signaling pathways were analyzed, and it was found that the mutant genes in subtype2 were less than those in the other two subtypes (Additional file 6: Fig. S3C). Tumor mutation burden (TMB) is a biomarker to predict clinical response to immunotherapy, it was found that the subtype3 had the lowest TMB, and significantly lower than subtype1 (Additional file 6: Fig. S3D). The fraction of genome altered of subtype2 had lower levels (Additional file 6: Fig. S3E). We further delineate the significantly CNVs of each subtype (Additional file 6: Fig. S3F). In each subtype, functional enrichment of genes which the expression levels were affected by CNV was applied. We found that subtype1 was enriched in cytokine receptor binding and subtype2 was enriched in phospholipid binding and oxidoreductase activity, and subtype3 has no specific enrichment function.

Analysis of ligand-receptor subtype at the scRNA-seq and cell level
As mentioned earlier, we identified a ligand-receptor subtype that showed upregulation in immune-inflammation-related pairs in the bulk RNA-seq datasets. Next, cluster analysis was performed on two scRNA-seq datasets, and the results showed that the LRI patterns of KUL-01, SMC-14 and SMC-20 were consistent with the characteristics of subtype2 (Fig. 4A). The cell interactions of subtype2 CRC were explored by spatial reconstruction (Fig. 4C). We obtained a significant difference from the previous tumor model of all specimens, where T cells were more closely linked to themselves and other cells, and the density plot illustrated the appearance of another interaction center (Fig. 4D). In addition, we found that tumor cells of subtype2 was significantly enriched in PI3K-AKT, MAPK, TNF, and TGFB pathways, (Fig. 4), and the effect of inhibitors of these pathways on CRC cells vitality was verified by CCK8 experiments (Additional file 7: Fig. S4A).
In order to verify whether subtype2 affects the interaction with T cells through the above-mentioned pathways, and further affects the survival of CRC patients and the efficacy of immunotherapy, we conducted the co-culture of tumor cells and CD8 + T cells. We isolated T cells from the peripheral blood of healthy people and activated the main tumor-killing component, CD8 + T-cells (Fig. 4F). The activated CD8 + T-cells were co-cultured with tumor cells treated with different pathway inhibitors. The results of LDH assay showed that PI3K inhibitor and MAPK inhibitor could significant promoted the tumor-killing ability of CD8 + T-cell in SW480 and HT29 (Fig. 4G, J right side). After TNF receptor inhibitor treatment, LDH was significantly reduced and the tumor killing ability was weakened in HCT116, SW480 and CACO2 (Fig. 4H right side). After TGF-β receptor inhibitor treatment, the tumor killing ability was weakened in SW480 and CACO2, HCT116 and HT29 had the opposite phenotype ( Fig. 4I right side). By analyzing the transcription profile and mutation profile of these cell lines, it is found that significant differences exist in genetic background ( Fig. 4G-J left side, Additional file 7: Fig. S4B).

Effect of LRI on tumor immune infiltration in CRC
To investigate the relationship between TILs and LRIs, we analyzed the correlation between LRI intensity and the abundance of B cells, T cells, DCs, macrophages, and NK cells in three datasets (Fig. 5A). The results demonstrated that LRIs significantly influenced the infiltration of DCs and macrophages in all three datasets. The intersecting ligand-receptor pairs of the three datasets showed that TGFB1:SDC2, HEBP1:FPR3, and CD14:ITGB1 were significantly positively correlated with the abundance of DCs and macrophage infiltration (Fig. 5B). The other 9 ligand-receptor pairs were significantly associated with DC infiltration. We further explored the role of these ligand-receptor pairs on DCs using scRNA-seq data. We assessed the infiltration of various types of cells in the sample (Fig. 5C). We analyzed the correlation between the strength of the LRI and the proportion of DC infiltrates, as well as the significance of the difference in the strength of ligand-receptor pair interactions between DCs and other cells (Fig. 5D, E, Additional file 8: Fig. S5). Through correlation analysis and strength analysis of the LRIs, we found that through ICAM1:IL2RA, the specific interaction among mDCs, DCregs and regulatory T cells affects the infiltration of DCs in CRC tissues. In addition, ICAM1 expressed on mDCs and Dregs interacts with ITGAM on M1 macrophages, significantly and specifically affecting the infiltration of DCs. The clustering process for monocytes and DCs is shown in Additional file 9: Fig. S6.

Generation and validation of the LRI-based prognostic risk score model to predict CRC patient survival
As mentioned earlier, the different LRI subtypes are related to different prognostic outcomes. To better explore the influence of LRIs on the survival of CRC patients, we constructed a prognostic prediction model for CRC patients. Univariate Cox analysis was performed and identified 421 prognosis-associated ligand-receptor pairs in the GSE39582 training set (Fig. 6A, Additional  file 3). Bubble plots showed the ligand-receptor pairs that influence prognosis in CRC patients, and these pairs were mainly involved in cell proliferation, extracellular matrix, cell communication, and immune response (Fig. 6B). LASSO regression followed by multivariate Cox analysis was performed, and 30 survival-related pairs were identified (Additional file 10: Fig. S7A-C). 30 ligand pairs were further analyzed by functional enrichment analysis (Additional file 11: Fig. S8). The prognostic risk score was developed based on these pairs. Risk scores were calculated for all patients in the training set (GSE39582) and validation set (GSE144735 and TCGA-COREAD), and the patients were divided into either a high-risk (high score) group or a low-risk (low score) group using the median value of the risk score as the cutoff value (Fig. 6C). KM survival analysis demonstrated that patients in the high-risk group had significantly poorer OS than those in the low-risk group (Fig. 6D). We constructed a prognostic model including the risk score combined with clinicopathological data and drew a nomogram (Fig. 6G). Time-dependent ROC analysis indicated that this model showed excellent performance in predicting the 3-year and 5-year survival rates of CRC patients (Fig. 6E). The calibration plots showed that the predictions were in concordance with the actual observations for the 3-year and 5-year survival rates in the training set and validation set (Fig. 6F). The survival-related LRI network was delineated to further demonstrate the intercellular interactions of these pairs (Fig. 7). These findings suggested the appreciable reliability of the prognostic model, which can be applied in diverse CRC patients.

Discussion
In the tumor microenvironment, the communication between different cell types is associated with mechanisms of oncogenesis, tumor progression, therapeutic resistance, immune infiltration, and inflammation [2,11]. scRNA-seq is an effective method to analyze the LRIs that occur in CRC [3,5]. With the development of cell-cell interaction-based immunotherapy in colorectal cancer, it is crucial to understand the existing ligandreceptor pairs further. By analyzing scRNA-seq data from human and mouse CRC, Lei et al. also revealed that DCs and macrophages are key regulators of cell communication in the tumor microenvironment and further identified the specific target of anti-CSF1R treatment [12]. Our study of the existing ligand-receptor pairs can provide a holistic and realistic view of cell-to-cell communication in the tumor microenvironment. By analyzing LRI strength, we found that myeloid cells and stromal cells were the core of LRIs. As early as 2011, Piero Dalerba carried out single-cell research, exploring the multilineage differentiation processes and cellular heterogeneity of CRC [13]. Intratumoral heterogeneity plays an important role in the resistance to cancer therapy, and finding new intratumoral cell types can help to identify therapeutic targets and explore drug resistance [14]. We identified a subtype (subtype2) of malignant epithelial cells characterized by low cell adhesion, inhibition of neutrophil chemotaxis, and inhibition of antigen presentation, which has the potential for immune escape and distant metastasis. We further analyzed the differences between subtypes at the SNP and CNV levels, and found that subtype2 was characterized by lower oncogenic-related mutations and lower fraction genome altered. The spatial reconstruction of tumors based on LRIs, which was accomplished by using CSOmap, a computational tool for inferring cell-to-cell interactions, makes up for the lack of spatial information in single-cell data and provides validation for the ligand-receptor subtypes. In addition, we identified an LRI subtype in the transcriptome of CRC patients that shows upregulation in chemokines and intercellular adhesion-related ligand-receptor pairs. Chemokines play an important role in immune and inflammatory responses, are involved in CRC progression and metastasis and are also associated with poor prognosis, which is consistent with our findings [15,16]. In addition, analyzing the related pathways of tumor cells in the subtypes shows that pathways related to immune escape are significantly enriched, including PI3K/AKT, TGF-β, and MAPK pathways. Previous studies have suggested that inhibiting these pathways can reverse tumor immune escape [17][18][19]. Our research results also found that these pathways inhibitors can enhance the cytotoxicity of CD8 + T-cells when co-cultured with some CRC cell lines. However, the role of TNF in tumors is a double-edged sword. On the one hand, it kills tumor cells through cytotoxicity, and on the other hand, it also participates in the occurrence and development of tumors. TNFR inhibitors can significantly decreased anti-tumor activity of CD8 + T-cells coculture in vitro experiment. However, co-culture systems may have different responses to inhibitors due to differences in genetic background between different cell lines and the effect of inhibitors on tumor cell viability.
Through the analysis of the influencing factors of TILs, we revealed that the infiltration of monocytes in CRC tissues was significantly affected by LRIs. We also found that ICAM1:IL2RA and ICAM1:ITGAM could increase the infiltration of DCs by interacting with Tregs and monocytes/macrophages, respectively. ICAM-1 plays an important role in leukocyte-mediated inflammation and T cell activation and was found to interact with ITGAM as early as 1990 [20]. Previous studies have shown that ITGAM-encoded CD11b (MAC-1) is highly expressed on macrophages and mediates intercellular adhesion [21], which was also observed in the CRC scRNA-seq datasets. Few studies on the function of ICAM1:ITGAM [22] and no studies on this pair between DCs and macrophages have been found before. We found that the infiltration of mDCs in tumor tissues can be promoted by the expression of ICAM1 and interaction with ITGAM expressed in macrophages. Furthermore, several studies have indicated that Tregs can downregulate the maturation process of DCs and abrogate the antigen-presenting capacity of DCs [23][24][25]. Previous studies have shown that IL2 is crucial for Treg development and maturation, and the long-term survival of Tregs requires continuous IL2 signaling [26]. Here, we speculate that ICAM1 inhibits IL2RA through competitive binding with IL-2 and inhibits the inhibitory effect of Tregs on DCs.
In previous research, prognostic risk models of CRC based on LRI have not been constructed; however, many prognostic models have included ligand-receptor pairrelated genes [27][28][29][30], and several genes in this model have been shown to influence CRC development [31][32][33][34][35][36]. Our survival model based on LRIs has good predictive efficacy for 3-year and 5-year survival rates, providing a new tool for the early assessment of adjuvant intervention in patients with CRC.
The current study has some limitations. The singlecell data lacked spatial information, and we used spatial reconstruction, but further spatial transcriptome verification is needed. The factors affecting TILs in this study need to be further verified by in-depth experiments. This prognostic model has yet to be further validated for clinical effectiveness in a larger cohort of CRC patients.

Conclusions
We delineated the LRI network of CRC using scRNA-seq and bulk RNA-seq datasets and calculated each ligandreceptor pair's strength across different cells. We identified a malignant epithelial cell subtype with potential for distant metastasis and immune escape and a CRC subtype that showed upregulation in immunoinflammatory ligand-receptor pairs, which is also associated with CRC immunotherapy response and survival. And it was found that this CRC subtype induces the hypo-responsiveness of CD8 + T-cells by regulating PI3K/AKT, MAPK, TGF-Beta and TNFA pathways. We explored ligand-receptor pairs affecting TILs and found that ICAM1:IL2RA and ICAM1:ITGAM specifically affected the infiltration of DCs. We also constructed a prognostic prediction model for colorectal cancer based on LRIs to predict patients' prognostic risk and provide suggestions for further guiding treatment strategies.