Gene Regulatory Network Rewiring in the Immune Cells Associated with Cancer

The gene regulatory networks (GRNs) of immune cells not only indicate cell identity but also reveal the dynamic changes of immune cells when comparing their GRNs. Cancer immunotherapy has advanced in the past few years. Immune-checkpoint blockades (i.e., blocking PD-1, PD-L1, or CTLA-4) have shown durable clinical effects on some patients with various advanced cancers. However, major gaps in our knowledge of immunotherapy have been recognized. To fill these gaps, we conducted a systematic analysis of the GRNs of key immune cell subsets (i.e., B cell, CD4, CD8, CD8 naïve, CD8 Effector memory, CD8 Central Memory, regulatory T, Thelper1, Thelper2, Thelp17, and NK (Nature killer) and DC (Dendritic cell) cells associated with cancer immunologic therapies. We showed that most of the GRNs of these cells in blood share key important hub regulators, but their subnetworks for controlling cell type-specific receptors are different, suggesting that transformation between these immune cell subsets could be fast so that they can rapidly respond to environmental cues. To understand how cancer cells send molecular signals to immune cells to make them more cancer-cell friendly, we compared the GRNs of the tumor-infiltrating immune T cells and their corresponding immune cells in blood. We showed that the network size of the tumor-infiltrating immune T cells’ GRNs was reduced when compared to the GRNs of their corresponding immune cells in blood. These results suggest that the shutting down certain cellular activities of the immune cells by cancer cells is one of the key molecular mechanisms for helping cancer cells to escape the defense of the host immune system. These results highlight the possibility of genetic engineering of T cells for turning on the identified subnetworks that have been shut down by cancer cells to combat tumors.


Introduction
The human body has an intricate ability to combat genetic disorders such as cancer. The immune system could combat cancer cells by recognizing cancer cells and killing them via sophisticated molecular mechanisms [1]. The immune system has different types of immune cells, which have different functions [2,3]. For example, T cells play an active role in combating cancer cells, while B cells play an innate role in attacking cancer cells. The relationship between immune cells and tumors is evolving with time and cancer progression. Cancer cells bypass immunological suppression and proliferate uncontrollably [4]. In this process, cancer cells may send out molecular signals that can modify the gene regulatory networks of the major immune cells. Interestingly, this process could also provide infiltrating cells) and from healthy donors (i.e., for CD8 naïve, CD8 Effector memory, and CD8 Central Memory T cells) [21]. The data source above has been already filtered to have the standard peaks through the systematic and conventional protocol. These high-quality peaks were used for constructing GRNs. In order to construct the receptors' regulatory subnetworks of each immune cell subset, we obtained a list of receptors (i.e., an immune cell receptor list) that are regulated in the immune cells from a recent study [3]. Furthermore, we used UpSet [22] to visualize intersecting sets in order to effectively analyze the unique and common genes, TFs (Transcription factors), and receptors across nine immune cell subsets. This tool could accurately establish Venn diagram relationships across multiple datasets and systematically illustrate the connections between multiple datasets.
Finally, in order to incorporate gene expression profiles into GRNs derived from chromatin profiles, we obtained RNA-seq data from Tirosh et al. (GSE72056) [23] and Pulko et al. (GSE80306) [24]. Gene expression profiles of the T cells of healthy people and the tumor infiltrating T cells have been generated in these studies. Further, we normalized the RNA-seq data and used the seaborn package, which was coded in python to produce a heatmap of the RNA-seq data from above [25].

Regulatory Network Construction
Using the high-quality peaks from either the DNase-seq or the ATAC-seq files for immune cell subsets, we constructed the GRNs of the immune cells. This was conducted by: (1) annotating the downstream genes for each peak; (2) finding the transcription factors that can bind to the peaks; and (3) mapping the downstream genes to the binding transcription factors (TFs).
In order to annotate the downstream genes that are associated with each peak, we employed Hypergeometric Optimization of Motif EnRichment (HOMER) [26], which is an efficient tool for discovering up-and down-stream genes associated with a DNase-seq or the ATAC-seq peak. It is a cluster of C++ and Perl programs, which can be used in command line operating systems such as Linux. We set the cutoff at 2 kb ± of the peak position to filter the genes that are closely associated with a peak. Next, we applied the Multiple Em for Motif Elicitation (MEME) suit to obtain the list of TFs that putatively bind to the peaks [27]. MEME is an efficient tool for uncovering DNA-binding motifs in a class of related protein or DNA sequences. We specifically used the Find Individual Motif Occurrences (FIMO) program of the MEME suite to estimate the occurrence of particular DNA-binding sites and the corresponding proteins that bind to them [28]. FIMO compares the sequences of the peak with the DNA-binding motif databases to establish protein binding regions. The DNA-binding motif database is from the MEME website (http://meme-suite.org/db/motifs/). The TFs were filtered based on the P-value, and a default P-value of 1 × 10 −4 was taken as the threshold value.
Once the TFs and their regulated downstream genes have been determined, a shell script was used to link the TFs with their regulated genes; then the GRNs were constructed. Subsequently, to find the subnetworks that are associated with the cell receptors of an immune sunset, we mapped the network genes to the immune cell receptor list and identified their upstream TFs. The unique cell receptors of an immune cell subset and their regulators form a subnetwork of the cell receptors of that immune cell subset.

Analyzing of the TFsand Hub Regulators in the Networks
Once the networks have been constructed, we evaluated the most important TFs by using the network approach. We utilized network analysis algorithms to find out which TF had the most downstream genes under its control for each immune cell subset.
Further, to compare the results obtained from the analysis mentioned above, we employed a 'PageRank algorithm' [29] to identify hub regulators.

Network Visualization and Comparison
To conduct an analysis and comparison of the GRNs and regulatory circuits, we applied Gephi [30] and Cytoscape [31] the GRNs for analyses. Both tools are effective for analyzing and visualizing interactive networks and biological pathways. Moreover, they were used for integrating and annotating networks with various biological data such as gene expression profiles. The GRNs of each immune cell subset were analyzed using a shell script and then visualized in Gephi. Furthermore, Cytoscape was used to identify key regulators of the GRNs.
To elucidate the signaling pathways that are associated with uniquely regulated genes in the GRN of the T cells in PBMC (peripheral blood mononuclear cell) versus the TIC T cells, we conducted an analysis of their GRNs using Enrichr [32]. Enrichr is an extensive web server to perform gene set enrichment analyses with over 184 annotated gene sets from 102 gene set libraries [33].

Characterize Gene Regulatory Networks of Key Immune Cells Associated with Cancer Immunotherapy
The key immune cells associated with cancer immunotherapy have been reported. Among them, we selected nine key immune cell subsets in PBMC (Table 1), the chromatin profiles of which have been determined by either DNase-seq or ATAC-seq. To construct the GRN of an immune cell subset, we obtained the DNase-seq or ATAC-seq peak files of the immune cell subset to identify their chromatin open regions. Based on these regions, we scanned the DNA-binding motifs and their associated TFs (for details, see Methods), and then constructed the GRN. The networks of the nine immune cell subsets contain 8000 to 12,000 nodes (genes). However, they are not quite similar. Each of the immune cells has its own unique genes and genes that are common across the immune cell subsets. For instance, Thelper17 has around 3900 unique genes, while NK cell has only 39 unique genes in the GRNs. The characteristics of the key immune cell subsets' GRNs have been summarized in Table 1 and Figure 1a. Table 1 contains the genes annotated for each immune cell subset and the list of the genes that are unique to each immune cell subset. Furthermore, in order to find out the biological pathways that are associated with the uniquely regulated genes of the immune cells, we conducted enrichment analyses using Enchr on the unique genes of each immune cell's GRN. Table 2 listed the most enriched signaling pathways for each immune cell subset. For example, both the IL-2 receptor beta chain for T cell activation pathway and the PDGF signaling pathway are enriched in B cells. The Ras pathway, the ErbB1 downstream signaling pathway, and the mTOR signaling pathway are enriched in CD4 T cells. In the NK and regulatory T cells, the ATR (ataxia telangiectasia and Rad3-related protein) signaling pathway is enriched. These results suggest that the immune cell subsets have specific signaling pathways for cell-cell communications. annotating networks with various biological data such as gene expression profiles. The GRNs of each immune cell subset were analyzed using a shell script and then visualized in Gephi. Furthermore, Cytoscape was used to identify key regulators of the GRNs.
To elucidate the signaling pathways that are associated with uniquely regulated genes in the GRN of the T cells in PBMC (peripheral blood mononuclear cell) versus the TIC T cells, we conducted an analysis of their GRNs using Enrichr [32]. Enrichr is an extensive web server to perform gene set enrichment analyses with over 184 annotated gene sets from 102 gene set libraries [33].

Characterize Gene Regulatory Networks of Key Immune Cells Associated with Cancer Immunotherapy
The key immune cells associated with cancer immunotherapy have been reported. Among them, we selected nine key immune cell subsets in PBMC (Table 1), the chromatin profiles of which have been determined by either DNase-seq or ATAC-seq. To construct the GRN of an immune cell subset, we obtained the DNase-seq or ATAC-seq peak files of the immune cell subset to identify their chromatin open regions. Based on these regions, we scanned the DNA-binding motifs and their associated TFs (for details, see Methods), and then constructed the GRN. The networks of the nine immune cell subsets contain 8000 to 12,000 nodes (genes). However, they are not quite similar. Each of the immune cells has its own unique genes and genes that are common across the immune cell subsets. For instance, Thelper17 has around 3900 unique genes, while NK cell has only 39 unique genes in the GRNs. The characteristics of the key immune cell subsets' GRNs have been summarized in Table 1 and Figure 1a. Table 1 contains the genes annotated for each immune cell subset and the list of the genes that are unique to each immune cell subset. Furthermore, in order to find out the biological pathways that are associated with the uniquely regulated genes of the immune cells, we conducted enrichment analyses using Enchr on the unique genes of each immune cell's GRN. Table  2 listed the most enriched signaling pathways for each immune cell subset. For example, both the IL-2 receptor beta chain for T cell activation pathway and the PDGF signaling pathway are enriched in B cells. The Ras pathway, the ErbB1 downstream signaling pathway, and the mTOR signaling pathway are enriched in CD4 T cells. In the NK and regulatory T cells, the ATR (ataxia telangiectasia and Rad3-related protein) signaling pathway is enriched. These results suggest that the immune cell subsets have specific signaling pathways for cell-cell communications.
(a)    TFs play a prominent role in controlling the functionality of cells. GRNs encode the key global regulators and the hierarchy structure in the networks of the immune cells [34]. To find out the most important TFs for each immune cell subset, we analyzed the GRN of each immune cell subset and found the key regulators that control the networks. Table 1 illustrates the total number of TFs that can be found across each immune cell subset and the TFs that are unique to each of them. The distribution of TFs across the cell subsets is illustrated in Figure 1b. Further, we employed the network analysis algorithm to find the top 20 important TFs that influence the most downstream genes. Table 3 showed the top 20 TFs (i.e., based on their regulated gene numbers) in each immune cell subset. The TFs ELF2, KLF14, ZFX, ZN639, PLAG1, TBX1, KLF6, ZN148, KLF4, EGR4, TBX15, WT1, KLF16, THAP1, SP3, SP4, MAZ, SP2, TFDP1, EGR1, and SP1 are shared by immune cell subsets such as B, CD4, CD8, DC, NK, Regulatory T, Thelper1, and Thelper2 cells. However, Thelper17 has unique TFs such as PROX1, MNT, PURA, and PAX5. To substantiate the results obtained from above analysis, we used another computational procedure, the 'PageRank' algorithm [29]. Table 3 lists the top 20 TFs obtained from both methodologies. The results ( Table 3) showed that the TFs derived from the two methods are quite indistinguishable.
To examine the cell receptors of the immune cells, which, not only represent cell-cell communication, but also identify the GRNs of the immune cell subsets, we analyzed the receptor subnetwork for each immune cell subset. The results are shown in Table 1 and Figure 1c. The receptor subnetworks of the immune cell subsets are very different (Figure 2). For example, although DC and Thelper17 cells have 68 receptors, respectively, none of them are shared by the two immune cell subsets. On the other hand, Thelper2 and NK cells do not have any unique receptors. Interestingly, we found that the P2RX4 receptor is specific to CD4 cells, ESX1 is specific to regulatory T cells, and the IL6, ATXN10, and LAMB1 receptors are specific to CD8 cells. These specific receptors and their associated subnetworks could determine their cell identity or cellular state. Table 3. Key regulators in the gene regulatory networks of immune cells.

Cell Type Computational Method Employed Key_Regulators
B Network-analysis The first column depicts the various normal immune cell types used for our analysis, while the second column highlights the computational methods employed to scrutinize the most important TFs. Finally, the third column enlists the TFs that are found by the corresponding computational methods.
In summary, these results showed that most of the key regulators are similar in the immune cell GRNs, suggesting that the top regulators have been shared by many immune cell subsets. However, the regulated genes of the key TFs change a lot. Further, the signaling pathways that are associated with uniquely regulated genes for each immune cell subset are quite different, which means that the key regulators may control multiple different pathways. However, the specific cell receptors and the subnetworks of the immune cell subsets could operate in different immune cells to determine their identity.

Melanoma Cells Shut Down Many Network Activities of the CD8 T Cells
The CD8 T cell is one of the major immune cells that play a key role in killing cancer cells. However, it has been quite elusive how the cancer cell has evaded the immune system. To understand the underlying molecular mechanisms, we constructed and compared the GRNs for PD hi CD8 T tumor infiltrating cells (i.e., CD8 T cells with a high expression of PD-1 in melanoma) and three corresponding T cell subsets (i.e., CD8 naïve, CD8 Effector memory, and CD8 Central Memory T cells) in the PBMC of healthy people. Interestingly, we found that the GRN network size of the PD hi CD8 T cells changed compared to the GRNs of each PBMC's T cell subset. For example, the GRNs of the CD8 Central Memory and CD8 Effector memory have three million gene regulatory pairs, whereas the GRN of the PD hi CD8 T cells has 2.8 million gene regulatory pairs. These changes suggest that cancer cells could shut down specific T cell activities (i.e., regional subnetworks), which could be responsible for attacking cancer cells during tumorigenesis.
To carefully examine the differences in the GRNs of the PD hi CD8 T cells and the three corresponding T cell subsets of the healthy people, we focused on comparing their cell receptors and their associated subnetworks. To do so, we mapped the immune cell receptor list to each GRN and constructed cell receptor subnetworks (see Methods), and we found that specific receptors and their subnetworks are different between the PD hi CD8 T cell and its corresponding PBMC T cells.  Figure 3 showed the cell receptor subnetworks, which are differential receptor networks between the PD hi CD8 T cells and each of the corresponding PBMC T cells. Figure 3a illustrated the differential receptor regulatory network derived from the GRNs of the CD8 naïve cells and the PD hi CD8 T cells and showed that both receptors, P4HB and CD7, which are originally regulated in the 179 and 332 TFs, respectively, in the GRN of the CD8 naïve cells, disappeared in the GRN of the PD hi CD8 T cells. We found that the receptor CD7, which is activated and regulated by 493 and 488 TFs in the differential receptor regulatory networks derived from the GRNs of the CD8 effector memory cells and the PD hi CD8 T cells (Figure 3b) and from the GRNs of the CD8 central memory T cells and the PD hi CD8 T cells (Figure 3c), respectively, was inactivated in the GRN of the PD hi CD8 T cells. These results suggest that both P4HB and CD7 receptors could be associated with the process of tumor infiltration for reducing the immunological function of the T cells. Moreover, RNA-Seq data from Tirosh et al. and Pulko et al. clearly showed that CD7 and P4HB have higher expression levels in the T cells of the healthy people but lower expression levels in the tumor-infiltrating T cells ( Figure  4). This result validated the findings from our GRN analysis, mentioned above. Furthermore, an enrichment analysis of the differentially modulated genes between the T cells (Table 4) revealed that interleukin signaling pathways are highly enriched in the tumor-infiltrating T cells.  Figure 3 showed the cell receptor subnetworks, which are differential receptor networks between the PD hi CD8 T cells and each of the corresponding PBMC T cells. Figure 3a illustrated the differential receptor regulatory network derived from the GRNs of the CD8 naïve cells and the PD hi CD8 T cells and showed that both receptors, P4HB and CD7, which are originally regulated in the 179 and 332 TFs, respectively, in the GRN of the CD8 naïve cells, disappeared in the GRN of the PD hi CD8 T cells. We found that the receptor CD7, which is activated and regulated by 493 and 488 TFs in the differential receptor regulatory networks derived from the GRNs of the CD8 effector memory cells and the PD hi CD8 T cells (Figure 3b) and from the GRNs of the CD8 central memory T cells and the PD hi CD8 T cells (Figure 3c), respectively, was inactivated in the GRN of the PD hi CD8 T cells. These results suggest that both P4HB and CD7 receptors could be associated with the process of tumor infiltration for reducing the immunological function of the T cells. Moreover, RNA-Seq data from Tirosh et al. and Pulko et al. clearly showed that CD7 and P4HB have higher expression levels in the T cells of the healthy people but lower expression levels in the tumor-infiltrating T cells (Figure 4). This result validated the findings from our GRN analysis, mentioned above. Furthermore, an enrichment analysis of the differentially modulated genes between the T cells (Table 4) revealed that interleukin signaling pathways are highly enriched in the tumor-infiltrating T cells.

Discussion
Cancer immunologic therapies have been advanced in the past few years. Immune-checkpoint blockade (i.e., blocking PD-1, PD-L1, or CTLA-4) has shown durable clinical effects in some patients with various advanced cancers. Although amazing clinical responses have been observed with these therapies, the fact remains that only a relatively small subset of patients derives substantive clinical benefit from the therapy. There are major gaps in our knowledge of immunotherapy. One of the critical unanswered challenges is how immune cells become cancer-cell friendly and do not attack cancer cells. To uncover the underlying molecular mechanisms, we constructed and analyzed the GRNs of the key immune cell subsets associated with cancer immunologic therapies. We first analyzed the GRNs of the key PBMC's immune cell subsets, including B cell, CD4, CD8, CD8 naïve, CD8 Effector memory, CD8 Central Memory, regulatory T, Thelper1, Thelper2, Thelp17, and NK and DC cells to understand their activation profiles, regulatory mechanisms, and molecular pathways. It should be noted that this is the first study to systematical analyze the GRNs of immune cells.
We constructed GRNs using ATAC-seq and DNase-seq data. To check that the inferred regulatory relationships are functional, it is essential to examine the GRNs using gene expression profiles. Previously, we had used the ATAC-seq, DNase-seq, and gene expression data of the immune cells from the ENCODE database to check the parameters of the MEME tool for filtering TF-binding motifs. For example, by changing the parameters, we were looking for a GRN in which most of the both predicted TFs and their regulated genes are expressed (via checking the gene expression profiles of the same cell for that GRN) and most of the non-predicted TFs and their potentially regulated genes are not expressed (unpublished results). In this study, we applied the appropriate parameters (i.e., we used the p-value, 1 × 10 −4 , for filtering TF-binding sites) to construct the GRNs.
By comparative-network analysis of the GRNs of the key immune cells, we showed that most of the GRNs share key important hub regulators and that the differences between the immune cells are from the local regulatory subnetworks, which control the cell type-specific cell receptors. The cell receptors are the hallmark of each immune cell subset and allow us to pinpoint the immune cell subsets. In fact, these cell receptors could be used as cell markers to study immune cell subpopulations. The network structures and organizations of the immune cells where global regulators are shared in the GRNs and the subnetworks of the local, cell type-specific cell receptors have specificity, which suggests that the transformation of the immune cell subsets or cell types is relatively easily: the key hub regulators can be kept and the local regulatory networks can be changed between immune cell subsets. This nature of the GRNs of the immune cell subsets suggests that immune cell subsets have plasticity and could have a fast molecular mechanism for cell transformation between the immune cell subsets. This fast cell transformation mechanism allows rapid responses to distinct environmental signal cues.
To understand how cancer cells educate immune cells to make them more cancer-cell friendly, we constructed the GRNs of the TICs and compared their corresponding immune cells in PBMC in the same melanoma patient. In human melanoma, we found that the GRN of the tumor-infiltrating T cells has reduced its network size compared to those of the corresponding T cell subsets in PBMC. These results clearly showed that cancer cells send signals to shut down many cellular activities of the TIC T cells. Based on this result, we highlighted that one of the key underlying molecular mechanism for making T cells more cancer-cell friendly is shutting down a lot of a T cell's cellular activities, which could originally allow T cells to recognize and attack cancer cells during tumor infiltrating T cells. Carefully analyzing of the shutdown-subnetwork of the TIC's T cells suggests that CD7 and P4HB receptor subnetworks have been turned off by cancer cells.
LGALS1 binds to glycoproteins on the T cell surface receptor CD7. It has been shown that upregulated CD7 in T cells can trigger a pro-apoptotic signal during LGALS1-induced T cell death. P4HB is a hormone binding protein that can bind to Galectin-9 [36]. Both Galetin-9 and LGALS1 can trigger T cell death. It has been shown that the exhausted T cells express high levels of PD-1 by losing CD7 and P4HB expression. The will decrease the level to trigger T cell death, which may make the T cells get into a proliferative defect, exhausted state. Tumor-infiltrating T cells have two fractions: PD1 hi (T cells with high-expressed PD1) and PD1 lo T cells (T cells with low-expressed PD1) (see Ref 23). An enrichment analysis of them using RNA-seq data (Table 4) showed that both the glucocorticoid receptor regulatory network and the ATF-2 transcription factor network are enriched in the PD1 hi T cells. The pathways of IL12-mediated signaling, IL2 signaling events mediated by STAT5, and calcineurin-regulated NFAT-dependent transcription are enriched in most of the functional CD8 T cells (functional CD8 vs exhausted CD8 T cells). Philip found that a cell subset of PD1 hi T cells had higher levels of CD38 and CD101 and lower levels of CD5, which could be cell markers for the exhausted T cells [23]. Pauken et al. found that PD1 blockade can reprogram the epigenetic landscapes of exhausted T cells to functionally effect T cells at some level [37]. These results could suggest that losing CD7 and P4HB expression in the PD1 hi T cells could induce the exhausted state of the T cells.
In addition, CD7 is specifically expressed in T and NK cells. Its expression could activate T cells. Further, studies [38] have shown that a ligand for CD7, SECTM1 (secreted and transmembrane protein 1), is highly expressed in many tumors, including melanoma cells. The interaction of the SECTM1 and CD7 significantly increases monocyte migration by activation of the PI3K pathway. Thus, we suspect that tumor cells could send a certain signal to turn off or down-regulate CD7 expression in the TIC T cells so that it could reduce the migration of the T cells into the tumors. We envisioned that if we could genetically engineer the T cells so that we can constitutively turn on the subnetworks that have been shut down by cancer cells, we could let T cells recognize and attack cancer cells and pave a new way to treat tumors using T cell genetic engineering.