Systems and methods for predicting post transcriptional gene regulation using deep learning

Abstract

Provided are systems and methods for tissue specific modeling of post-transcriptional gene regulation (PTGR) using a neural network. The PTGR model operates to predict tissue- specific measures of mRNA stability such as miRNA binding and degradation. The systems and methods may use a dilated convolutional neural network for predicting mRNA stability based on miRNA binding and mRNA degradation in a sequence-to-sequence manner. Also provided are methods for generating PTGR models using miRNA binding data and mRNA degradome data.

Description

SYSTEMS AND METHODS FOR PREDICTING POST TRANSCRIPTIONAL GENE REGULATION USING DEEP LEARNING CROSS-REFERENCE TO RELATED APPLICATIONS

[0001] The present application claims priority to US Provisional Application 63 / 557,943 filed February 26th, 2024 the entire contents of which are incorporated herein by reference. FIELD

[0002] The present embodiments related generally to the field of post-transcriptional gene regulation and more specifically to systems and methods for predicting post- transcriptional gene regulation using deep learning. INTRODUCTION

[0003] MicroRNAs and RNA binding proteins are crucial elements of post- transcriptional gene regulation, which governs the fate of mRNA molecules in the cell. However, the landscape of these regulatory interactions, particularly across different mammalian cell types, remains underexplored.

[0004] In eukaryotes, gene expression is regulated at various stages including transcription, RNA processing (e.g., splicing and RNA editing), nuclear export, and translation. Post-transcriptional gene regulation (PTGR) is particularly critical and affects mRNA stability and localization, protein synthesis, and ultimately cell function. PTGR is mediated primarily by interactions between mRNAs and two main classes of regulatory factors: RNA-binding proteins (RBPs) and microRNAs (miRNAs).

[0005] miRNAs are short regulatory RNAs that guide the RNA-Induced Silencing Complex (RISC) to target mRNAs via partial sequence complementarity, leading to repression through transcript destabilization and translational inhibition. Interactions between miRNAs and target RNAs are either canonical, with perfect base pairing within the miRNA seed region (positions 2-7), or non-canonical, consisting of imperfect seed matches and occasional extended 3’ end pairing. Conversely, RBPs directly interact with both pre-mRNA and mature mRNA to influence PTGR through splicing, RNA editing, transport, localization, degradation, and translation.

[0006] Over 2,300 miRNAs and 1,500 RBPs have been identified in humans, and collectively engage mRNA via a combination of sequence-specific, structural, and cofactor- dependent interactions. The diversity of functions carried out by these molecules reflects the complexity of the post-transcriptional regulatory networks they participate in. Despite considerable progress in elucidating the functions of miRNAs and RBPs, the full complexity of their interactions and the associated regulatory outcomes across different cellular environments remains only partially understood.

[0007] Existing machine learning models for predicting miRNA-mediated target repression either rely on hand crafted features

[0013] or are trained on a small number of miRNAs and are limited to very short input sequences

[0014] . Thus, there may be additional sequence determinants of miRNA targeting that are not fully captured by existing models. Other deep learning approaches have relied on artificial sequences as negative sets for model training, restricting performance due to their limited ability to learn the endogenous biology governing miRNA targeting [17-20]. A substantial limitation of most existing approaches is that they do not take cell type into account, despite the fact that tissue and cell type specific miRNA expression and targeting are relevant to both development and disease [21-23].

[0008] Moreover, the training of PTGR models and the prediction of miRNA-mediated target repression require a massive amount of data processing due to the combinatoric explosion of various mRNA sequences and the number of potential miRNA binding candidates. The nature of training these models and their use in predicting miRNA binding to mRNA sequences means they are by necessity computer-based solutions to computer- based problems.

[0009] Efforts to incorporate (non polyadenylation-related) site-specific endonucleolytic cleavage data have been primarily focused on identifying miRNA targets rather than predicting general mRNA degradation [24-26], or have been more descriptive in nature [6, 27-31]. Regarding data used to train models, new technologies have made transcriptome-wide identification of miRNA-specific target sites possible [32-34]. However, these methods have been applied to only a limited number of cell types that often do not include primary cell lines most relevant to disease modeling and the development oftherapeutics. Furthermore, efforts to profile molecular consequences of miRNA and RBP binding, such as site-specific cleavage and degradation across mammalian transcriptomes, have been limited [27, 35]. SUMMARY

[0010] A key unmet need in genetic medicine is the ability to predict post- transcriptional gene regulation across multiple tissue types based on a given mRNA sequence. The present embodiments provide improved systems and methods for predicting microRNA (miRNA) binding as well as mRNA degradation. The use of deep learning models trained on miRNA binding and / or degradome data from multiple tissue types allows for the prediction of tissue-specific differences in post-transcriptional gene regulation including the identification of cis-elements that regulate mRNA stability. Furthermore, the use of degradome and miRNA binding data in combination for training deep learning models captures additional complexities and improves the predictive capacity of post-transcriptional gene regulation.

[0011] Present embodiments describe a deep learning model that may be used to predict cell-type-specific microRNA binding and mRNA degradation directly from RNA sequence. In one aspect, the embodiments described herein use a dilated convolutional neural network for predicting mRNA stability based on miRNA binding and mRNA degradation in a sequence-to-sequence manner. Notably, the model has been demonstrated to reveal biological insights such as identifying repressive non-canonical miRNA target sites and decoding the regulatory effects of sequence context and miRNA binding site multiplicity. The model also demonstrated improvements relative to other advanced methods and neural architectures on a comprehensive suite of seven orthogonal tasks, including identifying genetic variants that affect microRNA binding, predicting out-of-distribution data from massively parallel reporter assays, and predicting canonical and non-canonical miRNA mediated repression. The models described herein can be used to make predictions regarding mRNA stability associated with sequence variants, provide insights into novel biology and be used in the design of RNA therapeutics.

[0012] As set out in the Examples, a detailed cell-type specific map of transcriptome wide miRNA binding and degradation events occurring in the cell was constructed usingstate-of-the-art NGS technologies, shedding light on the transcriptome-wide picture of miRNA and RBP interactions. Deep learning was used to investigate and build a cell-type specific model of PTGR (also referred to herein as “REPRESS”), that can accurately predict miRNA binding and mRNA degradation directly from sequence at a single bp resolution. The model captures non-linear complex relationships of miRNA and RBPs in a cell-type specific manner, learning to accurately predict cell type specific miRNA binding and degradation at a single-bp resolution. The model has been demonstrated to accurately capture validated aspects of biology and exhibits robustness to synthetic and out-of-distribution sequences. The model can also be used to aid in diagnostics and target discovery, including but not limited to predicting causal variants implicated in PTGR-related diseases. The model may also be used for investigating various therapeutic modalities for gene regulation like oligonucleotides, RNA editing and synthetic (e.g. non-naturally occurring) mRNA design.

[0013] As set out in the Examples and as shown in FIGs. 6G, 7A-7E, the present PTGR model provided improves upon the training and predictive capabilities found in conventional machine learning models. The present embodiments thus provide a computer- based solution to a computer-based problem.

[0014] In a first aspect, there is provided a computer-implemented method for predicting tissue-specific microRNA (miRNA) binding to a messenger RNA (mRNA). In one embodiment, the method comprises: providing, in a memory, a post-transcriptional gene regulation (PTGR) model comprising a base and at least one head; receiving, at a processor in communication with the memory, an RNA input sequence of length L corresponding to the mRNA; and determining, at the processor, an miRNA binding prediction matrix from the at least one head of the PTGR model, the miRNA binding prediction matrix comprising an L x K matrix comprising miRNA binding predictions at each position of the L nucleotides in the RNA input sequence for K tissue types, the miRNA binding prediction matrix determined using the RNA input sequence as input at the base of the PTGR model.

[0015] In one or more embodiments, the methods described herein comprise predicting tissue-specific degradation of the mRNA. For example, in one embodiment the PTGR model comprises at least two heads: a first head providing a first miRNA binding prediction matrix binding for K tissue types, optionally K human tissue types; and a secondhead providing a first mRNA degradome prediction matrix for K tissue types, optionally K human tissue types, the first mRNA degradome prediction matrix determined using the RNA input sequence as input at the base of the PTGR model. In one or more embodiments, the first mRNA degradome prediction matrix comprises an L x K matrix comprising degradome read coverage predictions at each position of the L nucleotides in the RNA input sequence for K tissue types.

[0016] The systems and methods described herein may also be advantageously trained on data from one or more different organisms such as human and mouse. In one or more embodiments, the at least two heads may comprise at least four heads: a third head providing a second miRNA binding prediction matrix for tissue types for a second organism, optionally mouse; and a fourth head providing a second mRNA degradome prediction matrix for tissue types for the second organism, optionally mouse.

[0017] In one or more embodiments, the PTGR model may comprise a convolutional neural network. For example, in one or more embodiments, the base may comprise: at least one base convolution layer; and at least one base residual block receiving an output from the at least one base convolution layer. In one embodiment, the base residual block comprises gated convolutions.

[0018] In one or more embodiments, each head of the at least one head may comprise a first head convolutional layer receiving input from the base, at least one head residual block receiving an output of the first head convolutional layer, and a second head convolutional layer receiving an output of the at least one head residual block. In one embodiment, the at least one head residual block does not comprise gated convolutions.

[0019] In one or more embodiments, each base residual block may comprise: a normalization layer receiving an input to the base residual block; a first pair of convolutional layers receiving the output from the normalization layer, a first of the first pair having a tanh activation and a second of the first pair having a sigmoid activation; a first element-wise multiplication operation on outputs of the first pair of convolutional layers; a second pair of convolutional layers receiving the output from the first element-wise multiplication operation, a first of the second pair having a tanh activation and a second of the second pair having a sigmoid activation; a second element-wise multiplication operation on outputs of the secondpair of convolutional layers; and an output of the base residual block determined by adding an output of the second element-wise multiplication operation and the input to the base residual block in a skip connection.

[0020] In one or more embodiments, each of the at least one head residual block may comprise: a normalization layer receiving an input to the head residual block; a first convolutional layer receiving the output from the normalization layer, the first convolutional layer having a gelu activation; a second convolutional layer receiving the output from the first convolutional layer, the second convolutional layer having a gelu activation; an output of the head residual block determined by adding an output of the second convolutional layer and the input to the head residual block in a skip connection.

[0021] In one or more embodiments, the method may further comprise: outputting, at a display device in communication with the processor, at least one of the miRNA binding prediction matrix and the mRNA degradome prediction matrix.

[0022] In one or more embodiments, the method may further comprise: receiving, at an input device in communication with the processor, one or more user selected nucleotides in the RNA input sequence; generating a first visualization comprising the miRNA binding predictions at the one or more user selected nucleotides in the RNA input sequence; and outputting, at the display device in communication with the processor, the generated first visualization.

[0023] In one or more embodiments, the method may further comprise: generating a second visualization comprising the mRNA degradome prediction matrix at the one or more user selected nucleotide in the RNA input sequence; and outputting, at the display device, the generated second visualization.

[0024] In one or more embodiments, the method may further comprise: converting, at the processor, the RNA input sequence of length L to an L x 4 input matrix, wherein the L x 4 input matrix is one-hot encoded.

[0025] In one or more embodiments, the miRNA binding predictions may be indicative of non-specific miRNA binding to the mRNA.

[0026] In one or more embodiments, the K tissue types may comprise at least one of liver, CNS, heart, kidney, muscle, cancer cells, HEK cells, HeLa cells, iPSCs, primary human hepatocytes (PHHs) and peripheral blood mononuclear cells (PBMCs).

[0027] In a second aspect, there is provided a system for predicting tissue-specific microRNA (miRNA) binding to a messenger RNA (mRNA). In one embodiment, the system comprises: a memory comprising: a post-transcriptional gene regulation (PTGR) model comprising a base and at least one head; a processor in communication with the memory, the processor configured to: receive an RNA input sequence of length L corresponding to the mRNA; and determine an miRNA binding prediction matrix from the at least one head of the PTGR model, the miRNA binding prediction matrix comprising an L x K matrix comprising binding predictions at each position of the L nucleotides in the RNA input sequence for K tissue types, the miRNA binding prediction matrix determined using the RNA input sequence as input at the base of the PTGR model.

[0028] In one or more embodiments, the processor may be further configured to predict tissue-specific degradation of the mRNA and the at least one head comprises at least two heads: a first head providing a first miRNA binding prediction matrix binding for K tissue types, optionally K human tissue types; and a second head providing a first mRNA degradome prediction matrix for K tissue types, optionally K human tissue types, the first mRNA degradome prediction matrix determined using the RNA input sequence as input at the base of the PTGR model.

[0029] In one or more embodiments, the first mRNA degradome prediction matrix may comprise an L x K matrix comprising degradome read coverage predictions at each position of the L nucleotides in the RNA input sequence for K tissue types.

[0030] In one or more embodiments, the at least two heads may comprise at least four heads: a third head providing a second miRNA binding prediction matrix for tissue types for a second organism, optionally mouse; and a fourth head providing a second mRNA degradome prediction matrix for tissue types for the second organism, optionally mouse.

[0031] In one or more embodiments, the PTGR model may comprise a convolutional neural network.

[0032] In one or more embodiments, the base may comprise at least one base convolution layer; and at least one base residual block receiving an output from the at least one base convolution layer. In one embodiment, the base residual block comprises gated convolutions.

[0033] In one or more embodiments, each head of the at least one head may comprise a first head convolutional layer receiving input from the base, at least one head residual block receiving an output of the first head convolutional layer, and a second head convolutional layer receiving an output of the at least one head residual block.

[0034] In one or more embodiments, each base residual block may comprise: a normalization layer receiving an input to the base residual block; a first pair of convolutional layers receiving the output from the normalization layer, a first of the first pair having a tanh activation and a second of the first pair having a sigmoid activation; a first element-wise multiplication operation on outputs of the first pair of convolutional layers; a second pair of convolutional layers receiving the output from the first element-wise multiplication operation, a first of the second pair having a tanh activation and a second of the second pair having a sigmoid activation; a second element-wise multiplication operation on outputs of the second pair of convolutional layers; and an output of the base residual block determined by adding an output of the second element-wise multiplication operation and the input to the base residual block in a skip connection.

[0035] In one or more embodiments, each of the at least one head residual block may comprise: a normalization layer receiving an input to the head residual block; a first convolutional layer receiving the output from the normalization layer, the first convolutional layer having a gelu activation; a second convolutional layer receiving the output from the first convolutional layer, the second convolutional layer having a gelu activation; an output of the head residual block determined by adding an output of the second convolutional layer and the input to the head residual block in a skip connection.

[0036] In one or more embodiments, the system may further comprise: a display device in communication with the processor, and the processor may be further configured to: output to the display device at least one of the miRNA binding prediction matrix and the mRNA degradome prediction matrix.

[0037] In one or more embodiments, the system may further comprise: an input device in communication with the processor; and the processor may be further configured to: receive from the input device one or more user selected nucleotides in the RNA input sequence; generate a first visualization comprising the miRNA binding predictions at the one or more user selected nucleotides in the RNA input sequence; and output to the display device in communication with the processor, the generated first visualization.

[0038] In one or more embodiments, the processor may be further configured to: generate a second visualization comprising the mRNA degradome prediction matrix at the one or more user selected nucleotide in the RNA input sequence; and output to the display device, the generated second visualization.

[0039] In one or more embodiments, the processor may further configured to: convert the RNA input sequence of length L to an L x 4 input matrix, wherein the L x 4 input matrix is one-hot encoded.

[0040] In one or more embodiments, the miRNA binding predictions may be indicative of non-specific miRNA binding to the mRNA.

[0041] In one or more embodiments, the K tissue types may comprise at least one of liver, CNS, heart, kidney, muscle, cancer cells, HEK cells, HeLa cells, iPSCs, primary human hepatocytes (PHHs) and peripheral blood mononuclear cells (PBMCs).

[0042] In a third aspect, there is provided a method determining an effect of a variant mRNA sequence relative to a control mRNA sequence on post transcriptional gene regulation. In one embodiment, the method comprises: predicting tissue specific miRNA binding to the variant mRNA sequence, and optionally tissue specific degradation of the variant mRNA sequence, according to a method described herein; comparing the miRNA binding prediction matrix, and optionally the mRNA degradome prediction matrix, for the variant mRNA sequence to a miRNA binding prediction matrix, and optionally a mRNA degradome prediction matrix, for the control mRNA sequence; and determining the effect of the variant mRNA sequence on post transcriptional gene regulation based on any differences between the miRNA binding prediction matrix, and optionally the mRNA degradome prediction matrix, for the variant mRNA sequence, relative to the miRNA binding predictionmatrix, and optionally the mRNA degradome prediction matrix, for the control mRNA sequence. In one embodiment, an increase in the miRNA binding prediction values for a variant sequence averaged across a window of 1, 2, 5, 10, 15, 20 or 25 bp relative to the corresponding miRNA binding prediction values for a reference sequence is indicative of decreased mRNA stability.

[0043] In one or more embodiments, the method may further comprise determining the miRNA binding prediction matrix, and optionally the mRNA degradome prediction matrix, for the control mRNA sequence according to a method described herein. Alternatively, the miRNA binding prediction matrix, and optionally the mRNA degradome prediction matrix, for the control mRNA sequence may be a set of predetermined values.

[0044] In one or more embodiments, the variant mRNA sequence may have between 1 and 20 single nucleotide polymorphisms relative to the control mRNA sequence.

[0045] In one or more embodiments, the variant mRNA sequence may have a variant of unknown clinical significance, a putative disease-causing mutation, a masked sequence corresponding to a SBO binding site, an ADAR editing site or a synthetic (non-naturally occurring) mRNA sequence.

[0046] In one or more embodiments, the method may further comprise synthesizing the variant mRNA and testing the variant mRNA for gene expression.

[0047] In a fourth aspect there is provided a system for determining an effect of a variant mRNA sequence relative to a control mRNA sequence on post transcriptional gene regulation, optionally on mRNA stability. In one embodiment, the system comprises: a memory comprising a PTGR model; a processor in communication with the memory, the processor configured to: predict tissue specific miRNA binding to the variant mRNA sequence using the PTGR model, and optionally tissue specific degradation of the variant mRNA sequence, according to a method described herein; compare the miRNA binding prediction matrix, and optionally the mRNA degradome prediction matrix, for the variant mRNA sequence to a miRNA binding prediction matrix, and optionally a mRNA degradome prediction matrix, for the control mRNA sequence; and determine the effect of the variant mRNA sequence on post transcriptional gene regulation based on any differences betweenthe miRNA binding prediction matrix, and optionally the mRNA degradome prediction matrix, for the variant mRNA sequence, relative to the miRNA binding prediction matrix, and optionally the mRNA degradome prediction matrix, for the control mRNA sequence.

[0048] In one or more embodiments, the processor may be further configured to determine the miRNA binding prediction matrix, and optionally the mRNA degradome prediction matrix, for the control mRNA sequence according to a method as described herein. Alternatively, the miRNA binding prediction matrix, and optionally the mRNA degradome prediction matrix, for the control mRNA sequence may be a set of predetermined values.

[0049] In one or more embodiments, the variant mRNA sequence may have between 1 and 20 single nucleotide polymorphisms relative to the control mRNA sequence. In one embodiment, the control mRNA sequence is a wild-type sequence.

[0050] In one or more embodiments, the variant mRNA sequence may have a variant of unknown clinical significance, a putative disease-causing mutation, a masked sequence corresponding to a SBO binding site, an ADAR editing site or a synthetic (non-naturally occurring) mRNA sequence, or the effect of ADAR editing.

[0051] While the predictive methods and systems described herein may be used independently to investigate post-transcriptional gene regulation, they may also be integrated into broader experimental platforms comprising in vitro or in vivo testing or corresponding molecules such as mRNA or modifications thereto. For example, in one embodiment, method further comprises synthesizing the variant mRNA and testing the variant mRNA for gene expression and / or for another biological activity.

[0052] In a fifth aspect there is provided a computer-implemented method for generating a post-transcriptional gene regulation (PTGR) model. In one embodiment, the method comprises: providing in a memory, a machine learning model comprising a base and at least two heads; providing, in the memory, a first data set comprising miRNA binding data for a plurality of mRNA sequences for a K plurality of tissues; training the machine learning model based on the first data set, wherein a prediction output of the machine learning model is a matrix of size L × K, K denoting a size of the K plurality of tissues, a k-th row of the matrix corresponding to a predicted miRNA binding probability each base pair for the k-th cell line.

[0053] In one or more embodiments, the miRNA binding data comprises mRNA sequence segments experimentally associated with miRNA binding.

[0054] In one or more embodiments, the method may further comprise: providing, in the memory, a second data set comprising mRNA degradome data for the plurality of mRNA sequences for the K plurality of tissues, wherein the training the machine learning model may further comprise training the machine learning model based on the first data set and the second data set.

[0055] In one or more embodiments, each of the mRNA sequences may comprise a degradome annotation, and the degradome annotation comprises a degradome matrix corresponding to degradome-seq data for each of the K plurality of tissue types.

[0056] In one or more embodiments, the mRNA degradome data may further comprise a read coverage value at each position of the mRNA sequence.

[0057] In one or more embodiments, the base and the at least two heads may comprise at least one residual block.

[0058] In one or more embodiments, the at least two heads may comprise a first head for providing a miRNA binding prediction matrix and a second head for providing a mRNA degradome prediction matrix and training the machine learning model may comprise training the base and the first head on the first data set and training the base and the second head on the second data set.

[0059] In one or more embodiments, the plurality of mRNA sequences may comprise a plurality of non-overlapping windows.

[0060] In one or more embodiments, the plurality of non-overlapping windows may be at least 3000 nucleotides long.

[0061] In one or more embodiments, each of the plurality of non-overlapping windows may comprise an mRNA context sequence about 500 bp to 10 kb, optionally about 6.25 kb around each non-overlapping window.

[0062] In one or more embodiments, a total input sequence length of each non- overlapping window may be 500 bp to 25 kb, optionally about 15.5 kb.

[0063] In one or more embodiments, the method may further comprise: converting each of the plurality of mRNA sequences to a one-hot encoded sequence.

[0064] In one or more embodiments, a loss function for the training the machinelearning model may comprise+ , wherein a binary cross-entropy losscomprises a miRNA binding output and a Poisson loss comprises a degradation output. Alternatively, the loss function for training the ML model may comprise other loss functions known in the art such as a mean square loss function.

[0065] In one or more embodiments, the machine learning model may comprise an ensemble model.

[0066] In one or more embodiments, the ensemble model may comprise at least four models, each of the four models trained based on a training data split across different sets of chromosomes.

[0067] In one or more embodiments, the training the machine learning model may comprise training the machine learning model using stochastic gradient descent, RMSprop or an Adam optimizer.

[0068] In a sixth aspect there is provided a system for generating a post-transcriptional gene regulation (PTGR) model, the system comprising a memory and a processor configured to perform any one of the methods described herein. DRAWINGS

[0069] FIGS.1A-F provide an overview and analysis of the data sets used for training the PTGR model. FIG. 1A: REPRESS framework for modeling post-transcriptional gene regulation. FIG. 1B: Summary of datasets used to train REPRESS. Top: Schematic illustration of AGO2-CLIP, miR-eCLIP, and Degradome-Seq assays used to generate REPRESS’s training dataset. Bottom: Example of a post-transcriptional gene regulation map of miRNA binding and mRNA degradation at the FAM3C 3’ UTR in HepG2 cells. FIG.1C: Relative degradome read coverage along protein-coding transcripts. Each row represents a transcript, each column corresponds to a segment of the 5’ UTR, CDS, or 3’ UTR. Valuesare scaled to the maximum read coverage within each transcript so that all transcripts can be visualized together (See Examples: Methods). Vertical bars on the left show the normalized average 3’ UTR length, number of 3’ UTR miRNA targets, and half-life of transcripts for each of the five major clusters obtained by hierarchical clustering of the data, depicted by the dendrogram. FIG.1D: The proportion of A549 miR-eCLIP peaks containing miRNA-specific seed matches, shown by individual miRNAs. Columns are ranked by the proportion of seed matches (including offset-m8 and offset-A1). FIG. 1E: Metaplots of average relative degradome read coverage at RBP eCLIP peak loci for FUBP3 (n = 1,021 for HepG2, no data for K562) and TIA1 (n = 479 for HepG2, n = 1,021 for K562), targets of conserved miRNAs (n = 11,585 for HepG2, n = 26,798 for K562), and UUAUUUAUU ARE motifs (n = 692 for HepG2, n = 643 for K562). Dashed lines indicate degradome read counts around control sites selected randomly from the same 3’ UTRs. Gray vertical bars represent the estimated location of the miRNA target / RBP binding site. To enable inter-cell line and target-control comparisons, degradome values for each dataset are scaled such that the minimum degradome value in the + / - 100 base window is 1.0. FIG. 1F: LinearCoPartition analysis of the rate of miRNA-target base pairing along the entire mature miRNA. Each row represents the average number of times the base was predicted to be paired to an mRNA across all targets of a single miRNA. K-means clustering was done with Kmeans++ initialization.

[0070] FIGS.2A-G provide an overview and analysis of various embodiments of the PTGR model described herein, also known as REPRESS. FIG.2A: REPRESS’s ConvNeXt inspired architecture and prediction pipeline. REPRESS takes RNA sequence as input along with bases from a 12.5 kb context window around the query sequence. The one-hot encoded sequence is then passed into a dilated convolutional neural network architecture with residual blocks and skip connections. Information from selected base layers is pooled and passed to specific output layers for generating the miRNA and degradome predictions respectively. FIG.2B: Heatmap comparing the performance of REPRESS to other popular miRNA binding models over a comprehensive suite of miRNA related tasks including predicting miRNA binding sites from miRTarBase, HEAP, miRAW and DeepMirTar test set; predicting miRNA mediated repression of synthetic sequences from the McGeary and Slutskin MPRA andpredicting the consequences of miRNA binding altering variants. * indicates non-orthogonal test data, i.e., that test data came from the same experimental setup as the training data. FIG.2C: REPRESS’s miRNA prediction on the 3’ UTR of NOTCH1. The figure illustrates one miRNA binding site of hsa-miR-30c-5p in REPRESS’s training dataset, but the model predicts four additional peaks not present in the training dataset. Two of these corresponded to the binding sites of hsa-miR-34a-5p and hsa-miR-144-3p which are validated binding sites from miRTarBase. The sequence attribution methods on the miRNA prediction revealed the identity of the seed sequence of the miRNA bound resulting in the predicted peaks. FIG.2D: REPRESS’s degradome predictions on the 3’ UTR of ENO2. The strongest predicted peak corresponded to an ARE site identified in the ARED-Plus dataset. Sequence attribution methods confirmed that the repeating overlapping ATTTA motif is driving the prediction of the strongest degradation peak. FIG. 2E: REPRESS learns underlying biology governing canonical and non-canonical miRNA mediated repression for 68 engineered target sites in a fixed sequence context with a Pearson R=0.89. REPRESS score is defined as the average REPRESS miRNA prediction over the entire query sequence. FIG.2F: REPRESS decodes the effect of miRNA binding altering variants and accurately predicts miRNA mediated repression for mutated and novel out-of-distribution sequences from the Slutskin MPRA (n=12,545) with a Pearson R=-0.62. REPRESS score is defined as the average K562 miRNA prediction across the entire MPRA construct. FIG.2G: REPRESS can efficiently design up- regulating ASOs ranking a validated miRNA target site blocking ASO from Sengupta et al., 1 out of the 2,037 possible ASOs that could target the 3’ UTR of the UTRN gene.

[0071] FIGS.3A-F provide an overview and analysis of the predictive capabilities of the PTGR model described herein, also known as REPRESS. FIG. 3A: Left: Sequence attribution base heights show the relative impact on miRNA binding that REPRESS assigned to each position (represented in terms of the corresponding miRNA position / base). Right, top: Proportion of miR-eCLIP peaks containing offset-6mers for hsa-let-7a-5p and hsa-miR- 148a-3p. P value is the result of a two proportion Z-test comparing the proportion of offset 6mer sites between hsa-let-7a-5p and hsa-miR-148a-3p. ”Other canonical” includes 6mer, 7mer-A1, 7mer-m8, and 8mer sites. Bottom: Conservation is the proportion of seed bases with PhyloP 100-way > 3. Different letters indicate statistically significant groupings of sitetypes (P < 0.05) after pairwise two proportion Z tests and Benjamini-Hochberg adjustment for multiple hypothesis testing. FIG.3B: Left: miRNAs are ranked by the odds ratio for non- seed region conservation rates in 6mer targets relative to 8mer targets. P values and confidence intervals (CI) are the result of Benjamini-Hochberg adjusted two-sided Fisher Exact Tests. Right: REPRESS sequence attributions for 6mer and 8mer targets of highly expressed miRNAs with substantial differences in non-seed region conservation between 6mer and 8mer targets. FIG.3C: REPRESS learns the importance of the local AU content of the dinucleotides flanking 8mer targets. For each model in the line plot, the mean score across all sites is normalized to sites with no A or U bases in either of the dinucleotides flanking positions 1-8. FIG. 3D: REPRESS predicts the impact of miRNA target site multiplicity on repression. Experimental repression (log2 fold-change) values and corresponding REPRESS predictions (average score across sequence) are shown for 3’ UTRs designed to have a variable number of miRNA targets embedded in two different sequence contexts. Examples of two miRNAs for which increasing the number of target sites differentially impacts observed repression are shown. FIG.3E: CDF plots showing the effect of miRNA modulation on transcript fold-change (log2) in two experiments. Left: hsa-miR-122- 5p transfection. Right: mmu-miR-26a-5p knockout. In both panels, the top 100 ranked transcripts by REPRESS (blue line) and TargetScan (orange line) are compared to the baseline of all transcripts analyzed (dashed black line, n = 19,195 for hsa-miR-122 experiment, n = 20,810 for mmu-miR-26a experiment). P values are the result of two-sided KS tests between REPRESS and TargetScan sets for each experiment. FIG.3F: Boxplots for relative REPRESS degradation predictions (normalized to median value for ARED-Plus cluster type = -1 are shown above boxplots for experimentally determined relative stability values for sequences stratified by ARED-Plus cluster type. The clusters are based on the number of overlapping AUUUA pantamers a sequence contains, with a higher number generally corresponding to more degradation and decreased stability. Horizontal lines indicate (from bottom to top) the first quartile, second quartile (median), and third quartile. Whiskers extend 1.5 × IQR below Q1 and 1.5 × IQR above Q3, where ”IQR” = interquartile range. ”ARED”: AU-Rich Element Database.

[0072] FIGS.4A-F provide an overview and analysis of the performance of the PTGR model described herein, also known as REPRESS. FIGs.4A, 4B: REPRESS’s wild type and mutant predictions on verified miRNA binding site altering variants. The FIG. 4A variant rs1876439052 disrupts the binding of hsa-miR-29b-3p and the FIG. 4B variant rs1063320 creates a hsa-miR-148a-3p binding site. (c): REPRESS’s ROC curve for predicting verified miRNA altering variants (n=100) from background variants (n=174) with high allele frequencies in gnomAD. Shaded area represents a 95% confidence interval resulting from bootstrapping. The diagonal dotted line represents the performance of a random classifier. FIG.4D: REPRESS variant scores on miRNA mediated P / LP variants, other P / LP variants and putative benign variants. REPRESS score for a variant is defined as the absolute difference between the mean REPRESS miRNA prediction over the wild type and mutant sequence. Horizontal lines indicate (from bottom to top) the first quartile, second quartile (median), and third quartile. Whiskers extend 1.5×IQR below Q1 and 1.5 × IQR above Q3, where ”IQR” = interquartile range. P values are the result of Mann-Whitney Wilcoxon (”M.W.W”) tests. FIG.4E: Scatter plot illustrating the relationship between REPRESS scores and log2 fold repression from the McGeary MPRA (n=952). Each point represents a synthetic sequence, with REPRESS scores on the x-axis and MPRA log2 fold repression on the y- axis. REPRESS accurately predicts miRNA mediated repression for synthetic sequences from the McGeary MPRA with Pearson R=0.68. REPRESS score is defined as the average miRNA binding prediction over the entire MPRA sequence. FIG. 4F: Performance comparison of different miRNA models on predicting miRNA mediated repression from MPRA datasets. REPRESS significantly outperforms other baseline miRNA models for synthetic sequences from the McGeary and Slutskin MPRAs.

[0073] FIGS.5A-F provide an overview and analysis of the performance of the PTGR model described herein, also known as REPRESS, for therapeutic development. FIG.5A: RREPRESS aids the design of novel antisense oligonucleotides that sterically block miRNA binding sites. The binding of the ASO is simulated by replacing the wild-type target sequence of the ASO with N’s and the ASO’s activity is ranked by the predicted drop in miRNA binding. To evaluate the performance of REPRESS, we curated 27 miRNA-blocking ASOs that up- regulate their target gene expression by 1.5x and internally screened 485 ASOs targeting the3’ UTR of PON1. FIG.5B: REPRESS’s miRNA predictions over the 3’ UTR of CHD9. The validated upregulating ASO perfectly corresponds to the strongest predicted REPRESS peak and is ranked at the 99.93th percentile (rank 2) out of all possible ASOs (2,606) targeting the 3’ UTR of CHD9. FIG.5C:Fraction of literature-curated ASOs captured by each model as a function of the oligo screening budget (x-axis, representing the number of ASOs screened). For each screening budget, the y-axis indicates the proportion 27 validated hits for which the lead ASO is identified using the different predictors to rank the ASOs. FIG.5D: REPRESSoutperforms the next- -regulation of PON1)ASOs from the PON1 dataset. The plot illustrates the fraction of PON1 ASO hits captured after screening the top x% of ASOs nominated by REPRESS or TargetScan predictions. REPRESS requires 2.5×fewer ASOs to identify 80% of the hits identified during screening. FIG. 5E: Successive sequence edits nominated by REPRESS have higher stability as predicted by Saluki, which is trained on orthogonal mRNA half-life data. Edits were made within the native 3’ UTRs for 116 genes, with miRNA-nominated edits resulting in a 44.6% average stability increase and degradome-nominated edits yielding a 20.6% average stability increase. FIG.5F: REPRESS’s miRNA prediction scores drop by 97.5% after incorporating 10 edits nominated by the miRNA predictions for ULK1.

[0074] FIGs.6A – 6G shows experimental performance evaluation of the REPRESS model as described in the Examples. FIG. 6A shows a performance comparison of REPRESS and TargetScan in identifying miRNA binding sites in each of the cell lines from the validation set of our miRNA CLIP dataset. REPRESS outperforms TargetScan across all cell lines (shown as dots in scatter plot) with higher AUROCs. FIG.6B shows a validation set performance of REPRESS on the Degradome-seq datasets. Each point corresponds to the performance on the validation set of each fold done for the 4-fold cross validation. FIG.6C shows performance of REPRESS in predicting miRNA binding sites validated in miRTarBase and DianaTarBase. This shows that the model is able to generalize beyond the CLIP dataset it was trained on. FIG. 6D shows a plot illustrating the importance of context sequence in predicted miRNA binding score for let-7 miRNA binding site in the 3’ UTR of UTRN. Each line indicates the predicted binding track after each incorporated non-seed edit. The edits were found using ISM within the context region of the let-7c binding site. After 30 non-seededits the predicted miRNA binding score drops by >50%. FIG.6E: Track plot illustrating the REPRESS predicted miRNA binding track for the wild type sequence and the sequence after 30 non-seed edits, showing the location of the corresponding edits. The furthest edit is 140 bp away from the seed region of the miRNA binding site. FIG.6F: The performance of both miRNA binding and degradome tasks of REPRESS increase with the model size (number of parameters). This is consistent with scaling laws observed in deep learning architectures for sequence modeling. Additionally, jointly training on miRNA and degradome tasks increases generalization performance for both but only for the larger models. FIG. 6G: Performance comparison of other popular ML architectures when trained on the REPRESS dataset. All the architectures were allowed to have 34M parameters. The Examples showed the ConvNeXt inspired REPRESS architecture outperforms all existing popular ML architectures. The miRNA binding performance of REPRESS is 10% better than the second best Mamba State-Space architecture and 17% better than a dilated CNN with similar number of parameters and the same 12.5 kb context window.

[0075] FIGs. 7A-7E and FIG. 8 show experimental comparative model performance as demonstrated in the Examples. Model performance values used in the heatmap for each of the individual datasets / tasks used for validation FIG. 7A-7E. The * indicates that the dataset was used as the validation / test set for those corresponding models, hence, is from the same distribution as the data used to train that model. FIG.7A: Performance comparison on the REPRESS validation set using the average AUPRC across the 29 different cell lines. REPRESS significantly outperforms all baseline models. FIG.7B: Performance comparison on the test set from the miRAW model. The performance of REPRESS is close to miRAW and miTAR which were trained on data from an identical distribution. FIG.7C: Performance comparison on the test set from the DeepMirTar model. The performance of REPRESS is close to DeepMirTar and miTAR which were trained on data from an identical distribution. miRAW fails to generalize on the DeepMirTar test set and vice-versa. FIG. 7D: Average AUPRC across 8 cell lines for predicting experimentally validated miRNA binding sites from miRTarBase and DianaTarBase. FIG. 7E: REPRESS shows a significant performance improvement over all baseline models in predicting miRNA binding sites from non-binding sites identified from the HEAP assay. HEAP is an orthogonal assay different from the CLIPdata used to train the model which REPRESS accurately generalizes to. FIG. 8: HEAP identifies miRNA-mRNA interactions and quantifies their repressive effect via log2 fold change in target expression. REPRESS outperforms existing models in identifying repressive HEAP identified sites from non-repressive HEAP identified sites for different log2 fold change thresholds.

[0076] FIG. 9 shows the performance of the miRNA-specific model on the cell-line- specific validation set by combining the predictions from the top 50 miRNAs expressed in that cell line. The miRNA expression weighted average was taken across the predictions of the top 50 miRNAs to get the prediction track for that corresponding cell line. The miRNA specific model maintains 97% of the average AURPC as that of the cell-line-specific version of REPRESS.

[0077] FIG.10 shows a system diagram of a PTGR prediction system in accordance with one or more embodiments.

[0078] FIGs.11A – 11B show PTGR model inputs and outputs in accordance with one or more embodiments.

[0079] FIG.11C shows PTGR prediction output from the PTGR binding models from FIG.11A in accordance with one or more embodiments.

[0080] FIG. 12A – 12C show PTGR models in accordance with one or more embodiments.

[0081] FIG.13 shows a model architecture of the PTGR model in accordance with one or more embodiments.

[0082] FIG.14 shows a device drawing of the server 1006a of FIG.10 in accordance with one or more embodiments.

[0083] FIG.15 shows a method drawing of PTGR prediction in accordance with one or more embodiments.

[0084] FIG.16 shows a method drawing for determining an effect of a variant mRNA sequence relative to a control mRNA sequence on post transcriptional gene regulation.

[0085] FIG.17 shows a method drawing of generating a PTGR model in accordance with one or more embodiments. DESCRIPTION OF VARIOUS EMBODIMENTS

[0086] Various embodiments will now be described below to provide an example of the claimed subject matter. No example described below limits any claimed subject matter and any claimed subject matter may cover embodiments such as systems or methods that differ from those described below.

[0087] Furthermore, it will be appreciated that for simplicity and clarity of illustration, where considered appropriate, reference numerals may be repeated among the figures to indicate corresponding or analogous elements. In addition, numerous specific details are set forth in order to provide a thorough understanding of the examples described herein. However, it will be understood by those of ordinary skill in the art that the examples described herein may be practiced without these specific details. In other instances, well-known methods, procedures and components have not been described in detail so as not to obscure the examples described herein. Also, the description is not to be considered as limiting the scope of the examples described herein.

[0088] It should also be noted that, as used herein, the wording “and / or” is intended to represent an inclusive-or. That is, “X and / or Y” is intended to mean X or Y or both, for example. As a further example, “X, Y, and / or Z” is intended to mean X or Y or Z or any combination thereof.

[0089] It should be noted that terms of degree such as "substantially", "about" and "approximately" as used herein mean a reasonable amount of deviation of the modified term such that the end result is not significantly changed. These terms of degree may also be construed as including a deviation of the modified term if this deviation would not negate the meaning of the term it modifies.

[0090] Furthermore, the recitation of numerical ranges by endpoints herein includes all numbers and fractions subsumed within that range (e.g., 1 to 5 includes 1, 1.5, 2, 2.75, 3, 3.90, 4, and 5). It is also to be understood that all numbers and fractions thereof are presumed to be modified by the term "about" which means a variation of up to a certainamount of the number to which reference is being made if the end result is not significantly changed.

[0091] Some elements herein may be identified by a part number, which is composed of a base number followed by an alphabetical or subscript-numerical suffix (e.g., 112a, or 1121). Multiple elements herein may be identified by part numbers that share a base number in common and that differ by their suffixes (e.g., 1121, 1122, and 1123). All elements with a common base number may be referred to collectively or generically using the base number without a suffix (e.g., 112).

[0092] The example systems and methods described herein may be implemented in hardware or software, or a combination of both. In some cases, the examples described herein may be implemented, at least in part, by using one or more computer programs, executing on one or more programmable devices comprising at least one processing element, a data storage element (including volatile and non-volatile memory and / or storage elements), and at least one communication interface. These devices may also have at least one input device (e.g., a keyboard, a mouse, a touchscreen, and the like), and at least one output device (e.g., a display screen, a printer, a wireless radio, and the like) depending on the nature of the device. For example, and without limitation, the programmable devices (referred to below as computing devices) may be a server, network appliance, embedded device, computer expansion module, a personal computer, laptop, personal data assistant, cellular telephone, smart-phone device, tablet computer, a wireless device or any other computing device capable of being configured to carry out the methods described herein.

[0093] In some examples, the communication interface may be a network communication interface. In examples in which elements are combined, the communication interface may be a software communication interface, such as those for inter-process communication (IPC). In still other examples, there may be a combination of communication interfaces implemented as hardware, software, and a combination thereof.

[0094] Program code may be applied to input data to perform the functions described herein and to generate output information. The output information is applied to one or more output devices, in known fashion.

[0095] Each program may be implemented in a high-level procedural, declarative, functional or object-oriented programming and / or scripting language, or both, to communicate with a computer system. However, the programs may be implemented in assembly or machine language, if desired. In any case, the language may be a compiled or interpreted language. Each such computer program may be stored on a storage media or a device (e.g., ROM, magnetic disk, optical disc) readable by a general or special purpose programmable computer, for configuring and operating the computer when the storage media or device is read by the computer to perform the procedures described herein. Examples of the system may also be considered to be implemented as a non-transitory computer- readable storage medium, configured with a computer program, where the storage medium so configured causes a computer to operate in a specific and predefined manner to perform the functions described herein.

[0096] Furthermore, the example system, processes and methods are capable of being distributed in a computer program product comprising a computer readable medium that bears computer usable instructions for one or more processors. The medium may be provided in various forms, including one or more diskettes, compact disks, tapes, chips, wireline transmissions, satellite transmissions, internet transmission or downloads, magnetic and electronic storage media, digital and analog signals, and the like. The computer useable instructions may also be in various forms, including compiled and non-compiled code.

[0097] Various examples of systems, methods and computer programs products are described herein. Modifications and variations may be made to these examples without departing from the scope of the invention, which is limited only by the appended claims. Also, in the various user interfaces illustrated in the figures, it will be understood that the illustrated user interface text and controls are provided as examples only and are not meant to be limiting. Other suitable user interface elements may be used with alternative implementations of the systems and methods described herein.

[0098] Referring to FIG. 10 there is shown a system diagram 1000 of a PTGR prediction system, optionally a system for determining an effect of a variant mRNA sequence relative to a control mRNA sequence on post transcriptional gene regulation, or optionally a PTGR model generation system in accordance with one or more embodiments.

[0099] In one embodiment, the PTGR prediction system provides miRNA binding predictions and mRNA degradation predictions. In one embodiment, the PTGR prediction system provides mRNA stability predictions, optionally based on miRNA binding predictions and / or mRNA degradation predictions.

[0100] The PTGR prediction system 1000 includes one or more user devices 1002, a network 1004 and a computing device 1006.

[0101] The one or more user devices 1002 may be used by a user such as an administrator, geneticist or technician to access a software application (not shown) running on server 1006a at remote service 1006 over network 1004. In one embodiment, the one or more user devices 1002 may access a web application hosted at server 1006a using a browser for determining outputs such as cell-type specific miRNA binding and mRNA degradation for therapeutic design, diagnostics or other purposes . For example, the clinician user at device 1002 may provide a list of one or more mRNA sequences to the server 1006a and request a binding prediction and / or degradation prediction of all possible SBOs of a given length (e.g. ~20 bp) in the 3’ UTR of the one or more mRNAs. In an alternate embodiment, the one or more user devices 1002 may download an application for predicting miRNA binding and / or mRNA degradation and may connect to the server 1006a using an API.

[0102] In an alternate embodiment, the one or more user devices 1002 may operate the remote service 1006 (either using the web application or via a downloaded application) to generate an miRNA binding prediction or degradation prediction, optionally a prediction of mRNA stability. This may include generating a PTGR model as described herein, such as a method as shown in FIG.17.

[0103] The one or more user devices 1002 may be any two-way communication device with capabilities to communicate with other devices. A user device 1002 may be a desktop computer, mobile device, or laptop computer. A user device 1002 may be a mobile device such as mobile devices running the Google® Android® operating system or Apple® iOS® operating system. A user device 1002 may be the personal device of a user, or may be a device provided by an employer.

[0104] The one or more user devices 1002 may be used by an end user to access the software application running on server 1006a over network 1004. In one embodiment, the one or more user devices 1002 may access a web application hosted at server 1006a and may allow a user to review outputs such as an miRNA binding prediction or degradation prediction in a database at data store 1006b, including historical miRNA predictions and historical degradation predictions.

[0105] The user at the one or more user devices 1002 may send or submit a nucleic acid sequence such as an mRNA nucleic acid sequence to be analyzed to the PTGR system running on the remote service 1006. The target nucleic acid sequence may be sent in a miRNA binding prediction or similar request. The miRNA binding prediction request may be a web application request, an Application Programming Interface (API) request, or another request. The nucleic acid target sequence may be provided in a variety of formats in the miRNA binding prediction request. For example, the nucleic acid target sequence may be provided by the user, optionally in plain sequence format, FASTQ format, EMBL format, or FASTA format, or in other similar formats containing nucleic acid sequence information. Alternatively, the nucleic acid target sequence may be manually entered through user device 1002 or by entering or referencing an accession number or database entry corresponding to the sequence for a nucleic acid target sequence. Upon receipt of the miRNA binding prediction request, the miRNA binding prediction system on remote service 1006 may use the methods described herein in order to determine a miRNA binding prediction, and transmit the miRNA binding prediction to the user at user device 1002. Alternatively, or in addition to this, the PTGR system may use the methods described herein for determining an effect of a variant mRNA sequence relative to a control mRNA sequence on post transcriptional gene regulation or on mRNA stability. The miRNA binding prediction response may be provided to the one or more user devices 1002 in an email, as a notification, or text message.

[0106] The input to the miRNA binding prediction system may be a gene sequence such as an mRNA sequence of length L. The output of the miRNA binding prediction system may be a matrix of L x K binding probabilities where L is the length of the input mRNA sequence and K is a number of cell lines (or tissue types) and K total number of cell lines. Each output prediction for position (l,k) in the output matrix may be the probability of anmiRNA binding at position l for cell type k or a predicted degradome read coverage of each base pair for the kth cell line. The input to the miRNA binding prediction model may be a matrix encoded in a one-hot representation. For an input gene sequence with length L, its representation may become an L x 4 matrix. Each column in this matrix may denote the presence of nucleotides A, G, C, or T, with a value of 1 indicating the occurrence of the respective nucleotide and 0 elsewhere.

[0107] The software application running on the one or more user devices 1002 may display one or more user interfaces on a display device of the user device. This may include one or more prediction interfaces including an miRNA binding prediction matrix and an mRNA degradome prediction matrix. This may further include one or more visualizations, where one or more user selected nucleotides in the RNA input sequence are selected using a user input device of the user device 1002, and a first visualization is generated including the miRNA binding predictions at the one or more user selected nucleotide(s) in the RNA input sequence; and the interface is output to the display device in communication with the processor. Alternatively, or in addition to this first interface, a second visualization may be generated comprising the mRNA degradome prediction matrix at the one or more user selected nucleotides in the RNA input sequence; and the second interface output to the display device.

[0108] Network 1004 may be any network or network components capable of carrying data including the Internet, Ethernet, fiber optics, satellite, mobile, wireless (e.g. Wi-Fi, WiMAX), SS7 signaling network, fixed line, local area network (LAN), wide area network (WAN), a direct point-to-point connection, mobile data networks (e.g., Universal Mobile Telecommunications System (UMTS), 3GPP Long-Term Evolution Advanced (LTE Advanced), Worldwide Interoperability for Microwave Access (WiMAX), etc.) and others, including any combination of these.

[0109] The remote service 1006 may include a server 1006a and a database 1006b. The server 1006a may further be in communication with a database 1006b. The database 1006b and the server 1006a may be provided on the same server device, may be configured as virtual machines, or may be configured as containers. The database 1006b may be provided by server 1006a, another server (not shown), or a cloud-based database service such as Amazon Web Services ® (AWS). The remote service 1006 may itself be providedby a service such as AWS. The remote service 1006 is in network communication with the one or more user devices 1002.

[0110] The server 1006a may host a web application or an Application Programming Interface (API) endpoint that the one or more user devices 1002 may interact with via network 1004. The server 1006a may make calls to the database 1006b to query miRNA datasets or other data. The requests made to the API endpoint of server 1006a may be made in a variety of different formats, such as JavaScript Object Notation (JSON) or eXtensible Markup Language (XML).

[0111] The database 1006b may store information including miRNA and / or degradome datasets as described herein, user data, PTGR model data including pre-trained models. The database 1006b may be a Structured Query Language (SQL) such as PostgreSQL or MySQL or a not only SQL (NoSQL) database such as MongoDB.

[0112] In one embodiment, the server 1006a may perform model training of the PTGR models, such as according to a method described in FIG.17. The model training by server 1006a may include querying the database 1006b for miRNA data and mRNA degradome data for use in model training as described herein.

[0113] Referring next to FIGs. 11A – 11B there are shown PTGR model inputs and outputs in accordance with one or more embodiments.

[0114] FIG.11A shows a first PTGR model trained to predict cell-line (or tissue type) specific miRNA binding. Notably, the model may be trained to predict generalized or non- specific miRNA binding to a target mRNA, in contrast to models trained to predict the binding of specific identified miRNAs to one or more mRNAs.

[0115] As described herein, a nucleic acid sequence 1102 of length L is received at PTGR model 1104. The nucleic acid sequence 1102 may be an mRNA sequence. The input sequence 1102 may be one-hot encoded. The input sequence 1102 may be 200 to 20000 nucleotides in length or more.

[0116] The PTGR model 1104 may create as output K sequences of length L, each k output sequence of the K sequences may correspond to the probability of miRNA binding on the mRNA sequence 1102 in the kth tissue type. For example, the K output sequences mayinclude a first cell line 1106, a second cell line 1108, and a third cell line 1110. The PTGR model 1104 may be as described in further detail in FIG.13.

[0117] The first cell line 1106 sequence output may predict the binding of miRNA at each lth position of the L positions in the mRNA sequence for the first cell line 1106.

[0118] The second cell line 1108 sequence output may predict the binding of miRNA at each lth position of the L positions in the mRNA sequence for the second cell line 1108.

[0119] The third cell line 1110 sequence output may predict the binding of miRNA at each lth position of the L positions in the mRNA sequence for the third cell line 1110.

[0120] Each of the output cell lines 1106, 1108 and 1110 may correspond to a head of a machine learning model as described further in FIG.13.

[0121] FIG 11B shows an alternate second PTGR model trained to predict the binding of specific miRNAs to an mRNA. As set out in Example 6, the PTGR model may be configured to output binding data for a plurality of specific miRNAs rather than for a plurality of tissues.

[0122] The second PTGR model 1154 may receive the nucleic acid sequence 1102 of length L (see e.g. as described in FIG.11A).

[0123] The second PTGR model 1154 may create as output K sequences of length L, each k output sequence corresponding to the probability of miRNA binding on the mRNA sequence 1102 of the kth miRNA type in a plurality K of miRNA sequences. For example, the K output miRNA types may include the first miRNA type 1156, the second miRNA type 1158 and third miRNA type 1160.

[0124] The first miRNA type 1156 sequence output may predict the binding of the first miRNA type 1156 at each lth position of the L positions in the mRNA sequence.

[0125] The second miRNA type 1158 sequence output may predict the binding of the second miRNA type 1158 at each lth position of the L positions in the mRNA sequence.

[0126] The third miRNA type 1160 sequence output may predict the binding of the third miRNA type 1160 at each lth position of the L positions in the mRNA sequence.

[0127] Referring next to FIG. 11C shows miRNA prediction output from the PTGR models from FIG. 11B in accordance with one or more embodiments. Specifically, the probability of binding of hsa-miR-20a, hsa-miR-124, hsa-miR-30a, and hsa-miR-9 miRNAs to a given input mRNA sequence.

[0128] Referring next to FIG.12A – 12C together, method diagrams 1200, 1230 and 1260 show embodiments of models with different number of heads and corresponding output sequences or matrices. The PTGR model may be a neural network 1204.

[0129] The outputs of the methods and systems described herein may be used for a variety of applications in genetic medicine including but not limited to ascertaining mRNA stability, diagnostics, to identify drug targets, as well as to predict the effect of modifying DNA or RNA (such as through ADAR dependent RNA editing) on mRNA stability and ultimately gene expression.

[0130] In one embodiment the present systems comprise structured computational architectures referred to herein as neural networks (NNs). NNs, also called artificial neural networks (ANNs), are a powerful class of architectures for applying a series of computations to an input to determine an output. The input to the NN is used to determine the outputs of a set of feature detectors, which are then used to determine the outputs of other feature detectors, and so on, layer by layer, until the output is determined. An NN architecture can be thought of as a configurable set of processors configured to perform a complex computation. The configuration is normally done in a phase called training, wherein the parameters of the NN are configured to maximize the computation’s performance on determining a desired output such as miRNA binding or, equivalently, to minimize the errors made on that task. Because the NN gets better at a given task throughout training, the NN is said to be learning the task as training proceeds. NNs can be trained using machine learning methods. Once configured, a NN can be deployed for use in the task for which it was trained and herein for predicting miRNA binding and / or mRNA degradation as described below.

[0131] In one embodiment, the models described herein are concurrently trained on a first data set comprising miRNA binding data and a second data set comprising degradome data. In one embodiment, the first data set is generated using cross-linking and precipitation (CLIP) techniques, optionally argonaute-CLIP (AGO-CLIP) techniques such as PAR-CLIP,HITS-CLIP, eCLIP or miR-eCLIP, including but not limited to the methods described in Example 1 and Figure 1A. In one embodiment, the second data set may be generated using PARE-SEQ techniques known in the art, including but not limited to the methods described in Example 1 and Figure 1A.

[0132] Neural networks are comprised of layers of feature detectors. The layers are ordered. The first layer is an input layer into which the inputs to the neural network are loaded. For example, the input layer may obtain a sequence such as an mRNA sequence represented as a vector sequence and optionally additional information. The last layer is the output layer, for example, the probability of miRNA binding at each position across the mRNA sequence in one or more tissue or cell types.

[0133] The head, base, and residual blocks of the neural network 1204 are described in further detail in FIG.13.

[0134] In one embodiment, the systems and methods described herein make use of NNs that are configured as a class of neural networks called convolutional neural networks (CNNs).

[0135] CNNs may be constructed to account for the complex relationships between biological sequences (such as mRNAs) and molecular phenotypes (such as tissue specific miRNA binding or mRNA degradation) that they may influence. Machine learning methods may be used to construct these computational models by extracting information from a dataset comprising measured or experimental data such as data describing miRNA binding to mRNA sequences in a tissue specific manner along with tissue specific degradation profiles of the corresponding mRNA sequences.

[0136] CNNs operate by: applying a set of convolutional filters (arranged as one or more convolutional layers) to the input sequence and applying non-linear activation functions to the outputs of the convolutional filters. These steps may be applied, recursively, to the feature map, by replacing the input sequence with the feature map, to obtain deeper feature maps. This may be repeated to obtain even deeper feature maps, and so on. At some point the output is obtained by applying a non-convolutional neural network to the deepest feature map.

[0137] In one embodiment, the convolutional filters in CNNs are shared across sequence positions and act as sequence feature detectors. The non-linear activation functions identify significant filter responses while repressing spurious responses caused by insufficient and often idiosyncratic matches between the filters and the input sequences.

[0138] It will be appreciated that the systems and methods described herein may make use of different variations of convolutional neural networks, including extensions such as recursive neural networks.

[0139] As described herein, a nucleic acid sequence 1202 of length L is received at a neural network 1204. The nucleic acid sequence 1202 may be an mRNA sequence, either a naturally occurring mRNA sequence or a synthetic mRNA sequence not typically found in nature. The input sequence 1202 may be one-hot encoded. The input sequence 1202 may be 200 to 20000 nucleotides in length or more. The outputs may be generated by one or more heads of a neural network 1204. The neural network 1204 may be a convolutional neural network (CNN).

[0140] As shown in FIG.12A, the neural network 1204 may have a single head 1206 producing tissue-specific probability of miRNA binding at each position of the input mRNA sequence 1202.

[0141] As shown in FIG.12B, the neural network 1204 may have two heads including a first head 1236 producing a tissue-specific probability of miRNA binding at each position of the mRNA sequence 1202 and a second head 1238 producing a tissue-specific indicator of degradation at each position of the input mRNA sequence 1202.

[0142] As shown in FIG.12C, the neural network 1204 may have four heads including a first head 1260 producing a human tissue-specific probability of miRNA binding at each position of the mRNA sequence 1202, a second head 1262 producing a mouse tissue- specific probability of miRNA binding at each position of the mRNA sequence 1202, a third head 1264 producing a human tissue-specific indicator of degradation at each position of the input mRNA sequence 1202, and a fourth head 1266 producing a mouse tissue-specific indicator of degradation at each position of the input mRNA sequence 1202. The heads 1260and 1264 may function to provide predictions for human tissue types and the heads 1262 and 1266 may function to provide predictions for mouse tissue types.

[0143] Referring next to FIG. 13 there is shown a model architecture of the PTGR model 1300 in accordance with one or more embodiments. The PTGR model 1300 may be, for example, REPRESS. The PTGR model 1300 may run on a server such as server 1006a (see e.g. FIG.10) and server 1400 (see e.g. FIG.14). The PTGR model 1300 may provide features including designing mRNAs with desired properties (such as increased or decreased mRNA stability), steric blocking oligonucleotides for manipulation of gene expression, and for predicting the impact of specific RNA edits. The PTGR model 1300 may be proficient in identifying disease-impacting variants and predicting the efficacy of steric blocking oligonucleotides (SBOs) in disrupting miRNA binding. This proficiency may highlight the PTGR model’s 1300 understanding of the underlying causal dynamics at play. Moreover, the PTGR model 1300 may have insights extending to the design of synthetic mRNA sequences, including optimizing them for increased half-life and stability, further evidencing its ability to navigate and manipulate the PTGR landscape. The PTGR model 1300 may generalize beyond naturally occurring sequences, applying its learned causal mechanisms to predict outcomes in novel contexts.

[0144] The PTGR model 1300 may consist of base layers 1304 and output layers (for example the first pair of heads 1314a and a second pair of heads 1314b), the core components of which are residual blocks 1308 composed of dilated convolution layers. All convolution layers may have the same number of filters, denoted by N, except for the last one in the output layers, which may be determined by the number of cell types to predict. Padding may be used in convolution layers so that the output sequence has the same length as the input sequence after convolution layers.

[0145] The PTGR model 1300 may be, for example, the neural network 1204 (see e.g. FIGs.12A – 12C). The model 1300 may include an input layer 1302 that represents an mRNA sequence, a base 1304, a concatenation layer 1310 which creates pooled head input 1312, a first pair of heads 1306a and a second pair of heads 1306b. The input layer 1302 may be processed through one or more convolutional layers before it is received at the base 1304. The input layer 1302 may receive the mRNA sequence in a 1-hot encoded format.

[0146] The input 1302 may first processed with an initial convolution layer and then may go through 16 residual blocks 1308 (RB) in base 1304 with gradually increasing convolutional dilation rate (receptive field). The dilation for each of the residuals blocks may be any arbitrary ordering of positive non-zero numbers. In a preferred embodiment, the dilation rate for each of the 16 residual blocks may be 1, 1, 1, 1, 4, 4, 4, 4, 10, 10, 10, 10, 25, 25, 25, 25.

[0147] The base 1304 may comprise a plurality of residual blocks. For example, 16 residual blocks may be used. At every fourth residual block 1308 in the base 1304, as well as the initial convolution layer, the base 1304 may have an extra convolution layer that will be concatenated 1310 as an input 1312 to the different heads 1306.

[0148] A residual block 1308 (or “RB” as indicated) is a stack of layers set in such a way that the output of a layer is taken and added to another layer deeper in the block. The non-linearity is then applied after adding it together with the output of the corresponding layer in the main path. This by-pass connection is known as the shortcut or the skip-connection.

[0149] The residual blocks of the model 1300 are described in further detail at 1308 and 1316. The residual blocks 1308 in the base 1304 and the residual blocks in the pairs of heads 1306a and 1306b may be configured differently. The residual blocks 1308 may have parameterized configurations of the embodied convolution layers including parameters related to the number of filters, the filter width, and the dilation rate.

[0150] In the residual block 1308, a layer normalization layer may first be applied, and then the output received by two convolution layers separately, the output of which goes through an element-wise multiplication operation as shown. The two convolution layers used have sigmoid and tanh activations respectively. After the element-wise multiplication is another two same layers of convolution followed by another element-wise multiplication. Finally, the input to the residual block 1308 is added to the output as a skip connection. The residual blocks in the base 1304 may have the same number of filters, and each of the four residual blocks may have filter width=11, 11, 21 and 41, and dilation=1, 4, 10, and 25.

[0151] The residual block 1308 configuration in the heads 1314 may be different than the base 1304, and may include a normalization layer, followed by two convolution layerswith gelu activation, the output of which is added to the input of the residual block 1308 as a skip connection. The four head residual blocks may have the same number of filters, which is the same as the ones in the base. However, they have filter width=11, 11, 21 and 41, and dilation=1, 4, 10, and 25 respectively.

[0152] In each residual block 1308, layer normalization may be applied first, the output of which may be passed through two separate convolution layers with sigmoid and tanh activation functions. Then the two outputs may go through an element-wise multiplication operation. The convolution and element-wise multiplication may be repeated twice. Finally, the input may be added to the output as a skip connection. All the residual blocks 1308 may have the same number of filters, and each set of four residual blocks has filter widths of 11, 11, 21 and 41. The output of the initial convolution layer, and the 4th, 8th, 12th and 16th of the residual blocks may then passed through a convolution layer separately. These five outputs may be accumulated, concatenated 1310 and passed into the output layers of the model.

[0153] In one embodiment, the PTGR residual blocks 1308 and 1316 (RB and RB-H respectively) may be changed to include replacing batch normalization layers with layer normalization and reducing the total number of normalization layers to only 1 per residual block instead of 2. As well, the GELU activation may be used in the residual blocks (RB-H) of the output layers instead of the commonly used ReLU activation function.

[0154] A residual block 1316 may be provided for the output layers 1306. In each residual block 1316 (RB-H) of the output module, a layer normalization layer may be applied first, followed by two convolution layers with GELU activation. Then the input to the residual block may be added to the output as a skip connection. Four residual blocks 1316 may be used in each of the heads in the first pair of heads 1306a and the each of the heads in the second pair of heads 1306b and may have the same number of filters as the ones in the base layers, but may have filter widths of 11, 11, 21 and 41, with dilations of 1, 4, 10, and 25, respectively. Before the final output convolution layer, the previous output may be cropped at both ends by half the context length (receptive field of the model) so that the final outputs 1314a and 1314b may have the same length as the initial input to the model 1300 without the added context sequence. For the last convolution layer, the number of filters may bebased on the number of cell lines in for the corresponding output, and the activation functions used may be sigmoid for miRNA binding prediction and softplus for degradation prediction.

[0155] A first pair of heads 1314a provide output matrices 1314a for human cell line miRNA binding prediction and mRNA degradome prediction.

[0156] A second pair of heads 1314b provide output matrices 1314b for mouse cell line miRNA binding prediction and mRNA degradome prediction.

[0157] The model 1300 has four heads based on output species and data type. That is, there is a first pair of heads 1314a that includes miRNA binding (ago-human) and miRNA binding (ago-mouse) and a second pair of heads 1314b that includes degradome-human and degradome-mouse. Each head has a similar structure, made up of a convolution layer and four residual blocks 1308, except for the last output convolution layer. Before the final output convolution layer, the previous output may be cropped at the two ends by half context length (receptive field of the model), so that the final output will have the same length L as the initial input to the model. For the last convolution layer, the number of filters may be based on the number of cell lines in for the corresponding output, and the activation function used may be sigmoid for miRNA binding (ago) and softplus for degradome.

[0158] The model 1300 may be a convolutional neural network and may comprise a layer of input values 1302 that represents an mRNA sequence (which may be referred to as an “input layer”), at least one set of convolutional layers and pooling layers 1310, the output 1312 of which are used as input to the one or more heads 1306 which generates output sequences or matrices 1314 that represent the computed relevance scores (which may be referred to as an “output layer”).

[0159] The PTGR model 1300 may be a dilated Convolutional Neural Network (CNN) model designed to predict species-specific miRNA binding across different cell types and RNA degradation across different cell types at single base pair resolution. Training of the PTGR model 1300 may involve randomly sampled regions from the human and mouse transcriptome spanning 15.5kb, with the model analyzing information within a 12.5kb context for each base pair prediction. This expansive contextual understanding allows the PTGRmodel 1300 to adeptly capture long-range dependencies and incorporate information across the various regions of a transcript.

[0160] The particular CNN 1300 shown in FIG. 13 is an example architecture; theparticular links between the convolutional feature detectors may differ in various embodiments, which are not all depicted in the figures. A person of skill in the art would appreciate that such embodiments are contemplated herein.

[0161] As shown in the system depicted in FIG. 13, the input to the CNN 1300comprises a biological sequence encoded by an encoder as a vector sequence or matrix.

[0162] One method that may be applied by the encoder is to encode the sequence ofsymbols in a sequence of numerical vectors, a vector sequence, using, for example, one-hot encoding.

[0163] For a RNA sequence of length L, the method may transform it into a one-hot-encoded sequence (A = [1,0,0,0], C = [0,1,0,0], G = [0,0,1,0], T = [0,0,0,1]) to obtain an inputmatrix of size L × 4. N = [0,0,0,0] may be used for the part that extends outside the transcript.This matrix may be used as the input to the PTGR model. The output of the model may be amatrix of size L × K, where K denotes the total number of cell lines. The k-th row of the matrixmay represent the predicted miRNA binding probability or predicted degradome read coverage of each base pair for the k-th cell line.

[0164] The symbol si is encoded in a numerical vector xi of length m: xi = (xi,1,… , xi,m)where xi,j = [si j] and [ ] is defined such that [True] = 1 and [False] = 0 (so called Iverson’snotation). One-hot encoding of all of the biological sequence elements produces an m × nmatrix X . For example, a DNA sequence CAAGTTT of length n = 7 and with an alphabet= (A, C, G, T), such that m = 4, would produce the following vector sequence:for representing biological sequences such asmRNA as numeric inputs to the neural network. It will be appreciated that other encodings ofmay be computed from linear or non-linear transformations of a one-hot encoding, so long as the transformed values are still distinct.

[0167] The CNN examples described above may all be implemented by the same orpossibly different CNN structures; that is, the number, composition and parameters of the filters, layers may or may not differ.

[0168] In one embodiment, the CNN may be trained by operating the CNN in amodified back-propagation mode using a dataset of examples, wherein each examplecomprises a biological sequence (mRNA) and corresponding miRNA binding data, andoptionally corresponding mRNA degradome data. For each example, the CNN is operated in the forward-propagation mode to ascertain the outputs of the CNN. Then, the CNN is operated in a modified back-propagation mode to determine the gradients for the parameters. These gradients are collected over examples, such as batches or minibatches, and are used to update the parameters. It will be appreciated that for all of the embodiments described above, the filter output can be differentiated with respect to the parameters of the filters.

[0169] Referring next to FIG. 14 there is shown a device drawing of the server 1006aof FIG.10 in accordance with one or more embodiments.

[0170] The device 1400 provides functionality described herein to determine a miRNAbinding prediction or a mRNA degradation prediction, and optionally for generating a model1426 for predicting miRNA binding or degradation, in accordance with one or moreembodiments. The server 1400 includes a communication unit 1404, a display unit 1406, a processor unit 1408, a memory unit 1410, an I / O unit 1412, a user interface engine 1414, and a power unit 1416.

[0171] The communication unit 1404 can include wired or wireless connectioncapabilities. The communication unit 1404 can include a radio that communicates using standards such as IEEE 802.11a, 802.11b, 802.11g, or 802.11n. The communication unit 1404 can be used by the server 1400 to communicate with other devices or computers.

[0172] Communication unit 1404 may communicate with a network, such as network2004 (see FIG.20).

[0173] The display 1406 may be an LED or LCD based display, and may be a touch sensitive user input device that supports gestures.

[0174] The processor unit 1408 controls the operation of the server 1400. The processor unit 1408 can be any suitable processor, controller or digital signal processor that can provide sufficient processing power depending on the configuration, purposes and requirements of the server 1400 as is known by those skilled in the art. For example, the processor unit 1408 may be a high-performance general processor. In alternative embodiments, the processor unit 1408 can include more than one processor with each processor being configured to perform different dedicated tasks. The processor unit 1408 may include a standard processor, such as an Intel® processor or an AMD® processor. In one embodiment, processor unit 1408 may comprise one or more Graphic Processing Units (GPUs), such as but not limited to a Compute Unified Device Architecture (CUDA)- compatible GPU with 4 GB of VRAM or higher.

[0175] The processor unit 1408 can also execute a user interface (UI) engine 1414 that is used to generate various UIs for delivery via a web application.

[0176] The memory unit 1410 comprises software code for implementing an operating system 1420, programs 1422, database 1424, PTGR model 1426, training unit 1428, and prediction unit 1430.

[0177] The memory unit 1410 can include RAM, ROM, one or more hard drives, one or more flash drives or some other suitable data storage elements such as disk drives, etc. The memory unit 1410 is used to store an operating system 1420 and programs 1422 as is commonly known by those skilled in the art.

[0178] The I / O unit 1412 can include at least one of a mouse, a keyboard, a touch screen, a thumbwheel, a track-pad, a track-ball, a card-reader, an audio source, a microphone, voice recognition software and the like again depending on the particular implementation of the server 1400. In some cases, some of these components can be integrated with one another.

[0179] The user interface engine 1414 is configured to generate interfaces for users to request model training, request miRNA binding predictions based on a target nucleic acidsequence, review miRNA interaction data, or other activities associated with training and miRNA binding prediction. The various interfaces generated by the user interface engine 1414 may be transmitted to a user device via the communication unit 1404.

[0180] The user interface engine 1414 may provide one or more visualizations of the output from the model (see e.g. PTGR model 1300 in FIG. 13). This may include a visualization of the input mRNA sequence, the output miRNA binding prediction matrix, and the output mRNA degradome prediction matrix. The visualization may be user selectable, for example, by a user clicking on or other selecting a nucleotide within the mRNA input sequence.

[0181] In one example, the visualization generated by the user interface engine 1414 may be the visualization in FIG.2D, which is a track visualization of the PTGR model’s output showing iNeuron miRNA and degradome predictions showing that the model predicts novel miRNA binding sites and predicts whether the site is derogatory.

[0182] In another example, the visualization generated by the user interface engine 1414 may be the visualization in FIG.2F, which is a track visualization of the PTGR model’s output showing identification of AREs.

[0183] In another example, the visualization generated by the user interface engine 1414 may be the visualization in FIG.3G, which is a sequence attribution that can be used to identify the identity of the predicted miRNA binding; the PTGR model miRNA scores may be a reflection of the miRNA expression of the corresponding miRNAs in those cell lines.

[0184] In another example, the visualization generated by the user interface engine 1414 may be the visualization in FIG.4A, which is a track visualization of the PTGR model identifying miRNA binding disrupting and creating variants.

[0185] In another example, the visualization generated by the user interface engine 1414 may be the visualization in FIG.4G which is an illustration of small sequence changes in the slutskin MPRA results in large changes in expression and predicted the PTGR model scores, shows and shows the model’s sensitivity to small sequence changes.

[0186] In another example, the visualization generated by the user interface engine 1414 may be the visualization in FIG.5C and FIG.5D which show examples of the PTGR model efficiently identifying miRNA blocking SBOs in UTRN and CHD9.

[0187] The power unit 1416 can be any suitable power source that provides power to the server 1400 such as a power adaptor or a rechargeable battery pack depending on the implementation of the server 1400 as is known by those skilled in the art.

[0188] The operating system 1420 may provide various basic operational processes for the server 1400. For example, the operating system 1420 may be a server operating system such as Ubuntu® Linux, Microsoft® Windows Server® operating system, or another operating system.

[0189] The programs 1422 include various user programs. They may include several hosted applications delivering services to users over the network, for example, a web application and an API application, and other applications as known.

[0190] In one or more embodiments, the programs 1422 may provide an application such as a web-based application for miRNA binding prediction, or client-server based application via an API. The application may provide functionality for a user to submit a miRNA prediction request including a target nucleic sequence that may initiate the miRNA binding prediction methods described in FIG.25. The application may provide functionality for a user to submit a model generation request for a PTGR model 1426 that may initiate the PTGR model training methods described in FIGs.17. The application may further provide access for queries and data analysis of the miRNA interaction data in database 1424, PTGR models in database 1424, miRNA library data, or other data.

[0191] The database 1424 may be a database for storing PTGR model data and miRNA interaction data, such as but not limited to miRNA binding data and degradome data created using the models described herein.

[0192] Database 1424 may include binding data from a plurality of different miRNA datasets comprising binding data associating with various miRNAs, degradomes, or mRNA sequences.

[0193] The PTGR model 1426 may include a degradome model and an miRNA binding prediction model. The PTGR model 1426 may provide at least one of degradome predictions and miRNA binding predictions. The PTGR model 1426 may be trained and / or evaluated on empirical data, such as data derived from methods set out in Example 1 or otherwise known in the art describing the binding of miRNA to mRNA and / or degradome data.

[0194] In one embodiment, PTGR model 1426 is a machine learning model such as a neural network. The PTGR model 1426 may be the model in FIG.13, including an input layer, a base, one or more heads, and one or more residual blocks.

[0195] The training module 1428 is for generating the machine learning model (such as a neural network) for operation by the PTGR model 1426, and for storage in database 1424. The training module 1428 may perform methods as described herein in FIGs.17. Once the model is generated by the training module 1428, it may be stored in the database 1424 or used by the PTGR model 1426.

[0196] The prediction module 1430 may determine miRNA binding predictions using the PTGR model 1426. The miRNA binding predictions may be generated by the prediction module 1430 and stored in database 1424, or transmitted to a user at a user computing device. A user may then use the miRNA prediction to design or synthesize a therapeutic application of the miRNA, optionally by reprogramming existing transcription factors or other biomolecules to target a particular nucleic acid sequence.

[0197] Referring next to FIG. 15 there is shown a method drawing 1500 of miRNA binding prediction in accordance with one or more embodiments. The method is a computer- implemented method for predicting tissue-specific micro-RNA (miRNA) binding to a messenger-RNA (mRNA).

[0198] At 1502, a PTGR model is provided at a memory, the PTGR model comprising a base and at least one head.

[0199] At 1504, an RNA input sequence is received at a processor in communication with the memory, the RNA input sequence of length L corresponding to the mRNA.

[0200] At 1506, an miRNA binding prediction matrix is determined at the processor, the miRNA binding prediction matrix determined from the at least one head of the PTGR model, the miRNA binding prediction matrix comprising an L x K matrix comprising miRNA binding predictions at each position of the L nucleotides in the RNA input sequence for K tissue types, the miRNA binding prediction matrix determined using the RNA input sequence as input at the base of the PTGR model.

[0201] In one or more embodiments, the method may further comprises predicting tissue-specific degradation of the mRNA and the at least one head comprises at least two heads: a first head providing a first miRNA binding prediction matrix binding for K tissue types, optionally K human tissue types; and a second head providing a first mRNA degradome prediction matrix for K tissue types, optionally K human tissue types, the first mRNA degradome prediction matrix determined using the RNA input sequence as input at the base of the PTGR model.

[0202] In one or more embodiments, the first mRNA degradome prediction matrix may comprise an L x K matrix comprising degradome read coverage predictions at each position of the L nucleotides in the RNA input sequence for K tissue types.

[0203] In one or more embodiments, the at least two heads may comprise at least four heads: a third head providing a second miRNA binding prediction matrix for tissue types for a second organism, optionally mouse; and a fourth head providing a second mRNA degradome prediction matrix for tissue types for the second organism, optionally mouse.

[0204] In one or more embodiments, the PTGR model may comprise a convolutional neural network.

[0205] In one or more embodiments, the base may comprise: at least one base convolution layer; and at least one base residual block receiving an output from the at least one base convolution layer.

[0206] In one or more embodiments, each head of the at least one head may comprise a first head convolutional layer receiving input from the base, at least one head residual block receiving an output of the first head convolutional layer, and a second head convolutional layer receiving an output of the at least one head residual block.

[0207] In one or more embodiments, each residual block may comprise: a normalization layer receiving an input to the residual block; a first pair of convolutional layers receiving the output from the normalization layer, a first of the first pair having a tanh activation and a second of the first pair having a sigmoid activation; a first element-wise multiplication operation on outputs of the first pair of convolutional layers; a second pair of convolutional layers receiving the output from the first element-wise multiplication operation, a first of the second pair having a tanh activation and a second of the second pair having a sigmoid activation; a second element-wise multiplication operation on outputs of the second pair of convolutional layers; and an output of the residual block determined by adding an output of the second element-wise multiplication operation and the input to the residual block in a skip connection.

[0208] In one or more embodiments, the method may further comprise: outputting, at a display device in communication with the processor, at least one of the miRNA binding prediction matrix and the mRNA degradome prediction matrix.

[0209] In one or more embodiments, the method may further comprise: receiving, at an input device in communication with the processor, one or more user selected nucleotides in the RNA input sequence; generating a first visualization comprising the miRNA binding predictions at the one or more user selected nucleotides in the RNA input sequence; and outputting, at the display device in communication with the processor, the generated first visualization.

[0210] In one or more embodiments, the method may further comprise: generating a second visualization comprising the mRNA degradome prediction matrix at one or more user selected nucleotides in the RNA input sequence; and outputting, at the display device, the generated second visualization.

[0211] In one or more embodiments, the method may further comprise: converting, at the processor, the RNA input sequence of length L to an L x 4 input matrix, wherein the L x 4 input matrix is one-hot encoded.

[0212] In one or more embodiments, the miRNA binding predictions may be indicative of non-specific miRNA binding to the mRNA.

[0213] In one or more embodiments, the K tissue types may comprise at least one of liver, CNS, heart, kidney, muscle, cancer cells, HeLa cells, iPSCs, primary human hepatocytes (PHHs) peripheral blood mononuclear cells (PBMCs) and other tissue types or cell lines known in the art.

[0214] Referring next to FIG. 16 there is shown a method diagram 1600 for determining an effect of a variant mRNA sequence relative to a control mRNA sequence on post transcriptional gene regulation. Understanding the intricacies of PTGR with cellular environments may necessitate also understanding the profound impact genetic variants can have on this process. In diseases ranging from cancer to genetic disorders, variants act as pivotal turning points, either by initiating a pathogenic process or by exacerbating an existing condition. While coding variants, which occur within the coding region of genes, can directly alter amino acid sequences and potentially disrupt protein function, non-coding variants pose a subtler, yet equally significant challenge. These variants, located in introns, UTRs, or intergenic spaces, can influence gene expression and protein production without altering the protein sequence, affecting the delicate balance of cellular processes crucial for maintaining health.

[0215] At 1602, tissue specific miRNA binding to the variant mRNA sequence is predicted, and optionally tissue specific degradation of the variant mRNA sequence, according to the method of described in FIG.15.

[0216] At 1604, comparing the miRNA binding prediction matrix, and optionally the mRNA degradome prediction matrix, for the variant mRNA sequence to a miRNA binding prediction matrix, and optionally a mRNA degradome prediction matrix, for the control mRNA sequence.

[0217] At 1606, the effect of the variant mRNA sequence on post transcriptional gene regulation is determined based on any differences between the miRNA binding prediction matrix, and optionally the mRNA degradome prediction matrix, for the variant mRNA sequence, relative to the miRNA binding prediction matrix, and optionally the mRNA degradome prediction matrix, for the control mRNA sequence.

[0218] In one or more embodiments, the method may further comprise determining the miRNA binding prediction matrix, and optionally the mRNA degradome prediction matrix, for the control mRNA sequence according to the method of FIG.15.

[0219] In one or more embodiments, the variant mRNA sequence may have between 1 and 20 single nucleotide polymorphisms relative to the control mRNA sequence.

[0220] In one or more embodiments, the variant mRNA sequence may have a variant of unknown clinical significance, a putative disease-causing mutation, a masked sequence corresponding to a SBO binding site, an ADAR editing site or a synthetic mRNA sequence, or the effect of ADAR editing.

[0221] In one or more embodiments, the method may further comprise synthesizing the variant mRNA and testing the variant mRNA for gene expression.

[0222] Referring next to FIG.17 there is shown a method drawing 1700 for generating a PTGR model in accordance with one or more embodiments.

[0223] The PTGR model may be trained using annotations from the National Cener for Biotechnology Information (NCBI), such as refseq.v109 for the human cell lines and ncbi refseq.m38.v106 for the mouse cell lines. The model training may use a first dataset including miRNA binding data for various mRNA segments of nucleotides (or also referred to herein as transcripts). The training method 1700 may use only protein-coding transcripts to train the model. For each gene, the most principal transcript (such as according to APPRIS) may be used. Each transcript from both the human and mouse transcriptomes may be divided into nonoverlapping windows of 3000 nucleotides across its mature mRNA sequence. A batch of 10 windows may be randomly sampled during training, and 6.25 kb of transcript context sequence around the sampled window may be appended to each side, and may result in a total input sequence length of 15.5 kb used to train the model. The PTGR model may be trained on transcript data including all exonic CLIP peaks, and no filtering based on miRNA identity, conservation, or site type was performed.

[0224] For a given sequence, the true labels may correspond to the miRNA outputs, including a binary matrix of shape (3000, 29) corresponding to the length of the sequence and the total number of human and mouse cell lines / tissues with miRNA binding data andpredictions. The binary matrix may have a true label of 1 for every location corresponding to a CLIP identified peak and 0 everywhere else. Similarly, the true labels for the degradome outputs may be a real valued matrix of shape (3000, 10) corresponding to the degradome- seq read coverages for the given input sequence and the total number of human and mouse cell lines / tissues with Degradome-Seq data / predictions. The degradome-seq labels may be processed with squashed transformation (x0.375) to reduce the skewness.

[0225] The PTGR model’s miRNA binding outputs may be in the range (0, 1) while its degradation outputs may be in the range (0, infinity). Binary cross-entropy loss may be used for the miRNA binding output and Poisson loss may be used for the degradationoutput. The total loss is the weighted sum of the two may be:= + (Equation 1)

[0226] When a sequence from a human transcript is sampled, gradients from the corresponding human miRNA and degradome output heads may be used to update the network, and vice-versa when a sequence from a mouse transcript is sampled. The final model may be an ensemble of models trained on the four splits with different sets of chromosomes for training and validation. The four splits may use chromosomes (1, 2), (6, 7, 8, 9), (10, 11, 12, 13), and (14, 15, 16, 17, 18) for validation (4-fold cross validation) respectively and the remaining for training. The splits may be chosen as described to ensure that 75% of identified miRNA CLIP peaks are used for training and 25% for validation.

[0227] For model evaluation, predictions and labels of the validation intervals may first be processed. For each cell line, a pre-calculated window length based on median peak width for the CLIP data may be used, and a fixed length of 50 base pairs for the degradation data may be used.

[0228] First, an average pooling for both predictions may be applied and labels using the window length as the pool size may be used. Then, the predictions and labels may be divided into multiple chunks with the defined window length. For the miRNA predictions, the value 1 may be assigned assigned otherwise. Then, for each cell line, area under the curve (AUC) and area under the precision-recall curve (AUPRC) may be calculated across all intervals for miRNA binding predictions, andSpearman may be calculated for degradation predictions. The mean value of each metric across all cell lines as the final metrics may be reported for evaluation purposes.

[0229] The training method code may be implemented in TensorFlow 2.11.0. The model weights may be initialized using the Glorot Uniform distribution. The Adam optimizer may be used for model training1 2= 0.999, . Trained may be performed with a batch size of 5 on a single NVIDIA RTX A6000 GPU. The learning rate may be tuned and the loss weights based on the evaluation metrics on the validation set, and also conducted early stopping in 60 training epochs based on these metrics.

[0230] The optimal learning rate may be 0.00005 for the 133M parameter model, with 1 2 = 0.1. ReduceLROnPlateau may be used for learning rate scheduling, monitoring on validation loss with factor=0.5, patience=5, and cooldown=1.

[0231] At 1702, providing in a memory, a machine learning model comprising a base and at least two heads.

[0232] At 1704, providing, in the memory, a first data set comprising miRNA binding data for a plurality of mRNA sequences for a K plurality of tissues;

[0233] At 1706, training the machine learning model based on the first data set, wherein a prediction output of the machine learning model is a matrix of size L × K, K denoting a size of the K plurality of tissues, a k-th row of the matrix corresponding to a predicted miRNA binding probability each base pair for the k-th cell line.

[0234] In one or more embodiments, the miRNA binding data comprises mRNA sequence segments experimentally associated with miRNA binding.

[0235] In one or more embodiments, the method may further comprise: providing, in the memory, a second data set comprising mRNA degradome data for the plurality of mRNA sequences for the K plurality of tissues, wherein the training the machine learning model may further comprise training the machine learning model based on the first data set and the second data set.

[0236] In one or more embodiments, each of the mRNA sequences may comprise a degradome annotation, and the degradome annotation comprises a degradome matrix corresponding to degradome-seq data for each of the K plurality of tissue types.

[0237] In one or more embodiments, the mRNA degradome data may further comprise a read coverage value at each position of the mRNA sequence.

[0238] In one or more embodiments, the base and the at least two heads may comprise at least one residual block.

[0239] In one or more embodiments, the at least two heads may comprise a first head for providing a miRNA binding prediction matrix and a second head for providing a mRNA degradome prediction matrix and training the machine learning model may comprise training the base and the first head on the first data set and training the base and the second head on the second data set.

[0240] In one or more embodiments, the plurality of mRNA sequences may comprise a plurality of non-overlapping windows.

[0241] In one or more embodiments, the plurality of non-overlapping windows may be at least 1500, 2000 or 3000 nucleotides long.

[0242] In one or more embodiments, each of the plurality of non-overlapping windows may comprise an mRNA context sequence of about 500 bp to 10 kb, optionally about 6.25 kb around each non-overlapping window.

[0243] In one or more embodiments, a total input sequence length of each non- overlapping window may be about 500 bp to 25 kb, optionally about 15.5 kb.

[0244] In one or more embodiments, the method may further comprise: converting each of the plurality of mRNA sequences to a one-hot encoded sequence.

[0245] In one or more embodiments, a loss function for the training the machinelearning model may comprise = + , wherein a binary cross-entropy losscomprises a miRNA binding output and a Poisson loss comprises a degradation output. Alternatively, the loss function for training the ML model may comprise other loss functions known in the art such as a mean square loss function.

[0246] In one or more embodiments, the machine learning model may comprise an ensemble model.

[0247] In one or more embodiments, the ensemble model may comprise at least four models, each of the four models trained based on a training data split across different sets of chromosomes.

[0248] In one or more embodiments, the training the machine learning model may comprise training the machine learning model using an optimizer, such as stochastic gradient descent (SGD), RMSprop, or an Adam optimizer.

[0249] As set out in the following Examples, Developing REPRESS necessitated the generation of a large scale dataset of miRNA binding and mRNA degradation across 29 cell types and tissues in both human and mouse, and the creation of a novel ConvNeXt-inspired neural architecture. REPRESS’s neural architecture has a capacity that far exceeds existing methods (133 million parameters vs tens of thousands), and exhibits a strong scaling law of performance gains with model size. Interestingly, the REPRESS-specific architecture also outperformed alternative implementations utilizing alternative architectures including transformers.

[0250] REPRESS can overcome deficiencies in the underlying training data. For example, we demonstrated that REPRESS can correctly identify validated miRNA binding sites not identified in the specific miR-eCLIP datasets used for training. This increased sensitivity is particularly important when modeling repressive phenomena that may push targets of strong interactions below an assay’s detection threshold. The ability of REPRESS to overcome the limitations of assays and ‘denoise the training data’ is advantageous. By jointly modeling miRNA binding and mRNA degradation, REPRESS achieves synergies during training and when making predictions. We found that training on the joint dataset yielded better test performance on the individual miRNA binding and mRNA degradation tasks, compared to training on task specific data. Notably, we found that this performance increase became larger when we increased the model size. Although REPRESS cannot generalize to cell lines not present in the training dataset, we propose a prototype of a ”miRNA-specific” variation of REPRESS and provide preliminary evidence showing that itcan predict miRNA targets in unseen cell lines and tissue types based solely on miRNA expression, at least in part alleviating the need to generate additional training data.

[0251] In support of the development of RNA based therapeutics across multiple modalities, REPRESS enables in silico experimentation at scales that are often not experimentally tractable. For instance, REPRESS can nominate genetic variants or therapeutic molecules (e.g. ASOs) that dampen repression and increase overall gene expression for therapeutic targets exhibiting a loss-of-function phenotype. These prioritized variants and / or therapeutic molecules can then be subsequently screened using experimental assays, such as gene editing or screening a prioritized set of ASOs in disease relevant models. REPRESS also enables the design of novel UTR sequences optimized to maximize tissue-specific expression for both mRNA and DNA encoded payloads including gene replacement therapies. EXAMPLES Example 1: Transcriptome-wide regulatory maps of miRNA binding and mRNA degradation

[0252] To establish a transcriptome-wide map of miRNA targeting, first 18 existing AGO2-CLIP binding site datasets from the literature were curated. While abundant, AGO2- CLIP data does not directly identify the specific miRNA(s) binding at a particular target. However, publicly available datasets from miRNA specific binding assays such as miR-eCLIP

[0034] were not very abundant, especially for cell lines and tissues of therapeutic interest. Thus, three publicly available datasets were curated, and access was purchased to an additional two datasets, and generated novel miR-eCLIP data from four human cell lines and two mouse tissues, for a total of 11 miRNA-specific CLIP datasets. To establish a transcriptome-wide map of site-specific mRNA degradation, Degradome-Seq was conducted, which identified uncapped 5’ transcript fragment ends resulting from endonucleolytic cleavage, in six human cell lines and four mouse tissues, as existing Degradome-Seq datasets have primarily been generated using non-mammalian tissues [36, 37]. Collectively, these efforts resulted in the high quality compendium of datasets needed to train REPRESS, a deep learning model thatseparately predicts cell type specific miRNA targeting and mRNA degradation for any input sequence to enable biological insights (FIG.1A and FIG.1B).

[0253] Peak calling on the miR-eCLIP datasets revealed precise miRNA-specific binding at more than one million target sites across 24,000 protein-coding genes across the human and mouse transcriptomes. A high proportion (average of 44% across cell lines and tissues) of these miRNA peaks contained a seed match of the targeting miRNA identified (FIG.1D) similar to, or exceeding, rates obtained by cross-linking ligation and sequencing of hybrids (CLASH) and miR-eCLIP [32, 34]. Furthermore, RNA co-fold analysis using LinearCoPartition

[0038] revealed clusters of miRNAs with distinct miRNA-target base pairing frequencies within and beyond the miRNA seed region, highlighting the molecular detail captured by miR-eCLIP (FIG. 1F). Together, these results demonstrate the ability of miR- eCLIP to identify precise miRNA-specific target loci. Overall, a diverse set of miRNA targets were detected across cell lines and tissues; miRNA-seq revealed highly tissue-specific expression profiles across miRNAs, and a corresponding strong positive correlation between miRNA expression levels and the number of miR-eCLIP peaks identified for it was observed While most miRNAs demonstrated preferences for targeting near the beginning and, to a lesser extent, the end of 3’ UTRs, several cell types had small clusters of miRNAs with distinct targeting patterns, including increased targeting toward the middle of 3’ UTRs.

[0254] While CLIP-based approaches offer powerful detection of targets of miRNAs (and RBPs), they do not in isolation foretell the functional consequences of any particular binding event. To assess one functional aspect of PTGR, Degradome-Seq was leveraged to sequence uncapped 5’ transcript ends and identify site-specific RNA degradation events resulting from either miRNA or RBP binding, in a mechanism-agnostic fashion. This resulted in over 500,000 unique Degradome-Seq peaks. Degradation levels at individual loci demonstrated some cell type specificity (average inter-sample Pearson R = 0.53; P = 4.8 × was performed on the relative Degradome-Seq signal at positions along protein-coding mRNAs, and consistently identified distinct clusters of transcripts, including those with: i) shorter 3’ UTRs with few miRNA targets, longer half-lives,and degradation predominantly in the CDS; and ii) longer 3’ UTRs with many miRNA targets, shorter half-lives, and highly localized 3’ UTR degradation (FIG.1C).

[0255] To further explore potential associations between different classes of RNA regulatory elements and sites of degradation captured by Degradome-Seq, ENCODE eCLIP data for RBP binding

[0039] and AU-Rich Element (ARE) motifs was leveraged, in addition to our miR-eCLIP peaks. Notably, overall elevated degradation at each of these classes of regulatory elements was detected, highlighting the utility of Degradome-Seq data to uncover diverse functional post-transcriptional regulatory elements (FIG.1E). Among all RBPs with available ENCODE eCLIP binding data, elevated Degradome-Seq read coverage was detected most prominently for FUBP3, PUM2, TIA1, and UPF1. This finding is consistent with UPF1 recently being identified as having the most enriched binding (among all RBPs analyzed) near capped 5’ ends in human 3’ UTRs

[0040] , while also suggesting that other RBPs such as FUBP3, PUM2, and TIA1 may have underappreciated roles in regulating mRNA decay through endonucleolytic cleavage. Consistent with observations of occasional miRNA- directed target cleavage [35, 41], a more modest increase in Degradome-Seq read coverage was also detected around miRNA targets overall, and site-specific degradation for individually validated instances of miRNA directed endonucleolytic mRNA cleavage

[0035] . Lastly, robust elevation in Degradome-Seq read coverage was detected around AREs, suggesting these elements may direct degradation by endonucleases. Intriguingly, these findings also offer insight into potential mechanistic differences between mRNA cleavage events associated with miRNA or RBP binding, as inferred cleavage sites occurred most frequently at miRNA targets, and downstream of AREs and RBP binding sites. In summary, our results suggest site- and tissue-specific mRNA cleavage is an integral part of PTGR mediated by diverse molecular interactions. Together, these miRNA binding and mRNA degradation datasets provide valuable biological insights and facilitate the training of a robust deep learning model. Example 2: A single-base resolution, cell type specific model of miRNA binding and mRNA degradation

[0256] Building on the foundation described in Example 1, a dilated Convolutional Neural Network (CNN) model of post-transcriptional gene regulation termed “REPRESS”was developed to predict species-specific miRNA binding across 29 cell types and RNA degradation across 10 cell types at single base pair resolution.

[0257] FIG. 13 shows the model architecture for REPRESS including a base for receiving an input mRNA sequence as well as separate heads that output tissue-specific miRNA binding and degradation predictions at each position along the mRNA sequence.

[0258] Training involved randomly sampled regions from the human and mouse transcriptome spanning 15.5kb, with the model analyzing information within a 12.5kb context for each base pair prediction. This expansive contextual understanding allows REPRESS to adeptly capture long-range dependencies and incorporate information across the various regions of a transcript.

[0259] REPRESS makes use of a novel 133M parameter neural network that takes an input nucleotide sequence and separately predicts miRNA binding and mRNA degradation (degradome-seq read coverage) across 29 and 10 cell lines / tissue types, respectively, at single base-pair resolution (FIG.2A, FIG.13). REPRESS uses a novel neural architecture inspired by ConvNeXt

[0042] , including residual blocks with gradually increasing receptive fields that capture long range dependencies by incorporating 12.5 kb of context around the query sequence, which is sufficient to model information across the entire mature mRNA length for 99.2% of all protein-coding mRNAs in the human and mouse genomes.

[0260] To evaluate REPRESS’s performance on miRNA target detection, we applied it to a validation set (4-fold cross validation) containing 25% of transcripts from the human and mouse transcriptome. REPRESS achieves an average 4-fold cross validation Area Under Receiver Operating Characteristic (AUROC) curve of 0.88 (FIG.6A) and an average Spearman correlation of 0.60 when predicting degradation events (FIG.6B). To further test REPRESS’s performance on orthogonally generated sets of experimentally validated miRNA targets we leveraged data from miRTarBase

[0043] and DianaTarBase

[0044] ( see Examples: Methods). REPRESS accurately predicted cell-type specific miRNA targets with an average AUROC curve of 0.78 including those that were not present in the training dataset (FIG.5C).

[0261] REPRESS’s ConvNeXt-inspired dilated CNN architecture exhibits strong scaling laws and outperforms other state of the art architectures that we trained on the samedata. REPRESS’s accuracy, measured by performance on the validation set of both miRNA binding and degradome tasks, scales smoothly and strongly with the number of parameters in the model (FIG. 6F). This is consistent with the scaling laws observed in deep learning models used for natural language processing

[0045] and biological sequence modeling

[0046] . When the model capacity is kept fixed at a large size, we also found that models jointly trained on both miRNA and Degradome-Seq data performed better than models trained on either individual dataset. We evaluated the generalization performance of REPRESS’s ConvNeXt- inspired architecture against other widely used deep learning architectures that we trained on the same data, including a vanilla CNN, Transformers

[0047] , Mamba (a state-space model)

[0048] , and a vanilla dilated CNN (without the ConvNeXt inspired residual block), all utilizing the same 12.5 kb context sequence. REPRESS’s ConvNeXt-inspired architecture demonstrated superior generalization, outperforming the Mamba model, the next highest performing architecture, by 10%, and the vanilla dilated CNN by 17% (FIG.6G).

[0262] To evaluate and benchmark REPRESS’s predictive capabilities against existing models of miRNA targeting, we compiled a comprehensive suite of diverse tasks related to the classification and functional impact of miRNA binding sites. We evaluated performance on all tasks for the following models: TargetScan

[0013] , Biochemical model

[0014] , miRAW

[0049] , DeepMirTar

[0017] , miTAR

[0018] , miRBind

[0019] and DMISO

[0020] . Tasks included: 1. Identifying experimentally validated miRNA binding sites from i) REPRESS’s test set, ii) the miRAW test set

[0049] , iii) the DeepMirTar test set

[0017] , iv) miRTarBase, and v) HEAP binding data

[0050] ; 2. Predicting miRNA-mediated repression of synthetic sequences from two Massively Parallel Reporter Assays (MPRAs) designed to interrogate miRNA targeting [51, 52]; and 3. Predicting the impact of miRNA binding-altering variants curated from the literature. REPRESS substantially outperformed the other models on all tasks, with the exception of tasks where the evaluation data was from the same distribution used to train the corresponding model, as indicated by asterisks (*) (FIG. 2B, FIGs. 7A-7E, FIG. 8). Specifically, while miRAW performed best on its own test set, it performed poorly on the DeepMirTar test set, and vice versa, suggesting overfitting to their respective datasets. For both the miRAW and DeepMirTar tasks, REPRESS was still the best performingindependently trained model, demonstrating its ability to generalize across diverse miRNA datasets.

[0263] For any arbitrary sequence of interest, REPRESS separately predicts the likelihood of miRNA binding and Degradome-seq read coverage. In addition to these predictions, REPRESS can be used to query all possible single nucleotide variants along a sequence to predict precisely which nucleotides are important for miRNA binding or mRNA degradation. As an example of miRNA binding prediction, REPRESS accurately predicted a hsa-miR-30c-5p binding site in the 3’ UTR of human NOTCH1 present in the training data, but also two additional upstream binding sites (for hsa-miR-34a-5p and hsa-miR-144-3p) not present in the training data but independently experimentally validated

[0053] , with REPRESS attributing this binding to seed matches for the corresponding miRNAs (FIG. 2C). As an example of degradome read coverage prediction, the strongest REPRESS degradome peak in the ENO2 3’ UTR overlapped an AU-rich element (ARE) present in the ARED-Plus database

[0054] , with REPRESS accurately attributing this signal to an (ATTT) repeat characteristic of AREs (FIG.2D). More generally, REPRESS can be deployed as a powerful tool for biological discovery, decoding the impact of genetic variants, and designing genetic medicines. Highlights of these capabilities include demonstrated learning of diverse canonical and non-canonical [7] modes of binding used by endogenous miRNAs and accurate quantitative prediction of subsequent repression for a subset of sequences screened from the McGeary et al. MPRA

[0052] (FIG. 2E), strong correlations between predicted miRNA binding and observed repression for mutated derivatives of miRNA target sites, even at the level of protein (FIG. 2F), and extremely efficient design of validated ASOs to sterically block miRNA binding sites and increase expression of therapeutic targets (FIG.2G). Example 3: REPRESS enables discovery of post-transcriptional regulatory biology

[0264] The PTGR model developed in Example 2 was further investigated to determine whether it could uncover the fundamental properties driving miRNA binding and degradation. We next wanted to explore how REPRESS could be leveraged to discover both validated and novel aspects of endogenous biology. We first utilized REPRESS to generate sequence attributions for strongly predicted target sites of highly expressed miRNAs toassess to what extent each target site position / base contributed to miRNA binding. These sequence attributions were computed by running REPRESS on in silico mutated sequences (called in silico mutagenesis, ISM) and were represented by positional scores along the mature miRNA sequence for interpretability. As expected, REPRESS generally attributed miRNA binding as being most prominently driven by the seed region spanning positions 2-7 (FIG. 3A). Strikingly, REPRESS identified that hsa-miR-148a-3p had a shifted attribution pattern that spanned positions 3-8, instead of 2-7, corresponding to an offset-6mer site type. While it is well established that offset-6mers are generally less potent than their other canonical site type counterparts [14, 55], REPRESS offered a hypothesis that offset-6mer targets are of increased functional relevance for the binding of hsa-miR-148a-3p specifically. To test this hypothesis and explore the validity of REPRESS’s prediction, we analyzed both the underlying miR-eCLIP data and orthogonal cross-species conservation data. We found that relative to typical miRNAs (e.g. hsa-let-7a-5p), hsa-miR-148a-3p engaged offset 6mer target sites more frequently, and these targets were conserved across vertebrates at a significantly higher rate than 6mer targets (at a rate similar to that observed for 8mer targets) (FIG.3A).

[0265] REPRESS also predicted an influence of positions outside of the seed region for some miRNAs. This prompted us to generate REPRESS sequence attributions for different site types for each miRNA. When comparing generally weaker canonical 6mer sites to generally stronger canonical 8mer sites, REPRESS attributed increased importance to non-seed region positions for 6mer targets for many miRNAs. One possible interpretation of this finding is that the weaker seed region base pairing offered by 6mers is often augmented by additional extended miRNA-target interactions, whereas stronger 8mers are less dependent on such interactions for achieving higher affinity binding

[0056] . To obtain independent evidence for this hypothesis, we compared rates of conservation across vertebrates for non-seed regions of both 6mer and 8mer target sites. Nearly half (21 of 47 = 45%) of miRNAs eligible for analysis had significantly higher rates of conservation at non- seed positions in 6mer targets relative to 8mer targets (mean odds ratio = 1.64), despite the vast majority of miRNAs having lower rates of conservation within the seed region for 6mer targets relative to 8mer targets (FIG.3B). Strikingly, for all of the highly expressed miRNAswith substantially higher non-seed region conservation for 6mer targets, REPRESS also attributed increased importance to non-seed positions for 6mer targets. In contrast, for the miRNA with the highest non-seed region conservation for 8mer targets (hsa-miR-103a-3p), REPRESS correspondingly attributed increased importance of non-seed positions for 8mer targets (FIG.3B). In addition to these novel insights into miRNA biology, we also confirmed that REPRESS has learned well established general rules of miRNA targeting: The ranking of site types by mean REPRESS score matches the relative strength of each site type, 8mer > 7mer-m8 > 7mer-A1 > 6mer

[0055] , and there is a strong positive linear relationship between the AU content of the dinucleotides flanking the 8mer seed region and mean REPRESS score (FIG. 3C). The latter result is consistent with independent observations of high AU contexts permitting stronger binding which has been attributed to secondary structure- dependent target accessibility [14, 57]. Notably, among all of the other models tested, only the biochemical model

[0014] (mostly) demonstrated the expected linear relationship between AU content and target strength (FIG.3C).

[0266] The biology of miRNA-mediated regulation is complex and nearby binding sites can act additively, synergistically, or competitively [51, 57, 58]. To determine whether REPRESS’s large receptive field enables it to account for sequence context and binding site multiplicity, we examined data for a series of engineered sequences containing a variable number of target sites for one of five different miRNAs in two different background sequence contexts

[0051] . As the number of miRNA targets increased from zero to five, REPRESS predictions tracked experimental repression for all five miRNAs in both sequence contexts, and was more accurate than all other models tested. For hsa-miR-20a-5p and hsa-miR-21- 5p, REPRESS correctly predicted that one background sequence had higher repression thanthe other (P = 4.8 × 10 and 9.7×10 , respectively) and also predicted that for hsa-miR-20a-5p, repression strongly increased with multiplicity (P = 4.2 × 10 ), whereas for hsa-miR-21-5p, it did not (FIG.3D). Notably, no other method that we tested correctly predicted the observed context-dependent repression or the difference in the effect of multiplicity on repression for the two different miRNAs, likely because other models are trained on shorter sequences (less than 60 nt) or model binding sites individually.

[0267] Given that REPRESS was trained using unperturbed wild type cells and tissues, we next examined whether it could predict miRNA mediated repression in the context of miRNA overexpression and knockout / knockdown experiments [59-62]. For each experiment, we calculated the average miRNA binding predicted by REPRESS or TargetScan, and selected the top 100 targets for each model. We then compared the distribution of log2 fold-changes of the corresponding transcript sets relative to all transcripts analyzed as a baseline. For the representative hsa-miR-122-5p transfection, REPRESS successfully identified transcripts that were significantly more repressed after transfectionrelative to all transcripts (KS statistic = 0.3, P = 1.7 × 10 ) and those highly ranked byTargetScan (KS statistic = 0.23, P = 0.009).

[0268] Similarly, for the representative mmu-miR-26a-5p knockout, REPRESS identified transcripts that were significantly less repressed after knockout relative to alltranscripts (KS statistic = 0.36, P = 3.8 × 10 ) and those highly ranked by TargetScan (KSstatistic = 0.24, P = 0.0061) (FIG. 3E). Notably, for experiments involving transfection of miRNAs with very sparse binding in our miR-eCLIP datasets, transcripts identified by REPRESS were not significantly more repressed after transfection compared to both baseline and TargetScan. We also found that REPRESS’s degradome predictions were inversely correlated with RNA stability measurements from an MPRA of ARE-containing 3’ UTRs

[0063] . Specifically, UTRs with higher ARED-plus cluster numbers

[0054] and lower stability in the MPRA had corresponding higher REPRESS degradation predictions (FIG.3F). One notable example of the predictive capabilities of REPRESS’s degradome model was identifying an experimentally validated stabilizing ARE in the mouse Bcl2 gene whose deletion caused reduction in the measured RNA stability

[0064] . Collectively, these results demonstrate that in learning fundamental aspects of miRNA targeting and mRNA degradation, REPRESS can aid in the discovery of new biology. Example 4: REPRESS decodes the impact of genetic variants

[0269] Accurate prediction of the impact of variant alleles on gene expression, stability, and disease phenotypes is a foundational goal in the application of deep learning to genomics. Such predictions allow for an improved understanding of the regulatory mechanisms mediating variant effects, and can support fine-mapping efforts to identifycausal variants. To evaluate REPRESS’s capacity to determine the effect of sequence mutations on miRNA binding sites, we benchmarked its ability to discriminate a novel curated set of miRNA binding-altering variants against a negative set of background variants from gnomAD

[0065] . REPRESS successfully identified individual variants known to disrupt (rs1876439052)

[0066] (FIG.4A) and create (rs1063320)

[0067] (FIG.4B) miRNA binding sites, and accurately distinguished the miRNA binding altering variants from the background variants. REPRESS outperformed existing miRNA binding models with an AUROC of 0.73 (FIG. 4C) while the next best model TargetScan

[0013] attained an AUROC of 0.59. For ”pathogenic” and ”likely pathogenic” (P / LP) variants from ClinVar [68, 69], REPRESS predicted higher scores for P / LP variants expected to modulate miRNA targeting relative to other P / LP variants and putative benign variants (FIG.4D).

[0270] We further investigated REPRESS’s ability to identify the underlying causal mechanisms behind disease pathologies and its utility in supporting therapeutic target discovery. The rs712 variant associated with increased risk of non-small cell lung cancer (NSCLC) was experimentally shown to impact expression by disrupting a repressive let-7 miRNA binding site in the 3’ UTR of the oncogenic KRAS gene [70-72]. REPRESS successfully made a high-scoring prediction that this variant would decrease miRNA binding. We also evaluated REPRESS predictions against a set of 24,5743’ UTR variants of uncertain significance (VUSs) and found that rs886049197 was also predicted to disrupt let-7 binding, suggesting a similar mechanism through which KRAS upregulation may be conferred (Examples: Methods). The other methods that we tested had low precision and scored these variants more weakly than they scored known benign variants, and they also did not exhibit statistically significant enrichment when comparing P / LP to benign variants. To evaluate more broadly whether REPRESS’s predictions can generalize to novel sequences and variants, we used MPRA datasets from McGeary et al.

[0052] and Slutskin et al.

[0051] to evaluate REPRESS’s ability to decode the impact of variants on miRNA mediated repression. The McGeary MPRA consists of 952 synthetically designed constructs of 120 bp containing either a canonical or non-canonical binding site of let-7a. The non-canonical sites used in the MPRA library included imperfect seed matches with varying degrees of 3’ supplementary binding. REPRESS accurately predicted the fold repression for the synthetic sequences from theMPRA (R=0.68) (FIG.4E) and outperformed other miRNA binding models by a large margin (FIG.4F).

[0271] The Slutskin MPRA dataset measured changes in RNA and protein levels across multiple cell lines for 12,545 regulatory sequences that modulate miRNA binding by varying target complementarity, target multiplicity, and other sequence features. REPRESS successfully recovered the experimentally measured variant effects as shown by a log-linear relationship between predicted miRNA binding and observed fold-change (FIG. 2F), outperforming other competing models (FIG.4F). It is notable that the sequences evaluated from both MPRAs fall outside of the distribution of wild-type sequences that REPRESS was trained on, but even so, REPRESS was able to accurately generalize its predictions to these datasets in a ”zero-shot” manner. Example 5: REPRESS facilitates the efficient design of RNA therapeutics

[0272] By pinpointing post-transcriptional regulatory elements, REPRESS can be used to inform the design of RNA therapeutics that target these elements to modulate gene expression. For instance, the activity of antisense oligonucleotides (ASOs) can be simulated with REPRESS to identify those likely to increase gene expression upon blocking of a repressive regulatory element (FIG.5A) (see Examples: Methods). Two instances exemplify how REPRESS can predict ASO effects: In one study of the UTRN gene, the target region bound by the ASO conferring the highest increase in expression

[0073] corresponds to a strong miRNA binding site predicted by REPRESS (FIG.2G). Additionally, when ASO binding was simulated by masking the sequence in this region, REPRESS correctly predicted the loss of the miRNA target. Similarly, in a separate study focused on the CHD9 gene

[0074] , REPRESS predicted that the lead ASO blocks a strong miRNA binding event without affecting other miRNA events in the same 3’ UTR (FIG.5B). In both of these examples, the lead ASO was ranked by REPRESS as rank 1 and rank 2 out of all possible 2,037 and 2,606 ASOs, respectively.

[0273] To benchmark ASO predictions, we curated a total of 27 ASO hits from 22 studies reported to increase RNA or protein levels of their gene targets by at least 1.5-fold by blocking miRNA binding sites (FIG. 5A). For each possible per gene oligo screening budget, we determined the fraction of the 27 ASO hits that would have been found if we hadscreened every gene using only the budgeted number of ASOs with the highest model predicted scores (FIG. 5C). REPRESS ranked the ASO hits highest among all methods tested, with 18 of the 27 ASO hits scoring in the top 10% of their distributions. Relative to TargetScan, REPRESS required more than four times fewer ASOs on average. None of the existing miRNA models evaluated showed a significant improvement over a naive seed matching approach.

[0274] In a prospective study of PON1, which is involved in cholesterol metabolism and implicated in diseases including familial hypercholesterolemia [75, 76], we designed 485 ASOs targeting the 3’ UTR and tested them in primary human hepatocytes (PHHs), measuring the corresponding increase in PON1 expression using an AlphaL-ISA assay (PerkinElmer) [77, 78] (FIG.5A and Examples: Methods). Of the 485 ASOs, we found 23 hits that increased the expression of PON1 by more than a therapeutically relevant threshold of 1.5 fold. We used REPRESS to rank all possible ASOs that could target the 3’ UTR of PON1 and found that REPRESS identified 80% of the hits among its top 22% of predicted ASOs, while for the next-best model TargetScan this hit rate of 80% required 56% of predicted ASOs( 2.5× more) (FIG.5D). This highlights REPRESS’s ability to prioritize regulatory regions inthe 3’ UTR and efficiently identify ASOs that can upregulate expression of their target gene.

[0275] Finally, we turned to another therapeutic modality, that of synthetic mRNA and gene therapies. REPRESS’s sequence attribution methods can be used to identify bases that drive miRNA binding or mRNA degradation, enabling the iterative design of sequences predicted to be more stable. As a proof-of-concept study, we selected 116 natural transcripts and progressively introduced single nucleotide changes to their 3’ UTRs using ISM, resulting in a substantial decrease in predicted miRNA binding and degradation scores (FIG.5E). To validate the impact of these edits on mRNA stability, we scored each sequence variant using the Saluki model

[0079] , a predictor trained using orthogonal half-life data from transcriptional inhibitor and pulse-labeling based protocols to predict mRNA stability from sequence. Half- lives predicted by Saluki increased by an average of 45% after incorporating all 14 edits nominated by REPRESS’s miRNA predictions, and by 21% after incorporating all 14 edits suggested by its degradome predictions (FIG.5E). An illustrative example of this is shown inFIG.5F, which shows 10 edits that REPRESS predicted would decrease miRNA binding in the ULK13’ UTR by 97.5%. Example 6: Expanding the model to predict cell type specific miRNA binding for any tissue type

[0276] The reliance on CLIP-Seq data from specific cell lines to predict molecular phenotypes such as miRNA binding constrains its applicability to a limited set of cell types. To overcome this constraint and broaden REPRESS's utility, a training strategy was developed aimed at enhancing the model's ability to extend its miRNA binding predictions to any cell line, utilizing just miRNA-Seq data relevant to the target cell type. By adapting the model to make predictions based solely on miRNA-Seq data, unlocks the potential to tap into the vast reservoir of publicly available miRNA-Seq datasets. This approach not only circumvents the limitations imposed by the availability of miR-eCLIP data but also significantly expands REPRESS's predictive reach to encompass a wider array of cell lines across diverse cellular contexts.

[0277] Each output track of the model represents the binding score for a specific miRNA, rather than aggregating the binding scores of all miRNAs within a particular cell line. By analyzing the top 100 expressed miRNAs and their binding sites across our dataset of 10 miR-eCLIP cell lines, REPRESS is able to predict the binding behavior of approximately 600 unique miRNAs spanning both humans and mice. When applying REPRESS to a new cell line, the miRNA-Seq data is used to calculate a weighted average of the binding predictions for the top 100 miRNAs expressed in that cell line. This approach not only streamlines the prediction process but also ensures that REPRESS's insights are tailored to the specific miRNA expression profile of any given cell line. Remarkably, this approach maintains 98% of the predictive performance of the cell-type specific REPRESS. Through this miRNA- specific prediction framework, REPRESS significantly broadens its application scope, enabling precise and relevant miRNA binding predictions across a diverse range of cellular contexts. MethodsCell culture

[0278] HepG2 cells were cultured in Dulbecco’s Modified Eagle Medium (DMEM), high glucose, pyruvate [Gibco: 11995065] supplemented with 10% FBS [Gibco: 16140071], 1% penicillin / streptomycin [Gibco: 15140122] and grown at 37°C with 5% CO2. K562 were cultured in supplemented Iscove’s Modified Dulbecco’s Medium (IMDM) [Gibco: 12440053] while A549 cells were cultured in DMEM (Gibco: 11965092) supplemented with 10% FBS. Primary Human Hepatocytes (PHHs) were acquired from BioIVT for one male donor and one female donor. They were thawed in Cryopreserved Hepatocyte Recovery Medium (CHRM) [Gibco; CM7000], plated in InVitroGRO CP Hepatocyte Medium (BioIVT: Z99029) supplemented with ROCK inhibitor- Y-27632 [Tocris Small Molecules: 1254 / 1] for 24 hours, the media was then changed to Cellartis Power Primary HEP Medium (Takara Bioscience: Y20020) and the cells were maintained in culture until day 7.

[0279] TET-ON NGN2 transgenic iPSCs were grown on plates coated with Vitronectin XF™ [Stemcell: 07180] and grown in mTeSR™ Plus media [Stemcell: 5825] supplemented with Rock Inhibitor. For neural induction, TET-ON NGN2 transgenic iPSCs were cultured on matrigel (VWR: 354234) -coated flasks for 3 days in vitro (DIV) in doxycycline (Millipore Sigma: D9891) and N2 supplement (Gibco: 17502048). Following neural induction, neuronal cultures were maintained up to 14 DIV in neuronal maintenance medium consisting of NeuroBasal Plus (Gibco: A3582901), B27 Plus (Gibco: A3582801), GlutaMAX™ (Gibco: 35050061), MEM Non-Essential Amino Acids (Gibco: 11140050), supplemented with BDNF 10ng / ml, NT-3 CC095). miR-eCLIP

[0280] miR-eCLIP data for HEK293T and 8 weeks old C57 / BL6J mouse liver tissue were curated from

[0034] . Huh-7 data was curated from

[0080] . A549 and K562 data were directly purchased from Eclipsebio (San Diego, USA). miR-eCLIP data for all additional cell lines and tissues were generated by Eclipsebio, following a protocol developed by Van Nostrand et al. (2022) to identify miRNA-Ago2 chimeric pairs. All cell lines were cultured and processed by the authors, whereas whole tissue samples were prepared directly by Eclipsebio. The miR- eCLIP protocol was conducted on duplicate samples of HepG2, iPSC derived neurons, andprimary human hepatocytes; and on triplicate samples of Yecuris human liver cells (sourced from Yecuris Corporation), post natal 2-day old C57 / BL6 mouse cortex tissue (P2), and 8- weeks old C57BL / 6J mouse cortex (sourced from the Jackson Laboratory). Each replicate of HEPG2 consisted of 20 million cells, all other cell samples consisted of 40 million cells. The whole tissue samples were prepared directly by Eclipsebio. For all in-house cell preparations, cells plated in 10-cm dishes were covered with DPBS, no calcium, no magnesium [Gibco: 14190144] and UV-irradiated with 400 mJoules / cm2 using a 254-nM bulb. Samples were then scraped and collected, flash frozen, and packaged for sample transfer. Additional details on the miR-eCLIP protocol are available in

[0034] . Degradome-Seq

[0281] Triplicate total RNA samples from 3 million cells of A549, K562, HepG2, iPSC- derived neurons, day0 primary human hepatocytes, day7 primary human hepatocytes and -550 mercaptoethanol [Gibco: 21985023] using the RNeasy [Qiagen: 74104] kit. Triplicate total RNA samples of 20 mg of tissue derived from Mouse liver 8-weeks old, Mouse cortex 8- weeks old, Mouse CNS e18 day 4 and Mouse CNS e18 day 11 were homogenized and also extracted using the RNeasy kit. Briefly, the degradome libraries were generated as follows: Poly-T magnetic beads were used to purify polyadenylated RNAs, followed by RNA ligase- mediated adapter ligation to the uncapped 5’ ends of 3’ cleavage products. The libraries were size selected using AMPureXP beads and the cDNA was PCR amplified yielding a library of 200-400 bp fragments. Libraries were sent to LC Sciences (Houston, USA) for 50-bp single- end sequencing on the Illumina Hiseq2500

[0081] . Small RNA-Seq

[0282] Total RNA extracted by RNeasy [Qiagen: 74104] was collected from 1 million cells in triplicate for A549, K562, HepG2, iPSC-derived neurons, day0 primary human hepatocytes and day7 primary human hepatocytes. Total RNA was also processed from 40 cortices of post-natal 2-day old C57 / BL6 mouse tissue (P2) sourced from the Jackson laboratory, and from 5 Yecuris liver samples obtained from Yecuris Corporation. Further processing and sequencing was conducted by The Centre for Applied Genomics (TCAG), Sick Kids Hospital (Toronto, Canada). Briefly, the small-RNA libraries were prepared with theNextera XT DNA Library Preparation Kit (Illumina) and sequenced by the Illumina HiSeq 2500 platform generating at read lengths of 125 bp and depth of 60M reads per sample. Data generation and processing

[0283] Curation of publicly available AGO2-CLIP datasets was performed in order to create a comprehensive training dataset. Datasets were discovered by searching the Gene Expression Omnibus (GEO) Database using keywords including ”AGO2-CLIP” and ”PAR- CLIP”, ”HITS-CLIP”, ”eCLIP” in conjunction with ”miRNA”. We only included studies where AGO2-CLIP was performed in untreated cells, ensuring the data reflected miRNA binding in the cell’s unperturbed, wild-type state. Raw fastq files from each study were downloaded and processed separately. Only studies that yielded over 10,000 non-overlapping AGO2-CLIP peaks across all replicates after trimming, genome alignment, and peak calling were kept. Studies that did not yield sufficient peaks due to low sequencing depths or low genome alignment rates were not included in the final dataset to ensure only high quality data is used for model training. This yielded a total of 18 AGO2-CLIP datasets across human and mouse cell lines curated from the literature.

[0284] AGO2-CLIP data curated from the literature was processed in a protocol specific manner. For PAR-CLIP fastp

[0082] was used to trim the reads with the following arguments : --adapter sequence XXX -w 16 --trim poly x --cut tail 30. Once trimmed, the reads were aligned to the genome using the bowtie aligner

[0083] with the following parameters : -v 2 -m 1 --best --strata followed by peak calling using PARalyzer

[0084] .

[0285] For HITS-CLIP and eCLIP data, fastp was also used for trimming reads with the same parameters while bowtie2

[0085] was used to align the reads to the genome with following arguments : -D 20 -R 3 -N 0 -L 20 -i S,1,0.50 -k 4 -p 32. After genome alignment, the reads were de-duplicated using picard and peak calling was done using CLIPper

[0086] .

[0286] All miR-eCLIP data was generated and processed by EclipseBio. UMIs were first extracted using umi-tools

[0087] and the adapters were trimmed using Cutadapt

[0088] . The chimeric reads were then reverse-mapped to a database of mature miRNA sequences from miRBase

[0089] using bowtie to get the miRNA identity and then aligned to the genome using the STAR aligner

[0090] . Peaks were then called on the chimeric reads using CLIPper.

[0287] All Degradome-Seq data was processed using the same pipeline. The reads were first trimmed using the fastp trimmer with the following arguments : --adapter sequence XXX --length required 20. After trimming, the reads were then aligned to the genome using HISAT2

[0091] , a splice aware aligner. Peak calling was not performed on the degradome data and the model is trained to predict the Degradome-Seq read coverage as a function of the input RNA sequence.

[0288] All miRNA-Seq data was processed using the following pipeline. The small RNA reads were first trimmed using cutadapt with -m 15. Once trimmed, the reads were aligned to the genome using the STAR aligner

[0090] with the following parameters : -- alignEndsType EndToEnd --outFilterMismatchNmax 1 --outFilterMultimapScoreRange 0 -- quantMode TranscriptomeSAM GeneCounts --outReadsUnmapped Fastx --outSAMtype BAM SortedByCoordinate –outFilterMultimapNmax 10 --outSAMunmapped Within – outFilterScoreMinOverLread --outFilterMatchNminOverLread 0 --outFilterMatchNmin 16 -- alignSJDBoverhangMin 1000 --alignIntronMax 1 --outWigType wiggle --outWigStrand Stranded --outWigNorm RPM. The featureCounts package

[0092] was then used to count the number of reads mapping to different pre-miRNAs and mature miRNAs in the genome with the following parameters -M -O -s 1 --minOverlap 3. After this the counts were transformed to CPM values, the final CPM for a given miRNA was taken as the average across all the replicates. Model Architecture

[0289] For a RNA sequence of length L, we first transform it into a one-hot-encoded sequence (A = [1,0,0,0], C = [0,1,0,0], G = [0,0,1,0], T = [0,0,0,1]) to obtain an input matrix of size L × 4. N = [0,0,0,0] is used for the part that extends outside the transcript. This matrix is used as the input to the REPRESS model. The output of the model is a matrix of size L × K, where K denotes the total number of cell lines. The k-th row of the matrix represents the predicted miRNA binding probability or predicted degradome read coverage of each base pair for the k-th cell line.

[0290] The REPRESS architecture consists of base layers and output layers, the core components of which are residual blocks composed of dilated convolution layers. All convolution layers have the same number of filters, denoted by N, except for the last one inthe output layers, which is determined by the number of cell types to predict. Padding is used in convolution layers so that the output sequence has the same length as the input sequence after convolution layers.

[0291] The input is first processed with an initial convolution layer and then goes through 16 residual blocks (RB) with gradually increasing convolutional dilation rate (receptive field), the dilation rate for each of the 16 residual blocks are - 1, 1, 1, 1, 4, 4, 4, 4, 10, 10, 10, 10, 25, 25, 25, 25. In each residual block, layer normalization is applied first, the output of which is passed through two separate convolution layers with sigmoid and tanh activation functions. Then the two outputs will go through an element-wise multiplication operation. The convolution and element-wise multiplication is repeated twice. Finally, the input is added to the output as a skip connection. All the residual blocks have the same number of filters, and each set of four residual blocks has filter widths of 11, 11, 21 and 41. The output of the initial convolution layer, and the 4th, 8th, 12th and 16th of the residual blocks are then passed through a convolution layer separately. These five outputs are accumulated, concatenated and passed into the output layers of the model.

[0292] The output layers consist of four components based on the data species and data type: human-miRNA, human-degradome, mouse-miRNA, and mouse-degradome.

[0293] Each component has an identical structure, made up of an initial convolution layer and four residual blocks, except for the final output convolution layer. In each residual block (RB-H) of the output module, a layer normalization layer is applied first, followed by two convolution layers with GELU activation. Then the input to the residual block is added to the output as a skip connection. The four residual blocks have the same number of filters as the ones in the base layers, but have filter widths of 11, 11, 21 and 41, with dilations of 1, 4, 10, and 25, respectively. Before the final output convolution layer, the previous output is cropped at both ends by half the context length (receptive field of the model) so that the final output has the same length as the initial input to the model without the added context sequence. For the last convolution layer, the number of filters is based on the number of cell lines in for the corresponding output, and the activation functions used are sigmoid for miRNA binding prediction and softplus for degradation prediction.

[0294] REPRESS’s novel residual block (RB and RB-H) is inspired from changes to the residual block proposed in the ConvNeXt publication

[0042] which improved the performance on CNNs on ImageNet helping it achieve state-of-the-art performance on par with transformers. The changes involved replacing batch normalization layers with layer normalization and reducing the total number of normalization layers to only 1 per residual block instead of 2. Also, the GELU activation was used in the residual blocks (RB-H) of the output layers instead of the commonly used ReLU activation function.

[0295] These changes alone boosted the validation performance of REPRESS on predicting miRNA binding and mRNA degradation by 17% and 10% respectively (FIG.6G). Model training and evaluation

[0296] REPRESS was trained using annotations from ncbi refseq.v109 for the human cell lines and ncbi refseq.m38.v106 for the mouse cell lines from GenomeKit. Only protein- coding transcripts were used to train the model. For each gene, the most principal transcript (according to APPRIS) was used. For genes with no APPRIS annotations, the first transcript in GenomeKit’s gene.transcripts table was used. Each transcript from both the human and mouse transcriptomes were divided into nonoverlapping windows of 3000 nucleotides across its mature mRNA sequence. A batch of 10 windows was randomly sampled, and 6.25 kb of transcript context sequence around the sampled window was appended to each side, resulting in a total input sequence length of 15.5 kb used to train the model. REPRESS’s miRNA predictions were trained on all exonic CLIP peaks, and no filtering based on miRNA identity, conservation, or site type was performed.

[0297] For a given sequence, the true labels corresponding to the miRNA outputs are a binary matrix of shape (3000, 29) corresponding to the length of the sequence and the total number of human and mouse cell lines / tissues with miRNA binding data and predictions. The binary matrix has a true label of 1 for every location corresponding to a CLIP identified peak and 0 everywhere else. Similarly, the true labels for the degradome outputs are a real valued matrix of shape (3000, 10) corresponding to the degradome-seq read coverages for the given input sequence and the total number of human and mouse cell lines / tissues with Degradome- Seq data / predictions. The degradome-seq labels are processed with squashed transformation (x0.375) to reduce the skewness.

[0298] REPRESS’s miRNA binding outputs are in the range (0, 1) while its degradation outputs are in the range (0, inf). We use binary cross-entropy loss for the miRNA binding output and Poisson loss for the degradation output. The total loss is the weighted sum ofthe two:(Equation 1)

[0299] When a sequence from a human transcript is sampled, only gradients from the corresponding human miRNA and degradome output heads are used to update the network, and vice-versa when a sequence from a mouse transcript is sampled. The final model is an ensemble of models trained on the four splits with different sets of chromosomes for training and validation. The four splits use chromosomes (1, 2), (6, 7, 8, 9), (10, 11, 12, 13), and (14, 15, 16, 17, 18) for validation (4-fold cross validation) respectively and the remaining for training. The splits were chosen as described to ensure that 75% of identified miRNA CLIP peaks were used for training and 25% for validation.

[0300] For evaluation, we first process predictions and labels of the validation intervals. For each cell line, we have a pre-calculated window length based on median peak width for the CLIP data, and use a fixed length of 50 base pairs for the degradation data.

[0301] First, we apply average pooling for both predictions and labels using the window length as the pool size. Then, we divide the predictions and labels into multiple chunks with the defined window length. For the miRNA predictions, we assign the value 1 to line, we calculate area under the curve (AUC) and area under the precision-recall curve (AUPRC) across all intervals for miRNA binding predictions, and Spearman for degradation predictions. We report the mean value of each metric across all cell lines as the final metrics.

[0302] The code was implemented in TensorFlow 2.11.0. We initialized the model weights using the Glorot Uniform distribution. For model training, we used the Adam optimizer1 2. We trained with a batch size of 5 on a single NVIDIA RTX A6000 GPU. We tuned the learning rate and the loss weights based on the evaluation metrics on the validation set, and also conducted early stopping in 60 training epochs based on these metrics.

[0303] The optimal learning rate we found is 0.00005 for the 133M parameter model, with1 2= 0.1. ReduceLROnPlateau is used for learning rate scheduling, monitoring on validation loss with factor=0.5, patience=5, and cooldown=1. Sequence attribution interpretation of REPRESS In-silico mutagenesis

[0304] Given a genomic interval starting and ending at positions pstartand pend, respectively, and a reference genome sequence x of length L, the attribution scores s are computed as follows: For each position i such that pstartend, and for each possible nucleotide u in the set of nucleotides { A, C, G, T}, we generate a variant of x by substituting the nucleotide at position i with u, denoted as xi,u. We then compute the ISM scores s for each position i and nucleotide u using REPRESS M, as given by the equation:

[0305] To visualize the effects across the interval, the ISM scores are averaged over all possible nucleotide substitutions, yielding the visualization scores( ):

[0306] The visualization scores ( ) represent the average impact of mutating eachposition within the interval to any of the four nucleotides. Smoothgrad

[0307] Let X RL×Adenote the input tensor, where L is the sequence length and A represents any of { A, C, G, T}). The SmoothGrad saliency map S is computed by:

[0308] 1. Generating Noisy Samples: We create N noisy samples of X by adding noisedrawn from a normal distribution N( , 2):

[0309] 2. Computing the Saliency Map: We calculate the gradients of REPRESS’s output with respect to( ), optionally focused on a particular class index c when a specific output cell-type head:

[0310] where f is the function applied to the saliency map and M denotes the REPRESS score.

[0311] 3. Averaging the Gradients: The final SmoothGrad saliency map S is obtained by averaging the model gradients over all noisy samples:miRNA frequency of complementarity map

[0312] We ran LinearCoPartition

[0038] with the central 100nt of each CLIP peak and the sequence of the corresponding mature miRNA using default settings to get the MFE folded structure and accessibility values for each miRNA-mRNA pair, where the miRNA was between 21-23 nt. After folding, we rejected any structure pair where the miRNA showed intramolecular folding to itself, as this is unlikely to occur while loaded into the RISC complex, and any pair that was predicted to have fewer than four paired positions, as this may represent noise. miRNAs with fewer than 20 instances in each cell-type were also dropped. For each miRNA in each cell type, we counted the average number of times that a position in the miRNA was predicted to be paired and clustered the miRNAs using scikit-learn with n clusters=5, using a k-means++ center initialization. Each cluster was further organized using hierarchical clustering with method=”ward”, to be then plotted using seaborn.clustermap. Transcriptome wide degradation map

[0313] We used hierarchical clustering to construct a heatmap of the Degradome-Seq profiles along all human protein-coding genes. For each gene, the most principal transcript (according to APPRIS)

[0093] from RefSeq (v109) was used. For genes with no APPRISannotations, the first transcript in GenomeKit’s gene.transcripts table was used. Degradome- Seq read coverage for each replicate was first normalized to the maximum value anywhere along the transcript, then averaged across replicates, and finally re-normalized to the maximum normalized average value. Replicates with no read coverage for a given transcript, and transcripts with no read coverage in any replicates, were excluded from analysis. To enable comparison of position-specific Degradome-Seq profiles across transcripts of varying lengths, each transcript was split into 205 segments: 15 segments for the 5’ UTR; 100 segments for the CDS, and 90 segments for the 3’ UTR. The number of UTR segments was based on the average length of these regions divided by the average CDS length. These values were then multiplied by 100 (the arbitrarily selected number of CDS segments) and rounded to the nearest five. Transcripts with a 5’ UTR less than 15 bases, CDS less than 100 bases, or 3’ UTR less than 90 bases were excluded from analysis. Degradome-Seq read counts were then averaged across all positions within a segment for clustering using seaborn.clustermap(method = ‘ward’). miRNA target counts were obtained from cell line / tissue-matched miR-eCLIP peaks. Half-life data (for available cell lines) based on 4sU data [94, 95] was obtained from Table S1 from

[0079] . Degradome and RNA element overlap analysis

[0314] To profile Degradome-Seq signal at various potential regulatory sites, we analyzed RBP binding sites for all RBPs with ENCODE eCLIP data

[0039] , miR-eCLIP peaks, and an AU-rich element (ARE) motif. For each of HepG2 and K562 cells, we assigned each peak / site to the most highly expressed overlapping Gencode (v29) transcript based on ENCODE RNA-Seq data. Transcripts were required to have at least one position with at least five Degradome-Seq reads to be included for analysis. Sites in introns, within 100 bases of the transcriptional start site, or within 500 bases of the poly(A) site were excluded from analysis. Degradome-Seq read counts from 100 bases upstream to 100 bases downstream of each site were averaged across replicates. To enable comparison across many transcripts with different baseline levels of expression / degradation, read counts were normalized to the read count at position -100. For each site, a control site was randomly selected from the same region (5’ UTR, CDS, or 3’ UTR) of the same transcript and analyzed in an identicalfashion. While RBP binding sites were analyzed at the level of individual RBPs, all miRNAs were analyzed collectively due to data sparsity. Validation set analysis

[0315] REPRESS was benchmarked against TargetScan v8, Biochemical model, miRAW, DeepMirTar, miTAR, miRBind and DMISO for predicting miRNA binding on its validation set. To benchmark TargetScan, for each transcript in the validation set we run TargetScan on the entire transcript and obtain context++ scores for every identified miRNA binding site for the top 50 miRNAs expressed in that cell line. TargetScan scalar values for each binding site were parsed into prediction tracks, where the interval spanned by an identified binding site is assigned the context++ score as the predicted score. Overlapping context++ scores (caused by different miRNAs predicted to bind to the same region) were averaged position wise. We do not provide 3’UTR MSA alignment input to TargetScan due to the difficulty in locating MSA alignments for all 4k transcripts in the validation set. Instead, we perform a lookup to precomputed Pct

[0055] values for each 3’ UTR-binding site, obtained from the TargetScan website. AIRs (affected isoform ratio for a binding site) are also obtained from a lookup table. The precomputed AIRs data was obtained from the ”3P- seq tag info” file from the TargetScan website. If no value was found then the value is set to 1. AUC’s are calculated on a window level classification task similar to REPRESS. Similarly, the top 50 expressed miRNAs were used to calculate the transcript wide miRNA binding using the biochemical model.

[0316] Unlike REPRESS, miRNA binding models miRAW, DeepMirTar, miTAR, miRBind and DMISO predict miRNA binding as a function of both the target sequence and the corresponding miRNA. miRAW models target sequences of size 40 nt; DeepMirTar and miTAR model target sequences of size 53nt; miRBind and DMISO models target sequences of length 50 and 60 respectively. In order to benchmark these models on the REPRESS validation set, we queried the top 30 miRNAs in each cell line against each of the transcripts in the validation set. For each miRNA-transcript pair we slide a window corresponding to the target sequence size of the model across the entire transcript with a stride of 1. Each window- miRNA pair is now scored with the model as this score is used to generate the prediction track across the entire transcript. We further tested querying only those windows that containa seed region of the corresponding miRNA and assigned a 0 score to any window that did not have a seed region of the corresponding queried miRNA, this consistently yielded better performance for all models. Predicting experimentally validated miRNA binding sites

[0317] For this task, miRNA binding sites from miRTarBase and DianaTarBase were cross-referenced to generate the positive set of cell type specific experimentally validated miRNA binding sites. DianaTarBase contains information about miRNA-target interactions and the corresponding cell type but does not tell you the exact interval of miRNA binding, while miRTarBase contains information about exact binding interval and the corresponding miRNA without cell type information. Hence, binding sites that were present in both these datasets were used to generate a cell type specific set of experimentally validated miRNA binding sites. For the negative set, we took intervals corresponding to 8-mer seed region matches of miRNAs that are not expressed in those corresponding cell lines ( Hence, using both the positive and negative sets we were able to benchmark REPRESS against this task by computing the corresponding AUCs and AUPRCs. For each the baseline models TargetScan, biochemical, miRAW, DeepMirTar, miTAR, miRBind and DMISO were directly queried on each target sequence miRNA pair for both the positive and negative sets. The predicted output scores were used to compute the AUCs and AUPRCs. Predicting HEAP identified miRNA binding sites

[0318] We analyzed miRNA binding sites identified in the HEAP dataset (GSE139344) to predict miRNA-target interactions in exonic regions. Specifically, we focused on peaks detected in the P13 mouse cortex samples (GSE139344 P12 cortex peaks.csv), as this was a cell type that was available in REPRESS (Mouse cortex p2). To ensure relevance, our analysis was restricted to miRNA binding sites located in exonic regions, removing intergenic and non-coding RNA features. Genomic sequences extracted using mm10 (ncbi refseq.m38.v106) for the corresponding intervals. REPRESS predictions were the average Mouse cortex miRNA binding prediction across the entire binding interval identified by HEAP.

[0319] We benchmarked against TargetScan, Biochemical, miRAW, DMISO, miRBind Deep-MirTar and miTAR incorporating flanking sequence to meet the length requirements ofeach model’s input format (e.g., 40-60 nt for miRAW, DeepMiRTar, DMISO). Each model was applied to each given sequence and corresponding binding miRNA identified from HEAP. Performance of each model was assessed by stratifying the predictions based on bins of log2FC values using predefined thresholds (<2, 2-3, 3-4, 4-5, 5-6, 6-7, 7+) to evaluate model performance in identifying repressive vs non repressive sites across the range of log2FC reported. We report ROC and for each model across the range of log2FC thresholds (FIG.8).

[0320] In addition to repressive vs non-repressive peak identification analysis, we also performed an analysis for identifying binding sites identified from HEAP vs non-binding sites. Similar to the validation set analysis, we predicted transcriptome wide binding for Mouse cortex, binned the predictions into windows of 71, and computed the AUPRC in predicting HEAP identified sites from background sites (FIG.7E). REPRESS performance comparison heatmap

[0321] All the performance values presented in the heatmap are normalized as a percentage of the highest performing model for that particular task. The metrics used for each of the following tasks are as follows - The validation set analysis of REPRESS used the average AUPRC across the 29 human and mouse cell lines; unnormalized values of the model performances can be found in FIG.7A. The miRTarBase prediction analysis used the average AUROC across the 8 cell lines for which data was cross referenced from miRTarBase and DianaTarBase; unnormalized AUPRC values for each of the models can be found in FIG. 7D. The HEAP analysis used AUPRC values in from predicting HEAP identified binding sites from non-binding sites; unnormalized values for this analysis can be found in FIG.7E. Both the miRAW and DeepMirTar test sets report AUPRC values for all the miRNA models; unnormalized values can be found in FIGs. 7B-7C. The McGeary MPRA dataset uses the Pearson correlation for predicted scores and repression for the 952 sequences from the let-7a transfection MPRA. The Slutskin MPRA uses the average Pearson correlation for predicted scores and repression across 4 values - K562 protein readout, K562 RNA readout, HEK 293 RNA readout, HepG2 RNA readout; unnormalized values for all models can be found in FIG. 4F. The variant effect prediction task uses theAUPRC values of predicting validated miRNA binding altering variants from background variants. REPRESS performance comparison heatmap

[0322] For each of the top 10 miRNAs expressed in each of the cell lines in our dataset, we took the top 500 predicted miRNA binding sites with a seed region match and computed the ISM sequence attribution scores for each of the binding sites. The sequence attribution scores for each miRNA’s corresponding binding sites were averaged across all of its top 500 sites and the average score at each position was used to determine the impact of complementarity at that particular position of the miRNA. Atypical miRNA binding analysis

[0323] REPRESS sequence attributions were performed on the highest scoring 100 miReCLIP targets according to REPRESS for each of the top 20 most highly expressed miRNAs in each cell line / tissue, as described in Examples: Methods. Sequence attribution scores were aligned, averaged across all targets, scaled to the largest value, and presented in terms of the corresponding mature miRNA sequence (i.e. the reverse complement of a perfectly paired target). Average attributions < 0 were plotted as 0 for clarity. For each miRNA with at least 250 unique seed-containing miR-eCLIP targets across all human cell lines profiled with miR-eCLIP, the proportion of canonical targets (including offset 6mer, 6mer, 7mer-A1, 7mer-m8, and 8mer site types) with an offset 6mer site type was determined. A two proportion Z-test was used to compare the proportion of offset 6mer sites between hsa-let- 7a-5p and hsa-miR-148a-3p.

[0324] Conservation was calculated as the proportion of target bases opposite miRNA positions 2-7 with a Phylo-P 100-way score > 3. For hsa-let-7a-5p and hsa-miR-148a-3p, statistically significant groupings of site types were identified using two proportion Z-tests and Benjamini-Hochberg (BH) adjustment for multiple hypothesis testing. Comparison of 6mer targets and 8mer targets of miRNAs

[0325] All human miRNAs with at least 50 unique 6mer targets and 50 unique 8mer targets in our miR-eCLIP dataset were eligible for this analysis. For each target, we determined the number of conserved bases in the non-seed region–corresponding to bases-13 through -2, where base 0 is the most 5’ base of the 6mer seed match. Bases with PhyloP 100-way values > 3 were considered to be conserved. For each miRNA, the odds ratio for the odds of a base being conserved versus non-conserved in 6mer targets versus 8mer targets (and the associated BH-adjusted P values and confidence intervals) were calculated using the statsmodels Python package. Sequence attributions were generated separately for 6mer and 8mer targets of all highly expressed (top 20 in any cell line) miRNAs with absolute log2 odds ratios > 0.5 (n = 6), for comparison of sequence attributions and conservation rates in the non-seed regions of 6mer and 8mer targets. miRNA binding multiplicity analysis

[0326] The Slutskin et al. paper included 557 MPRA constructs specifically to analyze the effect of miRNA binding site multiplicity for 5 different miRNAs - hsa-miR-320a-3p, hsa- miR-21-5p, hsa-miR-92a-3p, hsa-miR-19b-3p and hsa-miR-20a-5p. Each MPRA construct was designed to have between zero and five miRNA binding sites for each the corresponding miRNAs The locations of the five possible binding sites were fixed across all sequences and in each construct, each of these locations each site contained either a miRNA binding site or one of two control context sequences.

[0327] For each of the five miRNAs, all sequences with the same number of miRNA target sites and the same context sequence were used to calculate average values for: i) observed experimental repression, and miRNA binding predictions using ii) REPRESS, iii) TargetScan, iv) Biochemical model, iv) miRAW, v) DeepMirTar, vi) miTAR, vii) DMISO, and viii) miRBind. Biochemical model fold-change predictions were multiplied by -1 to represent repression. All model predictions excluding REPRESS were scaled to the maximum value predicted by each model across the entire dataset. P values for the impact of context and multiplicity on REPRESS predictions were calculated for all sequences with < 5 miRNA targets using formula.api.ols() with ”score context + multiplicity” and stats.anova lm() with ”type=2” from the statsmodels Python package. miRNA transfection and knockout datasets

[0328] We systematically processed datasets from GSE127211, GSE197363, GSE97060 and GSE123311 to extract log2FC information pertaining to each miRNA-mRNA pair.

[0329] For each transcript, we first calculated REPRESS miRNA binding prediction across the entire 3’ UTR, then using the sum of REPRESS predictions across all identified miRNA seed sequences we combined these predictions into a single composite score corresponding to the transcript. With this composite score, we ranked the transcripts andselected the top 100 (N=100) with the highest predictive scores to construct the eCDFs.= ({ ( } ) (Equation 2)

[0330] where Sseedrepresents the set of scores for each identified seed sequence of the miRNA being tested. To obtain the TargetScan TargetScan was run on every transcript- miRNA pair. The repression score for the transcript is the sum of TargetScan context++ scores for all identified binding sites in the 3’UTR. Again we ranked the transcripts and selected the top 100 to construct the eCDFs. Positional preference of miRNA targeting analysis

[0331] For this analysis, we first intersected miR-eCLIP peaks with human 3’ UTRs from principal transcripts (according to APPRIS)

[0093] of protein coding genes in v109 of NCBI refseq. For genes with no APPRIS annotations, the first transcript in GenomeKit’s gene.transcripts table was used. Peaks only from miRNAs annotated with conservation > -1 in miR Family Info.txt and with > 500 peaks in the cell line of interest were considered for this analysis. Peaks passing these filters were then assigned a bin between 1 and 100 based on the relative position of their seed match (or middle of the peak for peaks with no seed matches) along their 3’ UTR, where the first base downstream of the stop codon corresponds to bin 1 and the last base before the poly(A) tail corresponds to bin 100. For each miRNA, the number of peaks in each bin was divided by the total number of peaks. Min-max normalization was applied to these values to enable comparison across miRNAs. These miRNA-specific values were clustered and plotted with seaborn.clustermap(). Predicting validated miRNA binding altering variants

[0332] We curated the literature to create a dataset of 100 experimentally validated miRNA binding sites altering (disrupting or creating) variants. This dataset was used as the positive set in order to benchmark REPRESS’s performance on predicting miRNA binding altering variants. For the negative set, common variants from gnomAD with an allele set were taken. The assumption here being, common variants are likely benign and do not disrupt any major underlying cellular processes. This dataset was then used to benchmark the performance of REPRESS on predicting miRNA altering variants.

[0333] In order to compute the score of a variant, we computed the mean prediction across the miRNA binding interval for the wild type and mutant prediction and took the absolute difference between these two values. We also averaged this difference across all of REPRESS’s human cell line predictions to generate a scalar value corresponding to the impact score of that variant. For all baseline models, we considered a union of the top 50 miRNAs expressed in all the cell lines and calculated the mean miRNA binding prediction across for the top 50 miRNAs for the wild type and mutant sequences.

[0334] The absolute difference between the two mean predictions were used to score the variants. Application of REPRESS to pathogenic and likely pathogenic (P / LP) Clinvar variants that are curated to impact miRNA binding

[0335] The dataset consisting of pathogenic or likely pathogenic (P / LP) 3’ UTR variants and putative benign variants was obtained from

[0096] . For each variant, one prediction was made for the wildtype sequence and one was made for the mutant sequence as follows: The region in the genome affected by the variant was expanded 20bp on either side. This interval was given to REPRESS for inference and the human cell line track predictions corresponding to miRNA binding across the entire interval were extracted.

[0336] For each cell line head, the average across the prediction interval was taken to get a single value per cell line. The mean of the resulting cell line scalar values was taken to get a single mean scalar value for all the human cell lines. The variant effect was computed as the difference between the wildtype and mutant prediction (mut - wt).

[0337] The dataset was stratified into P / LP variants that were curated to impact miRNA binding (n=7), P / LP variants curated to operate through another mechanism or have an unknown mechanism (n=19), and putative benign variants (n = 67). We plot the absolute value of the variant effect and perform a Mann-Whitney Wilcoxon 2-sided test between the variant effect scores for each group Applying REPRESS to VUS 3’ UTR variants in Clinvar Variants from ClinVar (last accessed Apr. 30, 2023) that were classified as Variants of Uncertain Significance (VUS) and in 3’ UTR of the transcripts that they were reported in were extracted. These variants were scored using REPRESS in the same manner as the section: ”Applying REPRESS to pathogenic and likely pathogenic (P / LP) Clinvar variants that are curated to impact miRNA binding”. From the dataset of variants that impact miRNA binding and common variants (Fig 4b), we split the miRNA binding altering variants into those that create miRNA binding sites and those that disrupt miRNA binding sites. In separate classification tasks for these two sets of variants, we derived thresholds for the 5% FPR and annotated this on a distribution histogram of VUS variant scores from REPRESS. Applying REPRESS to VUS 3’ UTR variants in Clinvar

[0338] The Slutskin MPRA dataset included experimental measurements of repression in four cell lines: K562, HEK293, HepG2, and MCF7, covering 12,545 synthetically designed 176 nt sequences. Since REPRESS was trained on CLIP data from each of these specific cell lines, predictions were made directly for the corresponding cell lines and correlated with the measured fold change. REPRESS was queried using the 176 nt synthetic constructs, and the average prediction across the length of the construct was used to calculate the REPRESS score for each sequence.

[0339] The McGeary MPRA library consisted of 952 distinct 120 nt 3’ UTR fragments designed to include both canonical and non-canonical targets of let-7a with varying levels of 3’ complementarity. Let-7a was transfected into DMEM cells and used to measure the corresponding repression on each of the MPRA constructs. Since let-7a was transfected into the cells, we picked the cell line corresponding to the highest endogenous expression of let- 7a to query from REPRESS, which was Yecuris liver. REPRESS predictions were averagedacross the entire 120 nt sequence to generate the REPRESS score which were then correlated with measured values of repression from the MPRA.

[0340] Unlike REPRESS, which makes a direct sequence-to-binding prediction and does not require specification of any miRNA(s), existing models require an individual miRNA to be specified for each query. Therefore, we queried existing models using the top 100 expressed miRNAs in each of the cell lines from the Slutskin MPRA, averaging predictions across miRNAs to generate scalar values representing predicted repression for each MPRA construct. For the McGeary MPRA, each model was queried with let-7a specifically. The MPRA constructs were divided into overlapping windows, each with a stride of 1 and a length equal to the maximum input length tolerated by each model. Predictions were then averaged across all windows and miRNAs to produce a final scalar value representing predicted repression for each model. For miRAW, DeepMir-Tar, miTAR, miRBind and DMISO we further tested querying only those windows that contain a seed region of the corresponding miRNA and assigned a 0 score to any window that did not have a seed region of the corresponding queried miRNA, this consistently yielded better performance for the Slutskin MPRA for all models. Curated miRNA blocking ASOs Dataset and ASO prediction

[0341] 177 unique ASOs targeting miRNA binding sites were curated from the literature. Information like assay used to measure fold change, ASO chemistry, target gene, in vitro / in vivo and corresponding cell type used to measure ASO efficacy was also tracked. 27 out of the 117 curated ASOs increased the expression of their target gene by at least 1.5X. In order to measure the efficacy of REPRESS in identifying efficacious miRNA blocking ASOs we measured the reduction in number of ASOs we would have to screen in order to identify the hits when screening ASOs in decreasing order of their predicted score. We simulate the effect of an ASO with REPRESS by masking out the interval corresponding to regions complementary to the ASO and zeroing out the one-hot encodings of that interval. We then take the difference between the mean prediction with and without masking to measure the impact of that ASO in that corresponding cell line. If the corresponding cell line the ASO was tested in was not available we took the average difference across all cell lines.

[0342] For all baseline models, a similar approach as the validation set analysis was used to generate track-like predictions from the model using the top 30 expressed miRNAs in that cell line. The average prediction over the ASO interval is taken as the score for that ASO in that particular cell line. For miRAW, DeepMirTar, miTAR, miRBind and DMISO we further tested querying only those windows that contain a seed region of the corresponding miRNA and assigned a 0 score to any window that did not have a seed region of the corresponding queried miRNA, this consistently yielded better performance for all the models. PON1 screening SBO Dataset

[0343] We screened 485 PS-MOE ASOs designed to modulate miRNA interactions within the 3’ UTR of the PON1 gene (RefSeq transcript NM 000446). These were screened using an AlphaLISA assay [PerkinElmer: AL389C] to ascertain PON1 protein expression within primary human hepatocyte (PHH) cells. 18,000 PHH cells per well were reverse transfected with ASO using RNAiMAX [Thermofisher: 13778-150] in a 384-well format in InVitroGRO CP Hepatocyte Medium (BioIVT: Z99029) supplemented with ROCK inhibitor- Y-27632 [Tocris Small Molecules: 1254 / 1] for 24 hours, the media was then changed to Cellartis Power Primary HEP Medium (Takara Bioscience: Y20020) until day 7. The cells were then lysed in Alphalisa lysis buffer [PerkinElmer: AL003F] supplemented with 1X HALT protease buffer [Thermofisher: 78439] and samples were analyzed using the PON1 (human) AlphaLISA Detection Kit [PerkinElmer: AL389C].

[0344] To compute fold changes, ASO-treated wells were normalized to the mean of three wells treated with a non-targeting ASO control, after both sets of measurements were fit to a standard curve. ASO fold changes were then reported as the mean fold change across three replicates. ASO hits were defined as ASOs that achieved a mean fold change of at least 1.5x. At all five lengths ranging from 18 nt to 22 nt, we scored using REPRESS and TargetScan the 1268 possible intervals falling within the PON13’ UTR of the given length, yielding 6,340 scores for each model. REPRESS scores were computed as the mean difference between wildtype and N-masked sequence intervals, when used to predict drop in miRNA binding in PHH for this transcript.

[0345] TargetScan scores were computed as the mean score across the 50 most highly expressed miRNAs in PHH, with the score defined as the difference between the wildtype and N-masked sequence. For the 23 ASO hits exceeding the 1.5x expression threshold, we then determined where their corresponding interval scores fell within the distribution of all 6,340 interval scores for each model. REPRESS for synthetic mRNA design

[0346] We performed ISM on 116 protein-coding genes with high endogenous miRNA binding predictions whose 3’ UTRs were between 1500 and 2000 nt for the human cell lines. The squashing transformation in the degradome tracks was inverted using x1.0 / 0.375, where x is the model score. Cell lines corresponding to miRNA and degradome were averaged independently, to produce the average ISM result across all cell lines and tissues for degradome and miRNA, respectively. Using the ISM matrix, we then selected the position with the minimum most value (i.e. the mutation that would most reduce the model prediction); importantly, wild-type positions were masked so as to ensure a mutational event at each iteration and prevent getting stuck into local minima. At each iteration, for 14 cycles, the variant resulting in the minimum-most score of the model was selected and used as the input for the next cycle. We also utilized Saluki as our oracle for assessing sequence stability and half-life; for each iteration and mutation, we scored the sequence using Saluki with its default parameters and model weights from the original publication to evaluate improvements in stability. To analyze and visualize the changes in both Saluki and REPRESS scores relative to the wild-type (wt), unedited sequence, we first applied minmax scaling to the scores of each gene across all edits to standardize them. We then subtracted the score of the wt sequence from each edited sequence’s score.

[0347] The present invention has been described here by way of example only. Various modification and variations may be made to these exemplary embodiments without departing from the spirit and scope of the invention, which is limited only by the appended claims.

[0348] All publications, patents and patent applications are herein incorporated by reference in their entirety to the same extent as if each individual publication, patent or patentapplication was specifically and individually indicated to be incorporated by reference in its entirety.REFERENCES [1] Dassi, E. Handshakes and fights: the regulatory interplay of rna-binding proteins. Frontiers in molecular biosciences 4, 67 (2017). [2] Halbeisen, R. E., Galgano, A., Scherrer, T. & Gerber, A. Post transcriptional gene regulation: from genome-wide studies to principles. Cellular and molecular life sciences 65, 798–813 (2008). [3] Lewis, B. P., Burge, C. B. & Bartel, D. P. Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are microrna targets. cell 120, 15–20 (2005). [4] Bartel, D. P. Micrornas: target recognition and regulatory functions. cell 136, 215–233 (2009). [5] MacFarlane, L.-A. & R Murphy, P. Microrna: biogenesis, function and role in cancer. Current genomics 11, 537–561 (2010). [6] Xu, K., Lin, J., Zandi, R., Roth, J. & Ji, L. Microrna-mediated target mrna cleavage and 3- uridylation in human cells. sci rep 6, 30242 (2016). [7] Sheu-Gruttadauria, J., Xiao, Y., Gebert, L. F. & MacRae, I. J. Beyond the seed: structural basis for supplementary micro rna targeting by human argonaute2. The EMBO journal 38, e101153 (2019). [8] Castello, A. et al. Insights into rna biology from an atlas of mammalian mrna-binding proteins. Cell 149, 1393–1406 (2012). [9] Baltz, A. G. et al. The mrna-bound proteome and its global occupancy profile on protein- coding transcripts. Molecular cell 46, 674–690 (2012).

[0010] Alles, J. et al. An estimate of the total number of true human mirnas. Nucleic acids research 47, 3353–3364 (2019).

[0011] Mukherjee, N. et al. Deciphering human ribonucleoprotein regulatory networks. Nucleic acids research 47, 570–581 (2019).

[0012] Gerstberger, S., Hafner, M. & Tuschl, T. A census of human rna-binding proteins. Nature Reviews Genetics 15, 829–845 (2014).

[0013] Agarwal, V., Bell, G. W., Nam, J.-W. & Bartel, D. P. Predicting effective microrna target sites in mammalian mrnas. elife 4, e05005 (2015).

[0014] McGeary, S. E. et al. The biochemical basis of microrna targeting efficacy. Science 366, eaav1741 (2019).

[0015] Hogan, D. J., Riordan, D. P., Gerber, A. P., Herschlag, D. & Brown, P. O. Diverse rna- binding proteins interact with functionally related sets of rnas, suggesting an extensive regulatory system. PLoS biology 6, e255 (2008).

[0016] Barreau, C., Paillard, L. & Osborne, H. B. Au-rich elements and associated factors: are there unifying principles? Nucleic acids research 33, 7138–7150 (2005).

[0017] Wen, M., Cong, P., Zhang, Z., Lu, H. & Li, T. Deepmirtar: a deep-learning approach for predicting human mirna targets. Bioinformatics 34, 3781–3787 (2018).

[0018] Gu, T., Zhao, X., Barbazuk, W. B. & Lee, J.-H. mitar: a hybrid deep learning-based approach for predicting mirna targets. BMC bioinformatics 22, 1–16 (2021).

[0019] Klimentova´, E. et al. mirbind: A deep learning method for mirna binding classification. Genes 13, 2323 (2022).

[0020] Talukder, A., Zhang, W., Li, X. & Hu, H. A deep learning method for mirna / isomir target detection. Scientific Reports 12, 10618 (2022).

[0021] Sood, P., Krek, A., Zavolan, M., Macino, G. & Rajewsky, N. Cell-type1178 specific signatures of micrornas on target mrna expression. Proceedings of the National Academy of Sciences 103, 2746–2751 (2006).

[0022] Nam, J.-W. et al. Global analyses of the effect of different cellular contexts on microrna targeting. Molecular cell 53, 1031–1043 (2014).

[0023] Nowakowski, T. J. et al. Regulation of cell-type-specific transcriptomes by microrna networks during human brain development. Nature neuroscience 21, 1784–1792 (2018).

[0024] Addo-Quaye, C., Miller, W. & Axtell, M. J. Cleaveland: a pipeline for using degradome data to find cleaved small rna targets. Bioinformatics 25, 130–131 (2009).

[0025] Folkes, L. et al. Paresnip: a tool for rapid genome-wide discovery of small rna / target interactions evidenced through degradome sequencing. Nucleic acids research 40, e103– e103 (2012).

[0026] Thody, J. et al. Paresnip2: a tool for high-throughput prediction of small rna targets from degradome sequencing data using configurable targeting rules. Nucleic acids research 46, 8730–8739 (2018).

[0027] Bracken, C. P. et al. Global analysis of the mammalian rna degradome reveals widespread mirna-dependent and mirna-independent endonucleolytic cleavage. Nucleic acids research 39, 5658–5668 (2011).

[0028] Zhang, Y. et al. Identifying cleaved and noncleaved targets of small interfering rnas and micrornas in mammalian cells by spyclip. Molecular Therapy-Nucleic Acids 22, 900–909 (2020).

[0029] Won, J.-I., Shin, J., Park, S. Y., Yoon, J. & Jeong, D.-H. Global analysis of the human rna degradome reveals widespread decapped and endonucleolytic cleaved transcripts. International Journal of Molecular Sciences 21, 6452 (2020).

[0030] Zhou, F. et al. Identification of micrornas and their endonucleolytic cleavaged target mrnas in colorectal cancer. BMC cancer 20, 1–15 (2020).

[0031] Schmidt, S. A. et al. Identification of smg6 cleavage sites and a preferred rna cleavage motif by global analysis of endogenous nmd targets in human cells. Nucleic acids research 43, 309–323 (2015).

[0032] Helwak, A., Kudla, G., Dudnakova, T. & Tollervey, D. Mapping the human mirna interactome by clash reveals frequent noncanonical binding. Cell 153, 654–665 (2013).

[0033] Hejret, V. et al. Analysis of chimeric reads characterises the diverse targetome of ago2- mediated regulation. Scientific Reports 13, 22895 (2023).

[0034] Manakov, S. A. et al. Scalable and deep profiling of mrna targets for individual micrornas with chimeric eclip. BioRxiv 2022–02 (2022).

[0035] Karginov, F. V. et al. Diverse endonucleolytic cleavage sites in the mammalian transcriptome depend upon micrornas, drosha, and additional nucleases. Molecular cell 38, 781–788 (2010).

[0036] German, M. A. et al. Global identification of microrna–target rna pairs by parallel analysis of rna ends. Nature biotechnology 26, 941–946 (2008).

[0037] Addo-Quaye, C., Eshoo, T. W., Bartel, D. P. & Axtell, M. J. Endogenous sirna and mirna targets identified by sequencing of the arabidopsis degradome. Current Biology 18, 758–762 (2008).

[0038] Zhang, H. et al. Linearcofold and linearcopartition: linear-time algorithms for secondary structure prediction of interacting rna molecules. Nucleic Acids Research 51, e94–e94 (2023).

[0039] Van Nostrand, E. L. et al. Principles of rna processing from analysis of enhanced clip maps for 150 rna binding proteins. Genome biology 21, 1–26 (2020).

[0040] Haberman, N. et al. Abundant capped rnas are derived from mrna cleavage at 3’utr g- quadruplexes. bioRxiv (2023). URL https: / / www.biorxiv.org / content / early / 2023 / 06 / 22 / 2023.04.27.538568.

[0041] Yekta, S., Shih, I.-h. & Bartel, D. P. Microrna-directed cleavage of hoxb8 mrna. Science 304, 594–596 (2004).

[0042] Liu, Z. et al. A convnet for the 2020s, 11976–11986 (2022).

[0043] Huang, H.-Y. et al. mirtarbase update 2022: an informative resource for experimentally validated mirna–target interactions. Nucleic acids research 50, D222–D230 (2022).

[0044] Karagkouni, D. et al. Diana-tarbase v8: a decade-long collection of experimentally supported mirna–gene interactions. Nucleic acids research 46, D239–D245 (2018).

[0045] Kaplan, J. et al. Scaling laws for neural language models. arXiv preprint arXiv:2001.08361 (2020).

[0046] Hesslow, D., Zanichelli, N., Notin, P., Poli, I. & Marks, D. Rita: a study on scaling up generative protein sequence models. arXiv preprint arXiv:2205.05789 (2022).

[0047] Vaswani, A. et al. Attention is all you need. Advances in neural information processing systems 30 (2017).

[0048] Gu, A. & Dao, T. Mamba: Linear-time sequence modeling with selective state spaces. arXiv preprint arXiv:2312.00752 (2023).

[0049] Pla, A., Zhong, X. & Rayner, S. miraw: A deep learning-based approach to predict microrna targets by analyzing whole microrna transcripts. PLoS computational biology 14, e1006185 (2018).

[0050] Li, X. et al. High-resolution in vivo identification of mirna targets by halo enhanced ago2 pull-down. Molecular cell 79, 167–179 (2020).

[0051] Vainberg Slutskin, I., Weingarten-Gabbay, S., Nir, R., Weinberger, A. & Segal, E. Unraveling the determinants of microrna mediated regulation using a massively parallel reporter assay. Nature Communications 9, 529 (2018).

[0052] McGeary, S. E., Bisaria, N., Pham, T. M., Wang, P. Y. & Bartel, D. P. Microrna 3- compensatory pairing occurs through two binding modes, with affinity shaped by nucleotide identity and position. Elife 11, e69803 (2022).

[0053] Xu, H. et al. Nfix circular rna promotes glioma progression by regulating mir-34a-5p via notch signaling pathway. Frontiers in molecular neuroscience 11, 225 (2018).

[0054] Bakheet, T., Hitti, E. & Khabar, K. S. A. Ared-plus: an updated and expanded database of au-rich element-containing mrnas and pre-mrnas. Nucleic acids research 46, D218–D220 (2018).

[0055] Friedman, R. C., Farh, K. K.-H., Burge, C. B. & Bartel, D. P. Most mammalian mrnas are conserved targets of micrornas. Genome research 19, 92–105 (2009).

[0056] Wang, X. Composition of seed sequence is a major determinant of microrna targeting patterns. Bioinformatics 30, 1377–1383 (2014).

[0057] Grimson, A. et al. Microrna targeting specificity in mammals: determinants beyond seed pairing. Molecular cell 27, 91–105 (2007).

[0058] Sætrom, P. et al. Distance constraints between microrna target sites dictate efficacy and cooperativity. Nucleic acids research 35, 2333–2342 (2007).

[0059] Li, P. et al. Differential inhibition of target gene expression by human micrornas. Cells 8, 791 (2019).

[0060] Awan, H. M. et al. Comparing two approaches of mir-34a target identification, biotinylated-mirna pulldown vs mirna overexpression. RNA biology 15, 55–61 (2018).

[0061] Luna, J. M. et al. Argonaute clip defines a deregulated mir-122-bound transcriptome that correlates with patient survival in human liver cancer. Molecular cell 67, 400–410 (2017).

[0062] Fujiwara, Y. et al. mir-23a / b clusters are not essential for the pathogenesis of osteoarthritis in mouse aging and post-traumatic models. Frontiers in Cell and Developmental Biology 10, 1043259 (2023).

[0063] Siegel, D. A., Le Tonqueze, O., Biton, A., Zaitlen, N. & Erle, D. J. Massively parallel analysis of human 3 utrs reveals that au-rich element length and registration predict mrna destabilization. G312, jkab404 (2022). -Mun˜oz, M. D., Bell, S. E. & Turner, M. Deletion of au-rich elements within the bcl2 3 utr reduces protein expression and b cell survival in vivo.PloS one 10, e0116899 (2015).

[0065] Karczewski, K. J. et al. The mutational constraint spectrum quantified from variation in 141,456 humans. Nature 581, 434–443 (2020).

[0066] Verdura, E. et al. Disruption of a mi r-29 binding site leading to col4a1 upregulation causes pontine autosomal dominant microangiopathy with leukoencephalopathy. Annals of neurology 80, 741–753 (2016).

[0067] Tan, Z. et al. Allele-specific targeting of micrornas to hla-g and risk of asthma. The American Journal of Human Genetics 81, 829–834 (2007).

[0068] Landrum, M. J. et al. Clinvar: public archive of relationships among sequence variation and human phenotype. Nucleic acids research 42, D980–D985 (2014).

[0069] Landrum, M. J. et al. Clinvar: improving access to variant interpretations and supporting evidence. Nucleic acids research 46, D1062–D1067 (2018).

[0070] Chin, L. J. et al. A snp in a let-7 microrna complementary site in the kras 3 untranslated region increases non–small cell lung cancer risk. Cancer research 68, 8535–8540 (2008).

[0071] Li, Z.-H. et al. A let-7 binding site polymorphism rs712 in the kras 3 utr is associated with an increased risk of gastric cancer. Tumor Biology 34, 3159–3163 (2013).

[0072] Hu, H., Zhang, L., Teng, G., Wu, Y. & Chen, Y. A variant in 3-untranslated region of kras compromises its interaction with hsa-let-7g and contributes to the development of lung cancer in patients with copd. International Journal of Chronic Obstructive Pulmonary Disease 1641–1649 (2015).

[0073] Sengupta, K. et al. Genome editing-mediated utrophin upregulation in duchenne muscular dystrophy stem cells. Molecular Therapy-Nucleic Acids 22, 500–509 (2020).

[0074] Song, D., Zhang, Q., Zhang, H., Zhan, L. & Sun, X. Mir-130b-3p promotes colorectal cancer progression by targeting chd9. Cell Cycle 21, 585–601 (2022).

[0075] Idrees, M. et al. Decreased serum pon1 arylesterase activity in familial hypercholesterolemia patients with a mutated ldlr gene. Genetics and Molecular Biology 41, 570–577 (2018).

[0076] Zhao, Y. et al. Association between pon1 activity and coronary heart disease risk: a meta-analysis based on 43 studies. Molecular genetics and metabolism 105, 141–148 (2012).

[0077] Beaudet, L. et al. Alphalisa immunoassays: the no-wash alternative to elisas for research and drug discovery (2008).

[0078] Bielefeld-Sevigny, M. Alphalisa immunoassay platform—the “no-wash” high-throughput alternative to elisa. Assay and drug development technologies 7, 90–92 (2009).

[0079] Agarwal, V. & Kelley, D. R. The genetic and biochemical determinants of mrna degradation rates in mammals. Genome Biology 23, 245 (2022).

[0080] Moore, M. J. et al. mirna–target chimeras reveal mirna 3-end pairing as a major determinant of argonaute target specificity. Nature communications 6, 8864 (2015).

[0081] Zhou, S. et al. Degradome sequencing reveals an integrative mirna-mediated gene interaction network regulating rice seed vigor. BMC plant biology 22, 269 (2022).

[0082] Chen, S., Zhou, Y., Chen, Y. & Gu, J. fastp: an ultra-fast all-in-one fastq preprocessor. Bioinformatics 34, i884–i890 (2018).

[0083] Langmead, B., Trapnell, C., Pop, M. & Salzberg, S. L. Ultrafast and memory-efficient alignment of short dna sequences to the human genome. Genome biology 10, 1–10 (2009).

[0084] Corcoran, D. L. et al. Paralyzer: definition of rna binding sites from par-clip short-read sequence data. Genome biology 12, 1–16 (2011).

[0085] Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with bowtie 2. Nature methods 9, 357–359 (2012).

[0086] Lovci, M. T. et al. Rbfox proteins regulate alternative mrna splicing through evolutionarily conserved rna bridges. Nature structural & molecular biology 20, 1434–1442 (2013).

[0087] Smith, T., Heger, A. & Sudbery, I. Umi-tools: modeling sequencing errors in unique molecular identifiers to improve quantification accuracy. Genome research 27, 491–499 (2017).

[0088] Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet. journal 17, 10–12 (2011).

[0089] Kozomara, A., Birgaoanu, M. & Griffiths-Jones, S. mirbase: from microrna sequences to function. Nucleic acids research 47, D155–D162 (2019).

[0090] Dobin, A. et al. Star: ultrafast universal rna-seq aligner. Bioinformatics 29, 15–21 (2013).

[0091] Kim, D., Paggi, J. M., Park, C., Bennett, C. & Salzberg, S. L. Graph-based genome alignment and genotyping with hisat2 and hisat-genotype. Nature biotechnology 37, 907– 915 (2019).

[0092] Liao, Y., Smyth, G. K. & Shi, W. featurecounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30, 923–930 (2014).

[0093] Rodriguez, J. M. et al. Appris: selecting functionally important isoforms. Nucleic acids research 50, D54–D59 (2022).

[0094] Cao, J., Zhou,W., Steemers, F., Trapnell, C. & Shendure, J. Sci-fate characterizes the dynamics of gene expression in single cells. Nature biotechnology 38, 980–988 (2020).

[0095] Schofield, J. A., Duffy, E. E., Kiefer, L., Sullivan, M. C. & Simon, M. D. Timelapse-seq: adding a temporal dimension to rna sequencing through nucleoside recoding. Nature methods 15, 221–225 (2018).

[0096] Bohn, E., Lau, T. T.,Wagih, O.,Masud, T. &Merico, D. A curated census of pathogenic and likely pathogenic utr variants and evaluation of deep learning models for variant effect prediction. Frontiers in Molecular Biosciences 10 (2023).

Claims

CLAIMS:

1. A computer-implemented method for predicting tissue-specific microRNA (miRNA) binding to a messenger RNA (mRNA), the method comprising: - providing, in a memory, a post-transcriptional gene regulation (PTGR) model comprising a base and at least one head; - receiving, at a processor in communication with the memory, an RNA input sequence of length L corresponding to the mRNA; and - determining, at the processor, an miRNA binding prediction matrix from the at least one head of the PTGR model, the miRNA binding prediction matrix comprising an L x K matrix comprising miRNA binding predictions at each position of the L nucleotides in the RNA input sequence for K tissue types, the miRNA binding prediction matrix determined using the RNA input sequence as input at the base of the PTGR model.

2. The method of claim 1, wherein the method further comprises predicting tissue- specific degradation of the mRNA and the at least one head comprises at least two heads: - a first head providing a first miRNA binding prediction matrix binding for K tissue types, optionally K human tissue types; and - a second head providing a first mRNA degradome prediction matrix for K tissue types, optionally K human tissue types, the first mRNA degradome prediction matrix determined using the RNA input sequence as input at the base of the PTGR model.

3. The method of claim 2, wherein the first mRNA degradome prediction matrix comprises an L x K matrix comprising degradome read coverage predictions at each position of the L nucleotides in the RNA input sequence for K tissue types.

4. The method of any one of claims 1 to 3, wherein the at least two heads comprise at least four heads:- a third head providing a second miRNA binding prediction matrix for tissue types for a second organism, optionally mouse; and - a fourth head providing a second mRNA degradome prediction matrix for tissue types for the second organism, optionally mouse.

5. The method of any one of claims 1 to 4, wherein the PTGR model comprises a convolutional neural network.

6. The method of any one of claims 1 to 5, wherein the base comprises: - at least one base convolution layer; and - at least one base residual block receiving an output from the at least one base convolution layer, optionally wherein the base residual block comprises gated convolutions.

7. The method of claim 6, wherein each base residual block of the at least one base residual block comprises: - a normalization layer receiving an input to the base residual block; - a first pair of convolutional layers receiving the output from the normalization layer, a first of the first pair having a tanh activation and a second of the first pair having a sigmoid activation; - a first element-wise multiplication operation on outputs of the first pair of convolutional layers; - a second pair of convolutional layers receiving the output from the first element-wise multiplication operation, a first of the second pair having a tanh activation and a second of the second pair having a sigmoid activation; - a second element-wise multiplication operation on outputs of the second pair of convolutional layers; and- an output of the base residual block determined by adding an output of the second element-wise multiplication operation and the input to the base residual block in a skip connection.

8. The method of any one of claims 1 to 7, wherein each head of the at least one head comprises a first head convolutional layer receiving input from the base, at least one head residual block receiving an output of the first head convolutional layer, and a second head convolutional layer receiving an output of the at least one head residual block.

9. The method of claim 8, wherein each of the at least one head residual block comprises: - a normalization layer receiving an input to the head residual block; - a first convolutional layer receiving the output from the normalization layer, the first convolutional layer having a gelu activation; - a second convolutional layer receiving the output from the first convolutional layer, the second convolutional layer having a gelu activation; - an output of the head residual block determined by adding an output of the second convolutional layer and the input to the head residual block in a skip connection.

10. The method of any one of claims 1 to 9, further comprising: - outputting, at a display device in communication with the processor, at least one of the miRNA binding prediction matrix and the mRNA degradome prediction matrix.

11. The method of any one of claims 1 to 10, further comprising: - receiving, at an input device in communication with the processor, one or more user selected nucleotides in the RNA input sequence; - generating a first visualization comprising the miRNA binding predictions at the one or more user selected nucleotides in the RNA input sequence; and- outputting, at the display device in communication with the processor, the generated first visualization.

12. The method of claim 11, further comprising: - generating a second visualization comprising the mRNA degradome prediction matrix at the one or more user selected nucleotides in the RNA input sequence; and - outputting, at the display device, the generated second visualization.

13. The method of any one of claims 1 to 12, further comprising: - converting, at the processor, the RNA input sequence of length L to an L x 4 input matrix, wherein the L x 4 input matrix is one-hot encoded.

14. The method of any one of claims 1 to 13, wherein the miRNA binding predictions are indicative of non-specific miRNA binding to the mRNA.

15. The method of any one of claims 1 to 14, wherein the K tissue types comprise at least one selected from the group of liver, CNS, heart, kidney, muscle, cancer cells, HEK cells, HeLa cells, iPSCs, primary human hepatocytes (PHHs) and peripheral blood mononuclear cells (PBMCs).

16. The method of any one of claims 1 to 15 wherein the PTGR model comprises an ensemble model.

17. The method of claim 16 wherein the ensemble model comprises at least four models, each of the four models trained based on a training data split across different sets of chromosomes.

18. A system for predicting tissue-specific microRNA (miRNA) binding to a messenger RNA (mRNA), the system comprising: - a memory comprising: - a post-transcriptional gene regulation (PTGR) model comprising a base and at least one head;- a processor in communication with the memory, the processor configured to: - receive an RNA input sequence of length L corresponding to the mRNA; and - determine an miRNA binding prediction matrix from the at least one head of the PTGR model, the miRNA binding prediction matrix comprising an L x K matrix comprising binding predictions at each position of the L nucleotides in the RNA input sequence for K tissue types, the miRNA binding prediction matrix determined using the RNA input sequence as input at the base of the PTGR model.

19. The system of claim 18, wherein the processor is further configured to predict tissue- specific degradation of the mRNA and the at least one head comprises at least two heads by: - a first head providing a first miRNA binding prediction matrix binding for K tissue types, optionally K human tissue types; and - a second head providing a first mRNA degradome prediction matrix for K tissue types, optionally K human tissue types, the first mRNA degradome prediction matrix determined using the RNA input sequence as input at the base of the PTGR model.

20. The system of claim 19, wherein the first mRNA degradome prediction matrix comprises an L x K matrix comprising degradome read coverage predictions at each position of the L nucleotides in the RNA input sequence for K tissue types.

21. The system of any one of claims 18 to 20, wherein the at least two heads comprise at least four heads: - a third head providing a second miRNA binding prediction matrix for tissue types for a second organism, optionally mouse; and - a fourth head providing a second mRNA degradome prediction matrix for tissue types for the second organism, optionally mouse.

22. The system of any one of claims 18 to 21, wherein the PTGR model comprises a convolutional neural network.

23. The system of any one of claims 18 to 22, wherein the base comprises: - at least one base convolution layer; and - at least one base residual block receiving an output from the at least one base convolution layer, optionally wherein the base residual block comprises gated convolutions.

24. The system of claim 23, wherein each base residual block comprises: - a normalization layer receiving an input to the base residual block; - a first pair of convolutional layers receiving the output from the normalization layer, a first of the first pair having a tanh activation and a second of the first pair having a sigmoid activation; - a first element-wise multiplication operation on outputs of the first pair of convolutional layers; - a second pair of convolutional layers receiving the output from the first element-wise multiplication operation, a first of the second pair having a tanh activation and a second of the second pair having a sigmoid activation; - a second element-wise multiplication operation on outputs of the second pair of convolutional layers; and - an output of the base residual block determined by adding an output of the second element-wise multiplication operation and the input to the base residual block in a skip connection.

25. The system of any one of claims 18 to 24, wherein each head of the at least one head comprises a first head convolutional layer receiving input from the base, at least one head residual block receiving an output of the first head convolutional layer, and a second head convolutional layer receiving an output of the at least one head residual block.

26. The system of claim 25 wherein each of the at least one head residual block comprises: - a normalization layer receiving an input to the head residual block; - a first convolutional layer receiving the output from the normalization layer, the first convolutional layer having a gelu activation; - a second convolutional layer receiving the output from the first convolutional layer, the second convolutional layer having a gelu activation; - an output of the head residual block determined by adding an output of the second convolutional layer and the input to the head residual block in a skip connection.

27. The system of any one of claims 18 to 26, further comprising: - a display device in communication with the processor, and - wherein the processor is further configured to: - output to the display device at least one of the miRNA binding prediction matrix and the mRNA degradome prediction matrix.

28. The system of claim 27, further comprising: - an input device in communication with the processor; and - wherein the processor is further configured to: - receive from the input device one or more user selected nucleotides in the RNA input sequence; - generate a first visualization comprising the miRNA binding predictions at the one or more user selected nucleotides in the RNA input sequence; and - output to the display device in communication with the processor, the generated first visualization.

29. The system of claim 28, wherein the processor is further configured to:- generate a second visualization comprising the mRNA degradome prediction matrix at the one or more user selected nucleotides in the RNA input sequence; and - output to the display device, the generated second visualization.

30. The system of any one of claims 18 to 29, wherein the processor is further configured to: - convert the RNA input sequence of length L to an L x 4 input matrix, wherein the L x 4 input matrix is one-hot encoded.

31. The system of any one of claims 18 to 29, wherein the miRNA binding predictions are indicative of non-specific miRNA binding to the mRNA.

32. The system of any one of claims 18 to 30, wherein the K tissue types comprise at least one of liver, CNS, heart, kidney, muscle, cancer cells, HEK cells, HeLa cells, iPSCs, primary human hepatocytes (PHHs) and peripheral blood mononuclear cells (PBMCs).

33. A method for determining an effect of a variant mRNA sequence relative to a control mRNA sequence on post transcriptional gene regulation, the method comprising: - predicting tissue specific miRNA binding to the variant mRNA sequence, and optionally tissue specific degradation of the variant mRNA sequence, according to the method of any one of claims 1 to 16; - comparing the miRNA binding prediction matrix, and optionally the mRNA degradome prediction matrix, for the variant mRNA sequence to a miRNA binding prediction matrix, and optionally a mRNA degradome prediction matrix, for the control mRNA sequence; and - determining the effect of the variant mRNA sequence on post transcriptional gene regulation based on any differences between the miRNA binding prediction matrix, and optionally the mRNA degradome prediction matrix, for the variant mRNA sequence, relative to the miRNA binding prediction matrix, and optionally the mRNA degradome prediction matrix, for the control mRNA sequence.

34. The method of claim 33, further comprising determining the miRNA binding prediction matrix, and optionally the mRNA degradome prediction matrix, for the control mRNA sequence according to the method of any one of claims 1 to 17.

35. The method of claim 33 or 34, wherein the variant mRNA sequence has between 1 and 20 single nucleotide polymorphisms relative to the control mRNA sequence, optionally wherein the control mRNA sequence is a wild-type sequence.

36. The method of claim 33 to 35, wherein the variant mRNA sequence has a variant of unknown clinical significance, a putative disease-causing mutation, a masked sequence corresponding to a SBO binding site, an ADAR editing site or a synthetic mRNA sequence.

37. The method of any one of claims 33 to 36, further comprising synthesizing the variant mRNA and testing the variant mRNA for gene expression.

38. A system for determining an effect of a variant mRNA sequence relative to a control mRNA sequence on post transcriptional gene regulation, the system comprising: - a memory comprising a PTGR model; - a processor in communication with the memory, the processor configured to: - predict tissue specific miRNA binding to the variant mRNA sequence using the PTGR model, and optionally tissue specific degradation of the variant mRNA sequence, according to the method of any one of claims 1 to 14; - compare the miRNA binding prediction matrix, and optionally the mRNA degradome prediction matrix, for the variant mRNA sequence to a miRNA binding prediction matrix, and optionally a mRNA degradome prediction matrix, for the control mRNA sequence; and - determine the effect of the variant mRNA sequence on post transcriptional gene regulation based on any differences between the miRNA binding prediction matrix, and optionally the mRNA degradome prediction matrix, for the variant mRNAsequence, relative to the miRNA binding prediction matrix, and optionally the mRNA degradome prediction matrix, for the control mRNA sequence.

39. The system of claim 38, wherein the processor is further configured to determine the miRNA binding prediction matrix, and optionally the mRNA degradome prediction matrix, for the control mRNA sequence according to the method of any one of claims 1 to 17.

40. The system of claim 38 or 39, wherein the variant mRNA sequence has between 1 and 20 single nucleotide polymorphisms relative to the control mRNA sequence, optionally wherein the control mRNA sequence is a wild-type mRNA sequence.

41. The system of claim 38 or 39, wherein the variant mRNA sequence comprises one or more of a variant of unknown clinical significance, a putative disease-causing mutation, a masked sequence corresponding to a SBO binding site, an ADAR editing site or a non- naturally occurring mRNA sequence.

42. A computer-implemented method for generating a post-transcriptional gene regulation (PTGR) model, the method comprising: - providing in a memory, a machine learning model comprising a base and at least two heads; - providing, in the memory, a first data set comprising miRNA binding data for a plurality of mRNA sequences for a K plurality of tissues; - training the machine learning model based on the first data set, wherein a prediction output of the machine learning model is a matrix of size L × K, K denoting a size of the K plurality of tissues, a k-th row of the matrix corresponding to a predicted miRNA binding probability each base pair for the k-th cell line.

43. The method of claim 42, wherein the miRNA binding data comprises mRNA sequence segments experimentally associated with miRNA binding.

44. The method of any one of claims 42 or 43, further comprising:- providing, in the memory, a second data set comprising mRNA degradome data for the plurality of mRNA sequences for the K plurality of tissues, - wherein the training the machine learning model further comprises training the machine learning model based on the first data set and the second data set.

45. The method of any one of claims 42 to 44, wherein each of the mRNA sequences comprises a degradome annotation, and the degradome annotation comprises a degradome matrix corresponding to degradome-seq data for each of the K plurality of tissue types.

46. The method of any one of claims 44 or 45, wherein the mRNA degradome data further comprises a read coverage value at each position of the mRNA sequence.

47. The method of any one of claims 42 to 46, wherein the base and the at least two heads comprise at least one residual block.

48. The method of any one of claims 42 to 47, wherein the at least two heads comprise a first head for providing a miRNA binding prediction matrix and a second head for providing a mRNA degradome prediction matrix and training the machine learning model comprises training the base and the first head on the first data set and training the base and the second head on the second data set.

49. The method of any one of claims 42 to 48, wherein the plurality of mRNA sequences comprises a plurality of non-overlapping windows.

50. The method of claim 49 wherein the plurality of non-overlapping windows are at least 3000 nucleotides long.

51. The method of any one of claims 42 to 50, further comprising: - converting each of the plurality of mRNA sequences to a one-hot encoded sequence.

52. The method of any one of claims 42 to 51 wherein a loss function for the training themachine learning model comprises+ , wherein a binary cross-entropy losscomprises a miRNA binding output and a Poisson loss comprises a degradation output.

53. The method of any one of claims 42 to 52 wherein the machine learning model comprises an ensemble model.

54. The method of claim 53 wherein the ensemble model comprises at least four models, each of the four models trained based on a training data split across different sets of chromosomes.

55. A system for generating a post-transcriptional gene regulation (PTGR) model, the system comprising a memory and a processor configured to perform any one of claims 42 to 54.