#microCLIP super learning framework uncovers functional transcriptome-wide miR

[naivelocus]

microCLIP super learning framework uncovers functional transcriptome-wide miRNA interactions

Dataset collection

6724 high confidence MREs were retrieved from direct experiments, including reporter gene assay techniques indexed in DIANA-TarBase repository15,16, miRNA-chimeras from CLASH (crosslinking, ligation, and sequencing of hybrids)13 and CLEAR-CLIP (covalent ligation of endogenous Argonaute-bound RNAs)14 experiments, as well as additional miRNA-target chimeric fragments derived from a previous meta-analysis of published AGO-CLIP datasets12. In order to quantify miRNA-induced mRNA expression changes and to identify functional binding sites, 101 miRNA perturbation experiments were analyzed (89 microarray and 12 RNA-Seq experiments, Supplementary Tables 4–6). This process enabled the formation of ~3900 and 4000 positive and negative miRNA-target pairs, respectively. A set of five ribosome profiling sequencing (RPF-Seq) libraries after miRNA overexpression, capturing differentially ribosome-bound transcripts, and six pSILAC (quantitative proteomics) experiments were an additional source for detecting more than 5900 miRNA effects at protein expression level (Supplementary Tables 7–8). The inclusion of AGO-IP and biotin pull-down high-throughput experiments upon specific miRNA perturbation yielded ~2600 miRNA-binding events (Supplementary Table 9). The aforementioned miRNA perturbation experiments enabled the detection of deregulated targets without specifying the exact binding sites15. miRNA-targeted regions were extracted from AGO-bound enriched regions present in at least 1 of 24 AGO-PAR-CLIP-sequencing libraries (Supplementary Table 1). Published background PAR-CLIP libraries32, stably expressing a commonly utilized non-RBP control (FLAG-GFP), were incorporated in our pipeline to identify non-specific AGO-bound transcripts and deduce more than 24,000 negative miRNA-binding sites. A compendium of 96 AGO-CLIP-Seq experiments was derived from DIANA-TarBase and used to further select background-derived MREs displaying no overlap with AGO-enriched regions (Fig. 2a).

Analysis of high-throughput miRNA perturbation experiments

High-throughput experiments were collected to measure gene expression alterations after specific miRNA transfection, silencing, or knockout. Log2 fold-change values as calculated from differential expression analyses of control versus post-treatment state enabled the formation of miRNA–mRNA positive and negative interactions.

A total of 44 microarray studies of distinct experimental conditions (Supplementary Table 4, 6) covering 43 human cell lines and 49 miRNAs were examined to deduce positive and negative miRNA-target interactions. In-house analysis was initiated from microarray raw data (Affymetrix.CEL files). Probe set summarization was implemented using Robust Multi-Array Average (RMA) with R packages affy33 or oligo34. Annotation of probe sets to Ensembl Gene IDs was accomplished using the chip-specific annotation R packages hgu133a2.db, hgu133plus2.db or hugene10sttranscriptcluster.db. miRNA-treated and control samples in each experiment were analyzed independently of other cell lines or miRNA treatments. Log2 fold-change ratios and p values were calculated with limma package35, following package instructions on Single-Channel Designs. Probe sets mapping to the same gene were averaged to calculate its fold-change. A log2 fold-change cutoff of ±1 ( > 1 or < −1, respectively), depending on the type of regulation, was applied to determine negative and positive interaction subsets. For GSE8501 experiment conducted in Rosetta–Merck microarrays, error-weighted log10 intensity ratios were retrieved and transformed to log2-scale.

Ribosome profiling sequencing (RPF-Seq) and RNA-Seq libraries treated with specific miRNA overexpression, 12 experimental conditions in total (Supplementary Tables 5–7), were retrieved from Eichhorn et al.36, Nam et al.11, Pellegrino et al.37, Zhang et al.38. To identify positive/negative miRNA interactions, a ±0.5 log2 fold-change threshold was applied to genes presenting > 10 RPKM expression.

Quantitative proteomics datasets (pSILAC) in HeLa cells following the individual overexpression of five human miRNAs (let-7b, miR-1, miR-16, miR-30a, and miR-155) or knockdown of let-7b (Supplementary Table 8) were derived from Selbach et al.20. Positive/negative miRNA interactions were deduced using a ±1 log2 fold-change threshold, respectively.

Analysis of AGO-PAR-CLIP and (s)RNA-Seq expression datasets

AGO-PAR-CLIP datasets from nine studies, corresponding to 13 cell lines in human species, were derived from GEO7,39 and DDBJ40 repositories (Supplementary Table 1). Fifteen small RNA-Seq and 9 RNA-Seq experiments of similar cell types with PAR-CLIP libraries were analyzed following methodologies as described by Vlachos et al.41 to infer expressed miRNAs and transcripts. (s)RNA-Seq datasets were derived from the ENCODE repository and from a series of studies (Supplementary Tables 10, 11). Whole transcriptome depleted from ribosomal RNAs and poly-A selected RNA-Seq libraries were analyzed.

Pre-processing and alignment of PAR-CLIP datasets was realized as described by Vlachos et al.15. Initially, libraries were quality checked using FastQC (www.bioinformatics.babraham.ac.uk/projects/fastqc/). Adapter sequences were retrieved from the original publication or GEO/SRA entries, when available. Contaminants were detected utilizing in-house developed algorithms and the Kraken suite42. Pre-processing was performed utilizing Cutadapt43. PAR-CLIP reads were aligned against human reference genome (GRCh37/hg19) with GMAP/GSNAP44 spliced aligner, appropriately parameterized to identify known and novel splice junctions. microRNA expression was quantified using miRDeep245. Ensembl v7546 and miRBase v1847 were used as annotation for genes and microRNAs, respectively. Top expressed miRNAs and AGO-PAR-CLIP data in each cell type, were jointly utilized as input to microCLIP in silico framework for miRNA-target identification. Specifications on the processed 36 datasets are provided in Supplementary Tables 1, 10, 11.

For the analyses presented in Figs. 3–5 and Supplementary Figs. 4, 7, enriched AGO-CLIP peaks covered with reads having at least 20% cross-linked sites in the same position are defined as T-to-C targeted regions.

Analysis of PARS experimental data

PARS sequencing data on total RNA isolated from lymphoblastoid cells were obtained from Wan et al. study18 (GEO accessions GSM1226157, GSM1226158). The identification of single or double stranded regions, across the human transcriptome, was derived from deeply sequenced RNA fragments generated from RNase S1 or V1 nuclease treatment of GM12878 cells, respectively.

Raw reads of 51nt length, accordingly pre-processed for quality control and contaminant removal, were aligned against human reference genome (GRCh37/hg19) with GSNAP spliced aligner. This analysis resulted in ~130 M uniquely mapped PE-sequenced fragments per sample. In order to derive structural signals in RNase S1 or V1 nuclease experiments at single base resolution, we calculated single(S1) and double(V1) stranded raw reads initiating on each nucleotide. The number of PARS tags per sample starting at each base were normalized by sequencing library depth. These base intensities were subsequently combined into the formula described by Wan et al. to compute PARS scores.

RSS were defined by estimated PARS scores in the vicinity of PAR-CLIP-derived miRNA-binding sites in four lymphoblastoid cell lines from the study of Skalsky et al.2. miRNA-mRNA interactions were identified in both T-to-C and non-T-to-C PAR-CLIP clusters, corresponding to transcripts with > 1 TPM expression in GM12878 cells. For expressed miRNAs ( ≥ 50 aligned reads per miRNA) in respective EFD3-AGO2, LCL-BAC, LCL-BAC-D1, and LCL-BAC-D3 EBV infected lymphoblastoid cells, we included collapsed miRNA-binding sites residing within the PAR-CLIP clusters. For the performed comparisons, we incorporated negative MREs extracted from different high-throughput miRNA perturbation experiments (more detailed description in Methods “Analysis of miRNA transfection/knockdown high-throughput experiments”). MREs utilized for the assessment of RSS signatures on AGO-bound clusters and the derivation of (non-)functional conformations of miRNA-target base pairings, were localized on coding and 3′UTR regions. The examined sites had to present S1 and V1 signals in at least half of their occupied bases.

sRNA-Seq and RNA-Seq datasets were retrieved from ENCODE consortium (GEO accession numbers GSM605625, GSM1020026, GSM1020027, GSM1020028, GSM1020029, and GSM1020030).

microCLIP in silico framework

Feature set description: A set of 131 descriptors (Supplementary Data 1) with non-zero variance was included in microCLIP. The extracted features were retrieved from positive/negative miRNA interactions, identified on AGO-bound locations in different PAR-CLIP datasets. They comprised PAR-CLIP-specific descriptors, such as substitution ratios and distance of conversions from the MRE start, as well as coverage metrics. Aggregate substitution ratios, positions, and distances independent of the transition type were also included. In order to estimate MRE and AGO-peak respective sequencing coverage, we calculated normalized RPKM values for miRNA-target sites and clusters.

Moreover, single base and dinucleotide contents for miRNA-binding and respective flanking regions, complexity features for the MRE and proximal upstream/downstream sequences were introduced to microCLIP model. BLAST’s DUST score48 and Shannon-Wiener Index49 constituted measurements for masking sequence complexity. Other descriptors were formed to represent energy-related variables for the duplex structure, while metrics capturing sequence content skewness/asymmetry (GC-skew, AT-skew, purine-skew, Ks-skew) and biases of codon usage were added. Entropy, enthalpy, free energy, and melting temperature (Tm) thermodynamic properties were calculated for MRE sequences in R.

miRNA-target hybrids were associated with different descriptors such as the binding type, duplex structure energy calculated with the Vienna package50, positions and nucleotide composition of (un)paired nucleotides. Distinct features have been established to model (mis)matches, bulges, loops, and wobble pairs for miRNA-MRE hybrid structure and sub-domains encountered in the duplex. The distinct domains for miRNA sequences, as defined by microCLIP, are: (i) seed region (2–8 positions), (ii) central region (9–12 positions), (iii) 3′ supplementary/compensatory region (13–16 positions), (iv) tail region (17-3′ miRNA end). Similar regions were designated on the MREs based on the miRNA-binding anchors upon duplex formation.

Our approach incorporates conservation of the MRE and upflank/downflank-MRE regions. phastCons pre-computed scores from genome-wide multiple alignments were downloaded from the UCSC repository51 in bigwig format and were utilized to deduce respective evolutionary rates. Conservation signals were computed as mean intensities of the phastCons base-wise scores on miRNA targeted regions, as well as their flanking regions. The conservation of the 5′ MRE binding nucleotides was independently modeled. microCLIP integrates additional features corresponding to the location of the MRE within the AGO-enriched cluster and binding length ratios of miRNA and/or target regions.

The applied super learning scheme benefits from the incorporation of the complete array of features, maximizing their contribution through their parallel use in different classification models in every node. The impact of weaker features and classifiers in optimal super learner design and behavior is shown in Supplementary Fig. 9, where microCLIP performance was compared to three different classification schemes using the independent validation set available in Supplementary Data 3.

Description of the algorithm: microCLIP operates on AGO-PAR-CLIP-sequencing reads, requiring a SAM/BAM alignment file and a list of miRNAs as minimum input. It initially seeks for AGO-enriched regions and resolves coverage and observed transitions. A sensitive pipeline is adopted to scan read clusters for putative targeted sites including a wide range of binding types. The algorithm supports an extended set of (non-)canonical matches including 6mer to 9mer, offset 6mer, 3′supplementary and compensatory sites as well as (im)perfect centered bindings. microCLIP extracts features for each candidate MRE and subsequently scores sites through a super learning scheme.

The adopted framework incorporates two distinct levels of classification. The first layer comprises a group of 9 different nodes (base classifiers), which are aggregated in the meta-classifier of the second layer. The learning procedure is decentralized through the distinct nodes and relevant base classifiers that specialize in different subsets of features (Fig. 2b). “Region Features” node comprises CLIP-Seq-derived features, such as RPKM coverage, substitution frequencies, and region-related descriptors, including nucleotide composition, conservation, sequence energy, complexity, content asymmetry, and biases of codon usage. A set of five additional base classifiers were designed for MRE-specific features. Binary binding vectors of miRNA/MRE hybrid were separately incorporated in a base classifier (“Binding Vectors”). Each vector element corresponds to one (un)paired position in the duplex. Matches per miRNA/MRE sub-domain were added to a distinct base classifier introducing a group of 13 features regarding total and consecutive matches in the miRNA-target structure as well as in MRE and miRNA relevant sub-domains. Another base model consists of miRNA-target duplex descriptors (“Duplex Features”) including miRNA-target duplex structure energy, bulges, internal loops, GU wobbles, and AU base pairing features for the specified miRNA and/or target and relevant sub-domains. The “Base pairing” node encompasses composition descriptors (A, T, G, C) of the (un)paired nucleotides. An extra base learner incorporates MRE general descriptors such as the degree of overlap with the respective cluster, conservation of MRE bound nucleotides, MRE location within the cluster, MRE binding type as well as metrics for duplex paired nucleotides content asymmetry/skewness. The latter five base models are dedicated to the determination of genuine miRNA-binding sites. Non-overlapping feature sets from the aforementioned base nodes are combined into three supplementary classifiers also incorporated into microCLIP framework.

Eight of the nine base nodes adopt a super learning scheme that assembles the output of seven individual Random Forest (RF), Generalized Linear Model (GLM), Gradient Boosting Model (GBM), Deep Learning (DL) classifiers (2 RF, 2 GBM, 2 DL, 1 GLM models). The “Region features” is analyzed by an RF classification scheme. The retrieved scores from each node are aggregated in a final GBM meta-classifier.

Model training: The DL models developed for the microCLIP framework adopt a feed-forward multi-layer architecture. The input layers match the respective feature space and values are subsequently propagated within three hidden layers. We utilized a rectifier activation function to retrieve weighted combinations of the inputs transmitted to interconnected neuron units. Dropout regularization was added to achieve model optimization and avoid overfitting. A cross entropy cost-function was selected to adapt weights during the learning process by minimizing the loss. Bernoulli distribution function was used along with cross entropy (log-loss) to model the response variables. The output layer at the end of the network applies a Softmax activation function so that each neuron (predicted class) results in a probabilistic interpretation. The DL network depth, width, and topology, as well as activation functions and learning parameters were modeled with a tuning-in grid search algorithm using H2O52 R package. The RF, GBM, GLM learning models were developed, parameterized, and tuned with the caret53 and H2O52 R packages.

Base classifiers were trained against a collection of 8693 positive and 21,789 negative miRNA interactions (Supplementary Table 3). The final GBM meta-learner that aggregates the base classifier outcomes was trained against an independent dataset comprising 3276 and 6702 positive and negative instances, respectively. Ten-fold cross-validation was performed on the training data to estimate each model’s accuracy and finalize the algorithm’s learning architecture. Distribution of base model scores on positive and negative instances and their respective performance, in terms of sensitivity and specificity in an independent test set of ~4000 instances (Supplementary Table 3), are depicted in Supplementary Fig. 10. The individual performance of internal classifiers (DL, RF, GBM, GLM) in microCLIP base models adopting a super learner approach is shown using the same set in Supplementary Figs. 11, 12. Additionally, the performance of the super learner scheme against Random Forest models was tested (Supplementary Fig. 9). microCLIP training and required computations for model optimizations were multi-threaded.

A microCLIP model adopting the same super learning scheme, including information only from T-to-C enriched sites, microCLIP T-to-C, was also deployed. Clusters from the training set incorporating adequate T-to-C transition sites were selected as input to re-train the super learning classifier. Additional support for the robustness of CLIP-guided super learner classification irrespective of non-T-to-C site inclusion is provided through evaluations of microCLIP T-to-C algorithm performance (Supplementary Fig. 8), as well as of its prediction efficacy (Supplementary Fig. 13) using the same test and datasets as in Fig. 6.

Expression correlation analysis on (s)RNA-Seq TCGA samples

271 TCGA ductal breast cancer RNA-Seq and sRNA-Seq samples were obtained from Firehose (http://gdac.broadinstitute.org/runs/stddata__2016_01_28). mRNA and miRNA pre-computed expression values were RPM and RPKM units, respectively. In downstream analyses, miRNAs/mRNAs with zero expression value in at least 70% of the samples were filtered out. miRNAs presenting > 10 RPM in > 10% samples were included, based on miRBase criteria for defining a high confidence set47. The mRNA set was specified by applying a threshold of more than 1 RPKM in at least 10% of the processed samples. Zero mRNA expression values were replaced by the lowest non-zero RPKM value per sample. This pipeline resulted in a set of 13,346 mRNAs and 322 expressed miRNAs. Pearson correlation coefficient was computed for each miRNA-target pair across samples.

Functional analysis of AGO-PAR-CLIP-derived miRNA targets

miRNA-target pairs were retrieved from the analysis of MCF7 PAR-CLIP library (Farazi et al.3) with microCLIP in silico framework. The 100 most highly expressed miRNAs and their targets in 3′ UTR regions were retained. Gene set enrichment analysis of AGO-PAR-CLIP-detected miRNA targets was performed for KEGG pathways54 using the R package limma35. Enrichment P values were corrected for multiple comparisons using Benjamini–Hochberg false discovery rate and a 0.01 p-value threshold was applied. R package Pathview55 was used to visualize targeted pathway members in KEGG pathways.

miRNA interactions from in silico implementations

In order to form a complete list of interactions for MIRZA4, microMUMMIE6, and PARma5 computational approaches, each algorithm was evaluated on a set of 7 PAR-CLIP HEK293 libraries obtained from Kishore et al.21 and Memczak et al.31 studies (GEO accessions GSM714644, GSM714645, GSM714646, GSM714647, GSM1065667, GSM1065668, GSM1065669, and GSM1065670). The proposed settings for each implementation were retrieved from the relevant publications.

The MIRZA biophysical model was executed in the “noupdate” mode. The algorithm provides an optional parameterization to introduce miRNA expression profiles. We realized two different runs of MIRZA, with and without cell type-specific miRNA expression values that were extracted from the CLIPZ web server (http://www.clipz.unibas.ch). MIRZA input data were 51nt AGO-bound sequences centered on T-to-C sites and mature miRNA sequences of 21nt length as reported in the model’s restrictions. The “target frequency” score was utilized to evaluate the quality of MIRZA-detected sites.

microMUMMIE algorithm was tested in both Viterbi and posterior decoding modes. Following microMUMMIE’s constraints, we utilized PARalyzer v1.57 to define the set of T-to-C AGO-enriched peaks. An extra pre-requisite annotation step to complement PARalyzer detected clusters was implemented with the PARpipe tool (https://github.com/ohlerlab/PARpipe). Derived files, comprising annotated AGO clusters with positions of T-to-C transitions, constituted the input of the microMUMMIE core algorithm. Predictions with signal-to-noise ratio (SNR, generally correlated with sensitivity) equal to 9.95 were retained, while posterior probabilities were utilized for the evaluation of microMUMMIE’s performance.

PARma was applied on AGO-PAR-CLIP aligned data that were prepared following the algorithm’s described format. The required input files contained genomic locations of aligned CLIP reads along with positions of observed conversion sites. PARma predictions are coupled with Cscore and MAscore scores for the cluster and miRNA seed family activity, respectively. The latter score was utilized for PARma-detected miRNA-target sites evaluation.

Precompiled (non)conserved miRNA site context++ scores for representative transcripts were downloaded from the Targetscan v7.2 site (http://www.targetscan.org/cgi-bin/targetscan/data_download.vert72.cgi). Targetscan v7 algorithm was additionally executed following the proposed settings in order to cover a greater transcript collection, as well as the whole spectrum of Targetscan-detected interactions including 6mer sites. Gene annotation files were retrieved from the Targetscan v7.2 official download page, and the miRNA seed sequence file that is a pre-requisite for the execution of the model was provided by Targetscan developers. The local Targetscan run complements the precompiled data with miRNA-target interactions on transcripts presenting the longest 3′UTR, in cases they are not deposited on the online repository.

Median fold changes

The comparison between microCLIP and existing implementations was performed using six gene expression profiling datasets following individual transfection of highly expressed miRNAs into HEK293 cells (GEO accessions GSE60426, GSE52531, GSE21901, GSE14537, GSE35621, Supplementary Table 6). Genes with unchanged expression levels (zero log2 fold change) following miRNA transfection and/or knockdown have been filtered out. Subsequent measurements were realized at the gene level. miRNA-gene interactions retrieved from each implementation were sorted according to their prediction scores. Each miRNA-target pair was characterized by the highest scored miRNA-binding site overlapping coding or 3′UTR exons, since utilized algorithms provided MRE-oriented prediction scores. In cases of multiple transcript-gene associations, the transcript with the longest 3′UTR was selected. Median expression log2 fold changes were estimated in consistence with the number of top predicted targets. Aggregated expression changes of genes were calculated by applying stepwise score thresholds. Paired comparisons required tested programs to have targets at every computational cutoff. Lower mean log2 fold changes correspond to stronger downregulation of the detected targets upon miRNA transfection. The statistical differences in the mean log2 fold-change values obtained by each implementation were assessed using two-tailed Wilcoxon signed-rank test. Identified targets by each algorithm were also juxtaposed against averaged log fold changes of 1000 randomly selected genes (without replacement). The mean log2 fold-change values of the randomly selected genes in different stepwise thresholds were taken and the median curve derived from these values was calculated. Genes with zero fold-change indication were filtered out from the random selection process.

Statistics

Enrichment analyses were performed using one-sided Fisher’s exact test. Correlations between quantitative parameters were assessed by calculating Pearson’s correlation coefficient. Comparisons between two or more groups were conducted using Wilcoxon’s rank-sum test and Kruskal-Wallis’ test, respectively. In the latter, Wilcoxon’s rank-sum test was performed as a post hoc test in order to assess between-group differences. The one-sided Kolmogorov–Smirnov test was used to test for greater functional efficacy. In cases of multiple hypothesis testing, Benjamini–Hochberg’s false discovery rate was applied to control family-wise error rate. P values < 0.05 were considered as statistically significant.

— Nature Communications