Screening of Target Genes and Regulatory Function of miRNAs as Prognostic Indicators for Prostate Cancer

Background

Prostate cancer, the second leading cause of cancer-related deaths in men, is a complex and multifactorial disease, including genetic variation, daily habits, and environmental factors. The presence of prostate cancer may be indicated by symptoms, a physical examination, serum prostate-specific antigen (PSA), or biopsy [1]. However, as the most common tool used to detect prostate cancer, PSA testing increases cancer detection but does not decrease prostate cancer-associated mortality rates [2]. Therefore, identification and development of new potential biomarkers could be beneficial for the diagnosis and prognosis of prostate cancer patients [3].

MicroRNAs, a class of small non-coding RNAs, regulate protein-coding gene expression by repressing translation or cleaving RNA transcripts in a sequence-specific manner [4]. miRNAs could play important roles in various processes such as cell development, differentiation, proliferation, cell-cycle control, apoptosis, and metabolism [5]. Aberrantly expressed miRNAs (i.e., up-regulated and down-regulated miRNAs) are known to disrupt the expression of several mRNAs and proteins, and can be involved in the pathogenesis of various human cancers [6–8]. A growing body of evidence suggests that microRNAs expression profiling of prostate cancer is increasing in importance due to its usefulness in diagnosis, staging, progression, prognosis, and response to treatment. In 2010, Szczyrba et al. [9] demonstrated that 33 miRNAs in prostate cancer tissues were either up-regulated or down-regulated by more than 1.5-fold compared to that found in normal tissues. They suggested that miR-143 and miR-145 could be involved in the development of prostate cancer. Recently, Song et al. [10] showed that 20 abnormal prostate cancer-related miRNAs, including hsa-miR-26a, hsa-miR-152, hsa-miR-19a, hsa-miR-30c, hsa-miR-19b, and hsa-miR-146b-5p, could be used to predict the clinical behavior of prostate cancer by using their high-throughput miRNA expression profiling. Consequently, miRNAs are being regarded as new oncogenes or tumor suppressor genes and new prognostic and predictive biomarkers for prostate cancer prognosis and treatment [11].

The main purpose of this work is to identify specific miRNAs that are closely associated with prostate cancer by analyzing significantly altered miRNAs in a large dataset, and to investigate the signal pathway and function of their target genes, which might provide evidence that these miRNAs can serve as prognostic and predictive biomarkers for early diagnosis and treatment for prostate cancer.

In this study, we identified differentially expressed miRNAs between normal and prostate cancer tissues, then ROC curves were generated to determine the accuracy for miRNAs in the diagnosis of cancer. Besides, StarBase, and TargetScan software were used to predict the miRNA’s target group and target genes, respectively. GO functional and KEGG pathway analyses were conducted on up/down-regulated expressed miRNAs. Eventually, we found that 6 miRNAs – hsa-mir-153-2, hsa-mir-92a-1, hsa-mir-182, hsa-mir-29a, hsa-mir-10a, and hsa-mir-221 – could function as prognostic and predictive markers for use in treatment of prostate cancer.

Material and Methods

TCGA MIRNA DATASET AND PATIENT INFORMATION:

Expression data of miRNAs and the corresponding clinical data for patients were downloaded from The Cancer Genome Atlas (TCGA, http://cancergenome.nih.gov/) database. Expression levels of miRNAs were detected using Illumina HiSeq system. As the third stage data in TCGA, standardized miRNAs data (consisting of 1046 annotation miRNAs expression values) were extracted with zero expression data removed. Because the interior standardization of each downloaded miRNAs data were completed, we used Limma package of R language to achieve standardization between samples to remove experimental error.

SCREENING OF MIRNAS DIFFERENTIALLY EXPRESSED:

We screened differentially expressed miRNAs between normal prostate tissue samples and prostate cancer tissue samples using SAMR [12] of R software, in a condition of FC=2 and FDR<0.05. To discriminate the normal prostate tissue samples and prostate cancer tissue samples, principal component analysis of the differentially expressed miRNAs was used.

RECEIVER OPERATING CHARACTERISTIC ANALYSES:

ROC analysis for 6 significantly differentially expressed miRNAs related to patients was performed using MedCalc software [13]. Corresponding expression levels of miRNAs were sorted, and then displayed the status case-by-case (sicken for 1, survival for 0). The AUC under binomial exact confidence interval was calculated to generate the ROC curve.

ANALYSIS OF TARGET GENES OF THE DIFFERENTIALLY EXPRESSED MIRNAS ON LNCRNA, CIRCRNA, PSEUDOGENE, AND SNCRNA:

The most differentially expressed miRNAs were illustrated, and target sites of miRNAs on lncRNA, circRNA, pseudogene, and sncRNA were predicted by using the StarBase website. One experiment was at set up to prove that the miRNAs could bind to the corresponding RNAs.

FUNCTIONAL ANALYSES OF THE TARGET MIRNAS:

Putative targets of miRNAs were predicted by TargetScan software. Then the Functional Gene Ontology (GO) biological processes terms of the putative targets of candidate miRNAs, as well as the KEGG pathway analysis, were performed by DAVID [14]. For each GO score, the P value of the function enrichment and P value of Benjamini correction were calculated. The role of the MAPK pathway in prostate cancer was analyzed by BioChart and KEGG [15].

NETWORK ANALYSIS OF CORRELATION BETWEEN MIRNAS AND TARGET GENES:

High aggregation is an internal property of the biological network; it reflects a highly modular gene network. One method that can be used to distinguish those module groups with different scales and specific function is to decompose into relatively independent modules before analyzing the network containing a large number of genes.

Network topological property of proteins was analyzed by Network Analyzer [16], a software package of Cytoscape. Then Clusterone was used to modularize the network functions. The first 5 modules with enrichment P value <1.0e-05 were selected for functional analysis.

IDENTIFICATION OF CANDIDATE PROGNOSTIC MARKERS BY KAPLAN-MEIER METHOD:

We randomly divided the dataset into training and testing datasets, and applied the Cox regression model to the training dataset to determine the survival association markers (p<0.01). We then used this model to predict the testing dataset and risk score for each patient. According to this model, a high score means weak survival ability and the low score means strong survival ability. Finally, we used Kaplan-Meier analysis to identify the survival time distribution, and used log-rank analysis to judge the significance.

Results

DIFFERENT MIRNA EXPRESSION LEVEL BETWEEN SAMPLES:

There were a total of 1041 miRNA expression values in the dataset of 551 samples, including 52 normal individuals and 499 prostate cancer patients (418 patients had clinical information). Expression levels of 494 human miRNAs were collected and discrepancies of miRNA levels between groups were tested by principal component analysis and cluster analysis (Figure 1).

DIFFERENTIAL EXPRESSION ANALYSIS OF MIRNAS:

A total of 162 miRNAs were differentially expressed between normal samples and prostate cancer samples, of which 128 were up-regulated and 38 were down-regulated. The up-regulated miRNAs – hsa-mir-153-2, hsa-mir-92a-1, and hsa-mir-182 – accounted for 77.1% of the differentially expressed miRNAs, while the down-regulated miRNAs – hsa-mir-29a, hsa-mir-10a, and hsa-mir-221 – accounted for 22.9%.

ESTIMATION FOR THE ACCURACY OF MIRNAS IN THE DIAGNOSIS OF PROSTATE CANCER AS BIOMARKER:

AUC analysis of the top 6 differentially expressed miRNAs showed that 3 miRNAs were up-regulated and the other 3 were down-regulated (Figure 2). miRNAs with AUC greater than 0.5 could serve as biomarkers in the diagnosis and prognosis of prostate cancer. The AUC of the first 3 up-regulated miRNAs – hsa-mir-153-2, hsa-mir-92a-1, and hsa-mir-182 – were 0.929, 0.968, and 1.000 and the AUC of the first 3 down-regulated miRNAs – hsa-mir-29a, hsa-mir-10a, and hsa-mir-221 – were 0.700, 0.609, and 0.858 (Table 1).

PREDICTION FOR THE PUTATIVE TARGETS OF THE DIFFERENTIALLY EXPRESSED MIRNAS:

The most significantly differentially expressed miRNAs – hsa-mir-153-2, hsa-mir-92a-1, hsa-mir-182, hsa-mir-204, hsa-mir-10a, and hsa-mir-221 – could serve as biomarkers in the prediction of prostate cancer. Therefore, the putative targets of these 6 miRNAs were further analyzed (Figure 3, Table 2). Through use of the StarBase website, target sites of miRNAs on lncRNA, circRNA, pseudogene, and sncRNA were predicted.

FUNCTIONAL ANALYSIS OF THE TARGET GENES OF MIRNAS OF INTEREST:

Putative targets of the differentially expressed miRNAs were created by software, of which 7992 were down-regulated miRNAs and 35483 were up-regulated miRNAs. In gene ontology and KEGG-enrichment analysis (Figure 4), 7 GO annotations and 7 KEGG annotations showed significant enrichment of the down-regulated miRNAs’ target genes. Positive ion combination and JAK-STAT pathway were the most obvious annotation of GO and KEGG, respectively. There were 12 GO annotations and 11 KEGG annotations that showed significant enrichment of up-regulated miRNAs’ target genes. The function of plasma membrane and MARK signaling pathway was the most obvious annotation of GO and KEGG, respectively. Obviously, GO and KEGG annotations showed that both target genes of up/down-regulated miRNAs were mainly enriched in function of plasma membrane and the pathways in cancer. The other KEGG enrichment annotations for the target genes of up-/down-regulated miRNAs were also associated with prostate cancer and MARK signaling pathway. In addition, the KEGG network analysis showed that 72 target genes among the down-regulated ones were enriched in the MAPK signal pathway and 111 target genes among the up-regulated ones were also enriched in the MAPK signal pathway (Figure 5).

THE NETWORK ANALYSIS OF THE INTERACTIONS BETWEEN THE MIRNAS AND THEIR TARGET GENES:

Network topology calculation is shown in the following figures: the degree distribution (Figure 6A), the shortest path length distribution (Figure 6B), the average cluster coefficient degree (Figure 6C), and the intimate centrality (Figure 6D). These parameters showed that the network node degree distribution for the interaction between proteins encoded by the differentially expressed genes was consistent with the power-law distribution of the network and had structural properties.

Applying the Cluster one in Cytoscape software, the sub-modules in the interactional network of the proteins, which were correlated with differentially expressed genes, were screened. Two modules under the condition of more than 5 nodes were screened and modular significance P value less than 0.05 were obtained. Figure 7 describes the analysis of the function module in network. Module 1 was enriched in pathway in cancer, indicating these genes were mainly involved in cancers. The abnormally expressed genes caused diseases when the miRNAs combined with them. Module 2 was enriched in the Wnt signal pathway, and abnormal activation of the Wnt signaling pathway played a very important role in the pathogenesis and development of prostate cancer. Androgen receptor (AR) was the key in the transformation between the occurrence and development of prostate cancer and androgen-independent prostate cancer (AIPC). The AR signal pathway has been a major target for prostate cancer research and the Wnt signal pathway played an important role in the target therapy of tumors. Therefore, studies on the Wnt signal pathway contributed importantly to the treatment of disease.

SURVIVAL ANALYSIS OF EACH MIRNA:

We used a training dataset to construct a Cox multiple variance regression model; the prognostic function is Prognostic score=(−4.722× expression level of hsa-mir-483)+( −4.449× expression level of hsa-mir-675)+( −2.745× expression level of hsa-mir-139)+(4.866× expression level of hsa-mir-625)+(5.602× expression level of hsa-mir-187). In the model, 3 miRNAs were risky, and 2 miRNAs were protective. After divide the patients into high- and low-risk groups, we used the log-rank test model to test the difference between the 2 groups (Figure 8). The training dataset showed a very large difference between high- and low-risk patients and the result of Kaplan-Meier analysis and heatMap showed miRNA expression in high- and low-risk patients (Figure 9).

Discussion

CONCLUSIONS:

In summary, we identified hsa-mir-153-2, hsa-mir-92a-1, hsa-mir-182, hsa-mir-29a, hsa-mir-10a, and hsa-mir-221 and determined interactions between their target genes as more specific and accurate prognostic and predictive indicators to clinical outcome of prostate cancer patients, implying their utility as diagnostic and prognostic tools and treatments. However, further research is needed to investigate the individual role of these abnormal prostate cancer-related candidate miRNAs in the process of malignant progression of prostate cancer.