The co-expression networks of differentially expressed RBPs with TFs and LncRNAs related to clinical TNM stages of cancers

Differential expression analysis

We downloaded the original data including count matrix, expression FPKM values and clinical information of 20 cancer types and paired normal tissues from TCGA data base using SangerBox tool (SangerBox, http://sangerbox.com/). Then, we applied edgeR (Robinson, Mccarthy & Smyth, 2010) and DESeq2 (Anders, 2009) to select differentially expressed genes (DEGs) from count matrix for each cancer type with parameter padj < 0.05, and |log2FC| > 1. Thirdly, according 1,542 RBP genes information (Stefanie, Markus & Thomas, 2014), we divided the DEGs into two types for each cancer types: RBPs and non-RBPs. The expression FPKM values were used to compute the degree of dysregulation of DEGs and construct the co-expression networks. Clinical TNM information for 5,093 patients were used to co-expressed gene module detection. Our datasets contained 1,542 RBPs, 1,570 TFs and 4,147 lncRNAs; data acquisition information was presented in the Data Availability section. The ‘pheatmap’ package in R was used to generate a heatmap of differentially expressed RBPs shared by 16 cancer types in six systems based on log2-fold change values. The Pearson coefficient (R) was used to judge the similarity of RBPs genes in different cancers occurring at related tissues have similar expression patterns. The gene cluster with R > 0.5 was considered as having the similar expression pattern. The gene ontology categorization analysis tool DAVID (Huang, Sherman & Lempicki, 2009a; Huang, Sherman & Lempicki, 2009b) was used to determine biological processes and molecular functions of DE RBPs.

Standard deviation and proportion of dysregulation

In order to investigate the degree of dysregulation of RBPs, we computed the standard deviation of differentially expressed genes listed in Table 1 using the ‘std’ function in ‘MATLAB’ version R2015b, respectively. In detail, we firstly computed the standard deviation of differentially expressed RBPs in 20 cancers and paired normal tissues (Fig. 2A). Similarly, we next computed the standard deviation of differentially expressed non-RBPs in 20 cancers and paired normal tissues (Fig. 2B).

The definition of the standard deviation is as follows: ${\displaystyle \sigma =\frac{1}{N}\sqrt{\sum _{i=1}^N{\left(x_i-\mu \right)}^2}.}$ Here, N is the number of genes, x_i isthe expression FPKM value of the ith gene, µ is the mean FPKM value of all N genes. Larger σ values represent gene expression values that deviate from the mean value to a greater degree, which is indicative of greater gene dysregulation. The standard deviation values for RBP expression data is presented in Table S1.

In order to more clearly see the expression differences of RBPs and non-RBPs in cancers and normal tissues, we further computed their relative deviations. In detail, we first computed the relative deviation of RBPs in cancers and paired normal tissues, and then we computed the relative deviation of non-RBPs in cancers and normal tissues (Fig. 2C). The relative proportion of deviation is computed using the following function $:p=\frac{{\sigma }_c}{{\sigma }_n}$. Here, σ_c and σ_n are the standard deviations of the gene expressions in cancers and normal tissues respectively.

Similarly, in order to investigate the expression differences of up- and down-regulated RBPs, we also computed the standard deviation of them in cancers and normal tissues, respectively (Figs. 2D, 2E). Finally, the mean expression values of up- and down-regulated RBPs in 20 cancers and paired normal tissues were computed and presented in Figs. 2F & 2G.

Weighted co-expression network analysis

Weighted gene co-expression network analysis (WGCNA) (Peter & Steve, 2008) is a comprehensive R package that summarizes and standardizes methods and functions for co-expression network analysis. Module detection function of WGCNA was used to detect the correlations between co-expression gene modules and the clinical TNM information for 13 types of cancers with default settings one by one. Then, the threshold of correlation coefficient R > 0.5 and statistical significance P < 0.05 was used to select the cancer type whose TNM information related to co-expressed gene module. Genes with same expression pattern were clustered into one module and marked with same color. And then, the network construction module was used to construct co-expression networks for RBPs with TFs and lncRNAs in the CHOL and KICH cancer types, respectively. The detailed steps of constructing networks are as follows. Firstly, WGCAN network construction tool was used to generate the nodes and edges of genes by computing correlations of expression values. The nodes corresponded to genes, and the edges were determined by the pairwise correlations between the expression levels of genes. The corresponding called function in the R package was ‘blockwiseModules’ and the parameters were set as follows: ‘powers = 10, minModuleSize = 30, mergeCutHeight = 0.25’, other parameters we set to the default setting. Secondly, nodes with correlation r < 0.5 and edges with weighted threshold <0.3 were removed. Finally, the Cytoscape (https://cytoscape.org/) tool was used to plot the interactions using the nodes and edges of conserved genes.

Clinical TNM information processing

Clinical TNM information for 5,093 patients exhibiting 13 types of cancers downloaded from the TCGA database presented in Table S4. Generally, in the TNM system, ‘T’ refers to a primary tumor, ‘T1 ∼ T4’ represents the severity of primary cancer according to the increase in tumor volume and the extent of involvement of adjacent tissues, and ‘T0’ indicates no primary tumor. ‘N’ represents the tumor spreading to regional lymph nodes. ‘N1 ∼ N3’ represents the degree of spreading according to the extent of lymph node involvement. ‘M’ refers to tumor metastasis. No distant metastasis is expressed by ‘M0’, and distant metastasis is expressed by ‘M1’. To investigate the correlation between TNM information and gene expression values using WGCNA, we converted the TNM information into a weighted matrix. For example, if TNM information for a patient was ‘T1-N3-M1’, the corresponding weighted array is [1 3 1].

Statistical methods

We used the Mann–Whitney U-test (function ‘ranksum’ in software‘MATLAB’ version R2015b) to examine whether there is statistical significance between given two samples, the default significance level is 0.05 (Lian et al., 2018).

The deviations of expression levels of RBPs between cancerous and normal tissues

To investigate the degrees of dysregulation of RBPs in many human cancers, we computed the standard deviations of gene expression levels of RBPs and non-RBPs in 20 types of cancerous and normal (control) tissues, respectively. These results indicate that, relative to normal tissue, RBPs show greater variation in expression in almost all types of cancer (P < 0.05, Mann–Whitney U test, Fig. 2A). In contrast, non-RBPs show considerably less variation in expression in cancers relative to normal tissues (P < 0.05, Mann–Whitney U test, Fig. 2B). These results indicate that RBPs show a greater degree of dysregulated expression than non-RBPs in almost 20 types of cancers. Furthermore, the degrees of dysregulation of RBPs in different cancer types are significantly different. In particular, RBPs show severely dysregulated expression in eight types of cancers, including BLCA, LUAD, STAD, READ, LUSC, GBM and COAD. For these cancers, the corresponding relative proportions of dysregulation relative to normal tissues are 365.5%, 415.8%, 446.4%, 456.3%, 494.5%, 540.9% and 759.1%, respectively (Fig. 2C, Table S1). Interestingly, up-regulated RBPs show a smaller standard deviation in expression in cancer tissue than in normal tissue (P < 0.05, Mann–Whitney U test, Fig. 2D); However, down-regulated RBPs show considerable differences in the standard deviation of their expression in cancer tissue relative to normal tissues, especially in KICH, KIRC, KIRP, COAD, LUSC, CESC, UCEC, and GBM (P < 0.05, Mann–Whitney U test, Fig. 2E). This may indicate that down-regulated RBPs are more dysregulated in cancers than are up-regulated RPBs. Finally, biological process enrichment analyzing indicates both up- and down-regulated RBPs were highly enriched in the processes, such as rRNA metabolism, nuclear-transcribed mRNA catabolism, and ncRNA processing (Fig. S1). In addition, up-regulated RBPs were also enriched in mitochondrial gene expression and in the regulation of mRNA metabolic processes, while down-regulated RBPs were enriched in ribosome biogenesis and rRNA processing.

Up- and down-regulated RBPs show opposite expression patterns in cancer and normal tissue

To investigate cancer-specific differences in RBP expression, we analyzed the standard deviations and mean expression values of up-regulated and down-regulated RBPs in 20 types of cancers and in normal tissues (see methods and materials). Our results for all 20 cancer types suggests that, compared to normal tissues, up- and down-regulated RBPs show opposite patterns of expression in almost all cancers; what’s more, down-regulated RBPs tend to show the larger expression deviations in cancers than up-regulated RBPs (P < 0.05, Mann–Whitney U test, Figs. 2D, 2D). In particular, in almost all types of cancers, down-regulated RBPs show larger expression values than up-regulated RBPs. Furthermore, the expression deviations of up-regulated RBPs in cancers are lower than in normal tissues. However, the expression deviations of down-regulated RBPs are considerably greater in cancer tissue than in normal tissue; this is especially true for BLCA, GBM, HNSC, LUAD, LUSC, STAD, BRCA, KICH, READ, and UCSC (P < 0.05, Mann–Whitney U test, Figs. 2E, 2G). These results suggest that down-regulated RBPs show a severer dysregulation in cancers than up-regulated RBPs. Furthermore, the expression pattern of up- and down-regulated RBPs in 20 types of cancers are the opposite of those found in normal tissues; what’s more, molecular functions enrichment analyzing indicates that both up- and down-regulated RBPs showed enrichment in: catalytic activity acting on RNA, mRNA and mRNA 3′-UTR binding, nuclease and ribonuclease activity, and single-stranded RNA binding. In addition, we found enrichment in translation factor activity and RNA binding for up-regulated RBPs, and in catalytic activity acting on tRNAs for down-regulated RBPs (Fig. S1). In addition, we also were able to identify which specific biological processes and functions were regulated by up- and down-regulated RBPs, which may reveal how the normal processes of cells can be altered in a way that leads to cell carcinomatosis. These results provide a new insight into understanding the roles of up- and down-regulated RBPs in the process of cell carcinogenesis.

RBPs of different cancers in same system have a similar expression profile

To gain clearer insight into RBP expression in cancerous tissues, we analyzed expression heatmaps of 16 types of cancers, except four types of glandular cancer (PAAD, THCA, PRAD, and HNSC). Because the number of DE RBPs in these four types of cancers are too small (Table 1). We divided DE RBPs into six scale systems according to type (i.e., organization of canceration) and analyzed the expression profiles of 801 DE RPBs shared by all 16 cancers by system. In each system, we divided RBPs into different types according to their expression values and Pearson coefficient R. Gene clusters with R > 0.5 were considered as having the similar expression pattern, which was shown in heatmaps. This was true for all cancers except two cancers of the nervous system. The corresponding gene lists of RBPs relevant for each system are presented in Table S2.

In the respiratory and liver and gall systems (Figs. 3A, 3B), two pairs of cancers in same systems, LUAD and LUSC, CHOL and LIHC, show very similar RBP expression profiles. Among these two cancer systems, C1 and C4 genes showed the co-expression patterns of RBP for high expression and low expression respectively (R > 0.5). The proportions of co-expressed RBPs in the respiratory and liver and gall systems were 88.8% and 76.8%, respectively. However, RBP expression in two cancers of the nervous system (PCPG and GBM) showed opposite expression patterns (Fig. 3C). The proportions of highly expressed RBPs in PCPG and GBM cancer tissue were 60.5% and 36.1%, respectively. Highly expressed RBPs in PCPG showed a lower degree of expression in GBM, while highly expressed RBPs in GBM showed lower expression in PCPG. The proportion of RBPs showing opposing patterns of expression in PCPG and GBM was 68.2% (R < 0.5). In systems containing three types of cancers (i.e., the alimentary and reproductive systems) (Figs. 3D, 3E), the proportion of co-expressed RBPs decreased slightly, reaching 70.4% and 71.5% (R > 0.5), respectively. In the urinary system, which was affected by four types of cancers, the proportion of co-expressed RBPs reached its minimum value of 35.3% (Fig. 3F). Of the cancers of the urinary system, two (KIRC and KIRP) showed the strongest degree of RBP co-expression reaching 85.6% (R >0.5). Taken together, our results suggest that RBP expression in different cancers in similar tissues have similar expression profiles.

Co-expressed gene regulatory networks correlate with clinical TNM stage

Next, we investigated whether expression patterns are closely related to developmental stages of different cancers and constructed the co-expression networks for DE RBPs and non-RBPs using module detection and network construction tools of WGCNA (Peter & Steve, 2008). Results showed that for two types of cancers—CHOL and KICH—clinical TNM stage information was closely related to patterns of gene expression (Fig. 4). The results for the other 11 cancer types did not satisfy the threshold (Methods and Material, Figs. S2–S7).

In terms of module detection, we identified three modules (orange, green and red) related to the cancer metastasis stage (M-stage), as well as two modules (royal blue and red modules) closely related to the regional lymph node stage (N-stage) for CHOL (Fig. 4A). The red module was consistently in both the M and N stages. The gene list for each corresponding module is presented in Table S3. Furthermore, we found 10 differentially expressed RBPs in these three modules. These included six RBPs (ACO1, PPARGC1A, PUS7, KHDC1, ELAVL3, and BICC1) in the green module, two RBPs (PANBP17 and HNRNPA1) in the royal blue module, and two RBPs (DCPS and C2orf15) in the red module. For KICH, we found that the royal blue module was most closely related to the M-stage, with a correlation coefficient of 0.96 (P < 4e−15), while the red and green modules were closely related to the N-stage (Fig. 4B). The gene list for each corresponding module is presented in Table S3. Furthermore, we found seven DE RBPs related to different developmental stages of cancer. These included four RBPs (TNRC6A, MECP2, ZCCHC14, and POLR2F) in the green module and three RBPs (TDRD1, TDRD9, and CELF4) in the red module. The regulatory networks of CHOL (Fig. 4C) revealed that (1) in each sub-network, one RBP interacts with almost all non-RBPs, suggesting that these RBPs are key regulators for each module; (2) the RBPs in different sub-networks interact with each other, indicating that they work together to regulate the corresponding developmental stage of the cancer. (3) we found that most key RBPs are down-regulated, the proportion is 80%. In addition, in the red module, two RBPs—including one down-regulated RBP (C2orf15) and one up-regulated RBP (CDPS) also play such a regulatory role. For the regulatory networks of KICH, we also identified three sub-networks that corresponded to the M- and N-stages (Fig. 4D). First, we found that RBPs interacted with almost all non-RBPs in each module, indicating that RBPs are key regulatory factors of the genes in these modules. Second, we identified seven key RBPs in the three modules, of which five were down-regulated and two were up-regulated. In the green module, we identified four key RBPs—TNRC6A, MECP2, ZCCHC14, and POLR2F—all of which were down-regulated. In the red module, we identified two key up-regulated RBPs (TDRD1 and TDRD9) and one key down-regulated RBP, CELF4. These results suggest that dysregulated RBPs play a key role in the regulation of the development of the CHOL and KICH M-stage, which may provide a new perspective for potential prognostic biomarkers and therapeutic targets for patients with cancers at M stages in two cancer types CHOL and KICH.

The networks of key RBPs and lncRNAs for CHOL and KICH

To infer the potential regulatory mechanisms of lncRNAs with key RBPs related to M stages, we constructed the co-expression networks of key RBPs and differentially expressed (DE) lncRNAs for CHOL and KICH respectively and performed Gene Ontology and functional enrichment analyses.

Ten key RBPs and 2,943 DE lncRNAs were used to construct the co-expression network for CHOL. The resulting co-expression network consisted of two key RBPs (down-regulation C2orf15 and up-regulation DCPS) and 75 lncRNAs, which were grouped into three clusters (Fig. 5A). There were 63 up-regulated lncRNAs, and the proportion of up-regulated lncRNAs in three clusters was 86%, 75%, 83.3%, which probably suggest that up-regulated lncRNAs have a greater interaction with RBPs in the process of metastasis of CHOL cells. Functional enrichment analyzing demonstrated that cluster 1 and cluster 3 had similar functions and mainly enriched in functional categories involved in gene silencing and negative regulation of translation, such as post transcriptional gene silencing, negative regulation of translation and cellular amide metabolic process, cellular response to dsRNA, miRNA metabolic process. Cluster 2 had some special functions and mainly enriched in such as positive regulation of mRNA catabolic process, cellular response to interleukin-1 and calcium ion, RNA destabilization.

Seven key RBPs and 1,204 DE lncRNAs were used to construct the co-expression network for KICH. The resulting co-expression network consisting of 3 key RBPs (up-regulated TDRD1 and TDRD9, down-regulated CELF4) and 227 lncRNAs. There are 177 down-regulated lncRNAs in the network and were grouped into four clusters (Fig. 5B). The corresponding proportion of down-regulated lncRNAs was 46.7%, 65.3%, 75% and 95%, which probably suggest that down-regulated lncRNAs play more important roles in interacting with key RBPs TDRD1, TDRD9 and CELF4 in the process of metastasis of KICH cells. Functional enrichment analysis demonstrated that cluster 1 shared by three key RBPs mainly enriched in dsRNA fragmentation and production of miRNAs involved in gene silencing by miRNA. Cluster 2 regulated by TDRD9 mainly enriched in endoribonuclease and exon-exon junction complex, cluster 3 regulated by TDRD1 mainly enriched in mRNA catabolic process and regulation of mRNA metabolic process, cluster 4 regulated by CELF4 mainly enriched in transporting of RNA, mRNA and nucleic acids. Besides, cluster 1 and cluster 2 have some similar functions, such as telomere maintenance, histone mRNA and miRNA metabolic process, dosage compensation. Cluster 3 and cluster 4 have some similar functions, such as regulation of RNA stability, RNA localization and regulation of mRNA catabolic process.

These results provide a new insight into the understanding of the interactions of key RBPs with lncRNAs in the metastasis stage (M stage) of cancer cells.

The co-expression networks of DEG RBPs and TFs for CHOL and KICH

To investigate the interactions of RBPs and TFs, we constructed the co-expression networks of DEG RBPs and TFs for CHOL and KICH, respectively. The key regulatory RBPs were those (ten for CHOL and seven for KICH) detected in above section.

The co-expression network of CHOL revealed two important insights. First, we found the five largest transcription factor families, they are C2H2-ZF, Homeodomain, Nuclear receptor, bHLH and bZIP, and the corresponding proportion is 37%, 17%, 9%, 9% and 9% (Figs. 6A, 6B), which interacted with almost all differentially expressed RBPs. This result indicates that these transcription factors tend to show a co-expression pattern with DEG RBPs, which further suggest that they play a major regulatory role in RBP post-regulatory levels for CHOL. Second,we also identified several special TFs related to up- or down-regulated RBPs for CHOL. For example, Grainyhead, MADF, HMG-Sox and SAND are specific transcription factors associated with down-regulated RBPs and GTF2I-like, Myb-SANT, MADS box and CENPB are specific transcription factors associated with up-regulated RBPs for CHOL. The co-expression network for KICH also revealed the following insights. First, we found that the four largest transcription factor families, they are bHLH, bZIP, C2H2-ZF and Nuclear receptor, the proportion is 54%, 9%, 7% and 7% (Figs. 6C, 6D). Notably, transcription factor bHLH interact with all RBPs and it accounts for more than half of all interacted transcription factors, which indicate that bHLH transcription factor probably involved in regulation of all differentially expressed RBPs for KICH. Second, the proportion of up- and down-regulated RBPs co-expressed with TFs is 57% and 43% respectively. But, the up-regulated RBPs tend to show more interactions with TFs (Fig. 6C). These results provide insights into understanding the mechanism of interaction between transcription factors and RBPs.