BACKGROUND
control of gene expression is a highly complex process involving the coordinated activity of multiple and heterogeneous biological factors.
an underlying and intriguing general phenomenon is that biological molecules may act in a variety of different combinations to modulate cellular activities and to specifically react to changes in the biological milieu.
to provide coordinated and multiple complex responses, several mechanisms are known to integrate a number of molecules in a combinatorial way.
combinatorial post-translational modifications, such as methylation, acetylation, phosphorylation, and/or variations in regulatory trans-factors amounts, can influence the global regulation of gene expression at different levels.
for example, the combinatorial epigenetic tagging of dna and histones may enforce or reverse chromatin remodeling, thus playing a fundamental role in a variety of physiological and diseased cellular states  <cit> .
such combinatorial patterns of epigenomic modifications were found to be predictive of mrna and ncrna gene expression changes  <cit> .
it is also well known that combinatorial control is important for transcription, where the interaction of transcription factors is critical for gene regulation  <cit> .
the study of transcription factors combinatorics has been approached both on genome-wide  <cit>  and element-specific scales  <cit> .
as expected, the combinatorial arrangement of transcriptional regulation highly increases the overall possibilities to fine-tune the cellular response under different conditions.
out of the  <dig>  genes encoded in the human genome, about 70% produce transcripts that are alternatively spliced.
this yields multiple protein isoforms for each single pre-mrna, increasing the variability of the proteome in eukaryotes.
cis-acting rna sequence elements and enhancer complexes on pre-mrna splicing are thought to form a combinatorial control network that allows exon recognition and splicing to occur  <cit> .
evidence of splicing combinatorial mechanisms can be obtained from the parallel analysis of several genes, but a complete picture of the combinatorial rules underlying the control of splicing is still lacking  <cit> .
this is also true for the post-transcriptional control of gene expression, which is exerted by both cis-acting elements on the target mrna and trans-factors, such as rna binding proteins  and non-coding rnas .
the combined effect of trans-factors on mrnas has been hypothesized to organize the so-called “post-transcriptional rna operons ”  <cit> , but the global interplay of rbps and ncrnas on the same set of transcripts remains largely unexplored, despite being of paramount interest.
rbps play an important role in all the processing stages of rna fate, from synthesis to degradation.
an essential step for functionally understanding rbps is to identify their rna substrates and the sites at which the interactions take place.
the recent development of cross-linking and immunoprecipitation  coupled to rna-seq and related techniques has made it now possible to identify direct protein-rna interactions in vivo at a very high base resolution  <cit> .
these techniques provide positional information about the binding sites along the rna sequence  <cit> .
the combination of clip and rna-seq provides an unparalleled capability to identify transcriptome-wide protein-rna interactions.
modifications of the original method are also sprouting: the recent developments include iclip  <cit> , par-clip, incorporating photoactivable bases in rna  <cit> , and crac, an affinity-tag protocol  <cit> .
in  <dig>  a technique called global par-clip  was introduced, which allows the whole mrna-bound proteome and its global occupancy profile to be identified  <cit> .
the great progress in mapping protein-rna interactions using genome-wide tools represents a fundamental source of information for post-transcriptional regulation of gene expression, but it does not specifically address possible combinatorial patterns embedded in rna-protein and rna-rna interactions.
even if still largely incomplete, this binding site information creates for the first time a volume of data large enough to use for investigating possible combinatorial patterns of interactions by using machine learning techniques.
consequently, we can now start exploring post-transcriptional combinatorial rules in a systematic way.
to date, several tools have been developed to investigate and predict the interactions of transcription factors, mrnas and mirnas.
many of these bioinformatics approaches focus on predicting transcriptional networks by: i) modeling the expression level of a gene in terms of the predicted transcription factors that control its transcription rate  <cit> ; ii) identifying clusters of co-regulated genes  <cit> ; or, more generally iii) inferring portions of regulatory networks .
developing automated approaches to identify rules of combinatorial regulation at the post-transcriptional level would be of paramount interest to biologists.
a few attempts have focused on analyzing mirna-mediated interactions, by identifying putative feed-forward loops , in which a transcription factor regulates a mirna, and they together regulate a set of target genes  <cit> , or mirna-transcription factor motifs  <cit> .
this relies on in silico target predictions coupled to gene expression data  <cit> .
at a more general level, krek et al.
<cit>  developed pictar, a combinatorial approach to predict the binding affinity of a pre-specified set of candidate mirnas on a target mrna by combining the output of individual mirna target predictors.
coronello and benos  <cit>  refined the combinatorial model by integrating mirna expression levels in the scoring function.
albeit limited to mirna-mrna interactions, these methods can be seen as initial attempts to uncover the combinatorial nature of post-transcriptional regulations at a genome-wide level.
however, both methods require pre-specifying the set of mirna to be checked, thus preventing their general applicability for mining novel unknown combinatorial patterns.
the mining phase in fact requires an efficient search procedure to explore the space of possible combinatorial interactions, as a simple exhaustive enumeration of all combinations is computationally infeasible for all but the smallest sets of regulators.
given the assumption that functionally related genes are more likely to be co-expressed, transcriptome data have been used to derive potential gene networks.
using a similar approach, inference of clusters of co-expressed genes and potential regulatory programs in post-translational controls have been developed.
building on the work by segal et al.
on learning module networks  <cit> , joshi et al.
<cit>  developed a probabilistic approach to infer module networks in yeast, using both transcriptional and post-transcriptional regulators.
the results obtained in this work generated interesting sets of hypotheses for regulatory pathways or processes in specific biological conditions .
however, the method requires specific translational profiling time series as input.
here we propose a method, called ptrcombiner, to study the combinatorial nature of post-transcriptional trans-factors.
the method takes as input a collection of binding interactions and returns groups of factors sharing a conspicuous number of mrna targets.
it also analyzes the factors belonging to the same cluster in term of structural similarity of their binding sites.
the identification of clusters is cast as a boolean matrix factorization problem over the interaction matrix between trans-factors and mrnas.
this allows the simultaneous identification of multiple and possibly overlapping groups of factors that jointly cover as many interactions as possible.
while boolean matrix factorization has been employed in the data mining community to identify pattern sets, its application to the bioinformatics domain and especially to interaction data analysis is completely novel.
analyzing the trans-factors binding compatibility is cast as a binary classification problem aimed at discriminating between pairs of trans-factors in terms of their binding site similarity.
the classifier employs state-of-the-art graph kernels  <cit>  that account for predicted rna secondary structure in addition to sequence information.
thus, despite the still incomplete and noisy map of rna-protein interactions, this python-based tool will prove valuable in elucidating complex post-transcriptional networks.
RESULTS
ptrcombiner  is a new tool composed of multiple modules that infer meaningful combinatorial relationships between mrnas and their regulatory elements , namely rbps and ncrnas.
this goal is achieved by extracting combinatorial information through a pattern-set mining approach and a meta-analysis of genome-wide data.
ptrcombiner is divided into two activity components.
the first, “mining combinatorial features” , represents interaction data with a mathematical model.
the model involves an approximate boolean matrix factorization  of the interaction matrix, which identifies groups of regulatory elements acting on common mrna utrs, which we call clusters.
the second, “analyzing combinatorial features” , evaluates the biological characteristics of the clusters identified by the pattern-set miner.
each cluster is evaluated in terms of global biological meaning by gene ontology  analysis over its mrna targets and the binding site compatibility between the individual trans-factors in the cluster.
description of the dataset
the tool has been applied to the set of post-transcriptional interactions in human contained in the atlas of utr regulatory activity  <dig>  database .
aura  <dig> is a manually curated and comprehensive catalog of mrna untranslated regions  and their regulatory annotations, including interactions with trans-factors .
to date, aura  <dig> is the largest dataset of utr-centered regulatory annotations taken from the whole range of existing experimental techniques, such as clip, rip, selex, and rnacompete.
a subset of these techniques, namely clip and its variants followed by high-throughput sequencing, allows the positional annotation of the binding sites along the utrs, while the others  can only detect the presence of an interaction between a transcript and a trans-factor without the positional information of the specific binding site.
as a working example for ptrcombiner, we considered the whole set of human interactions annotated in aura  <dig>  the available number of utrs bound by each trans-factor ranges from  <dig> to  <dig> , with median of  <dig> and a mean of  <dig> .
symmetrically, the available number of trans-factors bound to the same utr ranges from  <dig> to  <dig>  with median  <dig> and mean  <dig> .
using these data, we built a boolean matrix c with  <dig>  rows  and  <dig> columns , where cij =  <dig> if the jth trans-factor interacts with the ith utr, and cij =  <dig> otherwise .
the annotated interactions are collectively  <dig>  , thus the sparsity of the interaction matrix is  <dig> .
mining combinatorial features
after obtaining the interaction matrix, the first step of the analysis was to identify clusters of trans-factors  that bind the same set of utrs.
each of these clusters could be a candidate combinatorial member of a post-transcriptional regulatory code.
this step thus aims to identify multiple overlapping clusters, which collectively cover most of the known interactions between trans-factors and utrs.
boolean matrix factorization  <cit>  provides this requirement by decomposing the matrix of known interactions  in the product of two boolean matrices.
one of the matrices represents the clusters of the trans-factors, while the other, the utrs in terms of their interactions within the clusters.
the algorithm takes two arguments: the number of clusters to return , and a threshold .
the τ value controls the minimal amount of shared utrs inside a cluster.
the higher the threshold, the more target utrs should be shared among trans-factors in order for these to be considered as a cluster.
the algorithm returns a list of clusters ordered by coverage, i.e.
the number of interactions associated to each cluster .
to analyze the behavior of the algorithm when varying its parameters, we produced a surface parameter graph  where the average number of trans-factors belonging to each cluster  was calculated using various combinations of k and τ values .
given a certain τ, different values of k do not affect the average cluster size, which appears to be mainly τ dependent.
we chose the τ value giving an average cluster size equal to the average number of trans-factors bound to a single utr .
we thus considered a τ of  <dig> , resulting in clusters composed of an average of six trans-factors.
interestingly, this value points to a stable region of the k and τ surface, in which the number of trans-factors for each cluster does not change drastically in the surrounding area .
table  <dig> shows the clusters obtained with the selected threshold.
the first nine clusters are composed exclusively of rbps, as well as the clusters r <dig> to r <dig>  r <dig> and r <dig>  the first cluster displaying co-occurrence of rbps and mirnas was r <dig>  followed by clusters r <dig>  r <dig>  r <dig> and r <dig>  no cluster uniquely composed of mirnas was present.
we would like to stress that  <dig> out of  <dig> clusters do not represent real combinations, as they comprise only one trans-factor.
we refer to these one-element clusters as “singletons”.
recalling that the algorithm aims to cover as many interactions as possible with the set of clusters, a singleton can be extracted by the algorithm whenever a trans-factor has many interactions that are not shared with any other trans-factor for which experimental interaction data are available.
clusters were classified according to their composition: “rbp” for those composed exclusively of rbps; “mirna”, composed exclusively of mirnas; “rbp-mirna”, composed of both rbps and mirnas; and “singleton”, composed of only one trans-factor.
as detailed in the previous section, different trans-factors have highly different numbers of annotated interacting utrs.
being driven by coverage, the algorithm is inherently biased towards clusters of “widely interacting” trans-factors.
therefore, when we analyzed the composition of the clusters, we observed that some trans-factors are present in multiple clusters.
for example, the argonaute proteins ago <dig> and ago <dig>  and the well-known rbp elavl1/hur, occur in  <dig>   <dig> and  <dig> out of  <dig> clusters, respectively.
ago <dig> and ago <dig> are components of the rna-induced silencing complex , the protein complex which is responsible for down-regulating mrnas  <cit> .
these proteins bind different classes of small ncrnas, such as mirnas and small interfering rnas , leading the argonaute proteins to their specific targets through sequence complementarity, thus silencing their targets.
therefore, it is not surprising to find ago <dig> and ago <dig> in almost all of the clusters, given the widespread activity of these proteins.
moreover, in accordance with our results, ago <dig> and ago <dig> have been found to interact with elavl1/hur  <cit> .
despite this interaction, the two proteins avoid any binding overlap on the target mrnas: ago proteins preferentially bind the boundaries of utrs, while elavl1/hur binds uniformly along utrs and disappears toward the stop codon and the polyadenylation site  <cit> .
elavl1/hur, a member of the elav family, is known to be broadly expressed in tissues and to bind au-rich elements in the 3’ utrs of thousands of mrnas  <cit> .
moreover, it has been demonstrated that elavl1/hur displays competitive and cooperative interactions with mirnas/risc  <cit> , that these interactions may depend on the proximity between the protein and mirna binding sites  <cit>  and that it is part of a complex mrna network for coordinating gene expression  <cit> .
these results are in agreement with the cluster composition in table  <dig> and support our finding that the elements displaying the highest number of interactions  are those that most frequently occur across the clusters.
for this reason, we called these elements “recurrent trans-factors” and these clusters “rk” .
to explore trans-factors that have a narrower spectra of interactions and thus occur less frequently in the clusters, we removed all recurrent trans-factors  and ran a second iteration of the algorithm.
this second iteration focused on trans-factors that appeared in none or only one of the clusters of the previous analysis.
we termed these trans-factors “sporadic”.
again, we studied the behavior of the algorithm when the k and τ parameters were varied, in order to select the most appropriate combination of values .
in this case, the best choice of τ was  <dig>  .
this value returned clusters comprised of an average number of three trans-factors, which corresponds to the average number of sporadic elements bound to each utr.
we called these clusters “sk”  .
it is clear that the majority of the clusters  are singletons.
in contrast to results obtained when recurrent factors were included, we observed that four clusters were formed exclusively by mirnas .
interestingly, pum <dig>  which was also found as a member of the recurrent clusters, formed two distinct clusters with different sets of mirnas .
in line with this, some evidence has suggested that pum <dig> associates with mirnas.
pum <dig> is known to act as a translational repressor in several organisms, contributing to dendritic rna localization and silencing  <cit>  and regulating synaptic formation  <cit> .
moreover, an extensive interaction between pum <dig> and pum <dig> with the mirna regulatory system has been suggested  <cit> , indicating that interactions between the rbps and the mirnas in translational regulation may be more frequent than previously thought.
in addition, a recent computational analysis suggested that specific groups of mirna binding sites localize within  <dig> nt from pum <dig> binding sites, giving support to a possible cooperativity between pum <dig> and mirna in mrna degradation  <cit>  and to the biological meaning of our clusters where pum <dig> acts on mrnas in combination with different mirnas.
in particular, hsa-mir- <dig> and hsa-mir- <dig>  which are part of cluster s <dig> together with pum <dig>  appear to colocalize with this rbp  <cit> .
clusters were classified according to their composition, as for table  <dig>
despite that only a small fraction of known trans-factors occur in the clusters , the vast majority of the existing interactions is covered by the identified clusters .
the fact that the most of the trans-factors are not part of any cluster is due to the paucity of available information, i.e.
the overall number of experimental data obtained until now for many of the trans-factors is still far from being complete.
this is not surprising given the novelty of the experimental techniques involved.
it is clear that more data are needed to obtain a reliable and exhaustive description of all the possible combinatorial interactions in human.
the availability in the next future of high-throughput experimental data covering an increasingly larger amount of trans-factors will open up new possibilities to uncover new and more reliable clusters and to obtain new combinatorial information.
analyzing combinatorial features - step 1: biological annotation
the first step in analyzing the clusters obtained from the mining module was to measure how much they differ in terms of: a) target mrnas ; b) enriched ontological terms ; and c) similarity among ontological terms .
overall, we expected to obtain a large amount of regulated genes that form clusters in the presence of recurrent trans-factors from the mining module, because they display the highest number of annotated interactions .
in fact, the module target mrnas overlap revealed that several hundred genes were co-regulated by recurrent trans-factors .
the average number of hgnc genes regulated by the first five clusters  was  <dig> , ranging from  <dig> for cluster r <dig> to  <dig>  for cluster r <dig>  the number of target genes was markedly reduced in the clusters of sporadic trans-factors, as expected given their greater specificity .
the average number of hgnc genes regulated by the first five clusters was  <dig>  ranging from  <dig> for cluster s <dig> to  <dig> for cluster s <dig>
we then explored the possibility of finding the same mrna targets in different clusters, by calculating the overlap among the populations of genes grouped in different clusters using the jaccard similarity .
for the clusters including recurrent trans-factors, the overlap was 21% on average .
this result suggests that the method is able to group sets of specific target mrnas even if they share common trans-factors.
however, some leakage is present: for example, cluster r <dig> and cluster r <dig> share 43% of their targets.
this is not surprising as careful inspection of the elements forming these clusters revealed that r <dig> and r <dig> share seven rbps and are individually characterized only by fxr <dig> and fus, respectively.
the same is observed for clusters r <dig> and r <dig>  which share 60% of their targets.
in this case, cluster r <dig> shares with r <dig> almost all its own trans-factors .
focusing on clusters of sporadic elements revealed that the average overlap decreased to 7%, one-third of that from the previous analysis .
in this case, the maximum overlap was between clusters s <dig> and s <dig>  which share 27% of their targets.
both of these clusters comprise caprin <dig> and mov <dig>  while they do not share fus, fxr <dig>  taf <dig> and c22orf <dig> .
these results confirm the effectiveness of the repeated run of analysis in identifying distinct, small-sized sets of genes, uniquely regulated by specific sets of sporadic trans-factors.
the gene ontology enrichment module of the tool allowed us to study the biological relevance of the mined clusters.
the module attempts to identify parts of common and biologically coordinated mechanisms or processes that govern coherent cellular outcomes and that could be characterized by the identified groups of preferential interactions.
in this module, the analysis is expanded from single genes to more general biological annotations, allowing the inferred clusters of trans-factors to be compared by the gene ontology  enrichment analysis of their target mrnas.
an initial and effective way to compare enrichments is to display the top enriched go terms for each cluster.
such a comparison is shown in figure 4e for the top five non-singleton rk clusters, using the molecular function  branch of go as example .
the modular blocks of enriched terms scattered along the columns of the heatmap clearly indicate that the clusters have marked and distinct molecular function consensuses.
similarly to what we observed for the mining module results, the modularity of the enrichments was further reinforced in the top five non singleton sk clusters .
here, the occurrence of terms enriched in multiple clusters is an infrequent event.
clusters s <dig> and s <dig> share the most similar enrichment signature, mirroring the strong similarity observed between the two clusters .
after comparing the enrichment of different clusters, we next assessed how the ontological enrichment of genes regulated by one cluster differs from the ontological enrichments of genes regulated by the individual trans-factors forming the cluster.
this intra-cluster comparison allowed us to potentially identify “emergent enrichments”, i.e.
to detect go terms enriched exclusively in a set of genes regulated by a set of trans-factors forming a cluster.
an example of this analysis is shown for cluster s <dig> .
the target genes in this clusters were specifically enriched for the biological process  term “cell division”, the cellular component  term “nuclear speck”, and the molecular function  term “transcription corepressor activity”.
these results strongly suggest that the clusters have emergent and specific combinatorial properties, as expected for combinatorial mechanisms.
to further estimate the similarity between the ontological enrichments associated with each cluster, the tool is able to calculate the semantic similarity values, which account for semantic similarity relationships between non-identical ontological terms .
the module gene ontology similarity calculates the semantic similarity values for all the three branches of go and between both each pair of clusters  and the single trans-factors belonging to the same cluster .
in figure  <dig>  the top rows of each panel list the semantic similarity values between enriched terms associated to single trans-factors and enriched terms associated to cluster s <dig>  here, the stronger semantic similarities are observed among the cc enrichments, while the weakest values lie among the bp enrichments.
this analysis can also help to rank the ontological distance between a cluster and its trans-factors according to the semantic similarity values.
for example, fxr <dig> globally shows the highest similarity with the enrichment of cluster s <dig> .
analyzing combinatorial features - step 2: classifying rbp-binding sites
to suggest some potential binding mechanisms based on the available experimental data  <cit> , the second step of the analysis focuses on the rbps forming the previously identified clusters.
the underlying idea is that whenever two rbps  are characterized by similar binding sites over the mrna, then a concurrent binding, either competitive or cooperative, has occurred.
although experimental methods  can be used to gather information about the proximity of binding sites for specific rbps, the resulting information is corrupted by several types of noise sources: 1) a considerable fraction of binding sites can remain undetected  because the methods are applied on cells of a particular type and in specific growth conditions, in which not all the potential bound mrnas are expressed; 2) interactions that are transient can be mistakenly identified as stable ; and 3) post-processing analysis, such as mapping and peak detection, can increase the number of false negatives due to the difficulty of dealing with splice sites and the stringent thresholds needed for a confident detection, respectively.
computational models for rbp target detection are therefore valuable tools in dealing with the low signal-to-noise ratio of current experimental techniques.
such models can significantly increase the precision with which target sites are resolved and can uncover sites that would be otherwise missed by experimental protocols.
to determine if two rbps are are likely to interact in a cooperative or competitive fashion, we first compiled an equivalent in silico model of the preferred rbp target sites.
given a cluster of rbps and their binding sites, the module trains a machine learning algorithm to discriminate between binding sites of two different rbps, for all possible pairwise combinations of proteins belonging to the same cluster.
two proteins are likely to have different binding sites when the algorithm effectively distinguishes between their binding areas.
in contrast, compatible binding sites are more likely to lead to a difficult discrimination task.
discrimination is based on a kernel machine binary classifier  <cit>  capable of computing similarity between base sequences in terms of their respective putative secondary structures  <cit>  and thus of the spatial conformation of their binding sites .
the structural component of this discrimination has a strong biological significance, since rna interactions are not exclusively driven by sequence specificities.
we report the results of the classification analysis performed on the first two clusters of sporadic trans-factors , namely, the cluster s <dig> that comprises caprin <dig>  fus, fxr <dig>  mov <dig>  and taf <dig>  and the cluster s <dig> that comprises c22orf <dig>  caprin <dig>  and mov <dig> .
for each rbp, we randomly selected  <dig>  mrna utr stretches  from available binding coordinates of large-scale experiments  stored in the aura  <dig> database .
the classification performances for clusters s <dig> and s <dig> are displayed in figure 6a and figure 6b, respectively.
performance is computed with the aurocc and f1-score measures.
aurocc is an aggregate measure that evaluates the quality of a classifier when varying the threshold, to decide when a prediction should be considered positive.
an aurocc value of  <dig>  corresponds to a random predictor, while an aurocc value of  <dig> indicates perfect discrimination.
the f <dig> score is the harmonic mean between precision and sensitivity, trading off the two complementary measures .
considering the cluster s <dig>  the classifier discriminates the binding sites of only a subset of the rbps in the cluster .
the most specific binding sites were observed for caprin <dig>  and mov <dig> .
fus and taf <dig> seem to have more compatible binding sites.
an aurocc of  <dig>  is in fact very close to the performance of a random predictor, suggesting that these proteins share a similar if not identical set of binding sites.
indeed, these two proteins are known paralogues, both belonging to the fet family of rna-binding proteins  <cit> .
higher aurocc values were achieved when distinguishing the binding sites of fxr <dig> from those of fus and taf <dig> , and maximal aurocc values were found for mov <dig>  suggesting different binding sites.
classification performances of cluster s <dig>  displayed in figure 6b, were generally high, suggesting that the utr stretches that are bound by the rbps forming the cluster  are different.
figures 6c and 6d show the distribution of the distances between couples of binding sites lying on the same mrna and targeted by distinct rbps, for cluster s <dig> and s <dig> respectively.
the distances between the binding sites of fus and taf <dig> are much lower than those of the other distributions, indicating that the two proteins tend to bind to the same or very close regions.
the average distance between fxr <dig> and fus or taf <dig> is also low, but it has a much larger spread.
a large average distance is observed for the other cases and indicates a high discrimination capacity of the classifier.
these examples show how we can start to use in silico modeling of rbps interactions to investigate their collaborative and/or competitive effects.
this type of modeling is effective when the experimental data is affected by noise, since it recovers missed interactions  and filters out accidental interactions .
a related predictive system was used in  <cit> , where it was shown that a model trained on a set of ago <dig> hits-clip sites could effectively identify binding sites missed by the high-throughput experiment.
these findings have been verified comparing the predicted ago <dig> targets with changes in transcript expression levels upon ago <dig> knockdown.
conversely this type of trend was not observed for the original hits-clip-detected sites indicating a significant impact of false negatives in the experimental setting.
the upgrade to predictive models capable of this level of accuracy will allow more sophisticated investigations, such as those warranted by the combinatorial analysis presented in this work.
for instance, a competitive effect can be hypothesized if two rbps, found in the same cluster, exhibit a compatible preference for the same target region, even when the experimental data are incomplete and do not report overlapping interaction areas.
conversely, if the target regions are predicted to be sufficiently close but not overlapping, a cooperative effect can be hypothesized, even when experimental data cannot resolve the distinct areas and these are therefore reported as overlapping.
alternative workflows
<dig> - filtering interactions by experimental technique
ptrcombiner allows a selection of the source interaction data to be use, thereby taking into account the different experimental approaches currently available for mapping protein-rna interactions  <cit>  and producing technically homogeneous results.
this method-selection feature of ptrcombiner can be employed to exclude from the analysis the interactions obtained from techniques that are considered more unreliable or inopportune for specific analyses.
as an example of the method-selection tool, we filtered the human interactions annotated in aura  <dig> according to the evidence type, identifying three subsets of the original dataset.
the first subset considered only par-clip experiments and includes a total of  <dig> trans-factors  and  <dig>  utrs, with an average number of interactions per utr of  <dig>  .
the second subset considered all the other clip experiments, namely clip, clip-seq, hits-clip and iclip.
it includes a total of  <dig> trans-factors  and  <dig>  utrs, with an average number of interactions per utr of  <dig>  .
the third subset used only interactions found by rna immunoprecipitation rip.
to date, it includes a total of  <dig> trans-factors  and  <dig>  utrs, with an average number of interactions per utr of  <dig>  .
in the first subset , we observed the presence of clusters composed exclusively of rbp .
the clusters are similar to those observed when the whole set of interactions was analyzed, excluding trans-factors with none or few interactions based on par-clip .
considering only the rbp interactions obtained with par-clip experiments, the similarity with the clusters analyzed and described in the previous paragraphs is  <dig>  .
after filtering for the clip, clip-seq, hits-clip, and iclip  interaction data, too few rbps  were left in the dataset to extract a reasonable number of clusters, so that no comparative analysis with the clusters obtained considering the whole set of interactions could be performed.
finally, the clusters extracted from the third dataset  had, on average, very few associated genes .
we expect that this lack of large scale interaction data will be overcome once more datasets are made available.
overall, these results suggests that, to date, the main component driving the selection of clusters is the set of interactions measured by par-clip experiments.
<dig> - filtering targets and regulators by expression level
so far, we have used ptrcombiner with a dataset that comprised all the binding evidence from multiple experiments in different biological systems that shared only the “human” source.
the biological interpretation of results would greatly benefit from combining interaction data with other sources of information, such as expression data when available, in order to generate more specific hypotheses.
to integrate ptrcombiner with expression data, we introduced the possibility to filter interactions according to: a) the expression levels of mrna targets; and b) the expression levels of the regulators.
in this way, the interactions can be filtered according to the specificity of the biological system being studied.
to show this functionality, we chose a well-studied cell line  and filtered the aura dataset by selecting: a) the transcripts known to be highly expressed in hela ; and b) rna binding proteins and mirna known to be present in hela .
the filtered dataset now comprises  <dig>  utrs ,  <dig> rbps , and  <dig> mirnas .
we ran ptrcombiner on this filtered dataset, with the optimal τ parameter corresponding to  <dig> .
the resulting clusters are displayed in additional file 9: table s <dig>  the majority of the clusters are composed of rbps, with the exception of six singleton clusters and three mixed clusters of mirnas and rbps.
to analyze the effect of the cell-specific filter on the resulting clusters, we compared these clusters  with those obtained in the first analysis that lacked any gene expression-abased filter .
the comparison is shown in additional file 9: table s <dig>  where each hela-cluster is matched to the general-cluster with maximal jaccard similarity.
the average jaccard similarity was  <dig>  and ranged from nearly identical clusters  to clusters sharing half of their members .
interestingly, one-third of the general clusters  could not be associated with any hela cluster.
this is not surprising since these clusters are mainly composed of trans-factors or corresponding targets that are not expressed in hela cells.
hela clusters therefore represent a subset of the general clusters.
<dig> - balancing the trans-factor sample size
an additional alternative workflow supported by ptrcombiner regards the mining step of the algorithm and arises from the previously discussed bias of the algorithm towards “widely interacting” trans-factors .
introducing a balanced trans-factor association score is a possible way to reduce this bias.
the difference between the standard version of ptrcombiner  and the proposed alternative workflow  is the way in which the pool of possible clusters is selected.
the standard version of the algorithm incrementally extracts clusters of trans-factors by taking each candidate trans-factor as a “seed” and computing its association score with other trans-factors.
the association score is the number of shared targets between the two factors, normalized by the number of targets of the seed .
this implies that the seed is required to have a significant fraction of targets in common with another trans-factor for that trans-factor to be included in the seed’s cluster.
by construction, this score is asymmetric and tends to associate trans-factors with many interactions  together with those with fewer interactions  that share a significant fraction of targets with them.
here, we used a novel, alternative version of the association score normalized by the square root of the product of the target number for each trans-factor , thereby rebalancing the interaction data to cope with factors having a widely different number of targets.
note that the optimal threshold τ in this case is different from the unbalanced one , as the alternative association score strongly affects clusters sizes.
clusters obtained with the balanced score covered a larger number of trans-factors  than did the standard  version of the algorithm  .
although the average length of the clusters was the same, the balanced ptrcombiner produced more singleton clusters  and a few very large clusters.
intuitively, the larger the cluster, the fewer its associated genes .
for this reason, the number of genes associated to the non-singleton clusters in the balanced case is lower than that in the unbalanced case .
also, the jaccard similarity between clusters, calculated over their trans-factors, is lower in the balanced case .
clusters in the balanced case are thus more specific.
indeed, heat maps of the enriched go terms have less overlap with respect to the unbalanced case .
a possible drawback of the balanced workflow is a tendency to create clusters that pick elements with a similar number of interactions.
this excludes for instance clusters with mirnas and rbps that emerge from the standard workflow even after removing recurrent trans-factors .
the two procedures thus have different characteristics, allowing users to discover different types of interesting combinatorial associations and to alternatively analyze the data according to their specific needs.
comparison to related work
ptrcombiner is a novel approach mining combinatorial post-transcriptional regulation patterns from interaction data at the genome-wide level.
as mentioned in the background, previous attempts have been made in recent years to develop automated approaches to identify combinatorial aspects of post-transcriptional regulation.
several computational tools have been developed for mirna target prediction using a number of different approaches, such as targetscan  <cit> , miranda  <cit> , targetboost  <cit> , pita  <cit> , magia <dig>  <cit> , mirtar hunter  <cit> .
target site predictions for different mirnas have been combined in the tools pictar  <cit>  and comir  <cit>  to identify potential groups of interacting mirnas.
by using simple model systems such as s. cerevisiae, some published studies have attempted to develop computational and experimental methods that identify functional modules based on combinatorial  rna-protein interactions.
for example, joshi and coworkers  <cit>  developed a probabilistic method for inferring regulatory module networks from expression profiles.
in the following, we provide a more in-depth comparison with these approaches, highlighting the differences with respect to our proposed method and reporting comparative quantitative analyses on benchmarks for which these alternative approaches can be run.
<dig> - comparison to pictar
pictar is a probabilistic predictive method that computes the probability of multiple mirnas co-binding to the same target mrna by combining the binding scores for each candidate mirna.
albeit focusing on combinatorial interactions, the algorithm differs from ptrcombiner in many aspects.
first, it explores the combinatorial interaction of mirnas only.
second, it relies on predictive methods for detecting potential binding sites rather than using experimental data.
the predictive approach can be reasonably accurate for mirnas but are still far from satisfactory for rbps.
the third and main caveat of the method is that it lacks a mining procedure that allows to efficiently search the combinatorial space of possible clusters.
rather, it either requires the user to specify a set of mirna to be jointly evaluated, or it needs to try all possible combinations in order to identify the high scoring ones.
ptrcombiner takes a different perspective, implementing a mining approach to efficiently explore the combinatorial space of candidate clusters of mirna and/or rbp, guided by their coverage of observed interactions with the target mrnas.
as discussed in the conclusions, we plan to adapt the method to deal with weighted interactions, e.g., the score of a predicted binding or the confidence of a certain experimental technique.
in order to quantitatively compare the two approaches, we analyzed clusters of mirna found by ptrcombiner using pictar.
we focused on the clusters s <dig> and s <dig>  because they are composed of only four mirnas, and evaluating larger clusters with pictar is too computationally expensive.
for each cluster, we took the set of its target genes, i.e.
the genes which interact with all mirnas in the cluster, and computed the pictar interaction score of each of them with the cluster .
the score is computed as the maximum value of the product of the binding scores of the mirnas in the cluster, with the constraint that binding sites of different mirnas should not overlap.
if this constraint cannot be satisfied, or one or more mirnas do not have predicted binding sites on the target, the score was set to zero.
binding scores for mirnas were downloaded from the dorina database  <cit> .
we then compared these cluster-target scores with those obtained by running the same procedure on the entire set of genes  using the welch’s two samples t-test.
the results of the statistical tests are presented in additional file  <dig>  for both clusters, the difference between cluster-target scores and general scores was statistically significant, with a confidence of approximately 99%.
the fact that the average on the former set is significantly higher than the average on the latter set is not a trivial information.
in fact, it indicates that also pictar estimated the clusters of trans-factors to be relevant exactly to the specific set of genes used to cluster them together, confirming the relevance of the clusters mined by ptrcombiner.
<dig> - comparison to comir
comir  <cit>  is a web tool for combinatorial mirna target prediction.
it uses mirna expression levels in combination with thermodynamic modeling and machine learning techniques to make combinatorial predictions.
comir sums up the weighted  scores of the single mirnas, computed according to the four different scoring schemes of miranda  <cit> , pita  <cit> , targetscan  <cit> , and mirsvr  <cit> .
these scores are combined through a support vector machine , which outputs the likelihood that the set of mirnas targets a specific gene.
as for pictar, the main difference with ptrcombiner is the lack of a mining procedure that would identify the clusters of mirna to be evaluated.
using comir, we analyzed all clusters composed only of mirnas extracted by ptrcombiner, namely, s <dig>  s <dig>  s <dig>  and s <dig> .
in computing scores, we gave no expression level to comir, resulting in uniform level for all mirnas.
as for the pictar case, we compared cluster-target scores with scores for the entire set of genes, i.e.
all genes in comir output .
welch’s two sample tests confirmed the statistical significance of the difference between cluster-target and general scores for all clusters, with a confidence of approximately 100% .
similarly to the comparison to pictar, this result confirms the relevance of the clusters mined by ptrcombiner.
<dig> - comparison to lemone
lemone  <cit>  is a probabilistic method for inferring regulatory module networks from expression profiles.
the module network model was introduced in the work by segal et al.
<cit>  as a probabilistic bi-clustering approach, extracting co-clusters of genes and conditions from a matrix of gene expression levels under different experimental conditions.
a co-cluster contains a subset of genes and conditions, such that the genes have similar expression levels under the selected conditions.
the subset of genes represents a regulatory module, while each subset of conditions for the same module is an expression state for the module.
lemone introduces an ensemble averaging strategy to generate more coherent modules from multiple runs of the module network inference process.
the approach was later adapted to infer transcriptional and post-transcriptional modules from both transcriptome and translatome expression profiles in s. cerevisiae <cit> .
as already discussed, this approach detects putative regulatory modules that characterize specific biological conditions , while ptrcombiner targets more general purpose, genome-wide combinatorial interactions.
nonetheless, it is interesting to study the relationship between clusters detected by the two methods, when both are applicable .
to compare this method with ptrcombiner, we ran our computational tool on the yeast dataset employed in  <cit> .
the interaction dataset  <cit>  contains rip-chip experiments involving  <dig> rbps and  <dig>  genes.
we obtained a list of interacting rbps and their relative number of target genes .
the distribution of the number of distinct trans-factors bound to the same gene are shown in figure 8b.
in total, the dataset contains  <dig>  annotated interactions, with the interaction matrix sparsity value of  <dig>  .
from the interaction matrix, ptrcombiner  extracted the clusters composed exclusively of rbps .
these clusters were then compared with the sets of rbps whose targets were enriched in post-transcriptional modules generated by lemone.
the list of rbp sets was extracted from the supplementary materials of  <cit>  .
ptrcombiner clusters were then compared with lemone clusters, according to the jaccard similarity maximized by matching each ptrcombiner cluster with a lemone cluster .
reassuringly, the agreement between the two methods is high, with half of the top ten ptrcombiner clusters identical to lemone clusters.
CONCLUSIONS
in this paper we present a computational tool for the combinatorial analysis of post-transcriptional regulation patterns involving multiple trans-factors.
this tool, called ptrcombiner, was tested on two sets of experimental interactions between post-transcriptional trans-factors and target mrnas, in human and yeast.
ptrcombiner enables the user to: a) detect groups of regulators that share a conspicuous amount of mrna targets; b) characterize the clusters biologically; and c) identify concurrent binding sites of trans-factors belonging to the same cluster.
further, the underlying boolean matrix factorization approach allows the user to identify multiple overlapping clusters of trans-factors that jointly account for as many interactions as possible, casting the problem into an approximate interaction coverage task.
this method naturally addresses the limitations of most clustering and motif mining approaches, which typically return non-overlapping clusters , clusters covering non-overlapping sets , or long lists of highly redundant clusters.
our tool is an original and comprehensive attempt to provide a computational pipeline for elucidating complex post-transcriptional combinatorial rules on a genome-wide level.
by integrating expression profiles of both trans-factors and target mrnas, the tool can also be used to mine combinatorial patterns in specific experimental conditions.
we plan to extend the method to deal with uncertainty in the interaction information, such as putative interactions from predictive algorithms, and to incorporate the binding strength by properly weighting the contribution of each interaction.
importantly, ptrcombiner is a versatile tool that is not limited to post-transcriptional regulation analysis.
in fact, it can be easily adapted to mine transcriptional and combined transcriptional/post-transcriptional regulation patterns.
an inherent limitation in the analyses that can be conducted with the tool is given by the paucity of interaction data available to date.
however, given the fast rate at which high-throughput interaction detection experiments are conducted, this tool may have the potential to unveil the complex mosaic of interactions underlying post-transcriptional regulation in the near future.

 control of gene expression is a highly complex process involving the coordinated activity of multiple and heterogeneous biological factors. an underlying and intriguing general phenomenon is that biological molecules may act in a variety of different combinations to modulate cellular activities and to specifically react to changes in the biological milieu. to provide coordinated and multiple complex responses, several mechanisms are known to integrate a number of molecules in a combinatorial way. combinatorial post-translational modifications, such as methylation, acetylation, phosphorylation, and/or variations in regulatory trans-factors amounts, can influence the global regulation of gene expression at different levels. for example, the combinatorial epigenetic tagging of dna and histones may enforce or reverse chromatin remodeling, thus playing a fundamental role in a variety of physiological and diseased cellular states. such combinatorial patterns of epigenomic modifications were found to be predictive of mrna and ncrna gene expression changes. it is also well known that combinatorial control is important for transcription, where the interaction of transcription factors is critical for gene regulation. the study of transcription factors combinatorics has been approached both on genome-wide and element-specific scales. as expected, the combinatorial arrangement of transcriptional regulation highly increases the overall possibilities to fine-tune the cellular response under different conditions. out of the genes encoded in the human genome, about 70% produce transcripts that are alternatively spliced. this yields multiple protein isoforms for each single pre-mrna, increasing the variability of the proteome in eukaryotes. cis-acting rna sequence elements and enhancer complexes on pre-mrna splicing are thought to form a combinatorial control network that allows exon recognition and splicing to occur. evidence of splicing combinatorial mechanisms can be obtained from the parallel analysis of several genes, but a complete picture of the combinatorial rules underlying the control of splicing is still lacking. this is also true for the post-transcriptional control of gene expression, which is exerted by both cis-acting elements on the target mrna and trans-factors, such as rna binding proteins and non-coding rnas. the combined effect of trans-factors on mrnas has been hypothesized to organize the so-called “post-transcriptional rna operons ”, but the global interplay of rbps and ncrnas on the same set of transcripts remains largely unexplored, despite being of paramount interest. rbps play an important role in all the processing stages of rna fate, from synthesis to degradation. an essential step for functionally understanding rbps is to identify their rna substrates and the sites at which the interactions take place. the recent development of cross-linking and immunoprecipitation coupled to rna-seq and related techniques has made it now possible to identify direct protein-rna interactions in vivo at a very high base resolution. these techniques provide positional information about the binding sites along the rna sequence. the combination of clip and rna-seq provides an unparalleled capability to identify transcriptome-wide protein-rna interactions. modifications of the original method are also sprouting: the recent developments include iclip, par-clip, incorporating photoactivable bases in rna, and crac, an affinity-tag protocol. in a technique called global par-clip was introduced, which allows the whole mrna-bound proteome and its global occupancy profile to be identified. the great progress in mapping protein-rna interactions using genome-wide tools represents a fundamental source of information for post-transcriptional regulation of gene expression, but it does not specifically address possible combinatorial patterns embedded in rna-protein and rna-rna interactions. even if still largely incomplete, this binding site information creates for the first time a volume of data large enough to use for investigating possible combinatorial patterns of interactions by using machine learning techniques. consequently, we can now start exploring post-transcriptional combinatorial rules in a systematic way. to date, several tools have been developed to investigate and predict the interactions of transcription factors, mrnas and mirnas. many of these bioinformatics approaches focus on predicting transcriptional networks by: i) modeling the expression level of a gene in terms of the predicted transcription factors that control its transcription rate; ii) identifying clusters of co-regulated genes; or, more generally iii) inferring portions of regulatory networks. developing automated approaches to identify rules of combinatorial regulation at the post-transcriptional level would be of paramount interest to biologists. a few attempts have focused on analyzing mirna-mediated interactions, by identifying putative feed-forward loops, in which a transcription factor regulates a mirna, and they together regulate a set of target genes, or mirna-transcription factor motifs. this relies on in silico target predictions coupled to gene expression data. at a more general level, krek et al. developed pictar, a combinatorial approach to predict the binding affinity of a pre-specified set of candidate mirnas on a target mrna by combining the output of individual mirna target predictors. coronello and benos refined the combinatorial model by integrating mirna expression levels in the scoring function. albeit limited to mirna-mrna interactions, these methods can be seen as initial attempts to uncover the combinatorial nature of post-transcriptional regulations at a genome-wide level. however, both methods require pre-specifying the set of mirna to be checked, thus preventing their general applicability for mining novel unknown combinatorial patterns. the mining phase in fact requires an efficient search procedure to explore the space of possible combinatorial interactions, as a simple exhaustive enumeration of all combinations is computationally infeasible for all but the smallest sets of regulators. given the assumption that functionally related genes are more likely to be co-expressed, transcriptome data have been used to derive potential gene networks. using a similar approach, inference of clusters of co-expressed genes and potential regulatory programs in post-translational controls have been developed. building on the work by segal et al. on learning module networks, joshi et al. developed a probabilistic approach to infer module networks in yeast, using both transcriptional and post-transcriptional regulators. the results obtained in this work generated interesting sets of hypotheses for regulatory pathways or processes in specific biological conditions. however, the method requires specific translational profiling time series as input. here we propose a method, called ptrcombiner, to study the combinatorial nature of post-transcriptional trans-factors. the method takes as input a collection of binding interactions and returns groups of factors sharing a conspicuous number of mrna targets. it also analyzes the factors belonging to the same cluster in term of structural similarity of their binding sites. the identification of clusters is cast as a boolean matrix factorization problem over the interaction matrix between trans-factors and mrnas. this allows the simultaneous identification of multiple and possibly overlapping groups of factors that jointly cover as many interactions as possible. while boolean matrix factorization has been employed in the data mining community to identify pattern sets, its application to the bioinformatics domain and especially to interaction data analysis is completely novel. analyzing the trans-factors binding compatibility is cast as a binary classification problem aimed at discriminating between pairs of trans-factors in terms of their binding site similarity. the classifier employs state-of-the-art graph kernels that account for predicted rna secondary structure in addition to sequence information. thus, despite the still incomplete and noisy map of rna-protein interactions, this python-based tool will prove valuable in elucidating complex post-transcriptional networks.  ptrcombiner is a new tool composed of multiple modules that infer meaningful combinatorial relationships between mrnas and their regulatory elements, namely rbps and ncrnas. this goal is achieved by extracting combinatorial information through a pattern-set mining approach and a meta-analysis of genome-wide data. ptrcombiner is divided into two activity components. the first, “mining combinatorial features”, represents interaction data with a mathematical model. the model involves an approximate boolean matrix factorization of the interaction matrix, which identifies groups of regulatory elements acting on common mrna utrs, which we call clusters. the second, “analyzing combinatorial features”, evaluates the biological characteristics of the clusters identified by the pattern-set miner. each cluster is evaluated in terms of global biological meaning by gene ontology analysis over its mrna targets and the binding site compatibility between the individual trans-factors in the cluster. description of the dataset the tool has been applied to the set of post-transcriptional interactions in human contained in the atlas of utr regulatory activity database. aura is a manually curated and comprehensive catalog of mrna untranslated regions and their regulatory annotations, including interactions with trans-factors. to date, aura is the largest dataset of utr-centered regulatory annotations taken from the whole range of existing experimental techniques, such as clip, rip, selex, and rnacompete. a subset of these techniques, namely clip and its variants followed by high-throughput sequencing, allows the positional annotation of the binding sites along the utrs, while the others can only detect the presence of an interaction between a transcript and a trans-factor without the positional information of the specific binding site. as a working example for ptrcombiner, we considered the whole set of human interactions annotated in aura the available number of utrs bound by each trans-factor ranges from to, with median of and a mean of. symmetrically, the available number of trans-factors bound to the same utr ranges from to with median and mean. using these data, we built a boolean matrix c with rows and columns, where cij = if the jth trans-factor interacts with the ith utr, and cij = otherwise. the annotated interactions are collectively, thus the sparsity of the interaction matrix is. mining combinatorial features after obtaining the interaction matrix, the first step of the analysis was to identify clusters of trans-factors that bind the same set of utrs. each of these clusters could be a candidate combinatorial member of a post-transcriptional regulatory code. this step thus aims to identify multiple overlapping clusters, which collectively cover most of the known interactions between trans-factors and utrs. boolean matrix factorization provides this requirement by decomposing the matrix of known interactions in the product of two boolean matrices. one of the matrices represents the clusters of the trans-factors, while the other, the utrs in terms of their interactions within the clusters. the algorithm takes two arguments: the number of clusters to return, and a threshold. the τ value controls the minimal amount of shared utrs inside a cluster. the higher the threshold, the more target utrs should be shared among trans-factors in order for these to be considered as a cluster. the algorithm returns a list of clusters ordered by coverage, i.e. the number of interactions associated to each cluster. to analyze the behavior of the algorithm when varying its parameters, we produced a surface parameter graph where the average number of trans-factors belonging to each cluster was calculated using various combinations of k and τ values. given a certain τ, different values of k do not affect the average cluster size, which appears to be mainly τ dependent. we chose the τ value giving an average cluster size equal to the average number of trans-factors bound to a single utr. we thus considered a τ of, resulting in clusters composed of an average of six trans-factors. interestingly, this value points to a stable region of the k and τ surface, in which the number of trans-factors for each cluster does not change drastically in the surrounding area. table shows the clusters obtained with the selected threshold. the first nine clusters are composed exclusively of rbps, as well as the clusters r to r r and r the first cluster displaying co-occurrence of rbps and mirnas was r followed by clusters r r r and r no cluster uniquely composed of mirnas was present. we would like to stress that out of clusters do not represent real combinations, as they comprise only one trans-factor. we refer to these one-element clusters as “singletons”. recalling that the algorithm aims to cover as many interactions as possible with the set of clusters, a singleton can be extracted by the algorithm whenever a trans-factor has many interactions that are not shared with any other trans-factor for which experimental interaction data are available. clusters were classified according to their composition: “rbp” for those composed exclusively of rbps; “mirna”, composed exclusively of mirnas; “rbp-mirna”, composed of both rbps and mirnas; and “singleton”, composed of only one trans-factor. as detailed in the previous section, different trans-factors have highly different numbers of annotated interacting utrs. being driven by coverage, the algorithm is inherently biased towards clusters of “widely interacting” trans-factors. therefore, when we analyzed the composition of the clusters, we observed that some trans-factors are present in multiple clusters. for example, the argonaute proteins ago and ago and the well-known rbp elavl1/hur, occur in and out of clusters, respectively. ago and ago are components of the rna-induced silencing complex, the protein complex which is responsible for down-regulating mrnas. these proteins bind different classes of small ncrnas, such as mirnas and small interfering rnas, leading the argonaute proteins to their specific targets through sequence complementarity, thus silencing their targets. therefore, it is not surprising to find ago and ago in almost all of the clusters, given the widespread activity of these proteins. moreover, in accordance with our results, ago and ago have been found to interact with elavl1/hur. despite this interaction, the two proteins avoid any binding overlap on the target mrnas: ago proteins preferentially bind the boundaries of utrs, while elavl1/hur binds uniformly along utrs and disappears toward the stop codon and the polyadenylation site. elavl1/hur, a member of the elav family, is known to be broadly expressed in tissues and to bind au-rich elements in the 3’ utrs of thousands of mrnas. moreover, it has been demonstrated that elavl1/hur displays competitive and cooperative interactions with mirnas/risc, that these interactions may depend on the proximity between the protein and mirna binding sites and that it is part of a complex mrna network for coordinating gene expression. these results are in agreement with the cluster composition in table and support our finding that the elements displaying the highest number of interactions are those that most frequently occur across the clusters. for this reason, we called these elements “recurrent trans-factors” and these clusters “rk”. to explore trans-factors that have a narrower spectra of interactions and thus occur less frequently in the clusters, we removed all recurrent trans-factors and ran a second iteration of the algorithm. this second iteration focused on trans-factors that appeared in none or only one of the clusters of the previous analysis. we termed these trans-factors “sporadic”. again, we studied the behavior of the algorithm when the k and τ parameters were varied, in order to select the most appropriate combination of values. in this case, the best choice of τ was. this value returned clusters comprised of an average number of three trans-factors, which corresponds to the average number of sporadic elements bound to each utr. we called these clusters “sk”. it is clear that the majority of the clusters are singletons. in contrast to results obtained when recurrent factors were included, we observed that four clusters were formed exclusively by mirnas. interestingly, pum which was also found as a member of the recurrent clusters, formed two distinct clusters with different sets of mirnas. in line with this, some evidence has suggested that pum associates with mirnas. pum is known to act as a translational repressor in several organisms, contributing to dendritic rna localization and silencing and regulating synaptic formation. moreover, an extensive interaction between pum and pum with the mirna regulatory system has been suggested, indicating that interactions between the rbps and the mirnas in translational regulation may be more frequent than previously thought. in addition, a recent computational analysis suggested that specific groups of mirna binding sites localize within nt from pum binding sites, giving support to a possible cooperativity between pum and mirna in mrna degradation and to the biological meaning of our clusters where pum acts on mrnas in combination with different mirnas. in particular, hsa-mir- and hsa-mir- which are part of cluster s together with pum appear to colocalize with this rbp. clusters were classified according to their composition, as for table despite that only a small fraction of known trans-factors occur in the clusters, the vast majority of the existing interactions is covered by the identified clusters. the fact that the most of the trans-factors are not part of any cluster is due to the paucity of available information, i.e. the overall number of experimental data obtained until now for many of the trans-factors is still far from being complete. this is not surprising given the novelty of the experimental techniques involved. it is clear that more data are needed to obtain a reliable and exhaustive description of all the possible combinatorial interactions in human. the availability in the next future of high-throughput experimental data covering an increasingly larger amount of trans-factors will open up new possibilities to uncover new and more reliable clusters and to obtain new combinatorial information. analyzing combinatorial features - step 1: biological annotation the first step in analyzing the clusters obtained from the mining module was to measure how much they differ in terms of: a) target mrnas; b) enriched ontological terms; and c) similarity among ontological terms. overall, we expected to obtain a large amount of regulated genes that form clusters in the presence of recurrent trans-factors from the mining module, because they display the highest number of annotated interactions. in fact, the module target mrnas overlap revealed that several hundred genes were co-regulated by recurrent trans-factors. the average number of hgnc genes regulated by the first five clusters was, ranging from for cluster r to for cluster r the number of target genes was markedly reduced in the clusters of sporadic trans-factors, as expected given their greater specificity. the average number of hgnc genes regulated by the first five clusters was ranging from for cluster s to for cluster s we then explored the possibility of finding the same mrna targets in different clusters, by calculating the overlap among the populations of genes grouped in different clusters using the jaccard similarity. for the clusters including recurrent trans-factors, the overlap was 21% on average. this result suggests that the method is able to group sets of specific target mrnas even if they share common trans-factors. however, some leakage is present: for example, cluster r and cluster r share 43% of their targets. this is not surprising as careful inspection of the elements forming these clusters revealed that r and r share seven rbps and are individually characterized only by fxr and fus, respectively. the same is observed for clusters r and r which share 60% of their targets. in this case, cluster r shares with r almost all its own trans-factors. focusing on clusters of sporadic elements revealed that the average overlap decreased to 7%, one-third of that from the previous analysis. in this case, the maximum overlap was between clusters s and s which share 27% of their targets. both of these clusters comprise caprin and mov while they do not share fus, fxr taf and c22orf. these results confirm the effectiveness of the repeated run of analysis in identifying distinct, small-sized sets of genes, uniquely regulated by specific sets of sporadic trans-factors. the gene ontology enrichment module of the tool allowed us to study the biological relevance of the mined clusters. the module attempts to identify parts of common and biologically coordinated mechanisms or processes that govern coherent cellular outcomes and that could be characterized by the identified groups of preferential interactions. in this module, the analysis is expanded from single genes to more general biological annotations, allowing the inferred clusters of trans-factors to be compared by the gene ontology enrichment analysis of their target mrnas. an initial and effective way to compare enrichments is to display the top enriched go terms for each cluster. such a comparison is shown in figure 4e for the top five non-singleton rk clusters, using the molecular function branch of go as example. the modular blocks of enriched terms scattered along the columns of the heatmap clearly indicate that the clusters have marked and distinct molecular function consensuses. similarly to what we observed for the mining module results, the modularity of the enrichments was further reinforced in the top five non singleton sk clusters. here, the occurrence of terms enriched in multiple clusters is an infrequent event. clusters s and s share the most similar enrichment signature, mirroring the strong similarity observed between the two clusters. after comparing the enrichment of different clusters, we next assessed how the ontological enrichment of genes regulated by one cluster differs from the ontological enrichments of genes regulated by the individual trans-factors forming the cluster. this intra-cluster comparison allowed us to potentially identify “emergent enrichments”, i.e. to detect go terms enriched exclusively in a set of genes regulated by a set of trans-factors forming a cluster. an example of this analysis is shown for cluster s. the target genes in this clusters were specifically enriched for the biological process term “cell division”, the cellular component term “nuclear speck”, and the molecular function term “transcription corepressor activity”. these results strongly suggest that the clusters have emergent and specific combinatorial properties, as expected for combinatorial mechanisms. to further estimate the similarity between the ontological enrichments associated with each cluster, the tool is able to calculate the semantic similarity values, which account for semantic similarity relationships between non-identical ontological terms. the module gene ontology similarity calculates the semantic similarity values for all the three branches of go and between both each pair of clusters and the single trans-factors belonging to the same cluster. in figure the top rows of each panel list the semantic similarity values between enriched terms associated to single trans-factors and enriched terms associated to cluster s here, the stronger semantic similarities are observed among the cc enrichments, while the weakest values lie among the bp enrichments. this analysis can also help to rank the ontological distance between a cluster and its trans-factors according to the semantic similarity values. for example, fxr globally shows the highest similarity with the enrichment of cluster s. analyzing combinatorial features - step 2: classifying rbp-binding sites to suggest some potential binding mechanisms based on the available experimental data, the second step of the analysis focuses on the rbps forming the previously identified clusters. the underlying idea is that whenever two rbps are characterized by similar binding sites over the mrna, then a concurrent binding, either competitive or cooperative, has occurred. although experimental methods can be used to gather information about the proximity of binding sites for specific rbps, the resulting information is corrupted by several types of noise sources: 1) a considerable fraction of binding sites can remain undetected because the methods are applied on cells of a particular type and in specific growth conditions, in which not all the potential bound mrnas are expressed; 2) interactions that are transient can be mistakenly identified as stable; and 3) post-processing analysis, such as mapping and peak detection, can increase the number of false negatives due to the difficulty of dealing with splice sites and the stringent thresholds needed for a confident detection, respectively. computational models for rbp target detection are therefore valuable tools in dealing with the low signal-to-noise ratio of current experimental techniques. such models can significantly increase the precision with which target sites are resolved and can uncover sites that would be otherwise missed by experimental protocols. to determine if two rbps are are likely to interact in a cooperative or competitive fashion, we first compiled an equivalent in silico model of the preferred rbp target sites. given a cluster of rbps and their binding sites, the module trains a machine learning algorithm to discriminate between binding sites of two different rbps, for all possible pairwise combinations of proteins belonging to the same cluster. two proteins are likely to have different binding sites when the algorithm effectively distinguishes between their binding areas. in contrast, compatible binding sites are more likely to lead to a difficult discrimination task. discrimination is based on a kernel machine binary classifier capable of computing similarity between base sequences in terms of their respective putative secondary structures and thus of the spatial conformation of their binding sites. the structural component of this discrimination has a strong biological significance, since rna interactions are not exclusively driven by sequence specificities. we report the results of the classification analysis performed on the first two clusters of sporadic trans-factors, namely, the cluster s that comprises caprin fus, fxr mov and taf and the cluster s that comprises c22orf caprin and mov. for each rbp, we randomly selected mrna utr stretches from available binding coordinates of large-scale experiments stored in the aura database. the classification performances for clusters s and s are displayed in figure 6a and figure 6b, respectively. performance is computed with the aurocc and f1-score measures. aurocc is an aggregate measure that evaluates the quality of a classifier when varying the threshold, to decide when a prediction should be considered positive. an aurocc value of corresponds to a random predictor, while an aurocc value of indicates perfect discrimination. the f score is the harmonic mean between precision and sensitivity, trading off the two complementary measures. considering the cluster s the classifier discriminates the binding sites of only a subset of the rbps in the cluster. the most specific binding sites were observed for caprin and mov. fus and taf seem to have more compatible binding sites. an aurocc of is in fact very close to the performance of a random predictor, suggesting that these proteins share a similar if not identical set of binding sites. indeed, these two proteins are known paralogues, both belonging to the fet family of rna-binding proteins. higher aurocc values were achieved when distinguishing the binding sites of fxr from those of fus and taf, and maximal aurocc values were found for mov suggesting different binding sites. classification performances of cluster s displayed in figure 6b, were generally high, suggesting that the utr stretches that are bound by the rbps forming the cluster are different. figures 6c and 6d show the distribution of the distances between couples of binding sites lying on the same mrna and targeted by distinct rbps, for cluster s and s respectively. the distances between the binding sites of fus and taf are much lower than those of the other distributions, indicating that the two proteins tend to bind to the same or very close regions. the average distance between fxr and fus or taf is also low, but it has a much larger spread. a large average distance is observed for the other cases and indicates a high discrimination capacity of the classifier. these examples show how we can start to use in silico modeling of rbps interactions to investigate their collaborative and/or competitive effects. this type of modeling is effective when the experimental data is affected by noise, since it recovers missed interactions and filters out accidental interactions. a related predictive system was used in, where it was shown that a model trained on a set of ago hits-clip sites could effectively identify binding sites missed by the high-throughput experiment. these findings have been verified comparing the predicted ago targets with changes in transcript expression levels upon ago knockdown. conversely this type of trend was not observed for the original hits-clip-detected sites indicating a significant impact of false negatives in the experimental setting. the upgrade to predictive models capable of this level of accuracy will allow more sophisticated investigations, such as those warranted by the combinatorial analysis presented in this work. for instance, a competitive effect can be hypothesized if two rbps, found in the same cluster, exhibit a compatible preference for the same target region, even when the experimental data are incomplete and do not report overlapping interaction areas. conversely, if the target regions are predicted to be sufficiently close but not overlapping, a cooperative effect can be hypothesized, even when experimental data cannot resolve the distinct areas and these are therefore reported as overlapping. alternative workflows - filtering interactions by experimental technique ptrcombiner allows a selection of the source interaction data to be use, thereby taking into account the different experimental approaches currently available for mapping protein-rna interactions and producing technically homogeneous results. this method-selection feature of ptrcombiner can be employed to exclude from the analysis the interactions obtained from techniques that are considered more unreliable or inopportune for specific analyses. as an example of the method-selection tool, we filtered the human interactions annotated in aura according to the evidence type, identifying three subsets of the original dataset. the first subset considered only par-clip experiments and includes a total of trans-factors and utrs, with an average number of interactions per utr of. the second subset considered all the other clip experiments, namely clip, clip-seq, hits-clip and iclip. it includes a total of trans-factors and utrs, with an average number of interactions per utr of. the third subset used only interactions found by rna immunoprecipitation rip. to date, it includes a total of trans-factors and utrs, with an average number of interactions per utr of. in the first subset, we observed the presence of clusters composed exclusively of rbp. the clusters are similar to those observed when the whole set of interactions was analyzed, excluding trans-factors with none or few interactions based on par-clip. considering only the rbp interactions obtained with par-clip experiments, the similarity with the clusters analyzed and described in the previous paragraphs is. after filtering for the clip, clip-seq, hits-clip, and iclip interaction data, too few rbps were left in the dataset to extract a reasonable number of clusters, so that no comparative analysis with the clusters obtained considering the whole set of interactions could be performed. finally, the clusters extracted from the third dataset had, on average, very few associated genes. we expect that this lack of large scale interaction data will be overcome once more datasets are made available. overall, these results suggests that, to date, the main component driving the selection of clusters is the set of interactions measured by par-clip experiments. - filtering targets and regulators by expression level so far, we have used ptrcombiner with a dataset that comprised all the binding evidence from multiple experiments in different biological systems that shared only the “human” source. the biological interpretation of results would greatly benefit from combining interaction data with other sources of information, such as expression data when available, in order to generate more specific hypotheses. to integrate ptrcombiner with expression data, we introduced the possibility to filter interactions according to: a) the expression levels of mrna targets; and b) the expression levels of the regulators. in this way, the interactions can be filtered according to the specificity of the biological system being studied. to show this functionality, we chose a well-studied cell line and filtered the aura dataset by selecting: a) the transcripts known to be highly expressed in hela; and b) rna binding proteins and mirna known to be present in hela. the filtered dataset now comprises utrs, rbps, and mirnas. we ran ptrcombiner on this filtered dataset, with the optimal τ parameter corresponding to. the resulting clusters are displayed in additional file 9: table s the majority of the clusters are composed of rbps, with the exception of six singleton clusters and three mixed clusters of mirnas and rbps. to analyze the effect of the cell-specific filter on the resulting clusters, we compared these clusters with those obtained in the first analysis that lacked any gene expression-abased filter. the comparison is shown in additional file 9: table s where each hela-cluster is matched to the general-cluster with maximal jaccard similarity. the average jaccard similarity was and ranged from nearly identical clusters to clusters sharing half of their members. interestingly, one-third of the general clusters could not be associated with any hela cluster. this is not surprising since these clusters are mainly composed of trans-factors or corresponding targets that are not expressed in hela cells. hela clusters therefore represent a subset of the general clusters. - balancing the trans-factor sample size an additional alternative workflow supported by ptrcombiner regards the mining step of the algorithm and arises from the previously discussed bias of the algorithm towards “widely interacting” trans-factors. introducing a balanced trans-factor association score is a possible way to reduce this bias. the difference between the standard version of ptrcombiner and the proposed alternative workflow is the way in which the pool of possible clusters is selected. the standard version of the algorithm incrementally extracts clusters of trans-factors by taking each candidate trans-factor as a “seed” and computing its association score with other trans-factors. the association score is the number of shared targets between the two factors, normalized by the number of targets of the seed. this implies that the seed is required to have a significant fraction of targets in common with another trans-factor for that trans-factor to be included in the seed’s cluster. by construction, this score is asymmetric and tends to associate trans-factors with many interactions together with those with fewer interactions that share a significant fraction of targets with them. here, we used a novel, alternative version of the association score normalized by the square root of the product of the target number for each trans-factor, thereby rebalancing the interaction data to cope with factors having a widely different number of targets. note that the optimal threshold τ in this case is different from the unbalanced one, as the alternative association score strongly affects clusters sizes. clusters obtained with the balanced score covered a larger number of trans-factors than did the standard version of the algorithm. although the average length of the clusters was the same, the balanced ptrcombiner produced more singleton clusters and a few very large clusters. intuitively, the larger the cluster, the fewer its associated genes. for this reason, the number of genes associated to the non-singleton clusters in the balanced case is lower than that in the unbalanced case. also, the jaccard similarity between clusters, calculated over their trans-factors, is lower in the balanced case. clusters in the balanced case are thus more specific. indeed, heat maps of the enriched go terms have less overlap with respect to the unbalanced case. a possible drawback of the balanced workflow is a tendency to create clusters that pick elements with a similar number of interactions. this excludes for instance clusters with mirnas and rbps that emerge from the standard workflow even after removing recurrent trans-factors. the two procedures thus have different characteristics, allowing users to discover different types of interesting combinatorial associations and to alternatively analyze the data according to their specific needs. comparison to related work ptrcombiner is a novel approach mining combinatorial post-transcriptional regulation patterns from interaction data at the genome-wide level. as mentioned in the background, previous attempts have been made in recent years to develop automated approaches to identify combinatorial aspects of post-transcriptional regulation. several computational tools have been developed for mirna target prediction using a number of different approaches, such as targetscan, miranda, targetboost, pita, magia, mirtar hunter. target site predictions for different mirnas have been combined in the tools pictar and comir to identify potential groups of interacting mirnas. by using simple model systems such as s. cerevisiae, some published studies have attempted to develop computational and experimental methods that identify functional modules based on combinatorial rna-protein interactions. for example, joshi and coworkers developed a probabilistic method for inferring regulatory module networks from expression profiles. in the following, we provide a more in-depth comparison with these approaches, highlighting the differences with respect to our proposed method and reporting comparative quantitative analyses on benchmarks for which these alternative approaches can be run. - comparison to pictar pictar is a probabilistic predictive method that computes the probability of multiple mirnas co-binding to the same target mrna by combining the binding scores for each candidate mirna. albeit focusing on combinatorial interactions, the algorithm differs from ptrcombiner in many aspects. first, it explores the combinatorial interaction of mirnas only. second, it relies on predictive methods for detecting potential binding sites rather than using experimental data. the predictive approach can be reasonably accurate for mirnas but are still far from satisfactory for rbps. the third and main caveat of the method is that it lacks a mining procedure that allows to efficiently search the combinatorial space of possible clusters. rather, it either requires the user to specify a set of mirna to be jointly evaluated, or it needs to try all possible combinations in order to identify the high scoring ones. ptrcombiner takes a different perspective, implementing a mining approach to efficiently explore the combinatorial space of candidate clusters of mirna and/or rbp, guided by their coverage of observed interactions with the target mrnas. as discussed in the conclusions, we plan to adapt the method to deal with weighted interactions, e.g., the score of a predicted binding or the confidence of a certain experimental technique. in order to quantitatively compare the two approaches, we analyzed clusters of mirna found by ptrcombiner using pictar. we focused on the clusters s and s because they are composed of only four mirnas, and evaluating larger clusters with pictar is too computationally expensive. for each cluster, we took the set of its target genes, i.e. the genes which interact with all mirnas in the cluster, and computed the pictar interaction score of each of them with the cluster. the score is computed as the maximum value of the product of the binding scores of the mirnas in the cluster, with the constraint that binding sites of different mirnas should not overlap. if this constraint cannot be satisfied, or one or more mirnas do not have predicted binding sites on the target, the score was set to zero. binding scores for mirnas were downloaded from the dorina database. we then compared these cluster-target scores with those obtained by running the same procedure on the entire set of genes using the welch’s two samples t-test. the results of the statistical tests are presented in additional file for both clusters, the difference between cluster-target scores and general scores was statistically significant, with a confidence of approximately 99%. the fact that the average on the former set is significantly higher than the average on the latter set is not a trivial information. in fact, it indicates that also pictar estimated the clusters of trans-factors to be relevant exactly to the specific set of genes used to cluster them together, confirming the relevance of the clusters mined by ptrcombiner. - comparison to comir comir is a web tool for combinatorial mirna target prediction. it uses mirna expression levels in combination with thermodynamic modeling and machine learning techniques to make combinatorial predictions. comir sums up the weighted scores of the single mirnas, computed according to the four different scoring schemes of miranda, pita, targetscan, and mirsvr. these scores are combined through a support vector machine, which outputs the likelihood that the set of mirnas targets a specific gene. as for pictar, the main difference with ptrcombiner is the lack of a mining procedure that would identify the clusters of mirna to be evaluated. using comir, we analyzed all clusters composed only of mirnas extracted by ptrcombiner, namely, s s s and s. in computing scores, we gave no expression level to comir, resulting in uniform level for all mirnas. as for the pictar case, we compared cluster-target scores with scores for the entire set of genes, i.e. all genes in comir output. welch’s two sample tests confirmed the statistical significance of the difference between cluster-target and general scores for all clusters, with a confidence of approximately 100%. similarly to the comparison to pictar, this result confirms the relevance of the clusters mined by ptrcombiner. - comparison to lemone lemone is a probabilistic method for inferring regulatory module networks from expression profiles. the module network model was introduced in the work by segal et al. as a probabilistic bi-clustering approach, extracting co-clusters of genes and conditions from a matrix of gene expression levels under different experimental conditions. a co-cluster contains a subset of genes and conditions, such that the genes have similar expression levels under the selected conditions. the subset of genes represents a regulatory module, while each subset of conditions for the same module is an expression state for the module. lemone introduces an ensemble averaging strategy to generate more coherent modules from multiple runs of the module network inference process. the approach was later adapted to infer transcriptional and post-transcriptional modules from both transcriptome and translatome expression profiles in s. cerevisiae. as already discussed, this approach detects putative regulatory modules that characterize specific biological conditions, while ptrcombiner targets more general purpose, genome-wide combinatorial interactions. nonetheless, it is interesting to study the relationship between clusters detected by the two methods, when both are applicable. to compare this method with ptrcombiner, we ran our computational tool on the yeast dataset employed in. the interaction dataset contains rip-chip experiments involving rbps and genes. we obtained a list of interacting rbps and their relative number of target genes. the distribution of the number of distinct trans-factors bound to the same gene are shown in figure 8b. in total, the dataset contains annotated interactions, with the interaction matrix sparsity value of. from the interaction matrix, ptrcombiner extracted the clusters composed exclusively of rbps. these clusters were then compared with the sets of rbps whose targets were enriched in post-transcriptional modules generated by lemone. the list of rbp sets was extracted from the supplementary materials of. ptrcombiner clusters were then compared with lemone clusters, according to the jaccard similarity maximized by matching each ptrcombiner cluster with a lemone cluster. reassuringly, the agreement between the two methods is high, with half of the top ten ptrcombiner clusters identical to lemone clusters.  in this paper we present a computational tool for the combinatorial analysis of post-transcriptional regulation patterns involving multiple trans-factors. this tool, called ptrcombiner, was tested on two sets of experimental interactions between post-transcriptional trans-factors and target mrnas, in human and yeast. ptrcombiner enables the user to: a) detect groups of regulators that share a conspicuous amount of mrna targets; b) characterize the clusters biologically; and c) identify concurrent binding sites of trans-factors belonging to the same cluster. further, the underlying boolean matrix factorization approach allows the user to identify multiple overlapping clusters of trans-factors that jointly account for as many interactions as possible, casting the problem into an approximate interaction coverage task. this method naturally addresses the limitations of most clustering and motif mining approaches, which typically return non-overlapping clusters, clusters covering non-overlapping sets, or long lists of highly redundant clusters. our tool is an original and comprehensive attempt to provide a computational pipeline for elucidating complex post-transcriptional combinatorial rules on a genome-wide level. by integrating expression profiles of both trans-factors and target mrnas, the tool can also be used to mine combinatorial patterns in specific experimental conditions. we plan to extend the method to deal with uncertainty in the interaction information, such as putative interactions from predictive algorithms, and to incorporate the binding strength by properly weighting the contribution of each interaction. importantly, ptrcombiner is a versatile tool that is not limited to post-transcriptional regulation analysis. in fact, it can be easily adapted to mine transcriptional and combined transcriptional/post-transcriptional regulation patterns. an inherent limitation in the analyses that can be conducted with the tool is given by the paucity of interaction data available to date. however, given the fast rate at which high-throughput interaction detection experiments are conducted, this tool may have the potential to unveil the complex mosaic of interactions underlying post-transcriptional regulation in the near future. 
the progress in mapping rna-protein and rna-rna interactions at the transcriptome-wide level paves the way to decipher possible combinatorial patterns embedded in post-transcriptional regulation of gene expression. post-transcriptional regulationboolean matrix factorizationrna binding protein binding site classificationkernel machinesmirnatranslationclip in addition the tool allows to analyze the binding site similarity of co-regulators belonging to the same cluster, given their positional binding information.
 the progress in mapping rna-protein and rna-rna interactions at the transcriptome-wide level paves the way to decipher possible combinatorial patterns embedded in post-transcriptional regulation of gene expression.  here we propose an innovative computational tool to extract clusters of mrna trans-acting co-regulators from pairwise interaction annotations. in addition the tool allows to analyze the binding site similarity of co-regulators belonging to the same cluster, given their positional binding information. the tool has been tested on experimental collections of human and yeast interactions, identifying modules that coordinate functionally related messages.  this tool is an original attempt to uncover combinatorial patterns using all the post-transcriptional interaction data available so far. ptrcombiner is available at http://disi.unitn.it/~passerini/software/ptrcombiner/. post-transcriptional regulationboolean matrix factorizationrna binding protein binding site classificationkernel machinesmirnatranslationclip 