Abstract

Noncoding, endogenous microRNAs (miRNAs) are fairly well known for regulating gene expression rather than protein coding. Dysregulation of miRNA gene, either upregulated or downregulated, may lead to severe diseases or oncogenesis, especially when the miRNA disorder involves significant bioreactions or pathways. Thus, how miRNA genes are transcriptionally regulated has been highlighted as well as target recognition in recent years. In this study, a large-scale investigation of novel cis- and trans-elements was undertaken to further determine TF-miRNA regulatory relations, which are necessary to unravel the transcriptional regulation of miRNA genes. Based on miRNA and annotated gene expression profiles, the term “coTFBS” was introduced to detect common transcription factors and the corresponding binding sites within the promoter regions of each miRNA and its coexpressed annotated genes. The computational pipeline was successfully established to filter redundancy due to short sequence motifs for TFBS pattern search. Eventually, we identified more convinced TF-miRNA regulatory relations for 225 human miRNAs. This valuable information is helpful in understanding miRNA functions and provides knowledge to evaluate the therapeutic potential in clinical research. Once most expression profiles of miRNAs in the latest database are completed, TF candidates of more miRNAs can be explored by this filtering approach in the future.

1. Introduction

Among functional noncoding RNAs (ncRNAs), microRNAs (miRNAs) are tiny molecules (~21–23 nt) with giant roles. miRNAs participate in gene regulation by targeting messenger RNAs (mRNAs) and influencing their stability and the initiation of translation. It is implied that if the expression of miRNA is aberrant, miRNA-mediated gene circuitries will be disordered, resulting in homeostatic imbalance, pathogenesis, and oncogenesis [1, 2]. In recent years, elucidating transcriptional regulatory mechanisms of miRNA genes has been highlighted when studying miRNA function. Core promoters of miRNA genes were promptly identified for depicting full-length primary transcripts [35]. High-throughput sequencing datasets derived from epigenetic signature and TSS-relevant experiments unfold transcriptional start sites (TSSs) of miRNAs and offer a practical strategy to determine miRNA promoters [69].

To further understand upstream regulatory elements controlling miRNA expression, transcription factors (TFs) and their binding sites of miRNA promoters were deciphered by either literature survey or computational prediction [1013]. By integrating the information of TF-miRNA regulatory relations and miRNA target interactions, regulatory networks that revolved around miRNAs provide a biological insight into how miRNAs dominate functional processes in biochemical reactions or metabolic pathways [1417]. However, most of the foregoing studies regarded the upstream regions of pre-miRNAs as promoters (e.g., 10 kb upstream or −900~ +100 of the 5′-start of the pre-miRNA). More convinced miRNA promoters should be used for searching putative TF-binding sites. In addition, the experimentally verified TF-miRNA relations are insufficient for current miRNAs. According to the statistics in the latest version 1.2 of TransmiR, only 735 entries, which include ~201 transcriptional factors, ~209 miRNAs, and 16 organisms from 268 publications, were curated. A large-scale investigation of novel cis- and trans-elements is necessary to fulfill the unmet need for locating transcription factor binding sites (TFBSs) within miRNA promoter regions.

Since sequence-specific TFs possess DNA-binding domains (DBDs) to recognize specific motifs in miRNA promoter sequences, potential binding sites can be detected by sequence-based computational approaches, for example, position weight matrix (PWM). A position weight matrix (also called position specific scoring matrix, PSSM) infers a pattern of DNA segment and is widely applied in searching TFBSs [1821]. Two well-known databases collecting matrix information are TRANSFAC [22] and JASPAR [23]. To avoid excess false positives due to short sequence motifs used for TFBS pattern search, cis-regulatory analysis of coregulated gene sets was executed to determine overrepresented TFBSs [24, 25]. Although previously published web server Pscan [26] and oPOSSUM [27] can search transcription factor binding sites in coexpressed gene promoters to remit the impact of false positive issue, users have to define their coexpressed gene groups with valid gene IDs. In addition, miRNA promoter sequences are not included in these two web servers.

Therefore, the purpose of this study is to create a computational pipeline to undertake large-scale investigation of novel cis- and trans-elements for human miRNA genes based on coexpression strategy. First, we constructed 255 coexpressed gene groups of human miRNAs. Moreover, instead of grapping the upstream regions of pre-miRNAs as promoter sequences, we exploited more concrete human miRNA promoters by following the processes as previously described [7]. Through the detection of transcription factor binding sites within promoter regions of human miRNAs and their coexpressed genes by using matrix information from TRANSFAC, the common TFBSs were identified. We then filtered the redundancy by not only the occurrence of common TFBSs but also the expression correlation between TF-encoded genes and corresponding human miRNA gene groups. Finally, more reliable TF-miRNA regulatory relations of 225 human miRNAs were provided. Furthermore, the liver-specific hsa-miR-122 was also selected as the case study to demonstrate the usage of this filtering approach and its practicability.

2. Materials and Methods

2.1. Human miRNA Promoter Sequences

The latest genomic coordinates of human miRNAs were retrieved from miRBase release 19 [28]. Human miRNA TSSs were identified by reproducing the computational procedures described in miRStart resource [7]. The formula which determines tag-enriched loci was modified by weighted tag density (tag density multiply by tag number) in the SVM step to precisely reveal the effect of tag intensity. Then, 1 kb upstream sequences of human miRNA TSSs were acquired from UCSC Genome Browser [29] (hg19/GRCh37 assembly) in FASTA format. The updated intergenic miRNA TSSs in human genome are listed in supplementary Table (available online at http://dx.doi.org/10.1155/2014/623078) in additional file for reference.

2.2. Expression Profiles of Human miRNAs and Annotated Genes

In order to determine human genes that are coexpressed with specific miRNA genes, GDS596 record, the Affymetrix gene expression profiles from 79 physiologically human normal tissues, was downloaded from Gene Expression Omnibus (GEO) [30, 31]. Among annotated human genes in GDS596, genes were filtered out if their HGNC symbols are invalid or have been withdrawn. On the other hand, the expression data of 345 miRNAs in 40 normal human tissues generated by a new type of real time reverse transcription- (RT-) PCR-based miRNA assays were also collected [32]. Mature miRNAs with eliminated miRBase IDs (release 19) were discarded. In both expression datasets, only 17 tissues (adrenal, brain, heart, kidney, liver, lung, lymph node, ovary, pancreas, placenta, prostate, skeletal muscle, testicle, trachea, thymus, thyroid, and uterus) are in common and were considered to be integrated expression profiles of annotated genes and miRNA genes with 17 conditions. To standardize expression levels across extensive range of both datasets, the raw intensity was -score transformed [33].

2.3. Coexpressed Gene Groups of Human miRNAs

After the transformation process between two expression datasets, Pearson’s correlation coefficient (PCC) was calculated to estimate which annotated genes are coexpressed with specific miRNAs. The formula of Pearson’s correlation coefficient (usually using the letter ) is as follows: where and denote the expression level of condition (tissue) of two genes and that are calculated for  , whereas and represent the average of corresponding expression levels in total conditions. Since the value of Pearson’s correlation coefficient ranges from +1 (positively correlated) to −1 (negatively correlated), we defined that two genes of interest are coexpressed if their Pearson’s correlation coefficient is more than 0.8. Subsequently, coexpressed genes of human miRNAs were determined, and their promoter sequences (1 kb upstream from TSS) were obtained from BioMart of Ensembl release 69 [34].

2.4. CoTFBSs in Human miRNA Promoters

In this work, we defined coTFBS as a common TFBS located in promoter regions of coexpressed genes. The TRANSFAC database [22], comprising data on transcription factors, their target genes, and regulatory binding sites, has been widely used when studying eukaryotic transcriptional regulation. After acquiring 1 kb promoter sequences of miRNA genes and their coexpressed genes in FASTA format, the Match program [35] was executed for sequence motif search of transcription factor binding sites according to the matrix information provided by TRANSFAC (version 2011.4). A customized profile which specified human matrices in TRANSFAC library was used for Match to minimize both false positive and false negative rates with core similarity and matrix similarity cut-off values for each matrix. As defined in the publication of Match, the core of each matrix refers to the first five most conserved consecutive positions of a matrix. A putative TFBS with core similarity less than 1 was filtered out. The occurrence of each coTFBS and its expression correlation with miRNA in coexpressed gene groups was calculated finally. Only those whose encoded genes are coexpressed with corresponding miRNA were considered.

3. Results

3.1. Human miRNAs with Coexpressed Gene Group for Detecting TF-miRNA Relations

Figure 1 summarizes the workflow to investigate human TF-miRNA regulatory relations based on expression profiles of miRNA and annotated genes. Pearson’s correlation coefficient (PCC) was applied to measure the similarity of expression patterns across 17 human normal tissues, which represents the coexpressed level between a specific miRNA gene and an annotated gene of interest. Due to the reason that Pearson’s cannot be calculated if the expression levels in all 17 conditions are identical, 29 mature miRNAs with such expression profiles were excluded. After discarding the mature miRNAs whose miRBase IDs are eliminated, 289 human miRNAs were selected to estimate PCC values with annotated human genes in GDS596 record.

Among them, however, only 255 human miRNAs have more than one coexpressed genes (PCC > 0.8). Moreover, 25 out of 255 human miRNAs have no identified TSSs by following the previous strategy [7] and have to be filtered. In total, 230 human miRNAs were qualified to discover putative cis- and trans-elements in their promoter regions. The number of members of each miRNA coexpressed gene group ranges from several to thousand (see Table in additional file for details).

3.2. Putative TF-miRNA Relations Were Explored according to the Occurrence of CoTFBS

Theoretically, a group of genes that are coexpressed may be regulated by common transcription factors. Based on this concept, a specific transcription factor binding site located in the promoter regions of most genes in each coexpressed group implies that its corresponding TF is the most possible one controlling the expression of these genes. Here, we introduce the term “coTFBS” which represents the common TFBS of a coexpressed gene group to filter redundant TF candidate of miRNA genes. For example, the “pink rectangle” binding motifs were detected in the proposed miRNA and other four gene promoters in Figure 2. The occurrence of this coTFBS is five. By scanning TFBSs in promoter regions of each miRNA gene and its coexpressed gene group using Match based on TRANSFAC library (version 2011.4) and following the filtering process, the putative TFs that regulate a specific miRNA were determined. In this work, we successfully identified putative TF-miRNA relations of 225 human miRNAs. The full list of putative TFs for each human miRNA promoter can be accessed in Table S3.

3.3. Case Study: hsa-miR-122

hsa-miR-122 is one of the intergenic miRNAs whose TSS and promoter have been experimentally characterized [7]. Previous studies reported that this liver-specific miRNA is significantly downregulated in hepatocellular carcinoma and profoundly affects carcinogenesis [36]. Because of the explicit promoter and biological importance, hsa-miR-122 was selected as the case study to investigate which TFs may regulate its gene expression. 261 annotated genes coexpressed with miR-122 were identified according to their expression levels with Pearson’s correlation coefficient (PCC) more than 0.8. Figure 3 compares the expression patterns of 261 coexpressed genes (the pink line represents the average value of their expression) with hsa-miR-122 among 17 human normal tissues, indicating that remarkable peaks appear in liver for all the coexpressed genes. The expression image (see Figure in additional file) also reveals the similar trends between hsa-miR-122 and its coexpressed genes.

Then, the promoter sequences (1 kb upstream from TSS) of hsa-miR-122 gene and 261 coexpressed genes were collected to identify coTFBSs using Match. The occurrence of putative transcription factors of miR-122 is listed in Table 1. Among TF candidates regulating hsa-miR-122, the TF binding motif of HNF-4alpha can be found in 191 coexpressed gene promoters. In 2011, Li et al. reported that HNF-4alpha is a key regulator positively controlling the expression of miR-122 in liver [37], proving our computational finding. They performed not only the luciferase reporter gene assay to detect the trans-activation effect of HNF-4alpha in miR-122 promoter but also the ChIP and EMSA assays to determine HNF-4alpha binding of miR-122 promoter in vitro and in vivo. Moreover, other liver-enriched transcription factors including HNF-1alpha, HNF-3alpha, HNF-3beta, and HNF-6 showed a strong positive correlation with miR-122. The knockout of HNF1A, FOXA1, and FOXA2 by RNAi assay reduces the expression of miR-122, suggesting that these transcription factors may bind to miR-122 promoter and transcriptionally regulate miR-122 [38]. Importantly, HNF-1alpha, HNF-3alpha, and HNF-3beta were identified in our list of TF candidates, and a HNF-6 binding site was determined (−2720 from miR-122 TSS) if the 3 kb promoter sequence of miR-122 was used.

In addition, the TF binding site of NR1B2 can be found in 226 coexpressed gene promoters. A previous research article published in 1987 indicated that the inappropriate expression of HAP gene (the official HGNC symbol is RARB) may relate to the hepatocellular carcinogenesis, and hap protein may directly participate in the hepatocellular transformation [39]. It is implied that transcription factor NR1B2 may bind to miR-122 promoter and regulate its expression.

4. Discussion

For large-scale investigation of human TF-miRNA relations in this study, the expression profiles of human miRNAs and over ten thousands genes from normal human tissues facilitate the acquirement of miRNA coexpressed gene groups. However, the latest usable RT-qPCR miRNA expression data for human normal tissues were published in 2007 [32]. The currently available data are almost derived from cancer cell lines or tumor tissues. Totally 230 coexpressed gene groups limit the coTFBS analysis of most miRNAs in existence. Besides, based on the cut-off 0.8 PCC value applied to define coexpression between a miRNA and an annotated gene, 34 miRNAs have no coexpressed genes and 50 miRNAs have less than ten coexpressed genes. Although each miRNA has sufficient coexpressed genes to analyze its putative cis- and trans-elements if the lower PCC cut-off was used, it may cause the trade-off of specificity when determining miRNA coexpressed genes. For example, hsa-let-7a-1 has no coexpressed gene when using 0.8 PCC value but has 29 coexpressed genes if 0.6 PCC value was applied.

It is noteworthy that the normalization process between two raw expression datasets was only -score transformed, without using log10 transformation before -score standardization. According to the definition of formula, Pearson’s correlation coefficient calculates linear correlation between two variables and has an invariant property in statistics. In fact, the log10 transformation tends to alter the original expression patterns and sharply reduces the scale of raw intensity, resulting in unexpected affection of PCC values. Figure illustrates an example of hsa-miR-122 expression patterns with and without log10 transformation. In addition to the high expression level in liver, two obvious peaks of brain and thymus appear in the log10 transformed expression profile of hsa-miR-122. The members of hsa-miR-122 coexpression group also reflect the variation. Unlike more than two hundred coexpression genes by using -score transformed expression data, merely five coexpressed genes left based on the cut-off 0.8 PCC value by -score and log10 transformed data.

Another limitation of TFBS/TF detection depends on present transcription factor binding matrices collected in TRANSFAC. As is well known, sequence motif search of transcription factor binding sites is executed by the Match program according to the matrix information in TRANSFAC library. A specific TFBS in miRNA promoter will not be detected if the corresponding matrix has not been obtained by in vitro selection studies. Here is an example. Unlike hsa-miR-122, miR-224 is upregulated in HCC through epigenetic mechanisms and controls several crucial cellular processes [40]. Wang et al. indicated that miR-224 expression is reciprocally regulated by HDAC1, HDAC3, and EP300. Our result shows that P300 (encoded by EP300) is the TF candidate of miR-224, and its TFBS can be found in 227 coexpressed gene promoters. Because matrices of histone deacetylases 1 and 3 are not available in TRANSFAC, no such transcription factors were predicted.

In conclusion, rather than traditional PWM search characterizing putative cis-elements in promoter regions, the coTFBS strategy was developed to determine more confident TFBSs for human miRNAs. The investigation was restricted by the incomplete expression profiles of present human miRNAs. Once most expression profiles of miRNAs in the latest database are available, TF candidates of more miRNAs can be explored in the future. Furthermore, although more and more ChIP-seq data were generated and were useful to identify transcription factor binding sites [4143], the corresponding ChIP-seq data of specific TFs are still the minority. It is expected that large-scale ChIP analysis of general TFs contributes more confident TFBSs by observing the aggregated peaks within promoter regions of human miRNAs.

5. Conclusions

In organisms, not all of ribonucleic acids (RNAs) are translated to proteins. miRNAs are such noncoding RNAs which play critical roles in gene regulation, even if it is generally believed that proteins convey vital information from genes and execute biological functions to maintain life processes. Although target prediction has been the mainstream when studying miRNA functions for a while, researchers start to explore TF-miRNA interactions and study the transcriptional regulation of miRNAs, which are necessary to depict how miRNAs participate in diverse biological processes. To determine putative TFs and TFBSs located in human miRNA promoters, we created a computational pipeline which not only allows large-scale investigation as long as the expression profiles of miRNAs are available, but also filters the redundancy when searching short sequence. This valuable information is helpful in understanding miRNA functions and provides knowledge to evaluate the therapeutic potential in clinical research.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Authors’ Contribution

Chia-Hung Chien, Shun-Long Weng, Wen-Chi Chang, Ann-Ping Tsou, and Hsien-Da Huang conceived and designed the experiments. Chia-Hung Chien and Yi-Fan Chiang-Hsieh analyzed the data and performed the experiments. Chia-Hung Chien wrote the paper.

Acknowledgments

Thanks are due to Pei-Wen Jiang for the assistance of programming and system maintenance and to Hsi-Yuan Huang for experience sharing in expression profile analysis. This research was supported by a grant from National Science Council of the Republic of China for financially supporting this research under Contract nos. NSC 102-2313-B-006-004, NSC 101-2311-B-009-003-MY3, NSC 100-2627-B-009-002, NSC 102-2911-I-009-101, and NSC 99-2628-B-006-016-MY3.

Supplementary Materials

The additional file provides the expression profile of hsa-miR-122 and its 261 co-expressed genes among 17 human normal tissues (Figure S1 and S2), and full list of TSSs, the number of co-expressed gene groups, and putative TFs of human miRNAs (Table S1-S3).

  1. Supplementary Materials