Integrated Analysis of MicroRNA (miRNA) and mRNA Profiles Reveals Reduced Correlation between MicroRNA and Target Gene in Cancer

Background Accumulating evidences demonstrated that microRNA-target gene pairs were closely related to tumorigenesis and development. However, the correlation between miRNA and target gene was insufficiently understood, especially its changes between tumor and normal tissues. Objectives The aim of this study was to evaluate the changes of correlation of miRNAs-target pairs between normal and tumor. Materials and Methods 5680 mRNA and 5740 miRNA expression profiles of 11 major human cancers were downloaded from the Cancer Genome Atlas (TCGA). The 11 cancer types were bladder urothelial carcinoma, breast invasive carcinoma, head and neck squamous cell carcinoma, kidney chromophobe, kidney renal clear cell carcinoma, kidney renal papillary cell carcinoma, liver hepatocellular carcinoma, lung adenocarcinoma, lung squamous cell carcinoma, stomach adenocarcinoma, and thyroid carcinoma. For each cancer type, we firstly obtained differentially expressed miRNAs (DEMs) and genes (DEGs) in tumor and then acquired critical miRNA-target gene pairs by combining DEMs, DEGs and two experimentally validated miRNA-target interaction databases, miRTarBase and miRecords. We collected samples with both miRNA and mRNA expression values and performed a correlation analysis by Pearson method for miRNA-target pairs in normal and tumor, respectively. Results We totally got 4743 critical miRNA-target pairs across 11 cancer types, and 4572 of them showed weaker correlation in tumor than in normal. The average correlation coefficients of miRNA-target pairs were different greatly between normal (-0.38 ~ -0.61) and tumor (-0.04 ~ -0.26) for 11 cancer type. The pan-cancer network, which consisted of 108 edges connecting 35 miRNAs and 89 target genes, showed the interactions of pairs appeared in multicancers. Conclusions This comprehensive analysis revealed that correlation between miRNAs and target genes was greatly reduced in tumor and these critical pairs we got were involved in cellular adhesion, proliferation, and migration. Our research could provide opportunities for investigating cancer molecular regulatory mechanism and seeking therapeutic targets.


Introduction
Cancer ranks among the most frequent causes of human death worldwide, and its incidence and mortality have risen rapidly recent years [1]. As a heterogeneous disease, cancer involves uncontrolled growth of malignant cells caused by alters of genome, methylation, mRNA, and miRNA expression [2][3][4]. Despite several advances being improved in diagnosis and therapy of tumor, approximate 5-year survival rate is still low [5]. Therefore, it is urgent for exploration of underlying molecules mechanisms, which would help to improve the early detection and therapy for cancer.
MicroRNAs (miRNAs) are defined as small regulatory RNA molecules with 19-24 nucleotides in length [6]. They repress gene expression by binding to complementary sequences in the 3' untranslated region (3' UTR) of mRNAs to target them for degradation and thereby prevent their translation [7]. Recent studies have demonstrated that 2 BioMed Research International miRNAs play an essential role in carcinogenesis, including cell migration and reproduction, by base pairing with target mRNAs of protein-coding genes [8]. Hence, screening critical miRNA-target gene pairs for cancer is a crucial work in tumor study. With the demands of miRNA-target research, the novel experimental technology like computational prediction, ChIP-seq, Microarray, and RNA sequencing had been widely used in miRNA-target identification [9], and the database for collecting experimentally validated miRNAtarget pairs had been established, like miRTarBase [10] and miRecords [11]. The wide application of coexpression analysis in miRNA-target identification began with the development of Microarray and RNA sequencing [12]. With the expression values of miRNA and mRNA in multisamples, researchers could calculate correlation coefficients of miRNA-target pairs by Kendall, Spearman or Pearson method [13]. As miRNA's negatively regulation on gene expression often resulted in an inverse correlation between miRNA and its target [14], negative correlation coefficient was proposed to be an effective way of miRNA-target pair identification [15].
With the application of correlation coefficients in miRNA-target identification, miRNA-target gene pairs have consistently been suggested to play a role in tumorigenesis and progression. For example, miR-182 could result in a sustained activation of HIF1 pathway via targeting PHD and FIH in prostate cancer, which might facilitate tumor cell adaption to hypoxic stress during prostate tumor progression [16]; miR-23b, miR-365-1, and miR-365-2 were identified as biomarkers for colorectal cancer detection by exploring colorectal cancer related miRNA-target gene pairs [17]; overexpression of miR-204 reduced the gastric cancer cell proliferation by targeting CKS B, CXCL , and GPRC A [18]. However, restricted by the number of samples, studies of correlation between miRNA and target genes in different clinical groups were still absent [19,20]. Without a clearly understanding of correlation coefficients of miRNA-target pairs, the application of it would be limited [21]. The popularization and application of low-cost RNA sequencing, including mRNA sequencing and miRNA sequencing, gave us opportunities to measure the changes in correlation coefficients of miRNA-target pairs between different clinical groups [22].
The Cancer Genome Atlas (TCGA, https://cancergenome .nih.gov/) project data portal provides a platform of RNA sequencing with mRNA and miRNA for 11 000 samples across 33 different human cancer types. It contains both tumor tissue samples and adjacent normal tissue samples. Here, using the comprehensive molecular data from TCGA, we performed miRNA-target coexpression analysis and characterized the correlation changes in 11 cancer types. The workflow in this study was shown in Figure 1 and the cancer types were chosen by considering their sample sizes of healthy controls and different staging tumors, which were bladder urothelial carcinoma (BLCA), breast invasive carcinoma (BRCA), head and neck squamous cell carcinoma (HNSC), kidney chromophobe (KICH), kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), liver hepatocellular carcinoma (LIHC), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), stomach adenocarcinoma (STAD), and thyroid carcinoma (THCA). Our study could provide opportunities for investigating cancer molecular regulatory mechanism. The miRNA-target gene pairs would be useful as a point in further research for finding novel diagnostic biomarkers and therapeutic targets of cancer.

Materials and Methods
. . Data Resources and Preprocessing. Raw read count data of 5680 RNA sequencing samples and 5740 miRNA sequencing samples and corresponding clinical annotation files for human cancer types were downloaded from TCGA by the officially recommended GDC data transfer tool. Samples with both miRNA and mRNA expression profiles data were obtained and, further, patients who donated their tumor tissues and adjacent nontumor tissues were collected as paired samples. As a sufficient number of paired samples was important to verify the correlation analysis results, we have finally chosen 11 cancer types in this study, which were BLCA, BRCA, HNSC, KICH, KIRC, KIRP, LIHC, LUAD, LUSC, STAD, and THCA. Samples were classified into 6 groups, normal tissues, tumor tissues from stage I to IV, and tumors without stage information, according to their pathologic stages on clinical annotation files [23]. We defined normal tissues as healthy controls, stage I and stage II as early stage (Early), and stage III and stage IV as advanced stage (Advanced), whereas tumors without stage information were as TNS. The detailed sample sizes for every cancer type were presented in Table S1 and clinical information, containing sample ID, patient ID, sample stages, and groups information, was given in Table S2.
For each type of cancer, only protein coding RNA and miRNA expression data were preserved, count data were normalized by DESeq2 [24] and NormExpression were used to evaluated the normalized data quality [25]. Since the data was obtained from TCGA, further approval by an ethics committee was not required. The experimentally validated miRNA-target interactions were collected from two databases, miRTarBase (2017, Version7.0) [10] and miRecords (April 27, 2013) [11]. We integrated these two databases and merged duplicated pairs for further analysis.
. . Differential Analysis. For each cancer type, we compared expression values of miRNAs and mRNAs in early stage and advanced stage tumors as well as healthy controls by DESeq2 [24]. The thresholds were set as |log2FC| > 1 and padj < 0.05 to select differentially expressed miRNAs or mRNAs. miRNAs and mRNAs with log2FC > 1 were defined as upregulated, whereas those with log2FC < −1 were defined as downregulated. We removed downregulated miRNAs or mRNAs whose average expression value in healthy controls was less than 5 and upregulated miRNAs or genes whose average expression value in early or advanced was less than 5. After filtering, the low abundant miRNAs and genes were filtered out and the peaks of log values density of miRNAs and mRNAs were increased ( Figure S1). Finally, we selected miRNAs and mRNAs that showed same status of regulation (upregulated or downregulated) in both early and advanced stage tumors for further analysis. . . Critical miRNA-Target Gene Pairs Screening and Correlation Analysis. To identify critical miRNA-target gene pairs affecting tumorigenesis, we combined DEGs, DEMs and two experimentally validated miRNA-target gene databases, miRTarBase and miRecords. miRNA-target gene pairs were collected from experimentally validated interactions when both miRNA and gene were differentially expressed in tumor tissues. Samples with both miRNAs and mRNAs expression data were obtained and divided into healthy controls, early stage and advanced stage tumors. With the normalized and log2 transformed expression values, Pearson correlation coefficients of miRNA-target gene pairs were calculated for healthy controls. The pairs with Pearson correlation coefficients < -0.2 and corresponding p-value < 0.05 in the healthy controls were considered as tissues-specific pairs. We performed the correlation analysis for these pairs in early stage and advanced stage separately and compared the results with that in healthy controls.

. . Functional Annotation and Interaction Visualization.
To gain insight into the functions of miRNA target genes, we performed Gene Ontology (GO) classification and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis based on the online tools of Database for Annotation, Visualization and Integrated Discovery (DAVID) [26,27]. Three categories were included in GO such as biological process involves pathways and cellular processes, molecular function refers to gene products activities, and cellular component to the location where gene products are active. Biological process, molecular function and cellular component. GO terms and KEGG pathways were selected using P < 0.05 as the cut-off. To present the regulation between miRNA and gene, critical miRNA-target interactions were used to construct miRNA-gene network by Cytoscape software [28]. Figure Plotting. Statistics and calculations in this study were conducted using the software R 3.3.0 (https://www.r-project.org/) with the Bioconductor packages [29]. | log 2FC| > 1 and padj < 0.05 was considered as a statistically significant difference between values. "ggplot2" were used to visualized the figures [30].

. . Differentially Expressed miRNAs and Genes in Tumor
Tissues. We totally identified 62 ∼ 176 DEMs and 1 483 ∼ 4 875 DEGs for the 11 cancer types ( Figure 2). As shown in this figure, the numbers of upregulated and downregulated miRNAs varied widely in 11 cancer types. For example, 143 miRNAs were upregulated and 18 miRNAs were downregulated in LIHC; 147 miRNAs were upregulated and 65 miRNAs were downregulated in BLCA; and 129 miRNAs were upregulated and 60 miRNAs were downregulated in LUSC ( Figure 2(a)). Similar phenomenon was observed for the upregulated and downregulated genes in different cancers. For example, 2 501 genes were upregulated and 1 002 genes were downregulated in LIHC, whereas, 2 754 genes were upregulated and 1 717 genes were downregulated in KIRC (Figure 2(b). These results suggested global expression differences between normal and tumor tissues. The DEMs and DEGs along with their detailed statistic items for each type of cancer were given in Tables S3 and S4 respectively.
To find out multicancer related DEGs and DEMs, we summarized DEGs and DEMs from 11 cancer types (Table  S5). Totally, 12 037 DEGs and 496 DEMs were found, and 51.6% of DEMs and 77.8% of DEGs were differentially expressed in multicancers. Table 1 listed mostly frequent DEMs and DEGs across 11 types of cancer, including genes like ADGRB , ASF B, CDKN A, CDT , DPT, E F , E F , MYBL , TCEAL and miRNAs like miR-183, miR-96, miR-1, miR-1258, miR-133a, miR-135b, miR-144, miR-182, miR-195, miR-21. These universal genes and miRNAs were mostly associated with enhanced proliferation, migration or invasion of various cancers. For example, MYBL which was targeted by miR-143-3p could regulate breast cancer cell proliferation and apoptosis [31]; CDT was a critical regulator for normal genome replication [32], knockdown of TGFBR in T24 cells could result in decreased cell growth, motility and invasion [33], and miRNAs like miR-182, miR-183, miR-21, miR-195 and miR-96 were all associated with tumorigenesis [16,[34][35][36][37]. Additionally, we identified 4 737 DEGs and 143 DEMs, which showed discordantly regulated status in different cancer types. For example, MYBL was downregulated in THCA but upregulated in the remaining 10 cancer types; FXYD was upregulated in KIRC but downregulated in other 10 cancer types. These results suggested the great commonness and heterogeneities among different types of tumors. The commonness could help us in finding biomarkers for multicancer, and the heterogeneities could give us guidance in individualized treatment.
. . Critical miRNA-Target Gene Pairs and Correlation Coefficients. We totally found 53 ∼ 1 151 critical miRNA-target gene pairs for each cancer type. The critical pairs followed with their Pearson correlation coefficients and corresponding p values in healthy controls, early and advanced stage tumor samples were given in Table S6. The distribution of Pearson correlation coefficients of these critical pairs in healthy controls, early and advanced stage tumors were shown in Figure 3. Obviously, the medians of Pearson correlation coefficient values of these critical pairs were higher in healthy controls than that of in early and advanced stage tumors for all 11 cancers studied. Additionally, we performed a correlation analysis for these critical pairs in samples stage I to IV separately and the results ( Figure S2) still exhibited the higher Pearson correlation coefficient values in healthy controls. Additionally, we divided miRNA-target gene pairs into two groups according their up-or downregulation gene, while seldom differences were observed.
The discrepancy in the above analysis was that the sample sizes of early and advanced stage tumors were different from those of healthy controls. Generally, the sample sizes of early or advanced were larger than those of normal, which might cause the biased results. To eliminate the impact of sample size, we found 476 patients, who not only donated their tumor tissues but also donated adjacent normal tissues, from 5650 samples. After reanalyzing, we got the correlation coefficients of miRNA-target genes which were still much higher in normal (the average Pearson correlation coefficient of -0.38 ∼ -0.61) than those of in tumor (the average Pearson correlation coefficient of -0.04 ∼ -0.26). The distribution of Pearson correlation coefficients of these critical pairs in normal and tumor tissues was shown in Figure 4. The above results proved that the correlation coefficient of miRNA-target pair had a difference between normal and tumor samples, and most of critical pairs showed much weaker correlation in tumor than in normal.
In order to find out universal miRNA-target gene pairs that might be associated with various cancer types, we summarized the pairs in all 11 cancer types and counted the frequency of each critical pair in 11 cancer type (Table S8). We totally got 4 743 critical miRNA-target pairs across 11 cancer types, and 4 572 pairs of them showed much weaker correlation in tumor than in normal. The most frequent pairs that appeared in at least 4 types of cancer were shown in Table 2. These common critical pairs could contribute to aberrant cell growth, survival and cell motility in various cancer types [38][39][40]. MiR-21 was the most common miRNA that was found in many critical miRNA-target gene pairs [37,41].
. . Functional Annotation and miRNA-mRNA Network. We firstly identified the functions by literatures for those outstanding critical pairs in each cancer types, because these pairs showed an obviously strong correlation in normal samples but a greatly weaker correlation in tumor samples. These pairs were given in Table 3  cancer were given in Table S7. After studying the function of these genes by literature, we found these genes of these outstanding miRNA-target gene pairs were mostly involved in cell junction and cell cycle regulation. The dysregulated cell junction genes including the genes which enhance migration and invasion of cancer cells and the dysregulated cell cycle genes might enable the cells out of normal cell cycle.
To gain global insights into the biological processes regulated by the critical miRNA-target gene pairs, we performed gene ontology (GO) analysis to identify the functional categories of the target genes. The top 10 items for the molecular function, cellular component, and biological process categories of GO annotation were shown in Table S9. These included protein binding (p = 3.42E-05,  Bonferroni correction) and cell adhesion (p = 1.66E-03, Bonferroni correction), which are both linked to migration of cells. The genes like GNE, TRIM , and PDE A greatly influenced cell junction [42][43][44]. Thus, dysregulated miRNAtarget gene pairs involved in these genes like miR-21/GNE, miR-369/TRIM , and miR-203a/PDE A might promote the migration of cancer cells. Moreover, 38 critical genes in BRCA and 15 critical genes in KICH were significantly enriched in cell division (p = 6.29E-08, Bonferroni correction) and DNA replication cell proliferation (p = 5.40E-04, Bonferroni correction). The dysregulated miRNA-target gene pairs involved in these genes like miR-215/SKA , miR-3653/EMP and miR-192/PIM , might enable the cells to evade normal constraints on cell proliferation [45][46][47].

DEGs or DEMs
We also performed KEGG pathway analysis to identify the pathways regulated by these critical pairs. Table S10 showed 199 KEGG pathways (p-value < 0.05) that were enriched in the 11 cancer types. These included known cancer-related pathways, such as Pathways in cancer (p = 3.78E-07, Bonferroni correction, BLCA), cell cycle pathway (p = 1.20E-05, Bonferroni correction, BRCA), p53 signaling pathway (p = 8.83E-06, Bonferroni correction, HNSC), and AMPK signaling pathway (p = 3.60E-04, Bonferroni correction, KIRC). The dysregulated miRNA-target gene pairs involved these pathways might modulate the development and progression of tumors. The most frequently dysregulated pathways were listed in Table 4  For a better visualization of these miRNA-target gene interactions, we constructed the miRNA-mRNA network by using miRNA-target gene interaction in Table S6. Common miRNA-target pairs were usually deserved more attention because their universality in multicancer types. We constructed a pan-cancer network with the miRNA-target gene pairs that appear in at least three cancer types ( Figure 5). The network consisted of 108 edges connecting 35 miRNAs and 89 target genes. Among these, miR-21 and miR-30a were the most common regulators, which were linked to 20 and 10 target genes, respectively.

Discussion
The correlation coefficient, like Pearson, Kendall, Spearman method correlation coefficient, was widely used to measure the intensity of association between two variables [48].
miRNA's negative regulation on gene expression could result in an inverse correlation between miRNA and its target, so negatively correlation coefficient was proposed to be an effective way of miRNA-target pair identification [49], while seldom study analyzed the changes of miRNA-mRNA interactions in tumor and normal [22]. Our data suggested that the negative correlation of several critical miRNA-target gene pairs had a global reduction in tumor tissues compared with in normal. And these correlation reduced critical pairs were mostly participant in tumorigeneses and tumor progression.
Since the expression changes of miRNA or mRNA regulated by various factors could cause the changes of correlation between them, the factors who result in the reduced correlation in tumor tissues were indistinct. We speculated that correlation regulatory networks of these critical pairs were rewired in cancer cells by following reasons. Firstly, decreased expression of core miRNA in tumor tissues might      diminish the negative regulation of its target genes. Only the most abundantly expressed miRNAs could affect their target mRNAs stability by effectively binding to the available mRNA target sites [50]. Secondly, the miRNA-target gene pairs had specific interact manner due to particular vivo environmental [51]. The condition-specific regulation of miRNA might cause different correlation of miRNA-target pair in tumor and normal tissues. Thirdly, changes in the intracellular environment, especially cellular metabolism, might reduce the correlation of miRNA-target gene pairs. Highly proliferating cancer cells showed profound metabolic changes in response to the increased demands of energy, reducing power and building blocks to sustain their elevated biosynthetic and physiological needs [52]. The Warburg effect indicated that cancer cells increase glycolysis and alter the intracellular acidbase balance in cancer cells [53], which might not conducive to miRNA binding its target mRNAs. Moreover, transcription factors [54] and endogenous long noncoding RNAs [55] might also affect the correlation of miRNA-target gene pair. Whole gene methylation and gene copy number were significantly different between tumor and normal, which have been confirmed to be closely related to gene and miRNA expression [56]. Cause of variation of pairs correlation was so complex, there was a need for systematic study. Several previous studies had measured the correlation between miRNA and mRNA expression in tumor samples using a multivariate linear model that eliminate the noise of copy number alters and methylation changes [57,58]. Based on their method, we estimated the effect on correlation of miRNA-target gene pairs by methylation and copy number altersing. The detailed processes were given in the Supplementary Material, and the results were shown in Figures S3-S6 and Table S11. From these results, we found that methylation and copy number altering really affected the correlation of miRNA-target pairs but had no significance. This conclusion was consistent with results in Jacobsen's and Andres-Leon's researches. For example, the correlation coefficient of miR-20b/RORC in Jacobsen's work just increased from -0.08 to -0.20. And the average correlation coefficients of pairs in Andres-Leon E et al. work were similar to ours (both close to -0.1), which were all much smaller than the average correlation coefficient values of the pairs in normal tissue (nearly to -0.5).

Conclusion
In summary, this integrative analysis had confirmed the changed correlation of miRNA-target gene pairs in different clinical status and the reduced correlation of several critical pairs in tumor compared with in normal which was not observed in previous studies. These critical pairs were found to be mostly participant in biological processes and signaling pathways, such as focal adhesion, p53, MAPK, and PI3K-Akt signaling pathways, who were hubs pathways in cancer progression, invasion, and metastasis. This study had presented the molecular regulatory mechanism and the critical pairs would be useful for finding novel diagnostic biomarkers and therapeutic targets of cancer. Figure S1. Expression density distribution of mRNAs and miRNAs before and after filtering. Figure S2. Distribution of the PCC values of critical miRNA-target gene pairs in healthy controls and 4 clinically staged tumors. Figure S3: Simple correlations coefficients (Cor) and partial correlation coefficients (Pcor) of critical pairs in normal and tumor samples of LIHC. Figure S4: Simple correlations coefficients (Cor) and partial correlation coefficients (Pcor) of critical pairs in normal and tumor samples of LUAD. Figure S5: Simple correlations coefficients (Cor) and partial correlation coefficients (Pcor) of critical pairs in tumor samples of LIHC. Figure S6: Simple correlations coefficients (Pearson correlation coefficient) and partial correlation coefficients (given CNAs effects) of critical pairs in cancer samples of LUAD. Table S1: Sample sizes for each cancer type. Table S2. Sample clinical and grouping information for each cancer type. Table S3. List of DEMs for each cancer type. Table  S4. List of DEGs for each cancer type. Table S5. Frequency of DEMs and DEGs in 11 cancer types. Table S6. List of critical miRNA-target gene pairs for each cancer type. Table  S7: Interactional references and gene functional references for some critical miRNA-target pairs. Table S8. Frequency of critical miRNA-target pairs in 11 cancer types. Table S9. List of the significantly GO annotation terms (p < 0.05) for the critical target genes in each cancer type. Table S10. List of the significantly KEGG pathway terms (p < 0.05) for the critical target genes in each cancer type. Table S11. Partial and simple correlation coefficients for each critical pair in normal and tumor samples of LUAD and LIHC when controlling methylation or CNAs effect. (Supplementary Materials)