{'fraction_non_alphanumeric': 0.026599860000736837, 'fraction_numerical': 0.010978889584791658, 'mean_word_length': 5.003981420039814, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Functional Representation of Large-Scale Heterogeneous RNA Sequences with Integration of Diverse Multi-omics, Interactions, and Annotations Data

Long non-coding RNA, microRNA, and messenger RNA enable key regulations of various biological processes through a variety of diverse interaction mechanisms. Identifying the interactions and cross-talk between these heterogeneous RNA classes is essential in order to uncover the functional role of individual RNA transcripts, especially for unannotated and newly-discovered RNA sequences with no known interactions. Recently, sequence-based deep learning and network embedding methods are becoming promising approaches that can either predict RNA-RNA interactions from a sequence or infer missing interactions from patterns that may exist in the network topology. However, the majority of these methods have several limitations, eg, the inability to perform inductive predictions, to distinguish the directionality of interactions, or to integrate various sequence, interaction, and annotation biological datasets. We proposed a novel deep learning-based framework, rna2rna, which learns from RNA sequences to produce a low-dimensional embedding that preserves the proximities in both the interactions topology and the functional affinity topology. In this proposed embedding space, we have designated a two-part"source and target contexts"to capture the targeting and receptive fields of each RNA transcript, while encapsulating the heterogenous cross-talk interactions between lncRNAs and miRNAs. From experimental results, our method exhibits superior performance in AUPR rates compared to state-of-art approaches at predicting missing interactions in different RNA-RNA interaction databases and was shown to accurately perform link predictions to novel RNA sequences not seen at training time, even without any prior information. Additional results suggest that our proposed framework can capture a manifold for heterogeneous RNA sequences to discover novel functional annotations.


Introduction
Regulatory long non-coding RNAs (lncRNAs) and microRNAs (miRNAs) that influences gene expression post-transcriptionally by interacting to target messenger RNAs (mRNA) form a complex network of transcriptomic interactions. These heterogeneous families of noncoding RNAs are associated with nearly all cellular processes, and recently, lncRNAs are gaining considerable attention as the largest and most diverse class of non-coding RNA, encompassing nearly 30,000 discovered transcripts in human. Among many of their known functional interaction mechanisms, lncRNAs are known to act as miRNA decoys, derepress gene expression by competing with miRNAs for shared mRNA targets, or directly regulate gene expression. 1 Determining the biological functions of the individual lncRNAs remains a challenge as most of these RNA transcripts are currently unannotated, and their known interactions are sparse. Recent advances in RNA sequencing (RNA-Seq), deep sequencing (CLIP-seq, LIGR-Seq), and computational methods allow for an unprecedented analysis of such transcripts and have enabled researchers to generate large-scale interaction and functional annotation databases. However, the interaction networks generated from such data are often scant and incomplete in the number of lncRNAs covered. Furthermore, although a large number of lncRNAs have been identified, only a few hundreds have had functional and molecular mechanisms determined to date, as observed in as lncRNAdb. 2 In other avenues, a growing number of lncRNAs are being assigned biological functions based on their cell-specific expressions 3 or disease associations. 4 Still, a vast majority of other lncRNAs only has basic genomic information such as locus biotype and transcript sequence assigned. These transcripts might support important biological cell functions and could potentially serve as targets for genomic, diagnostic, or therapeutic studies. Thus, in the effort to functionally characterize these "hypothetical" lncRNAs, the essential tasks are integrating the various -omics and annotation attributes between the known RNAs, while inferring the functions and interactions to the novel RNAs directly from sequence.
To address these challenges, this paper a machine learning methodology called "network embedding," which learn from the interaction topology and attributes of transcriptome-wide RNAs to accurately infer a functional representation that can be used to predict RNA-RNA functional interactions. Many graph-theoretic methods have been applied to biological networks with the intuition that RNAs close together in the interaction topology are more likely to be involved in many of the same functions. 5 Due to the extreme sparsity of the known interaction network among lncRNAs, it is pertinent to unravel the functional association between lncRNAs by considering its gene/transcript genomic annotations, tissue-specific expressions, gene-disease associations, and sequence similarity. 4 In recent studies, Kishan et al. 6 uncovered the secondorder proximity relationship between interacting genes by integration of the gene regulatory network and gene expression as side information. On the other hand, several structure-free sequence-based methods 7 have also been proposed for prediction of protein binding sites, family classification, structure prediction, or interaction prediction from sequence. Motivated by these works, we propose an integrative method to extract a functional representation from RNA sequences in order to accurately predict new interactions and functional annotations, while simultaneously integrates various existing biological multi-omics and annotations data.

Defining the Heterogeneous lncRNA-miRNA-mRNA Network
We formally define the heterogeneous network of lncRNA, miRNA, and mRNA interactions and functional similarity as two networks of directed and undirected edges. We denote the two networks, G 1 (V, E d ) and G 2 (V, E u ) having the same set of nodes V and two set of edges E d and E u . The set of nodes V = {v 1 , . . . , v n } can also be expressed as V = {L, M, N } s.t. |L| + |M | + |N | = n, where L, M, N are the sets containing the lncRNA, microRNA, and mRNA heterogeneous nodes, respectively. The set of directed edges E d = {e d ij } n i,j=1 represent directed regulatory interactions that each specify a source and a target. The undirected edges E u = {e u ij } n i,j=1 s.t. e u ji = e u ji represents the undirected functional affinity associated with the heterogeneous RNA nodes. Each edge e ij is associated with a weight such that 0 ≤ e ij ≤ 1, indicating the strength of the connection between RNA i and RNA j. If e ij > 0, we consider the edge a positive interaction/affinity, and if e ij = 0, we consider the edge a negative (non)interaction/affinity. In this paper, we consider weights e ij to be binary, indicating whether RNA i and RNA j has an interaction or affinity.

Materials and Data
The directed edges represent the directed regulatory interactions between lncRNAs, miRNAs, and mRNAs. We consider the directionality in the regulatory interactions by encoding directed edges, i.e., encoding miRNA-lncRNA interactions (e.g., miRNA inducing lncRNA decay) to be separate from lncRNA-miRNA interactions (e.g., lncRNAs acting as miRNA decoys). In this study, the different types of interaction collected from various experimentally-verified interaction databases in the lncRNA-miRNA-mRNA interactome considered are: which contained a total of 65,655 interactions between 1037 lncRNAs and 28,866 mRNAs. • mRNA-mRNA gene regulatory interactions obtained from the BioGRID v3.5 database, 12 which included more than 347,246 matched interactions among 19,607 mRNAs.
These heterogeneous interactions are combined into an integrated network, and the associated set of edges is E d , where the binary edge weight e d ij ∈ {0, 1} indicates whether a regulatory interaction from RNA node v i to RNA node v j has been observed in the literature.
We also outline the various multi-omics, sequence, and annotation databases utilized to provide functional attributes to individual RNAs. To maximize the number of transcripts and genes matched between the different databases, the lncRNA, miRNA, and mRNA transcripts are indexed by standard gene symbols.
• Tissue-specific median expressions obtained from the GTEx Portal, 13 which provides tissue-specific RNA-Seq transcript and gene expression for 15,598 non-diseased samples of various human tissue types. From the GTEx Analysis v7 dataset, the median TPM values by tissue were calculated for a total of 18653 mRNA genes, 1181 miRNA and 12706 lncRNA transcripts over 44 tissue types. • Genomic annotations obtained from the GENCODE Release 29 14 which contains the transcript biotype annotation and genomic location for lncRNAs, miRNAs, and mRNAs.
• Transcript sequences comprised of genome-wide human lncRNA and mRNA primary sequences from GENCODE Release 29 14 and miRNA hairpin sequences from miRBase. 15 • Functional annotations for mRNA were obtained from the HUGO Gene Names database. 16 In addition, GO terms for 162 matched lncRNAs were obtained from RNAcentral 17 which aggregated data from NONCODE 18 and lncipedia. 19 For miRNAs, family classification obtained from the TargetScan Release 7.2 (March 2018), 20 RNA structure family annotation from Rfam 13.0, 21 and GO terms from RNAcentral were included. • Disease associations from DisGeNet 22 covered 7577 mRNA genes. The HMDD miRNAdisease database 23 provided disease associations for 553 miRNAs, and the LncRNADisease database 24 provided disease associations for 150 lncRNAs.
Each of the RNA functional attribute data type, excluding transcript sequences, is encoded to a data matrix, A k ∈ |v|×m , where m is the dimensionality of the annotation type k, which can either be real-valued or categorical. If the attribute k is categorically-valued (e.g. disease association), A k can be an indicator matrix where entries are binary-valued, and m is the size of all possible diseases. Considering that the functional annotation to the non-coding RNAs is sparse, each RNA node can have up to K attribute types associated with it.

Undirected RNA-RNA Functional Affinity Edges
We aim to capture the functional similarity between two RNA nodes of the same class by calculating an affinity score as a similarity measure, suggesting a resemblance in RNA function or structure. For any categorical text-valued annotation attribute (e.g., disease association, transcript biotype, RNA structure family, or GO terms) that two RNA nodes v i and v j both have been annotated, the attributes in this annotation are first transformed to binary feature representation. For a categorical annotation denoted as k, the binary 1-D feature vector associated with RNA node v i is denoted as a k i . In this vector contains m binary entries that indicate whether or not the RNA node v i has been associated with each of the m total possible attributes of this annotation. Using the Sørensen-Dice coefficient score, 25 a similarity score between two binary vectors for node v i and node v j for feature k can be obtained by . This similarity measure ranges [0, 1] and gives higher weight to the common attributes present in both RNAs than by the attributes present in only one RNA. Since most RNAs have null annotations, the Dice coefficient score is only computed between pairs of RNA nodes that have both been annotated, while other pairs remain a null value. For the real-valued attribute data, i.e., tissue-specific expressions, we calculate the absolute value of the correlation between two the real-valued expression vectors a k i and a k j . Since a large number RNAs of have associated tissue-specific expression data, the co-expression correlation measure can be calculated for most of the RNAs.
To obtain an aggregate affinity score between a pair of RNA nodes across all K similarity values, we utilized a modified version of the Gower's Similarity score. 26 For each RNA-RNA pair, Gower's similarity aggregates similarity scores across all the annotation features and perform a weighted average. Typically, a similarity score for a pair of RNA nodes that do not have any associated annotation would be considered a 0, but in this study, we remove these null pairwise similarity from consideration. Thus, Gower's similarity will only aggregate the available pairwise similarity scores from annotations that exist between both nodes to compute the average. The Gower's similarity score is computed between all pairs of RNA nodes, and the resulting pairwise affinity matrix is S, where entries S ij = Kij k s k ij /|K ij | with K ij being the set of attribute types present in both RNA nodes v i and v j . The entries S ij will then be selected as edge weights e u ij that represent the functional similarity edges between node. Since our model currently only considers unweighted binary edges, we selected undirected edges to have an affinity score close to 1.0, or higher than a chosen hard-threshold, to be considered as a positive edge. In our experiments, the hard-threshold was arbitrarily chosen where the number of positive affinity edges covers no more than 0.1% sparsity of the entire affinity matrix. We also utilize the set of undirected affinity edges with a weight close to 0, indicating a negative edge that suggests functional dissimilarity between a pair of RNA nodes.

Network Embedding with Source-Target Contexts
A network embedding is mapping from each RNA node to a low-dimensional representation, denoted as a mapping function f : v i → y i ∈ d , ∀v i ∈ V , where d is the dimensionality of the embedding such that d n. The embedding y i associated with each node v i is learned such that, in this embedding space, nodes preserve some meaningful proximities to other nodes according to the given topology in the networks G d and G u . We propose the embedding space to have two components: source context and target context. That is, each embedding vector y i = [s i , t i ] is represented as a concatenation of the "source context", s i ∈ d/2 , and "target context", t i ∈ d/2 . This embedding representation can simultaneously capture directed and undirected edges by our proposed first-order directed proximity and second-order undirected proximity. The first-order directed proximity represents the directed regulatory interaction between node i's source context and node j's target context, defined as can take on a different value than d 1 (v j , v i ). The second-order undirected proximity represents the functional affinity between node i and node j, defined as The value of these proximities is the Euclidean distance, where the embeddings are trained such that if two nodes have a positive (directed or undirected) edge, its respective embeddings will be more similar, i.e., having a smaller distance. Otherwise, if two nodes have a negative or non-interactions, their embeddings should be more dissimilar, incurring a greater distance.

Representation Learning for RNA Sequences to Reconstruct the Interactions and Functional Topology
Aside from the interaction topology data, each RNA v i also has an associated transcript sequence with one-hot vector representation denoted by x i ∈ {1, 2, 3, 4} li , where l i is the variable length of an RNA sequence. We propose the network embedding function f : x i → y i ∈ d to be a deep neural network with Siamese architecture that maps RNA sequence input x i to an embedding y i of a fixed dimension d, see Fig.1  Siamese network is an architecture where a pair of objects can each be encoded, where its resulting embeddings can determine if the two objects are similar or dissimilar. More specifically, the network learns to output embeddings for a pair of RNA sequences, guided by edge weight e u ij as the label indicating whether the pair is functionally similar or dissimilar. Additionally, for an interacting pair of RNAs, the directed edge e d ij would indicate whether RNA i interacts with RNA j using the corresponding directed proximity. In order for the output embeddings to preserve the proximities across all edges in both G d and G u network topologies, we utilize the binary cross-entropy loss function, 28 defined as, After training is complete, given two RNA sequence inputs x i and x j , the learned model can output the embeddings y i and y j , which is used to predict whether a relationship exists between them by computing the respective proximity function. We use the proximity score d 1 (v i , v j ) or d 2 (v i , v j ) to either predict the existence of an interaction or functional similarity, respectively.
is applied to yield a pairwise affinity.

Model Optimization with Batch Sampling Strategy
The network weights in f (x) are trained with Stochastic Gradient Descent (SGD) using the standard back-propagation algorithm. At each SGD iteration, a batch of RNA nodes are sampled along with its associated sets of positive and negative, directed and undirected edges, and the neural network trains on the RNA sequences associated with these edges. We implement a biased random node sampling strategy where we first randomly select a set of nodes, then train on the set of edges induced by this sub-graph. The probability of selecting a node is a function of its degree, utilizing the function proposed by Riad et al.. 29 The sampling compression function is chosen to be the square-root function, which retains the power-law degree distribution while keeping the frequency ranking of each node. When a batch of nodes S is sampled without replacement from this distribution, each node i has a set of positive edges, P i . To obtain the negative directed edges, we then sample a number of nodes in S by adopting the approach of negative sampling as proposed in, 30 where the ratio of negative edges to that of positive edges incident to each node is between 2.0 and 5.0, as a free parameter. Given S, the sampled batch of nodes, and E d S , E u S , the set of directed and undirected edges containing both positive and negative interactions incident to S, we train the loss function with batch optimization with L(S, E d where λ is the coefficient parameter to control the effect of the second-order undirected proximity.

Large-scale Data Integration of lncRNA-miRNA-mRNA Sequences, Interactions, Expressions, and Annotations
We integrated various experimentally verified interaction databases to build a large-scale lncRNA-miRNA-mRNA interaction network. Additionally, various functional annotations, expressions, sequences, and disease associations were also integrated to enable extraction of the undirected attribute affinity edges. In total, there are 12725 lncRNAs, 1870 microRNAs, and 20284 mRNAs considered in this study, comprised of a comprehensive integration of the various databases illustrated in Table 1. To accomplish the primary task of predicting novel interactions not seen at training time, we propose an experimental setup using prospective evaluation. All models were trained exclusively using the prior version of each interactions databases. Then, we validate the link prediction model by using the set of new interactions from the latest database version update. This type of evaluation, rarely done in the literature, is extremely important as it allows us to mimic a realistic scenario where the task is to discover novel RNA-RNA interactions, based on our current knowledge. The training sets are comprised from the interaction network of database versions released before 2015, while the validation sets are comprised of updates from the most recently released database versions. After integration of the training databases, self-interactions and redundant interactions edges are removed, and only interactions between RNAs with an associated transcript sequence will be considered. In the validation databases, we selected only interactions that do not overlap with interactions from the training set.
After integration of various multi-omics and annotation attributes, RNA-RNA pairwise functional affinities were computed, and a number of undirected affinity edges were then added to the undirected interactions training set. After computing the affinities A for all lncRNA, miRNA, and mRNA pairs and filtering second-order undirected affinities at a 0.8 threshold, 65864, 405, and 362362 undirected edges were added, respectively. With the negative sampling ratio set at 5.0 per positive edge, a total of 329320, 2025, and 724724 negative edges were added to the undirected edges training set respectively for lncRNAs, miRNAs, and mRNAs.

Novel Link Predictions.
Our experiments include comparative analysis across different evaluation tasks with existing state-of-the-art network embedding methods. The methods considered are node2vec, 31 LINE, 32 HOPE, 33 and SDNE. 34 In the following experiments, each method was assessed by learning a 128-d embedding representation from the training network. All other free parameters are set according to the default values in the proposed methods. We compose our training set for link prediction task by a union of all ground-truth interactions set from the miRTarBase 6.0, lncBase v2, NPInter v2.0, lncRNA2Target v1.0, and BioGRID v3.4 databases. After the models have been trained, its estimated interaction adjacency matrix is computed and evaluated on the novel interactions from miRTarBase v7.0, lncBase predicted, NPInter v3.0, lncRNA2Target v2.0, and BioGRID v3.5 databases separately. For a test to differentiate between positive interactions and random noise interactions, we also uniformly sample a number of interactions from the set of all possible pairwise interactions to consider as negative interaction using the random node sampling distribution specified in Riad et al . 29 This set is denoted as E n , and the number of negative interactions is sampled such that the ratio of negative to positive interactions is 1.0. At evaluation time, the set of ground truth validation edges E d and random noise E n edges is used to calculate the precision and recall rates. All methods were evaluated on the same set of positive and sampled negative interactions. The area under the precision-recall curve (AUPR) scores shown in Fig. 2 highlights the comparison analysis across five different interaction databases. It is important to note that rna2rna can perform inductive link prediction to novel RNA sequences since the validation set contains a number of novel lncRNA, miRNA, and mRNA not seen at training time. While other transductive methods were evaluated only on the interactions incident to the RNA nodes seen at training time, rna2rna were evaluated on a greater number of interactions, yet still achieved the top performance.

Inferring Functional Similarity From Embeddings
In comparison analysis, we first obtained the embeddings from each of the methods and performed K-Means clustering only on the nodes that have an associated functional annotation. The number of clusters in K-Means is the same as the total number of unique labels in a particular annotation. The clustering result of different methods is compared over the RNA family and RNA type annotations in Table 2 and Table 3. The result shows that although there is a greater number of RNA nodes to assign to clusters, rna2rna embeddings can achieve the highest NMI score over the RNA functional family annotations. Since rna2rna embeddings have demonstrated functionally similarity in the experiments above, an important next step is to assign putative biological functions to novel lncRNAs. To do this, we perform gene set enrichment analysis on K-mean clusters of RNAs, select the cluster with the highest enriched functional term, then associate the lncRNAs belonging in this cluster with this term. In this experiment, the embeddings are trained from both training set and validation set over 32,741 different RNAs, and K-means cluster (k = 2000) is performed. We then performed enrichment analysis on these 2000 clusters using Enrichr 35 over the KEGG Human 2019 36 terms, which includes both functional and disease pathways. Among the 2000 clusters, 559 have an adjusted P-value of less than 0.01, and 139 have an adjusted P-value of less than 0.001. Interestingly, the highest-scoring gene sets often contain some lncRNAs not previously associated with these functional terms, which are shown in Table 4. It warrants additional experimental studies to verify the functional associations of these lncRNAs.

Discussions and Conclusion
With the framework we have developed, heterogeneous functional attributes and interaction data are integrated to enable characterization of RNA sequences using an embedding representation. While the method of integrating various functional attributes is simple, its purpose is to allow for characterizing the functional affinity for an extensive number of RNAs, even among sparsely annotated ones. While very few lncRNAs have been annotated for all of its attributes, especially functional annotation or disease association, most have already been annotated with tissue-specific expressions, transcript biotype, and sequence. Thus, the deep network can be trained to recognize a greater number of RNA sequences to achieve better performance.  Additionally, since our method was able to map the functional affinity between RNA nodes belonging in disconnected components in the interaction topology, we hypothesize rna2rna can effectively map any RNA transcript sequence to a functional manifold in the embedding representation. It is observed in the t-SNE visualization in Fig. 3 that the embeddings can preserve the local structure of the interactions and functional annotations, as well as exhibit good separation based on their transcript biotype classification. There is a clear separation between miRNAs to lncRNAs and mRNAs. While the protein-coding mRNAs and lncRNAs may have some overlap, is expected since the sequence structure of these two RNA classes is similar. Note that although no negative undirected edges between RNAs of different locus types (e.g., lncRNAs v.s. miRNAs) were sampled to explicitly indicate different RNA types to have dissimilar embeddings, the network can still make a distinction between their functional roles. This shows that rna2rna's source-target embedding is an effective representation that can encode an RNA's biological function only by its given directed interactions.
Our main contribution proposes a highly versatile architecture aimed at predicting interactions between heterogeneous RNA transcripts while characterizing the functional landscape of non-coding RNAs. In conclusion, we intend this method to be the groundwork for further down-stream analysis tasks, where various other downstream genomic prediction tasks such as prediction of gene annotation, gene-disease association, and discovery of unknown gene families can be readily applicable. Further works to this framework can provide invaluable tools to support significant discoveries in systems biology, especially for newly identified lncRNAs.