TargetScore used to reveal potential targets of miRNA203 and miRNA-146a in psoriasis by integrating microRNA overexpression and microarray data.

OBJECTIVE: Systematic tracking of microRNA (miRNA) targets remains a challenge. In our work, we aimed to use TargetScore to investigate the potential targets of miRNA203 and miRNA-146a in psoriasis by integrating miRNA overexpression information and sequence data, and to further uncover the functions of miRNA203 and miRNA-146a in psoriasis.
METHODS: This was a case-control bioinformatics analysis using already published microarray data of psoriasis. We calculated targetScores by combining log fold-change and sequence scores obtained from TargetScan context score, probabilities of conserved targeting, and derived the distribution of targetScores. The scoring cutoff was chosen based on the different targetScore distributions for the nonvalidated and validated targets. The potential target genes for miRNA-203 and miRNA-146a were predicted based on the targetScore threshold. To reveal the functions of miRNA-203 and miRNA-146a, we implemented pathway enrichment analyses for the targets of miRNA-203 and miRNA-146a.
RESULTS: TargetScore >0.4 was selected as the threshold to filter out less confidence targets because we observed little overlap between the 2 distribution at targetScore = 0.4. Based on the targetScore >0.4, 49 target genes for miRNA-203 and 17 targets for miRNA-146a were identified. Pathway enrichment results showed that the target genes of miRNA-203 (including KIR2DL1, HLA-DQA1, KIR3DL1) only participated in antigen processing and presentation. The target genes of miRNA-146a (covering ADORA3, CYSLTR2, HRH4) were only involved in neuroactive ligand-receptor interaction.
CONCLUSION: MiRNA203 and miRNA-146a played important roles in psoriasis progression, partially through regulating the pathways of antigen processing and presentation, and neuroactive ligand-receptor interaction, respectively.

Introduction

Psoriasis is a common, complicated, and chronic inflammatory disorder, affecting 2% to 3% of the general population, which is related to a high risk of developing systemic comorbidities, including cardiovascular disease, obesity, and diabetes mellitus.[ Psoriasis skin lesions are typically characterized by hyperproliferation and abnormal differentiation of keratinocyte, infiltration of inflammatory cells into epidermis, and increased vascularity in the dermis.[ Psoriasis is a lifelong autoimmune disease, which severely reduces the quality of patients’ life.[ However, the mechanisms of psoriasis remain still unclear and a complete cure is lacking.

In the span of a decade, substantial advances have been made in revealing the pathogenesis of psoriasis. Based on microarray technology, several differentially expressed genes have been detected in psoriatic skin, which are defined as the molecular signature of this disease.[ In addition, studies have indicated a key role for microRNAs (miRNAs) in the biology of psoriasis.[ Several studies have reported that abnormal miRNAs and the predicted target genes exert crucial functions in psoriasis pathogenesis, for example, miRNA-99a and target gene insulin-like growth factor 1 receptor, and miRNA-424 and target gene cyclin E1. These miRNAs and their corresponding target genes mainly participate in the cycle of cell proliferation and apoptosis in the epidermis.[ Significantly, miRNA-203 was the first skin and keratinocyte-specific miRNA, which was found to regulate the inflammatory pathways, and keratinocyte proliferation as well as differentiation, and was observed to be significantly upregulated in psoriasis patients.[ The miRNA-203 has been demonstrated to play a potential role in disease progression via inhibiting suppressor of cytokine signaling 3 (SOCS3).[ In addition, miR-203 has been recently observed to target SOCS6.[ In a nutshell, these findings indicate that overexpression of miRNA-203 facilitates repression of immunosuppressive genes and exerts potential functions in disease progression. Nevertheless, there is no formal proof of this assumption, and miRNA-203 has additional targets in a complex network of immune regulated interactions in which the miRNA-203 under certain circumstances serves as an activator and under other conditions as a suppressor of immune response. In addition, miRNA-146a has been reported to be overexpressed in psoriatic skin lesions.[ Moreover, miRNA-146a has an important function in the differentiation or maintenance of CD4+ T regulatory cells, emphasizing its role as a negative regulator of inflammation signaling.[ However, the functions of miRNA-146a dysregulation in psoriasis remains unclear. Thus, in our study, we decided to investigate the functions of miRNA-203 and miRNA-146a in psoriatic.

Of note, identification of miRNAs functional characterizations relies on accurate detection of their targets. Nevertheless, it is difficult to experimentally extract miRNA-mRNA interactions. Excitedly, computational prediction methods provide a rapid alternative tool to uncover putative miRNA targets. Most of these prediction methods are conducted on the basis of sequence complementarity, evolutionary conservation, and target site accessibility.[ Unfortunately, precise prediction of miRNA targets remains a challenge with <50% specificity and having poor overlap among them. Of note, miRNA overexpression data combined with mRNA expression profiling has been indicated to be a promising approach.[ However, expression data are noisy, and alterations in expression can be induced by indirect regulatory effect by miRNAs.[ To date, few studies are specifically developed for transfection-based miRNA target prediction. Of note, Sylamer is the earliest work generated to extract the seed regions among the top sorted targets according to P values.[ Worriedly, Sylamer method does not analyze the distribution of fold-changes.

Thus, TargetScore, a Bayesian probabilistic scoring method taking into account of the fold-change due to miRNA overexpression and sequence-based information, was used in our study to identify the potential targets of miRNA-203 and miRNA-146a in psoriasis.

Materials and methods

This was a case-control bioinformatics analysis using already published microarray data of psoriasis. The study protocol was approved by the People's Hospital of Rizhao. And all the data collection processes were in compliance with all laws and regulations.

In our analysis, miRNA-overexpression data corresponding to 44 psoriasis patients at baseline, 30 healthy controls, and 2 distinct miRNAs (miRNA-203 and miRNA-146a) were firstly downloaded. Then, calculation of log fold-change (logFC) for all genes was conducted. Next, we calculated TargetScores by combining logFC and sequence scores obtained from TargetScan context score (TSCS), probabilities of conserved targeting (PCT), and derived the distribution of TargetScores. The potential target genes for miRNA-203 and miRNA-146a were predicted based on the TargetScore values. In order to reveal the functions of miRNA-203 and miRNA-146a, we implemented pathway enrichment analyses for the targets of miRNA-203 and miRNA-146a, respectively.

MiRNA overexpression data collection

Firstly, miRNA overexpression data corresponding to 1 EMBL-EBI set, 1 platform, 81 human samples, and 2 distinct miRNAs (miRNA-203 and miRNA-146a) were collected. The microarray profile of E-GEOD-55201[ including 81 samples (44 psoriasis patients at baseline, 7 psoriasis samples treated with ixekizumab, and 30 healthy controls) were downloaded from the EMBL-EBI based on the platform of AFFY-44 - Affymetrix GeneChip Human Genome U133 Plus 2.0 [HG-U133_Plus_2]. In order to investigate the molecular mechanism of psoriasis, we only selected the psoriasis patients at baseline and healthy controls for further analysis.

Data preprocessing and calculation of logFC for all genes

The Affy package[ was used to preprocess the gene expression profile of E-GEOD-55201. Specifically, data were corrected using background adjustment, following by quartile data normalization,[ and probe summarization. Then, probes were mapped to the gene symbols.

As all know, logFC stands for the changes in the expression levels of genes. Thus, we determined the logFC values of each gene between psoriasis and control groups. For mRNAs interrogated by multiple probes in the current study, we took the average of the fold-changes. Finally, we obtained the logFC values for 20,514 genes.

Brief outline of proposed method

This proposed method was a probabilistic approach for miRNA target prediction algorithm via combining miRNA-overexpression information and sequence-based scores obtained from other prediction approaches. In brief, every score feature is considered an independent variable as input to a Variational Bayesian-Gaussian Mixture Model (VB-GMM). A Bayesian was selected over a maximum likelihood method to avert overfitting. Concretely, with regard to the expression fold-change, 3-component VB-GMM was used to infer downregulated targets accounting for genes with little or positive fold-change (because of off-target effects).[ If not, we applied 2-component VB-GMM to unsigned sequence scores. Variational Bayesian Expectation-Maximization (VB-EM) algorithm was employed to optimize the parameters of the VB-GMM. The mixture component with the maximum absolute averages of negative fold-change or sequence value was connected with the targets of miRNA and determined as “target component.” The other components correspond to the “background component.” As a result, inferring miRNA-mRNA pairs was equal to speculating the posterior distribution of the miRNA targets. TargetScore was calculated as the transformed fold-change, which was weighted by the mean posteriors of target components over logFC, TSCS, and PCT.

Bayesian mixture model

Supposing there were M genes, we defined x = (x1, …, xM) as the expression of logFC (xf) or sequence scores (xs). Therefore, for L sets of sequence scores, x ∊ {xf, x1, …, xL}. With the goal of simplifying the following formulas, we used x to stand for one of the independent variables without loss of generality. In order to deduce the target genes for a miRNA given x, the posterior distribution p(z|x) of the latent variable z ∊ {z1, …, zK} was needed to obtain, in which K = 3 (K = 2) for modeling signed (unsigned) scores including logarithmic fold-changes (sequence scores).

The standard Bayesian-GMM was followed according to Bishop with minor modifications. In detail, the latent variables z were sampled at probabilities (mixing coefficient), which followed a Dirichlet prior Dir (π|α0) with α0 = (α0,1,………α0,). In order to interpret the relative frequency of targets and nontargets for any given miRNA, we set the α0,1 (connected with the target component) to aM and other r0, = (1 − a) − M/(K − 1), in which a = 0.01 (default). Supposing x followed a Gaussian distribution M (x|μ, Λ−1), where Λ was the inverse covariance matrix, p(μ, Λ) together follow a Gaussian-Wishart prior (μ|m0, (β0Λ)−1)W(Λ|W0,ν0), where the hyperparameters {m0, β0, W0,ν0} = {μ, 1, I, D + 1}.

VB-EM analysis

Let θ was equal to {z, π, μ, Λ}. The marginal log likelihood was written in terms of lower bound L(q) as well as Kullback-Leibler divergence κL(q||p):

In this formula, q(θ) denotes a proposed distribution for p(θ|x). Because ln p(x) was a constant, maximizing L(q) demonstrated minimizing KL(q||p). The optimal solution ln was the expectation of variable j w.r.t with respect to other variables, Ei/ = j [ln p(x, θ)]. Particularly, we determined that q (z, π, μ, Λ) was equal to q(z)q(π)q(μ, Λ). The expectations for the 3 terms (at log scale), namely ln q∗(z), ln q∗(π), ln q∗(μ), had the same forms as the initial distributions due to the conjugacy of the priors. However, they require evaluation of the parameters {z, π, μ, Λ}, which in turn all were dependent on the expectations of z or the posterior of interest. The interdependence of the expectations and model parameters fell naturally into an EM framework, named as VB-EM.

TargetScore calculation and prediction of potential miRNA targets

In detail, TSCS is a sequence-based score for single target site computed by TargetScan,[ and PCT is the probability of conserved targeting for single target site.[ TSCS and PCT were available from TargetScan Web site (http://www.targetscan.org/). Target score is a measure of the mean effect of all neighbors serving as the targets, which ranges from 0 to 1.[ The higher the TargetScore was, the greater the accuracy in identifying known targets was. Thus, based on logFC, TSCS, and PCT values, we counted the TargetScore for the genes to further extract the potential targets for miRNA-203 and miRNA-146a. In detail, the TargetScore was defined as an integrative probabilistic score of a gene being the targets of an miRNA. The TargetScore was calculated using the following formula:

where p(s|X) was posterior. The distribution of TargetScore for validated and nonvalidated targets of all the miRNA-mRNA interactions each owning at least 1 validated targets was analyzed. We defined the predefined κ as the cutoff criteria when little overlap between the 2 distribution at TargetScore = κ existed.

Pathway analysis for the predicted targets of miRNA-203 and miRNA-146a

Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/pathway) is a knowledge base to systematically study gene functions, and KEGG provides a reference knowledge base for understanding cellular processes via the process of pathway aligning, which is to map genes to KEGG reference pathways to deduce systemic behaviors of the cell.[ In our study, to reveal the functions of miRNA-203 and miRNA-146a, we implemented pathway enrichment analyses for the targets of miRNA-203 and miRNA-146a. Specifically, all KEGG reference pathways were derived from the KEGG database. Next, the targets of miRNA-203 and miRNA-146a obtained in our work were mapped to the KEGG reference pathways to extract the potential pathways which were enriched by the targets of miRNA-203 and miRNA-146a. The pathways enriched by the targets of miRNA-203 and miRNA-146a were considered as the biological functions of miRNA-203 and miRNA-146a.

Results

Calculation of log FC values

After quality control, normalization, and combining multiple probes, a total of 20,514 genes were remained for further analysis. A volcano plot exhibiting the log FC distribution of the 20,514 analyzed genes was generated, as shown in Figure 1. From Figure 1, we found that most genes distributed between −0.5 and 0.5.

Figure 1

Volcano plot exhibiting log false change (log FC) values of 20,514 genes in 44 psoriasis patients at baseline and 30 healthy controls. x axis was the log FC values of 20,514 genes. y axis represented the log transformed odds.

Identification of the potential targets for miRNA-203 and miRNA-146a

We calculated TargetScores by combining logFC and sequence scores obtained from TSCS and PCT, and derived the distribution of TargetScores. Based on the distribution of TargetScore, we observed little overlap between the 2 distribution at TargetScore = 0.4. Accordingly, we selected TargetScore >0.4 as the threshold to filter out less confidence targets. Based on the TargetScore value >0.4, a total of 49 target genes for miRNA-203 were identified, and specific information is shown in Table 1. Among these interactions, 2 targets (SOCS3, and BCL2L2) were verified, and the remaining 47 were predicted with TargetScore >0.4. Moreover, there were 17 target genes for miRNA-146a, as listed in Table 2. The targeted gene CCNA2 was validated, but others were predicted with TargetScore >0.4, which called for experimental verification.

Table 1

The target genes of microRNA -203 based on TargetScore >0.4.

Table 2

The target genes of microRNA -146a based on TargetScore >0.4.

Pathway enrichment analyses of the predicted targets of miRNA-203 and miRNA-146a

Pathway enrichment results showed that the target genes (KIR2DL1, HLA-DQA1, KIR3DL1) of miRNA-203 only participated in the pathway of antigen processing and presentation. Moreover, the target genes (ADORA3, CYSLTR2, HRH4) of miRNA-146a were only involved in the pathway of neuroactive ligand-receptor interaction.

Discussion

MiRNAs, as a class of noncoding RNA molecules, have significant power to regulate biological processes.[ Given these, miRNAs offer a clue to elaborate the complex mechanisms of the inflammation disease, including psoriasis. In our analysis, to enhance the understanding of functions of miRNA, we proposed a TargetScore, a Bayesian probabilistic scoring method taking into account of the fold-change due to miRNA overexpression and sequence-based information, to identify the potential targets of miRNA203 and miRNA-146a. Significantly, based on the TargetScore value >0.4, a total of 49 target genes for miRNA-203 and 17 target genes for miRNA-146a were identified. Interestingly, 2 targets of miRNA-203 and 1 targets of miRNA-146a were validated; in contrast, 47 targets of miRNA-203 and 16 targets of miRNA-146a were predicted with high TargetScore but needed experimental verification. Pathway enrichment results showed that the target genes of miRNA-203 only participated in the pathway of antigen processing and presentation. Moreover, the target genes of miRNA-146a were only involved in the pathway of neuroactive ligand-receptor interaction.

Antigen presentation is controlled by major histocompatibility complex (MHC) class II molecules.[ Moreover, MHC class II molecules have been suggested to be coded by human leukocyte antigen (HLA)-DR, HLA-DQ, and HLA-DP.[ Furthermore, HLA-DQ is the subclass of presenting HLA, and HLA genes have been demonstrated to be susceptibility genes of psoriasis via triggering CD8+ T cell–specific responses.[ To our knowledge, antigen processing and presentation is a fundamental part of the autoimmune response, and psoriasis is a common, chronic, inflammatory skin diseases. Significantly, antigen presentation has been indicated to be related with psoriasis.[ KIR2DL1 and KIR3DL1 were 2 members of killer-cell immunoglobulin-like receptors (KIR) which interacted with HLA molecules and regulated natural killer (NK) cell activity.[ NK cells are large granular lymphocytes, which exert important functions in the innate immune response.[ Moreover, genetic associations between the KIR genes and psoriasis have been reported.[

In our study, 3 target genes (KIR2DL1, HLA-DQA1, and KIR3DL1) of miRNA-203a were significantly enriched in the pathway of antigen processing and presentation. Accordingly, miRNA-203 might exert crucial functions in psoriasis progression, probably by regulating the antigen processing and presentation pathway.

Significantly, the target genes (ADORA3, CYSLTR2, and HRH4) of miRNA-146a were involved in the pathway of neuroactive ligand-receptor interaction. Neuroactive steroids regulate the neurotransmitter receptors to further mediate neuronal activity.[ The effect of neuroactive steroid reveals a ligand-receptor interaction. Remarkably, neuroactive steroid influences the modulation of gamma-aminobutyric acid (GABA) receptor.[ Of note, former studies have indicated that GABA receptors control cell proliferation.[ Psoriasis lesions are typically characterized by hyperproliferation of keratinocyte.[ So far, neuroactive ligand-receptor interaction is not directly involved in psoriasis. Demonstrated here, we speculate that miRNA-146a might play a crucial role in psoriasis, partially by regulating the neuroactive ligand-receptor interaction pathway.

Among the 3 target genes enriched in the pathway of neuroactive ligand-receptor interaction, HRH4 owned the higher TargetScore (value = 0.556848). HRH4, as a recently observed histamine receptor, plays a role in immunological and inflammatory processes.[ Moreover, growing evidence has demonstrated that HRH4 exerts important functions in cell proliferation.[ Moreover, another study has indicated that HRH4 affects the migration of plasmacytoid dendritic cells (pDCs) and cytokine production.[ Furthermore, pDC have been shown to infiltrate lesions of inflammatory skin diseases, such as psoriasis.[ Gschwandtner et al[ also have suggested that HRH4 alone might be a promising therapeutic target in psoriasis. Accordingly, we supposed that miRNA-146a exerted significant functions in psoriasis, partially through regulating the expression of HRH4 to a certain degree.

However, some limitations must be taken into consideration in our analysis. Firstly, sample size was small. Secondly, the dataset used in our work were downloaded from the EMBL-EBI database, not produced by ourselves. Thus, there is a need to conduct a meta-analysis for the relevant datasets of psoriasis in the future. In addition, the results obtained from our work were bioinformatics-based prediction, further experiments need to be carried out to validate our results.

In a nutshell, the advantage of this proposed method emphasizes the importance of integrating individual informative predictors into a unified framework. Based on the results, we suggested that miRNA203 and miRNA-146a played important roles in psoriasis progression, partially through regulating the pathway antigen processing and presentation pathway, and neuroactive ligand-receptor interaction pathway, respectively. However, the validation relying on other datasets is conducted to demonstrate that these functions are useful in classifying psoriasis and normal sample.

Author contributions

Data curation: Hai-Xia Chen, Yan-Song Liu.

Methodology: Xiu-Jie Zhang.

Software: Xiu-Jie Zhang.

Writing – original draft: Hai-Xia Chen.

Writing – review and editing: Yan-Song Liu.