LncRNAs dominantly co-regulate immune-related functions across cancers.
To explore the function of lncRNAs, we proposed a multiple-step model to identify lncRNA-lncRNA functional cooperation in each cancer type (Fig. 1a). We first identified the cancer context-specific lncRNA-target pairs by integrating multi-omics data. The lncRNA-lncRNA cooperations were further identified based on the cancer-context lncRNA-target regulation across 20 cancer types (Supplementary Table S1)14. In each cancer type, one lncRNA pair was considered as cooperation, if and only if they could co-regulate at least a functional module, which was consisted of their shared targets. We identified 6449 lncRNA-lncRNA cooperation pairs among 2145 lncRNAs. The number of cooperative lncRNA accounted for only 16.86% of all the lncRNAs, but these lncRNAs were tightly connected and assembled into an lncRNA-lncRNA cooperative network in pan-cancer (Fig. 1b). Then, we found that approximately 45.13% (968 / 2145) of lncRNAs exhibit cooperation in at least two cancers, on the other hand, 94.73% (6109 / 6449) of lncRNA cooperative interactions are cancer-specific (Supplementary Figure S1). These results indicate that although near half of lncRNAs play cooperative roles in more than one cancer, the cooperative partners of these lncRNAs might be changed in different cancers. Furthermore, by revealing the topological structure with a power-law degree distribution and higher clustering coefficients than randomly linked networks (Supplementary Figure S2a and S2b), the cooperative network of pan-cancer exhibited the scale-free and modular characteristics, indicating that lncRNA-lncRNA cooperative interactions influence each other and effectively exchange regulated information both at a global and at a local scale. Through analysis of the genomic location and expression, we found that the majority of the cooperative pairs regulate the same function module in trans (Supplementary Figure S3, p = 0.003) and the cooperative lncRNAs tend to have similar expression pattern in cancer types with similar tissue origin (Supplementary Figure S4), such as lower-grade glioma (LGG) and glioblastoma multiforme (GBM).
To further explore the functional characteristics of cooperative lncRNAs, we ranked the functions co-regulated by lncRNAs in descending order according to the number of cooperative lncRNA pairs. We found that all the top functions are related with cancer development and progression, including cell cycle checkpoint, negative regulation of DNA replication and especially many immune-related functions. The proportion of immune-related functions increased with the order and is averagely as high as 58% in top10 to 40 (Fig. 1c). Particularly, the trend of lncRNAs co-regulating immune function is more obvious in lung squamous cell carcinoma (LUSC) and bladder carcinoma (BLCA), and respectively 100% and 90% of top 10 functions are related with immune. As shown in Fig. 1d, 12 (67%) of the top 18 functions are associated with immune, such as T cell activation, lymphocyte activation and interleukin 8 production. Up to 10.7% (696 / 6449) cooperative lncRNA pairs are involved in regulation of T cell activation. In LUSC, each of these immune functions is co-regulated by an average of 208 lncRNA pairs. The lncRNA MAGI2-AS3 could co-regulate cytokine biosynthetic and secretion, interleukin biosynthetic and T cell activation with 9 lncRNAs, and was overexpressed in LUSC (Fig. 1e). It has been proved that MAGI2-AS3 could up-regulate suppressor cytokine signalling 1 and suppress the proliferation of NSCLC cells15. The result suggests that co-regulating immune function with other lncRNAs may be another mechanism in cancer development for MAGI2-AS3. Moreover, we found that lncRNAs with more neighbours in the cooperative network tend to co-regulate more immune-related functions (Fig. 1f). These results imply that cooperative lncRNA pairs may contribute to carcinogenesis by regulating immune-related functions in multiple cancers. Our finding may provide a new entry that cooperative lncRNAs have contributions to modify tumor immune microenvironment.
Landscape of lncRNA immune-related co-regulated networks.
To further gain insight into the roles of cooperative lncRNAs in immunity, we extracted all the cooperative lncRNA pairs regulating immune-related functions in each cancer (Supplementary Table S1) and constructed LICNs in 17 cancers respectively based on these pairs(Supplementary Figure S5 and Figure S6). Furthermore, an immune lncRNA-lncRNA cooperative regulatory landscape in pan-cancer was constructed by integrating cancer-specific LICNs (Fig. 2a). This network involves 505 IC-lncRNAs and 1628 IC-lncRNA- IC-lncRNA cooperative interactions across 17 cancer types. Furthermore, we found that the expression of IC-lncRNAs is higher than the other lncRNAs in most of cancers (Supplementary Figure S7a). The IC-lncRNAs pairs also exhibited significantly stronger expression correlation than random in 10 cancers (Supplementary Figure S7b). Thus, we proposed that similar expression patterns might help IC-lncRNAs perform cooperative functions.
Next, we summarized these immune functions into 14 categories and found that these functions are widely distributed, involving multiple levels of immune processes (Fig. 2b). Overall, lymphocyte activation was co-regulated by the largest number of IC-lncRNAs. Moreover, we also found that some function categories with internal connections share more cooperative IC-lncRNAs, such as lymphocyte activation and T cell activation (Fig. 2c). The trend of the co-regulation for IC-lncRNAs to lymphocyte activation may offer new candidates for cancer immunotherapy. Next, we explored the co-regulation of IC-lncRNAs to immune function categories in each cancer. On average, each cancer type yielded 111 co-regulations among 63 lncRNAs participating in the immune process. The higher number of IClncRNA co-regulations were identified in the cancer types where immune checkpoint-blocking drugs are applicable in clinical (Fig. 2d)16. For example, there are 99 IC-lncRNAs and 737 co-regulation pairs in LUSC (Fig. 2e). These IC-lncRNAs accounted for ~ 23.5% (505 / 2145) of all cooperative lncRNAs (Fig. 2d). In addition, the LICN of SKCM contains 110 cooperations among 73 IC-lncRNAs (Fig. 2g). For example, IC-lncRNA LINC00324 which was highly expressed than others co-regulated 10 immune-related functions, such as lymphocyte activation, cytokine biosynthetic and secretion, and interleukin biosynthetic with 49 IC-lncRNAs. A previous study concluded that LINC00324 could regulate the expression of FasL, an apoptosis suppressor concentrated in immune cells and cancer cells, which has been shown to play a vital role in immune evasion17. Higher number of IC-lncRNAs regulated the T cell activation and lymphocyte activation across cancer types (Fig. 2d). Particularly, 489 IC-lncRNA pairs are involved in regulation of T cell activation in LUSC (Fig. 2f). Activation of T cells has become an important way to promote antitumor immune response in lung cancer18. Moreover, 51 IC-lncRNAs co-regulate T cell differentiation in SKCM (Fig. 2h). The dysregulation of T cell differentiation has been shown in melanoma progression 19. These results imply that the IC-lncRNAs might play critical roles in tumor immune microenvironment.
Conserved and rewired LICN hubs in tumor immunity.
A few nodes with a large number of neighbours as hubs hold the nodes together in each network. The presence of hubs seems to be a general feature of all biological networks, and these hubs fundamentally could determine the behaviour of networks20,21. To investigate the crucial nodes in the pan-cancer LICN, we first identified the top 10% of nodes with the highest connectivity as hub IC-lncRNAs. Then, we assessed the roles of these hub IC-lncRNAs among the different cancer LICNs and grouped these hubs into three categories: common hubs, cancer-specific hubs and other hubs. The common hubs signify the central roles of IC-lncRNAs in more than one cancer; the specific hubs identify IC-lncRNAs with specific central roles in a given cancer. In total, 39.2% (20 / 51) of hubs are common hubs, 39.2% (20 / 51) are specific hubs, and the remaining 21.6% are other hubs. We found that all the hub IC-lncRNAs are extensively involved in immune functions (Fig. 3a). The distribution of three categories of hubs is different across cancers; both LUSC and lung adenocarcinoma (LUAD) contain many common hubs. We further investigated whether three types of hubs have different roles in cancer immune. We found that most of these hubs co-regulated the majority of the immune function categories, such as T cell activation, cytokine biosynthetic and secretion and B cell activation (Fig. 3b). On the other hand, several function categories tend to be specifically co-regulated by different categories of hubs, such as inflammatory response, which are regulated only by the common hubs AC109826, RP11-25K19 and RP11-420G6, and interferon gamma biosynthetic are specifically regulated by the cancer-specific hubs (Fig. 3b). These results suggested that these hub IC-lncRNAs might play important roles in immune-related functions.
We further explored the contribution of the IC-lncRNAs to cancer development. First, we found that 18 out of 51 hub IC-lncRNAs could participate in cancer hallmarks (Fig. 3c). 9 hub IC-lncRNAs, particularly common hubs, regulate evading immune detection. At the same time, these IC-lncRNAs could regulate other non-immunological hallmarks. For example, other hub IC-lncRNA CTC-241F20 specifically regulated insensitivity to antigrowth signals. The results suggested that IC-lncRNAs may synergistically regulate immune-related functions and other cancer-related processes to involve in the development of cancer. We next found that 8 hub IC-lncRNAs are experimentally validated cancer-related lncRNAs, including 3 common hubs and 5 specific hubs (Fig. 3d). For example, the common hub MIR155HG with 17 partners to co-regulate B cell and T cell activation, cytokine biosynthetic and secretion has been associated with multiple cancer types, including CESC and KIRC. MIR155HG is also a hub IC-lncRNA respectively in CESC and KIRC LICNs. MIR155HG was known as an oncogenic lncRNA and dysregulated by PRDM1 in natural killer/T-cell lymphoma22. Then, we further investigated the expression of hub lncRNAs. We found that the average expression level of other hubs is higher than common hubs (Supplementary Figure S8, two-sided t test p = 0.011). These results suggest that IC-lncRNAs may provide a candidate list of cancer-related lncRNAs. Then, we investigated the association between hub IC-lncRNAs and anticancer drugs. We found that 66 out of 84 Food and Drug Administration (FDA) approved anticancer drugs could target hub IC-lncRNA target genes. The common and specific hubs have significantly more anticancer drug targets than other hubs, respectively with an average of 29 and 26 drug targets (Fig. 3e, two-sided t-test, p = 4.39e-05, p = 8.55e-04). For example, some target genes of the common IC-lncRNA RP11-445H22 have been used as anticancer drugs targets, such as CTLA4 as target of radotinib in cutaneous melanoma, EGFR and ERBB2 as target of afatinib in metastatic non-small cell lung cancer. The results suggest that target genes of hub IC-lncRNAs are more likely to be druggable and may be potential targets of anticancer drugs.
Tumor immune microenvironment is very complex involving a range of cell types and molecular mechanisms, especially tumor-immune cell interactions23. We surprisingly found that some hub IC-lncRNAs (RP11-445H22, RP11-25K19, and LINC00324) could regulate the genes mediating the direct interaction between tumor cells and immune cells such as, T cells, NK cells, DC cells and B cells (Fig. 3f). The common hub RP11-445H22 has the most immune cooperative partners and could regulate CRTAM, CADM1, CXCR3, CXCL10, PDCD1 and CD274 to mediate tumor - T cell interactions as shown in Fig. 3f, which play a role in promoting tumors and suppressing immunity. Among these receptor and ligand genes, CXCL10 and CXCR3 received strong regulation from RP11-445H22 and the regression coefficient respectively was 0.63 and 0.55 in the above lncRNA target identification model. RP11-445H22 is a KCNK15 and WISP2 antisense RNA and it has been found as an oncogene in lung cancer24,25. The regulation to tumor-immune cell interactions for RP11-445H22 may be another mechanism to participate in tumorigenesis. LINC00324 is a specific hub in SKCM, and could be involving in tumor cell-NK cell interactions by regulating TNFSF18, TNFRSF18, CRTAM and CADM1. Furthermore, the three hub IC-lncRNAs could co-regulate many immune functions with 363, 220 and 50 partners respectively, such as lymphocyte activation, inflammatory response and interleukin biosynthetic (Fig. 3g). The regulation of IC-lncRNAs to these genes mediating interactions can be leveraged to understand the cross-talk between tumors and immune cells and provide a promise to develop novel drugs or therapeutic strategies.
IC-lncRNAs were correlated with immune cell (B cells and T cells) infiltration.
Tumor immune microenvironment is broadly populated with immune cells26. Therefore, we reasoned that if these IC-lncRNAs participate in tumor immune microenvironment regulation, then they would be more likely to be highly expressed in immune cells and to be correlated with immune cell infiltration in tumors. Firstly, we found that a significantly higher proportion of IC-lncRNAs (90.7%) was expressed in immune cells by analysing the immune cells RNA-seq datasets (Supplementary methods and Fig. 4a, two-sided fisher’s exact test, p = 8.08e-92). In particular, 93.5% of the IC-lncRNAs co-regulating B cell-related functions were expressed in B cells. Moreover, we also found that the IC-lncRNAs co-regulating T cell-related functions were significantly more highly expressed than other lncRNAs in T cells and the B cell-related IC-lncRNAs exhibited significantly higher expression in B cells (Fig. 4b, two-sided wilcoxon test, B cells p = 1.32e-10, T cells p = 2.67e-33). These results suggest that IC-lncRNAs exhibit higher expression in immune cell populations.
We next estimated the associations between expression of cooperative lncRNAs and infiltration levels of six immune cells (B cells, CD4 T cells, CD8 T cells, macrophages, neutrophils, and dendritic cells) in each cancer. We found that a large number of cooperative lncRNAs were correlated with immune cells infiltration in most cancers, such as BRCA, LGG, LUSC, SKCM and THCA (Fig. 4c). In particular, 56.3% and 50.2% of cooperative lncRNAs were correlated with CD8 T cell infiltration in LUSC and in SKCM respectively. The correlations between the expression of IC-lncRNAs co-regulating B cell-related functions and B cell infiltration are significantly higher than others in LUSC (Fig. 4d, two-sided wilcoxon test,p = 1.16e-12). For example, the PCED1B-AS1 co-regulated B cell activation with 8 IC-lncRNAs and the correlation coefficient between the expression and B cell infiltration is as high as 0.69 (Fig. 4e). PCED1B-AS1 could govern aerobic glycolysis, and its overexpression was closely related to larger tumor size and shorter survival time27. In addition, the correlations between the expression of IC-lncRNAs co-regulating T cell-related functions and T cell infiltration are significantly higher than others in SKCM, such as RP11-861A13, RP11-25K19, LINC00324 (Fig. 5f, two-sided wilcoxon test, CD4 T cell p = 3.75e-05, CD8 T cell p = 2.22e-03). LINC00324 co-regulating T cell proliferation, activation and differentiation were correlated with CD4 and CD8 T cell infiltration with correlation coefficients 0.42 and 0.45 respectively (Fig. 4g). Indeed, LINC00324 has been found to be overexpressed in cancer and correlated with aggressive cancer progress 28.
We further used fisher’s exact test to investigate whether the IC-lncRNAs were likely to be associated with immune cell infiltration. The infiltration levels of CD8 T cells and CD4 T cells were used to divide the patients into two groups respectively (high immune cell infiltration group and low immune cell infiltration group), and it was found that a lot of IC-lncRNAs co-regulating the T cell functions are differently expressed in majority of cancer types (Fig. 4h, two-sided fisher's exact test). For example, the expressions of AP001055 and CTD-2587H19 which regulate T cell-related functions are significantly high in LUSC patients with high CD4 T cell infiltration (Two-sided t test, AP001055 FDR = 1.32e-05, FC = 2.05; CTD-2587H19 FDR = 1.04e-05, FC = 2.35) (Fig. 4I). In SKCM, RP11-861A13 has significant higher expression in the group with high infiltration of CD8 T cells (Two-sided t test, FDR = 1.29e-07, FC = 2.17). In summary, these results suggest that IC-lncRNAs exhibit higher expression in immune cells and are associated with immune cell infiltration, further validating the roles of the IC-lncRNAs in tumor immune microenvironment.
IC-lncRNA cooperation reveals similar regulation of cancer immune microenvironment.
Several lines of evidence have indicated that sharing molecular features can reveal similar carcinogenic mechanisms between cancer types29–31. Here, we hypothesized that if the cancer types exhibit more similar IC-lncRNA cooperation patterns, they are more likely to be with similar immune mechanism. Then, we computed a paired similarity score based on integrating the expression IC-lncRNAs and structure of the LICNs in each cancer (Materials and Methods). We found that some cancers showed greater similarity to each other than to other cancers, such as LUSC, BLCA, CESC, head and neck squamous cell carcinoma (HNSC), BRCA and SKCM (Fig. 5a and Supplementary Results), in which a general immune related trend has emerged32. Especially, LUSC, BLCA and CESC were closely clustered together. 41.2% (74/177) IC-lncRNAs were shared by at least two cancers (Supplementary Figure S9). We focused on IC-lncRNAs that were common for the three cancers (Fig. 5b). We found that these IC-lncRNAs were widely involved in co-regulation of immune-related functions, such as leukocyte activation, lymphocyte activation and T cell differentiation. We further discovered that the expression of shared IC-lncRNAs was significantly higher than other lncRNAs in three cancers (Fig. 5c, two-sided t test, all p < 9.43e-22). Both lung cancers and bladder cancer have shown relatively effectiveness under immunotherapy using checkpoint blockade33. These observations suggest that IC-lncRNA cooperations might operate in cancer types with similar immune mechanism.
Numerous studies elucidate how various components of the immune system control or contribute to cancer progression, thus revealing their prognostic value34. We next investigated whether these IC-lncRNAs are associated with survival of cancer patients. We defined risk scores for each patient based on the expression of shared IC-lncRNAs (Materials and Methods). Using the median risk scores as the threshold, cancer samples were classified into two groups with significantly different overall survival rates respectively in all three cancers (Fig. 5d-f). In LUSC, 21 IC-lncRNAs are protective factors which are highly expressed in the low-risk group and other 15 IC-lncRNAs acted as risky factors of which increased expression is associated with a poor survival outcome (Fig. 5g and h). For example, MIR155HG is significantly down-regulated in high-risk group as a protective factor, hazard ratio of which is 0.77 (Two-sided t test, p = 7.27e-05). It has been shown that the expression of MIR155HG could predict overall survival in multiple cancers35. The approach based on immune lncRNA-lncRNA cooperation could effectively identify prognostic biomarkers. Collectively, these results suggest that the cooperative pattern based on both the structure of LICNs and the expression of IC-lncRNAs could reveal new cancer clusters from immune view.
LICN contribute to subtype classification of SKCM.
The skin cutaneous melanoma is a malignancy of melanocytes, which accounts for the most number of deaths from skin cancer36. Therefore, we investigated to what extent the IC-lncRNAs can be applied to SKCM molecular subtyping. We first identified the top 10% IC-lncRNAs with larger expression fluctuation and further filtered by hubs in SKCM. Finally, 6 IC-lncRNAs (RP11-71G12, RP11-555F9, RP11-367G6, ITGB2-AS1, AP000233 and AL928768) were obtained. These 6 IC-lncRNAs participate in immune functions such as regulation of B cell activation, immune response, lymphocyte activation and immune system process. These IC-lncRNAs were correlated with immune cell infiltration and they were also more highly expressed in B cells (Supplementary Figure S10A and S10B, two-sided wilcoxon tests, B cell p = 0.0015). Next, we found that the SKCM patients can be classified into two subtypes based on the expression of the 6 IC-lncRNAs (Fig. 6a). These IC-lncRNAs are highly expressed in group 2 patients (Fig. 6b, two-sided t test, p = 2.29e-3). There are no significant differences in age, gender, and cancer stages observed between the two subtypes. We further compared six immune cell infiltration levels between two subtypes respectively. We found that group2 patients consistently have significantly higher infiltration levels than group1 patients for all immune cells (Fig. 6c, two-sided wilcoxon test, all p < 1.80e-09). Particularly, average infiltration level of CD8 T cells in group2 is 25%, while infiltration level is only 8% in group1. At the same time, the average infiltration level of dendritic cells (DC) of group2 patients was almost double than that in group1. It has been shown that the presence of mature DC within tumors is a positive prognostic factor in melanoma patients and melanoma-associated DCs is also related with immunotherapy37.
It has been demonstrated that some indicators are useful biomarkers for predicting the immune response, such as immune cytolytic activity (CYT) and major histocompatibility complex (MHC) scores and the immune score38,39. We further investigated the distributions of these scores among SKCM patients. We found that group2 patients had significantly higher CYT, MHC and immune scores than group1 (Fig. 6d, two-sided wilcoxon test, all p < 1.23e-12) Treatment with immune checkpoint blockade has transformed the outcome for patients with melanoma, such as programmed cell death protein 1 (PDCD1) inhibitors, CTLA-4 inhibitor40. Therefore we compared the expression of immune checkpoint genes between two subgroups. We found that PDCD1, CD274 and CTLA4 are significantly higher expressed in group2 (Fig. 6d, two-sided wilcoxon tests, all p < 3.25e-09). At the same time, the expression of all the 6 IC-lncRNAs are significantly positively correlated with PDCD1 and CD274 (Fig. 6e, pearson correlation, p < 0.05). Each IC-lncRNA was more closely related to PDCD2 than CD274 in expression. In particular, the correlation coefficient between IC-lncRNA ITGB2-AS1 and PDCD1 was as high as 0.58 (Fig. 6f). IC-lncRNA ITGB2-AS1 is antisense to ITGB2 and has been shown to be highly expressed and closely related to immunosuppression in AML41. ITGB2-AS1 may be as a potential immunotherapeutic target in cancers.
The degree of tumor infiltration by immune cells can predict a patient's clinical outcome in many cancer types42,43. Therefore, we explored the prognostic implications of molecular subgroups. We found that patients with group2 had significantly better overall survival than group1 (Fig. 6g, log-rank test p = 0.02). These results suggest that IC-lncRNAs identify different SKCM subtypes (group1 immunosuppressive subtype and group2 immunoactivated subtype) with remarkable immunology diversity (Fig. 6h), which is helpful to improve personalized cancer management.