Method for producing functional RNA sequence

Abstract

Provided is a method for producing a functional RNA sequence based on a deep generative model.　The method for generating a functional RNA sequence comprises: a) a step for obtaining input data X acquired by digitizing stored secondary structure information on a training RNA sequence; b) a step for inputting the input data X on a plurality of training RNA sequences into a learning model to obtain generated output data X'; and c) a step for presenting a plurality of generated RNA sequences on the basis of the output data X'.

Description

Methods for generating functional RNA sequences

[0001] The present invention relates to a method for generating functional RNA sequences.

[0002] Over the past few decades, leveraging the functionality of RNA has led to the invention of various synthetic molecular systems, which have had a profound impact on many fields, including basic research, biomanufacturing, and medical applications. Most of these experiments have focused on designing specific types of functional RNAs with desired properties, such as riboswitches, riboregulators, and ribozymes, based on manual manipulation, repetitive trial-and-error, or massively parallel arrays of random mutants. However, because the sequence space of RNA is vast and complexity increases exponentially with sequence length, experimental searches for desired RNA sequences remain costly and inefficient. Therefore, there is a strong need for versatile computational platforms to understand sequence space and efficiently design functional RNAs.

[0003] Computational attempts at functional RNA design have primarily focused on RNA reverse folding, an approach that finds sequences that fold into specific secondary structures based on RNA secondary structure prediction and individual optimization algorithms. Previous studies have used RNA reverse folding to design consensus secondary structures for desired functional RNAs, resulting in the acquisition of several functional RNAs. However, because RNA function is not characterized solely by structure, designed sequences must be selected ad hoc. Furthermore, RNA reverse folding lacks flexibility and versatility, making it difficult to generate sequences with acceptable variations for the desired function. The accuracy is inherently limited by the accuracy of the RNA secondary structure prediction and the optimization algorithm.

[0004] Unlike RNA reverse folding, the covariance model (CM), a statistical framework for RNA alignment and consensus secondary structure, allows for quantitative evaluation of sequence and structural variation without relying on RNA secondary structure prediction. Due to its flexibility, CM has remained the gold standard for RNA homology searches for decades, classifying most functional RNA species into thousands of "RNA families" (see, for example, Non-Patent Document 1). Although CM has not been used for functional RNA design to date, it may provide a promising framework for functional RNA design that circumvents the technical difficulties associated with RNA reverse folding.

[0005] The optimization difficulty associated with searching exponentially large sequence spaces using discrete algorithms can be alleviated by array representations (see, for example, Non-Patent Document 2). Array representations embed discrete sequences into continuous spaces based on machine learning. To efficiently design sequences with desirable properties, it is important to learn meaningful representations while maintaining continuity with the semantics of the sequences (e.g., sequence and structural similarity) (see, for example, Non-Patent Document 2).

[0006] Recently, deep generative models for learning semantically rich sequence representations have attracted much attention. In particular, deep generative models with an internal representation called a "latent space," such as variational autoencoders (VAEs), can generate synthetic data that mimic existing data while embedding the original data in the latent space (see, for example, Non-Patent Document 2). Because the latent space allows for continuous interpolation of data, VAEs can generate sequences with desired properties through continuous optimization, such as backpropagation and Bayesian optimization.

[0007] Eddy, S., Durbin, R. RNA sequence analysis using covariance models. Nucleic Acids Res. 22, 2079-2088 (1994). Hawkins-Hooker, A. et al. Generating functional protein variants with variational autoencoders. PLoS Comput. Biol. 17(2), e1008736 (2021).

[0008] A method for generating functional RNA sequences based on deep generative models is needed.

[0009] That is, the present invention includes the following aspects. [1] A method for generating functional RNA sequences, comprising the steps of: a) obtaining input data X by digitizing secondary structure information stored in training RNA sequences; b) inputting the input data X for multiple training RNA sequences into a training model to obtain generated output data X'; and c) presenting the generated multiple RNA sequences based on the output data X'. [2] A method for generating functional RNA sequences according to [1], wherein the step a) is a step of obtaining input data X by digitizing alignment information of the training RNA sequences in addition to secondary structure information of the training RNA sequences. [3] A method for generating functional RNA sequences according to [1], wherein the step a) is a step of obtaining the digitized input data X based on a covariance model (CM). [4] The method according to [3], wherein the step a) is a step of obtaining the digitized input data X by aligning the training RNA sequence to a CM and expressing a path on the CM in a numerical matrix, and separating non-terminal symbol portions, base pair portions, and non-base pair portions from the matrix of the numerical matrix expression, wherein the numerical matrix expression is a matrix with CM state indexes as columns and used transition / output rules as rows, assigning 1 to rules that have been used and 0 to others, and assigning the number of times the insertion transition rule has been used if multiple insertion transition rules are used. [5] The method according to [3], wherein the input data X is expressed as X=(tr, ss, bp), where tr represents non-terminal symbol transition information, ss represents single-strand output information, and bp represents base pair output information. [6] The method according to [3], wherein the step b) includes a step of transforming the plurality of pieces of input data X into latent variables Z using a first nonlinear function, and obtaining output data X' from the latent variables Z using a second nonlinear function, wherein the first and second nonlinear functions are parameterized, and the parameters are learned so as to minimize a first term that brings the input data X and the output data X' closer together, plus a second term that brings the distribution of the latent variables Z obtained by transforming the input data X closer to a prior distribution.[7] The method according to [3], wherein the learning model performs machine learning of the input data X using a variational autoencoder. [8] The method according to [6], wherein the first term represents the cross-entropy error between the input data X and the output data X', and the second term represents a term obtained by multiplying a weight by the Kullback-Leibler distance between the distribution of Z and a normal distribution with a mean of 0 and a variance of 1. [9] The method according to [6], wherein the first nonlinear function is a first neural network and the second nonlinear function is a second neural network.

[10] The method according to [3], wherein step c) includes: permuting parameters of the CM of step a) based on the output data X' to create an output CM; and obtaining an RNA sequence with a maximum probability from the output CM.

[11] The method according to any one of [1] to

[10] , wherein the generated functional RNA includes linear or circular mRNA or non-coding RNA.

[12] The method according to any one of [1] to

[10] , wherein the generated functional RNA exhibits activity equivalent to that of the training RNA sequence in vitro, in vivo, ex vivo, and / or in cellulo (e.g., in a test tube including in an aqueous solution or on paper, in a liposome, in a cell, in an organoid, in a tissue, and / or in a living organism).

[13] The method according to any one of [1] to

[10] , wherein the functionality includes RNA activity selected from RNA activity corresponding to an RNA family registered in the Rfam Database, RNA activity corresponding to all RNA families known from literature, etc., and any other RNA activity.

[14] The method according to any one of [1] to

[10] , wherein the functionality includes riboswitch activity, RNA interference activity, aptamer activity, tRNA activity, signal recognition activity, snRNA activity, snoRNA activity, guide RNA activity, Y RNA activity, long non-coding activity, eRNA activity, microRNA activity, piRNA activity, siRNA activity, rRNA activity, genome editing activity, RNA editing activity, transposon activity, ribosome activity, spliceosome activity, and / or ribozyme activity (including, for example, self-cleaving ribozyme activity and RNA methylation ribozyme activity).

[15] The method of

[13] or

[14] , wherein the step a) includes a step of obtaining input parameters X for training RNA sequences having a common sequence or structure.

[16] The method of [6], wherein in the step b), in a coordinate space visualized using t-SNE of a latent space represented by the latent variable Z, the relative distance between variables into which the multiple training RNA sequence information has been transformed represents the similarity in properties between the training RNA sequences.

[17] A method for generating functional RNA sequences, comprising: A) obtaining a new latent variable Z' for a trained latent variable Z, transforming the latent variable Z', and obtaining output data X'; and B) presenting multiple generated RNA sequences based on the output data X', wherein the latent variable Z has information on the sequence, structure, biological species, and / or activity level of input data X based on the training RNA sequences, and the latent variable Z and the new latent variable Z' have similar information.

[18] A computer program causing an information processing device having at least one processor to execute each step of the method of any one of [1] to

[17] .

[19] A non-transitory computer-readable medium comprising the computer program according to

[18] .

[20] The method according to any one of [1] to

[17] , executed by an information processing device having at least one processor.

[21] An information processing device having at least one processor, wherein the processor executes the method according to any one of [1] to

[17] .

[0010] The method for generating functional RNA sequences according to the present invention can generate a substantial number of novel functional RNA sequences from various datasets, providing broad applicability. In particular, the method according to the present invention can generate sequences from diverse RNA families, such as tRNAs, ribozymes, and riboswitches, encompassing a wide range of structures and lengths. Furthermore, the method according to the present invention can not only encapsulate sequence semantics within a latent space, but also predict variant effects, which is useful for designing sequences with desirable properties.

[0011] Figure 1A shows an overview of RfamGen, one embodiment of the present invention. RfamGen is a variational autoencoder trained on feature maps of sequence alignments on CM. RfamGen encodes each feature map into a latent space and decodes it from Z. Sequences are generated by taking the maximum probability from the reconstructed feature maps. The latent space aggregates the semantics of the training data, such as variant information and sequence activity. The sequences generated by RfamGen were evaluated through validation experiments, including massively parallel assays. Figure 1B shows the characterization in RfamGen. CM is an extended Markov model constructed from multiple sequence alignments and consensus secondary structures. CM includes states, transition rules between states, and exit rules from states. All states and rules are indexed by integers. A sequence ("gaac") is aligned to the CM (left) as a path on the CM, and its structure is inferred by the alignment ("(..)"). The paths on the CM can be decomposed into each rule for each state (center) and expressed as one-hot sequences, where rows and columns indicate the rule for the CM and the state index within the CM (right). Figure 2A shows the results of the RfamGen evaluation, showing the evaluation of the likelihood of RNA family membership ("bit scores") for sequences generated by RfamGen (alignment (+), secondary structure (+)) and three benchmark models (GCVAE, GVAE, and CVAE). Eighteen RNA families with multiple aligned sequences were used. Each benchmark lacked either alignment or secondary structure information: GCVAE (alignment (+), secondary structure (-)), GVAE (alignment (-), secondary structure (+)), and CVAE (alignment (-), secondary structure (-)). Scores were normalized by the GA threshold for each RNA family in the Rfam database to align the plots. Figure 2B shows the bit scores for sequences generated by RfamGen and GCVAE under various undersampling conditions. The line graph and shaded area show the mean and standard deviation of five trials. Figure 2C shows a scatter plot of the bit scores of the sequences of 628 RNA families generated by RfamGen and GCVAE.Figure 2D shows a scatter plot of the sequence information content of natural and generated sequences for 628 RNA families. Sequence information content was calculated from the position entropy of MSA. The ideal distribution is shown by the dashed line. Figure 2E shows a scatter plot of the secondary structure information content of natural and generated sequences for 628 RNA families. Secondary structure information content was calculated based on the consensus secondary structure defined by CM. The ideal distribution is shown by the dashed line. Figure 3A shows that the latent space has a semantically rich representation of sequences. tSNE visualization of the latent space trained for 5S ribosomal RNA (RF00001) and U1 spliceosomal RNA (RF00003), colored by phylogenetic information. Figure 3B shows the predicted activity of tRNA variants (y-axis) and predicted activity by linear regression of the latent space (x-axis). Figure 3C shows a comparison of phylogenetic distance (x-axis) and Euclidean distance (latent distance, y-axis) in the latent space trained by Twister Sister ribozyme (RF02681). Figure 3D shows a tSNE visualization of the latent space trained by Twister Sister ribozyme (RF02681) colored by structural variants (left). Examples of structural variants of Twister Sister ribozymes (SEQ ID NOs: 33 and 34), defined by the presence or absence of variable stem-loops indicated by black boxes, are shown (right). These were visualized using the forna server (http: / / nibiru.tbi.univie.ac.at / forna / ). Figure 4A shows a consensus diagram of the natural and generated sequences of the glmS ribozyme, with the ribozyme cleavage site indicated by a black arrow. Figure 4B shows a consensus diagram of the native and engineered sequences of the Pistol ribozyme (from Zasha Weinberg data), with the ribozyme cleavage site indicated by a black arrow. Figure 4C shows a consensus diagram of the native and engineered sequences of the Twister Sister ribozyme, with the ribozyme cleavage site indicated by a black arrow. Figure 4D shows the percent identity to the native sequence of the glmS ribozyme. Figure 4E shows the percent identity to the native sequence of the Pistol ribozyme. Figure 4F shows the percent identity to the native sequence of the Twister Sister ribozyme. Figure 4G shows a PAGE separation of the ribozyme cleavage products of the glmS ribozyme.In the graph, "mut" indicates a mutant with reduced activity. For the glmS ribozyme, cleavage was assessed with and without ligand, and arrows indicate the cleaved products. The percentage of cleaved products was quantified from gel images (n = 3). A paired t-test (one-tailed) was performed to compare the cleavage rate between each sequence and its corresponding control. Error bars indicate standard deviation. Figure 4H shows PAGE separation of ribozyme cleavage products from the Pistol ribozyme. Figure 4I shows PAGE separation of ribozyme cleavage products from the Twister Sister ribozyme. Figure 5A shows an overview of the massively parallel assay of sequences generated from RfamGen. The glmS ribozyme cleaves the 5' end of sequences in a ligand-dependent manner. Massively parallel assays of generated sequences, native sequences, and variants were performed as reported by Andreasson et al., 2020 (see Table 3). All sequences were pooled and autonomously cleaved upon addition of ligand. The cleavage reaction was initiated by the addition of ligand and terminated by the addition of quenching buffer after various incubation times (0, 1, 3, 5, 10, 30, and 60 minutes). Kinetic parameters (kcat) were estimated from the cleaved and uncleaved products sequenced by a massively parallel sequencer. Figure 5B shows the time-dependent cleavage rates of the generated sequences. Figure 5C shows the distribution of kcat values ​​for the native and generated sequences. The Rfam database logo was obtained from the Rfam database. Figure 6 illustrates the one-hot representation of the alignment to CM in one embodiment of step a) of the present invention. Figure 7 illustrates the division of the one-hot representation into three matrices in one embodiment of step a) of the present invention. Figure 8 illustrates the node information of CM in one embodiment of step a) of the present invention. Figure 9 illustrates the decoding of the alignment and generation of sequences in one embodiment of step c) of the present invention.

[0012] Hereinafter, embodiments of the present invention will be described, but the present invention is not limited to the embodiments described below.

[0013] 1. Method for generating functional RNA sequences (1) According to a first embodiment, the present invention provides a method for generating functional RNA sequences, comprising the following steps a) to c): a) obtaining input data X that digitizes secondary structure information stored in a training RNA sequence; b) inputting the input data X for multiple training RNA sequences into a training model to obtain generated output data X'; and c) presenting the generated multiple RNA sequences based on the output data X'.

[0014] As used herein, a functional RNA sequence refers to an RNA sequence that exhibits any activity in vitro, in vivo, ex vivo, and / or in cellulo. The activity is not particularly limited and may include all activities known in literature, etc., and those that will be discovered in the future. Examples of activities include, but are not limited to, riboswitch activity, RNA interference activity, aptamer activity, mRNA activity, long non-coding RNA (lncRNA) activity, tRNA activity, signal recognition activity, snRNA activity, snoRNA activity, guide RNA activity, Y RNA activity, long non-coding activity, eRNA activity, microRNA activity, piRNA activity, siRNA activity, rRNA activity, genome editing activity, RNA editing activity, transposon activity, ribosome activity, spliceosome activity, and / or ribozyme activity (e.g., including self-cleaving ribozyme activity, RNA methylation ribozyme activity, and transfer ribozyme activity).

[0015] The "RNA sequence" produced is a linear or circular single-stranded RNA, a double-stranded RNA, or a fusion thereof. The nucleotides constituting the produced RNA sequence include the natural bases A, G, C, and U, and may further include modified bases and DNA. The modified bases that can be included are not particularly limited and may be any modified base. For example, modifications may be deamination, methylation, pseudouridylation, halogenation, thiolation, or amination, but are not limited to these. More specifically, modified bases include, but are not limited to, inosine, 5-methylcytosine (5-me-C), 5-methyluridine (5-me-U), N6-methyladenosine (N6-me-A), 1-methyladenosine (1-me-A), pseudouridine, 1-methylpseudouridine (1-me-Pseudouridine), 5-bromouridine (Br-U), 5-iodouridine (5-IU), 6-thioguanosine (6-thio-G), 4-thiouridine (6-thio-U), 2-aminopurine, and 2,6-diaminopurine (2-amino-A). The type of RNA represented by the RNA sequence may be mRNA or non-coding RNA, and is not particularly limited. The "RNA sequence" generated may be of one type or of multiple types, and may differ depending on the input information.

[0016] As used herein, generating a functional RNA sequence refers to determining and outputting the primary sequence of the RNA. Due to the nature of nucleic acid bases, the primary sequence of the RNA may also contain structural information. In the method of the present invention, the specific "functional RNA sequence" generated may be referred to as the "target product." Each of the steps a) to c) is described below.

[0017] a) Step of Obtaining Input Data X In this step, the "training RNA sequence" refers to the RNA sequence that the training model is to learn from. The training RNA sequence may be an RNA sequence that shares a sequence and structure with the "functional RNA sequence" that is the target of the method of the present invention. Therefore, the training RNA sequence can be appropriately determined by a person skilled in the art depending on the target. It can be selected from known RNA sequences obtained from RNA databases or literature information, RNA sequences predicted from known RNA, RNA sequences experimentally obtained using a sequencer, RNA sequences experimentally obtained from libraries containing random sequences by in vitro selection (SELEX), and RNA sequences experimentally obtained from mutant libraries. RNA databases include, but are not limited to, the Rfam database, the Zasha Weinberg Database (ZWD), RNAcentral, the NT database, and the Ensembl database.

[0018] Like the target RNA, the training RNA sequence may be linear or circular, single-stranded RNA or single-stranded DNA. The length of the training RNA sequence is not particularly limited, as long as it has the intended "functionality." The length of the training RNA sequence may be, for example, 8 or more bases and approximately 2,000 or less bases. Furthermore, the nucleotides constituting the training RNA sequence may include the natural bases A, G, C, and U, and, like the target RNA, may further include DNA or modified bases.

[0019] The "conserved secondary structure information" of a training RNA sequence refers to evolutionarily conserved secondary structure information of RNA, and for known RNAs, this information is known from databases, etc. In one embodiment, the input data X may be quantified data of the secondary structure information of the training RNA sequence. Using quantified data of conserved secondary structure information is advantageous in that it allows for the elimination of unnecessary structural information generated by secondary structure prediction, etc., by considering only secondary structures important for functional expression. Conserved secondary structures, also known as consensus structures, can be derived by those skilled in the art using known methods. For example, RNAalifold, PETfold, TurboFold, RNAstructure, RNAz, Pfold, Construct, CentroidAlifold, RNAcofold, LinearAlifold, and other methods can be computationally determined using multiple sequence alignment-based, comparative genomic analysis, comutation analysis, energy minimization, or machine learning-based methods, and these methods are commonly used. Alternatively, biochemically determined structures or manual determination of structures are also possible.

[0020] In another embodiment, the input data X may be data obtained by digitizing the stored secondary structure information and alignment information of a training RNA sequence. The "alignment information" of a training RNA sequence refers to information about sequences common to multiple sequences. For known RNAs, alignment information is also known from databases or the like. Using data obtained by digitizing both the stored secondary structure information and alignment information can be advantageous in that it can generate RNA sequences with higher activity in terms of the desired functionality.

[0021] In this step, preferably, a plurality of training RNA sequences is used to obtain input data X. In this step, "a plurality" specifically means 2 or more, preferably 50 or more, and more preferably 200 or more, but is not limited thereto. In theory, there is no particular upper limit, but with training RNA sequence information of about 500, a target substance with sufficient activity can be obtained.

[0022] The quantification of the conserved secondary structure and alignment information of the multiple training RNA sequences can be performed, for example, based on the covariance model (CM). Alternatively, the quantification of the secondary structure and alignment information can be performed using other methods that input numerical representations of the secondary structure and alignment separately.

[0023] Next, we will explain in more detail how stored secondary structure information and alignment information are quantified based on CM. If the CM is known, an existing CM can be used. If a CM does not exist, a CM can be created using the cmbuild command. The sequence of a training RNA can be aligned to the CM. Alignment of an RNA sequence to a CM can be performed using software, for example, but not limited to, the cmalign command in Infernal (Nawrocki, EP & Eddy, SR Infernal 1.1: 100-fold faster RNA homology searches. Bioinformatics 29, 2933-2935 (2013)). Other methods for aligning an RNA sequence to a CM include the CYK algorithm and query-dependent banding. Specific examples and explanations of aligning an RNA sequence to a CM are described in detail in the Examples.

[0024] Next, the paths on the aligned CM are expressed as a numerical matrix. The numerical matrix representation is a matrix with the CM's state index as columns and the transition / output rules used as rows. The state index is a numerical value assigned to the state used in the CM and is determined by the order in which the non-terminal symbols of the context-free grammar appear when the syntax tree is created. The context-free grammar used is the G1 grammar (Eddy, S., Durbin, R. RNA sequence analysis using covariance models. Nucleic Acids Res. 22, 2079-2088 (1994)), but is not limited to this. Other context-free grammars that can parse bracket notation for RNA secondary structure, such as G2, G3, G4, and G5, can also be used. This process is illustrated in Figure 1B. A 1 is assigned to the transition / output rule that was used, a 0 is assigned to the others, and if multiple insertion transition rules are used, the number of times the insertion transition rule was used is assigned. This mimics one-hot representation. One-hot representation is a method of quantifying categorical data in machine learning. Categorical data refers to data that is expressed as a category rather than a numerical value, such as color (red, blue, green, etc.) or gender (male, female). In one-hot representation, a vector consisting of 0 and 1 is assigned to each category. For example, if it belongs to a certain category, it is set to 1, and if it does not, it is set to 0. The length of this vector is the same as the number of categories. In this case, the transition / output rules correspond to the categories, and the number of one-hot representations increases by the number of rule usage orders.

[0025] Next, the numerical matrix representation is further decomposed into a triplet of three matrices (tr, ss, bp). Here, tr represents the transition probability of the non-terminal symbol, ss represents the single-strand output probability, and bp represents the base pair output probability. Because the numerical matrix representation contains many zeros or missing values ​​that cause learning instability, we sought to eliminate these. ss and bp are transformed using a reindexing function, while tr is transformed as is, resulting in three matrices (trans, ss, bp). "Transformation" here refers to selecting specific rows and columns from second-order matrices such as tr, ss, and bp to create a new second-order matrix. For the transformation of ss, we selected only the columns corresponding to the state indexes that can output single-stranded RNA, and selected the rows corresponding to the output of single-stranded RNA to create a new second-order matrix. For the transformation of bp, we selected only the columns corresponding to the state indexes that can output double-stranded RNA, and selected the rows corresponding to the output of double-stranded RNA to create a new second-order matrix. Also, trans was transformed to tr by further using node indexing to obtain a representation without missing values.

[0026] The triplet X = (tr, ss, bp) obtained in this manner can be used as input information in step a). One input information triplet X is prepared for one training RNA. Therefore, in step a), the number of input information triplets X prepared corresponds to the number of training RNAs used. The triplet X created in step a) contains at least information on the conserved secondary structure.

[0027] b) Step of Obtaining Output Data X' In step b), the multiple pieces of input data X obtained in step a) are input into the learning model to obtain output data X'.

[0028] Step b) may include a step of converting a plurality of input data X=(tr, ss, bp) into a latent variable Z using a first nonlinear function, and obtaining output data X'=(tr', ss', bp') from the latent variable Z using a second nonlinear function. Both the first nonlinear function and the second nonlinear function may be neural networks. When the input data is generated based on CM in step a), the first nonlinear function may be a CM encoder, and the second nonlinear function may be a CM decoder. Note that the above output data X'=(tr', ss', bp') can be mathematically expressed as follows: is correct, but for convenience of format in this specification, this expression is replaced by X'=(tr', ss', bp').

[0029] The first and second nonlinear functions can be parameterized. The nonlinear function can be expressed as a loss function L including a first term and a second term. The first term of the loss function L is a term that brings the input data X = (tr, ss, bp) and the output data X' = (tr', ss', bp') closer together, and the second term is a term that brings the distribution of the latent variable Z obtained by transforming the input data X closer to the prior distribution. The parameters can be learned to minimize the loss function L. More specifically, the first term of the loss function L represents a term obtained by weighting the cross-entropy error between the input data X and the output data X' depending on the distribution of the training data X, and the second term of the loss function L represents a term obtained by multiplying the Kullback-Leibler distance between the distribution of Z and a normal distribution with a mean of 0 and a variance of 1 by a weight. Data weighting is often performed to reduce observation bias; the more similar data there are, the less weight there is of that data. In one embodiment, this process can be performed using a variational autoencoder (VAE). The details of the parameter formulas are shown in the Examples. The first and second nonlinear functions can be used independently. The first nonlinear function embeds data X = (tr, ss, bp) into latent variables for visualization and data exploration, and the second nonlinear function is used to generate data from the latent variables.

[0030] Let's explain VAE in more detail. A VAE consists of two functions: an encoder and a decoder. The encoder converts data X into a continuous representation Z, and the decoder converts the continuous representation Z back into data X. The general framework of a VAE was proposed by Kingma and Welling in 2013. After creating the data X and weights, this study follows this same framework for training. Training is completed by updating the parameters of both the encoder and decoder to minimize the loss function (cross entropy + Kullback-Leibler distance). The continuous representation Z does not necessarily correspond 1:1 to the data available at the time of training; rather, Z approximates the data X. This means that once the VAE has been trained, it is possible to acquire a variety of new data by selecting various values ​​of Z and inputting them into the decoder to convert them into data. Furthermore, the spread of the continuous representation can be specified during training. Typically, the continuous representation Z of the training data is specified to approximate a Gaussian distribution. This can be achieved by training the encoder to achieve this. Among multiple Zs, two or more neighboring Zs are converted into two or more neighboring data X', so by observing how the training data is spread, it is possible to search for Zs that contain the target data X' and Z's that are neighboring to it.

[0031] This step is performed by "learning" (also called training) using a learning model. "Learning" refers to updating the parameters of the first and second nonlinear functions so as to minimize the loss function L. "Learning" can generally be performed multiple times, for example, one or more times, preferably about five times, and more preferably about 50 times, but the number of times is not limited to a specific number. By performing the above-mentioned predetermined number of times of learning, output data X' = (tr', ss', bp') can be obtained. Note that the method of step b) itself is known as a conventional method, and any established method can be used. Therefore, the learning method is not limited to one using a variational autoencoder, and any learning method can be applied.

[0032] c) Step of presenting an RNA sequence based on the output data X' In step c), an alignment is reconstructed from the triplet of output data X' = (tr', ss', bp') to generate a sequence. Step c) can be performed by completely reversing the procedure of step a). More specifically, the parameters of the CM in step a) are replaced to create an output CM, and the sequence with the highest probability is obtained from the output CM, thereby generating the desired functional RNA sequence. More detailed information about the method of step c) is provided in Supplementary Note 2 of the Examples.

[0033] The functional RNA sequences obtained through steps a), b), and c) of the present invention have been confirmed to have sufficient activity with respect to the functionality focused on when determining the training RNA sequence, and in some cases, have higher activity than the training RNA sequence. According to the method of the present invention, highly active functional RNA sequences can be obtained using a relatively small number of training RNA sequences, and can be used to obtain desired functional RNA sequences in a wide range of fields, including basic research, biomanufacturing, and medical applications.

[0034] 2. Method for generating functional RNA sequences (2) According to a second embodiment, the present invention provides a method for generating functional RNA sequences, which includes the following steps A) and B): A) obtaining a new latent variable Z' for a trained latent variable Z, and converting the new latent variable Z' to obtain output data X'. B) presenting a plurality of generated RNA sequences based on the output data X'. In step A), the latent variable Z has information on the sequence, structure, biological species, and / or activity intensity of the input data X based on the training RNA sequence. The latent variable Z and the new latent variable Z' have similar information.

[0035] In the second embodiment, the definitions of functional RNA, training RNA, input data X, and output data X' are substantially the same as in the first embodiment. Therefore, machine learning using the training RNA and input data X based on it to obtain a latent variable Z can be performed in the same manner as in the first embodiment. This relates to a method for selecting a new latent variable Z' from a latent variable Z obtained using a first nonlinear function having parameters already trained on input data X of multiple training RNAs to obtain output data X'. This is based on the fact that in step b) of the first embodiment, in a coordinate space visualized using t-distributed Stochastic Neighbor Embedding (t-SNE) to visualize the latent space represented by the latent variable Z, the relative distance between variables converted from multiple training RNA sequence information represents the similarity in characteristics between each training RNA sequence. Therefore, by obtaining a new latent variable Z' that is different from the latent variable Z obtained by learning but follows a Gaussian distribution, desired output data X' can be obtained that is different from the output data converted from the latent variable Z. The new latent variable Z' can also be referred to as a latent variable Z' that is close to the latent variable Z. A latent variable Z' close to a latent variable Z can refer to proximity in the coordinate space visualized using t-SNE.

[0036] According to the method for generating a functional RNA sequence of the second embodiment, it is possible to obtain RNA having a primary sequence different from existing RNA sequences, and it is possible to obtain RNA having specific sequence or structural characteristics, RNA with altered function, or RNA having a function not present in existing RNA.

[0037] 3. Computer Program, Computer-Readable Medium, and Information Processing Device According to a third embodiment, the present invention relates to a computer program for implementing the method of the first or second embodiment, a computer-readable medium storing the method, and an information processing device that implements the method.

[0038] The computer program may cause a processor of an information processing device to execute the steps of the method of the first or second embodiment. Such a computer program may be stored on a computer-readable medium. The computer-readable medium may be, for example, but is not limited to, a RAM memory, a magnetic or optical drive, a floppy disk, or any other device. The computer program is stored on the computer-readable medium or stored in an information processing device such as a computer, and is configured to specifically execute the method of the first or second embodiment by the information processing device, which is a hardware resource.

[0039] Furthermore, the information processing device according to this embodiment includes at least one processor, a memory, an output device, and an input device. The processor may include a control device and an arithmetic device and may be capable of executing the instructions of the computer program. The memory may be capable of storing the instructions of the computer program for implementing the method of the first or second embodiment. The output device and input device may be any general-purpose device depending on the form and purpose of the information processing device, such as a keyboard or display of a general-purpose computer. By performing the steps of the method of the first or second embodiment in the information processing device, a novel functional RNA sequence can be obtained as output.

[0040] According to the third embodiment of the present invention, the steps of the method of the first or second embodiment can be configured to be executable by an information processing device such as a computer, and can be widely used.

[0041] The present invention will be described in more detail below with reference to examples, which are not intended to limit the scope of the present invention.

[0042] [Experimental Method] We performed RfamGen, an example of a method for generating functional RNAs developed by the present inventors. Sequences and CMs were collected from the Rfam database (version 14.7). Identical sequences and sequences containing nucleotides other than AUGC were removed. Sequences were aligned to CMs using the cmalign command in Infernal (version 1.1.3), and the alignment was converted to a one-hot representation in CM grammar. Details of CM grammar were as follows: Eddy, S.R. & Durbin, R., RNA sequence analysis using covariance models. Nucleic Acids Res. 22, 2079-2088 (1994).

[0043] CM is a context-free grammar that simultaneously models multiple sequence alignment (MSA) and consensus secondary structure. Each sequence in an MSA is aligned to a CM and can be represented as a single path on the CM. Because CMs have a repetitive structure, we represent them as a probability matrix, where each row represents the probability of a CM rule and each column represents the index of a CM state (see "Sequence Generation from CM" below). Similarly, sequences aligned to a CM are represented as one-hot sequences. Next, the matrix was further decomposed into three triplets: X = (tr, ss, bp). Each matrix represents the transition, single-strand output, and base pair output probabilities of non-terminal symbols, respectively. RfamGen was configured to accept the triplet X = (tr, ss, bp) as input and generate the reconstructed triplet X' as output. A more detailed process is described in "Supplementary Note 1" below.

[0044] To reduce the observation bias in the Rfam database, each sequence was reweighted by the number of sequences outside its neighborhood, as proposed in a previous study (Hopf, TA et al. Mutation effects predicted from sequence co-variation. Nat. Biotechnol. 35, 128-135 (2017)). To define the neighborhood of a sequence s, we use the distance D between two sequences s and t. H (s, t) is defined (Equation (1)). Next, D H Using (s, t), we defined the weight πs of the array s with a threshold θ (Equation (2)). where I is the index function. In previous studies, we set θ = 0.2 for protein MSA, but we found that θ = 0.1 gave slightly better results for many RNA families, so we set θ = 0.1 in this study. Detailed parameters are listed in Table 1.

[0045] RfamGen Architecture and Training As mentioned above, a sequence is represented by a triplet X = (tr, ss, bp). We used a variational autoencoder (VAE) to generate a synthetic triplet X' = (tr', ss', bp') that mimics the input triplet X. RfamGen consists of a CM encoder q_φ(Z|X) and a CM decoder p_θ(X|Z). The CM encoder converts the triplet X into a latent variable Z using three one-dimensional convolutional layers followed by two fully connected layers. The CM decoder recovers the triplet X from the latent variable Z using two fully connected layers followed by three one-dimensional deconvolutional layers. (The dimensionality of the fully connected layers is shown in Table 1 for the hidden state number. The dimensionality of the latent variable Z is shown in Table 1 for the latent dimension.) To train RfamGen, we defined the loss function L as follows (Equation (3)): During the ceremony, and KL(q|p) is the Kullback-Leibler divergence (KLD) between distributions p and q, and β is for relaxing KL vanishing for KL(q φ (Z|X)|p(Z)), where p(Z) is a Gaussian distribution. The parameters of the Gaussian distribution p(Z) used are referred to sampling in Table 1. For training RfamGen, the Adam optimizer (learning rate (refer to Table 1 for the learning rate), a periodic annealing schedule for β, and early stopping) was used. β was determined empirically by examining the KLD. The detailed parameters of this experiment are described in Table 1.

[0046]

[0047] Array Generation from CM Arrays were generated by obtaining the most likely array from the reconstructed CM using a similar approach to Kusner, M. J., Paige, B. & Hernandez-Lobato, J. Grammar Variational Autoencoder. Preprint at https: / / arxiv.org / abs / 1703.01925 (2017). The CM was reconstructed by replacing the values of the CM registered in the Rfam database, and the argmax rule was obtained from each state of the CM to generate the array. Formally, the CM is summarized as follows. Each state and rule is indexed by an integer, and the transition T i→j and the branch B i→(k,l) must have i < j and i < k < l. Given a non-terminal state S i , for each possible rule for S i : p(r i ) may be used, and is. To generate an array from the reconstructed CM, from the i-th rule r i the most likely rule The alignment was constructed on the CM by stepwise selection of the CMs. The detailed procedure is described in Supplementary Note 1.

[0048] In silico evaluation of RfamGen's generative ability. To evaluate RfamGen's generative ability computationally, we generated sequences from the RfamGen, CVAE, GVAE, and GCVAE models. Details of the benchmark models other than RfamGen are provided in Supplementary Note 2. We sampled 1,000 points in the latent space from a Gaussian distribution, decoded them into CMs, and generated argmax sequences from each CM. Next, we performed alignments with "ground truth" CMs registered in the Rfam database using Infernal (version 1.1.3) and compared the bit scores of the sampled sequences from the models. In Figure 2A, the bit scores were normalized by the GA threshold, which is the threshold used in the Rfam database to align the graph baseline.

[0049] For information content analysis, we used the strategy described by Wilburn, GW & Eddy, SR. Remote homology search with hidden Potts models. PLoS Comput. Biol. 16, e1008085 (2020). For information content at 1°, we calculated position-wise entropy. For information content at 2°, we calculated pairwise information (mutual information) constrained by the nested consensus secondary structure. We used two different consensus secondary structures: (1) a covariance model consensus secondary structure and (2) a consensus secondary structure determined by the Cacofold algorithm in R-scape software (http: / / eddylab.org / R-scape / ). We removed pseudoknots and extracted the nested secondary structure.

[0050] Design and Analysis of Diverse Self-Cleaving Ribozymes. The Pistol ribozyme sequence was curated from the ZWD, while the glmS ribozyme and Twister Sister ribozymes were curated from the Rfam database (version 14.7). Because the ZWD does not contain a CM for the Pistol ribozyme, a CM was constructed using the cmbuild command in Infernal (version 1.1.3). To prepare the dataset for training and evaluation of RfamGen, the dataset was split into training / validation / test in a ratio of 7:1.5:1.5. RfamGen was trained 20 times, and the optimal model was selected using the test dataset. Detailed parameters for this experiment are listed in Table 1. To evaluate the identity of the generated sequences, 1,000 sequences were randomly generated from a Gaussian distribution in the latent space and the minimum Levenshtein distance (global alignment) to the training sequence was calculated.

[0051] For experimental validation, three sequences were randomly selected from the sampled sequences. In addition to the three generated sequences, a native sequence and a loss-of-function mutant reported in a previous study were prepared as controls. ssDNA was purchased from Eurofins genomics, and RNA was prepared using the MEGAshortscript™ T7 Transcription Kit (Thermo Fisher Scientific, AM1354). The glmS ribozyme was cleaved in a buffer solution (10 mM GlcN6P, 50 mM HEPES pH 7.8, 150 mM KCl, 4 mM Mg(OAc)2, 100 μM Na2EDTA, 1 mM DTT) with or without 10 mM GlcN6P. The reaction was stopped by adding an equal volume of gel loading buffer II (included in Thermo Fisher Scientific, AM1354). For the pistol and twister sister ribozymes, spontaneous cleavage occurs during transcription, so purified RNA was used directly for electrophoresis. Reaction products were separated using denaturing 10% polyacrylamide / 7 M urea / 1×TBE and imaged using a Typhoon FLA 7000 (Cytiva). The consensus diagrams in Figures 4A, 4B, and 4C were visualized using R2R software (version 1.0.6). The sequences used in this assay are listed in the Sequence Listing, and their descriptions are provided in Table 2 below.

[0052]

[0053] Design of glmS Ribozymes for In Vitro Validation and Massively Parallel Cleavage Assays. Ribozymes were generated by sampling from a Gaussian distribution in the latent space and decoding to sequences from the best model used in the previous section. The decoding process was described in the Experimental Methods section. Oligo pools were designed by flanking the reverse-complementary sequences with the reverse-complementary T7 promoter sequence (CTATAGTGAGTCGTATTAGC (SEQ ID NO: 29)) and the RT primer binding site (CGTCGGCAGCGTCTTCCGATCT (SEQ ID NO: 30)). The designed oligo pools were purchased from TwistBioscience. To amplify and prepare the template for in vitro transcription (IVT), PCR was performed using primers (fwd: CGTCGGCAGCGTCTTC (SEQ ID NO: 31), rev: GCTAATACGACTCACTATAG (SEQ ID NO: 32)). IVT was then performed using the MEGAshortscript™ T7 Transcription Kit (ThermoFisher Scientific, AM1354). The cleavage reaction was carried out in the following buffer conditions: 10 mM GlcN6P, 50 mM HEPES pH 7.8, 150 mM KCl, 4 mM Mg(OAc)2, 100 μM Na2EDTA, and 1 mM DTT. The reaction was stopped by adding an equal volume of gel loading buffer II (contained in Thermo Fisher Scientific, AM1354). After the cleavage reaction, the RNA was purified using a denaturing 10% polyacrylamide / 7M urea / 1XTBE gel.Reverse transcription (RT) was performed in a final volume of 10 μL using RNA and RT primers ( / 5Phos / AGATCGGAACCGAGCCCACGAGAC (SEQ ID NO: 25) / iSp18 / TCGTCGGCAGCGTCTTCCGATCT (SEQ ID NO: 26), Integrated DNA Technologies) combined with TGIRT-III enzyme (InGex, IGX-TGIRT50-50), SUPERase In (Thermo Fisher Scientific, AM2694), 5 mM DTT (included in Thermo Fisher Scientific, 18090010), and custom RT buffer (50 mM HEPES-KOH, 75 mM KCl, 3 mM MgCl). Because Tris induces cleavage of the glmS ribozyme, the custom buffer contained HEPES instead of Tris. For reverse transcription, samples were incubated at 60 °C for 16 h. The cDNA was purified and circularized, and one-tenth of the circularized cDNA was used directly for barcoding PCR using Platinum SuperFi v2 (Thermo Fisher Scientific, 12361010) and the Nextera XT Index Kit v2 (Illumina, FC-131-2001). The PCR products were directly separated on a 10% polyacrylamide / 1XTBE gel, and the region from 130 bp to 400 bp was excised and purified by dissolving in water. The library was quantified using the Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific, Q32851) and sequenced for 150 PEs on the Illumina MiSeq platform (Core Laboratory, Kyoto University iPSC) using custom read primers (Read 1: TCGTCGGCAGCGTCTTCCGAT (SEQ ID NO: 27), Read 2: GTCTCGTGGGCTCGGTTCCGAT (SEQ ID NO: 28)).

[0054] NGS data analysis for reaction rate determination. FASTQ files of demultiplexed paired-end Illumina reads were merged using PEAR (version 0.9.6) after alignment with BWA (version 0.7.17). Cleaved / uncleaved reads were counted from the fastq files. From the count data, Bayesian regression was used to estimate the cleavage kinetics (kcat). Read counts of cleavage products observed by NGS at each time point t were calculated using NGS. cleaved (t) is the cutting rate R cleaved (t) and depth are binomial distributions of observations, and R cleaved (t) was assumed to follow the simple rate equation:

[0055]

[0056] where Binominal and Exponential represent the binomial and exponential distributions, respectively. We exploited the property that the Michaelis-Menten equation asymptotically approaches equation (6) in the presence of high concentrations of ligand. Using stochastic variational inference, we find that the prior distribution is exponential, K cat , we approximated the posterior distribution of t0. In this study, we set a = 0.1 and b = 1. Bayesian inference was performed using Pyro (version 1.7.0).

[0057] Prediction of variant effects. We manually curated deep mutation scan data of publicly available functional RNAs registered in the Rfam database with over 100 aligned sequences. The dataset was prepared by splitting training data, valid data, and test data in a ratio of 7:1.5:1.5. RfamGen was trained 10 times using the parameters listed in Table 1. The optimal model was selected using the test dataset for prediction. For fitness comparison, the loss function was averaged over five runs and compared with the fitness.

[0058] Results: Our method (RfamGen) is a VAE that leverages the use of the CM architecture to improve the efficiency of RNA family sequence generation (Figure 1A). CM is a context-free grammar that formalizes multiple sequence alignments (MSAs) under the constraint of consensus secondary structure. In this framework, each sequence in an MSA was aligned and represented as a path on the CM (Figure 1B). Because the CM explicitly describes the structural constraints of an MSA, we reasoned that coupling it with the CM architecture would enhance the generative capabilities of deep generative models. Using the VAE architecture, RfamGen is trained to recover alignments that mimic natural sequences while embedding them in a latent space. To achieve this goal, we vectorized the alignment features into a CM and used a VAE to generate paths on the CM instead of generating sequences directly. We disassembled the paths to the CM transition rules and output rules and converted them into one-hot representations, which were then further processed to stabilize the learning process (see "Experimental Methods" above and "Supplementary Note 1" below). Conversely, the CM was reconstructed, and the array was generated by taking the maximum argument (argmax) from the reconstructed CM and sending the array.

[0059] 1. Data Efficiency of RfamGen To examine the effect of explicitly utilizing CM on generative performance, we first compared RfamGen with three ablation models that do not consider information about secondary structure and / or alignment. To isolate the effects of data size and architecture (provided in the source data), we first evaluated 18 RNA families with complete alignments consisting of at least 10,000 sequences in the Rfam database. Each benchmark included GCVAE, which includes only alignments; GVAE, which includes only secondary structure; and CVAE, which includes neither (see Supplementary Note 2). Note that none of GCVAE, GVAE, or CVAE has been used for RNA generation before. For the purpose of comparison with RfamGen, we generated input data tailored to the purposes of GCVAE, GVAE, and CVAE using CM and trained these models (details are provided in Supplementary Note 2). Therefore, RNA generation using GCVAE, GVAE, and CVAE is also achieved for the first time in this study. The quality of the generated sequences was assessed by their likelihood as RNA families ("bit scores") calculated based on alignments with "ground truth" CMs obtained from the Rfam database, a standard approach for detecting RNA family sequences.

[0060] We quantified the average bit scores of 1,000 sequences randomly generated from the model for each RNA family (Figure 2A). The results showed that RfamGen generally performed best, except for RF00163, but RfamGen and GCVAE had superior capabilities for sequence generation, followed by GVAE and CVAE. These results indicate that incorporating alignment information significantly improves performance, while incorporating secondary structure information moderately improves it. Next, we evaluated the effectiveness and robustness of RfamGen and GCVAE for data of various sizes using undersampling. RfamGen maintained high generation capabilities across a wide range of sampling sizes, typically achieving near-peak performance with only approximately 500 input sequences for most RNA families (Figure 2B). The results showed that RfamGen mostly generated sequences with higher bit scores, demonstrating its efficient generation capability even with small datasets containing only a few hundred sequences (Figure 2B). These results make it possible to apply RfamGen to a larger number of RNA families with smaller data sets.

[0061] Finally, we trained RfamGen using 628 RNA families whose complete alignments consisted of at least 100 sequences in the Rfam database and compared it with GCVAE to confirm its broad applicability. We observed that the bit scores of sequences generated from RfamGen were higher than those from GCVAE, demonstrating the superiority of RfamGen in generating RNA family sequences (Figure 2C, p = 3.372*10 -30 , on-sided Wilcoxon signed-rank test).

[0062] Next, we investigated the ability of both models to learn the natural sequence information content. From natural and generated sequences, we calculated sequence entropy (sequence information content) and the mutual information of pairwise occurrence probability on the consensus secondary structure of CMs (secondary structure information content), as described in our previous study. In this analysis, we ignored higher-order information corresponding to pseudoknots, as it has previously been shown to be minor information in natural sequences. RfamGen was slightly less robust in learning the sequence information content of some RNA families compared to GCVAE (Figure 2D, p = 1.63*10). -7 , one-sided Wilcoxon signed-rank test). However, RfamGen significantly outperformed GCVAE in terms of secondary structure information reproduction (Fig. 2E, p = 4.70*10 -31 ( , one-sided Wilcoxon signed-rank test). Similar results were observed using consensus secondary structures obtained by R-scape software, a standard approach for statistically detecting conserved pairwise interactions from MSA (data not shown). In summary, even with a small dataset, RfamGen was able to generate novel sequences with a high probability of being RNA families and to reproduce the statistics of natural sequences.

[0063] 2. Learning a Semantically Rich Latent Space with RfamGen For sequence design, it is important that the latent space contains semantically rich information (Fig. 1A, center). Therefore, we investigated whether RfamGen's latent space contains meaningful sequence representations. A 3D t-SNE projection of the trained RfamGen's 16-dimensional latent space showed separated clusters of phylogenetic information (Fig. 3A). Similar trends were observed for other RNA families (data not shown). Furthermore, we investigated how RfamGen embeds sequence variants into the latent space using Deep Mutational Scanning (DMS) of tRNAs, as reported by Li, C., Qian, W., Maclean, CJ & Zhang, J. The fitness landscape of a tRNA gene. Science 352, 837-840 (2016). We performed linear regression to fit the activity landscapes in the latent space and found a moderate correlation between them (Figure 3B, Pearson's r = 0.2932). This result suggested that the latent space linearly represents tRNA activity. Another example, the twister sister ribozyme, showed a smooth continuum between phylogenetic and Euclidean distances in the latent space (Figure 3C, Spearman's r = 0.7404), separating clusters of structural variants with and without variable stem-loops (Figure 3D). Finally, we performed a computational search for sequences with the U1A protein-binding motif in a latent space trained on the natural type III hammerhead ribozyme (RF00008), but found no sequences without the binding motif. RfamGen successfully generated sequences with the U1A protein-binding motif and high bit scores, such as RF00008. tSNE embedding of the latent variables revealed organized clusters of sequences with the binding motif, demonstrating the utility of the latent space for discovering sequences with unnatural properties (data not shown). Overall, these results demonstrate that RfamGen can learn meaningful representations, which may be useful for designing sequences with desirable properties.

[0064] 3. Generation of Active Sequences of Diverse RNA Families with RfamGen In RNA engineering, self-cleaving ribozymes are one of the most extensively studied functional RNAs. We evaluated the applicability of RfamGen to generate functional sequences using a relatively small number of self-cleaving ribozymes (approximately 100–800) from the ZWD and Rfam databases. The sequences generated from RfamGen exhibited similar sequence and structural profiles to the native sequences (Figures 4A, 4B, and 4C), but differed in sequence identity (Figures 4D, 4E, and 4F). RfamGen restricts the evolutionarily conserved secondary structure regions of ribozymes, but can generate RNA families with high sequence diversity. To evaluate the functionality of the generated sequences, we randomly synthesized three generated sequences for each ribozyme and designed corresponding mutants that reportedly have reduced activity alongside the native sequences. The cleavage activity of the sequences was confirmed by electrophoresis (Figures 4G, 4H, and 4I). Unlike the sequences from RfamGen, sequences generated by random sampling from CM were found to be inactive (data not shown). Random sampling from CM theoretically had similar sequence / secondary structure priorities, but lacked higher structural constraints and many-body correlations such as base pair stacking. Thus, our results support the ability of RfamGen to learn many-body correlations from data. Collectively, these results demonstrated the broad applicability and high generative power of RfamGen.

[0065] 4. Predicting Mutation Effects in Diverse RNA Families with RfamGen RNA engineering often requires the design of sequences with tailored activity. To this end, high-throughput experiments have been conducted to screen mutants with optimal activity for engineering applications. Meanwhile, in protein engineering, it is empirically known that activity can be predicted by the generation probability modeled by a generative model trained on natural sequences. Several previous studies have made similar observations for tRNAs. However, it remains unclear whether the same approach can be generally applied to various RNA families. Therefore, we extended this analysis to more RNA families by manually curating the activities of sequence variants from previous DMS studies. Next, we investigated whether the activities of sequence variants can be predicted using the same strategy as for proteins.

[0066] Using a carefully selected dataset, we compared the activity of reported sequence variants with RfamGen's loss function. Because the VAE loss function is a lower bound on the log-probability ln p(x), we approximated the log-probability ln p(x) using RfamGen's loss function. We also evaluated EVMutation, a Potts model-based predictor of mutation effects (Table 3). Through this comparison, we confirmed that RfamGen predicts variant effects with comparable accuracy to EVmutation. These results support the potential use of RfamGen for sequence tuning.

[0067] *a: Andreasson, JO, Savinov, A., Block, SM & Greenleaf, WJ Comprehensive sequence-to-function mapping of cofactor-dependent RNA catalysis in the glmS ribozyme. Nat. Commun. 11, 1663 (2020). *b: Kobori, S. & Yokobayashi, Y. Analyzing and Tuning Ribozyme Activity by Deep Sequencing To Modulate Gene Expression Level in Mammalian Cells. ACS Synth. Biol. 7, 371-376 (2018). *c: Kobori, S. & Yokobayashi, Y. High-Throughput Mutational Analysis of a Twister Ribozyme. Angewandte Chemie Int. Ed. 55, 10354-10357 (2016).

[0068] 5. Massively Parallel Assay of Generated Ribozymes Massively parallel assays are a method for searching for sequences with desired functions using randomized libraries of functional RNA sequences. However, this approach involves a significant trade-off between mutational diversity and activity due to the stochastic loss of activity caused by random mutations. Furthermore, because this approach relies heavily on empirical sequence-activity correlations, sequence modifications are limited to RNAs with many known sequences. Similarly, one of the major technical challenges is the efficient in silico design of individual functional RNA sequences.

[0069] As an example of a self-cleaving ribozyme, we focused on the glucosamine-6-phosphate riboswitch ribozyme (glmS ribozyme), the only known ligand-dependent self-cleaving ribozyme in nature. This ribozyme is approximately 200 nucleotides long, contains a double pseudoknot structure in its core region, and controls self-cleavage by responding to a specific metabolite, glucosamine-6-phosphate (GlcN6P). Due to its unique nature, the glmS ribozyme is a natural model for gene expression regulators in RNA engineering (Andreasson et al.). Although the sequence characteristics of the glmS ribozyme have been studied for many years, our overall understanding of the relationship between sequence and activity remains limited, and the design of artificial sequences distinct from the natural sequence has yet to be achieved. We aimed to test whether RfamGen could address this technical challenge by verifying the functionality of in vitro-generated sequences.

[0070] First, we trained RfamGen with approximately 500 glmS ribozyme sequences and confirmed the ligand-dependent cleavage activity of all three randomly generated sequences by electrophoresis (Figure 4G). Next, we aimed to comprehensively evaluate the sequences using a massively parallel assay. We generated 1,000 sequences by sampling the RfamGen latent space without any additional selection process. In addition to the generated 1,000 sequences, we selected 100 sequences with reaction rates reported by Andreasson et al. as spike-ins, and added 761 sequences of medium length from the Rfam database, Rfam accession RF00234. For experimental validation, we established a high-resolution massively parallel assay (Figure 5A). This assay stably measured the kcat parameter over a dynamic range of approximately three orders of magnitude, even for sequences where kcat was difficult to measure due to the time resolution of the assay (data not shown).

[0071] Furthermore, we confirmed the reproducibility of the assay against external data from Andreasson et al. (data not shown). Similar kcat values ​​of the picked sequences were also observed in the gel-based assay (data not shown). These results enabled us to quantitatively compare and discuss the kcat values ​​of the sequences. Interestingly, we found that the kcat distribution of the generated sequences was similar to that of the natural sequences, and furthermore, we observed that the generated sequences exhibited higher kcat values ​​compared to the naturally occurring glmS ribozyme (Figure 5C, p = 5.223*10). -45 , Mann-Whitney U test). This indicates that RfamGen can generate functional RNA sequences with high efficiency.

[0072] 6. Results and Discussion We integrated the CM architecture into RfamGen to enhance its generative capabilities. CM only considers nested base pairs, which are known to constitute an important part of the information content of RNA, improving data efficiency. Furthermore, the deep generative model architecture allows RfamGen to compensate for higher-order information that cannot be explained by CM alone, further improving RfamGen's performance. This was indirectly supported by the lack of information from random sampling from CM (data not shown). Similar results were also confirmed by examining generated sequences randomly sampled from CM registered in the Rfam database (data not shown), confirming the effectiveness of argmax sampling. Indeed, although CM cannot consider pseudoknots, functional ribozymes with pseudoknots were successfully generated by RfamGen (Figures 4A-I and 5A-C). These results indicate that RfamGen implicitly learns higher-order information from natural sequences. Finally, we found a match between RfamGen's loss function and the performance of activity prediction from DMS data (data not shown). This suggests the usefulness of DMS data as external validation. Further accumulation of DMS data for functional RNAs may further accelerate the improvement of the RfamGen architecture.

[0073] Limitations of RfamGen may include data bias in RNA families, incorrect annotation of consensus secondary structures, and nested CMs. Specifically, the high conservation of related RNA families in the data space makes some parameters, such as θ and β, sensitive to data diversity, making training RfamGen challenging. For example, the diversity limitation in RF00163 could be strengthened by considering base pairs as a single object during processing in RfamGen, which could contribute to stronger regularization against reconstruction error and reduced quality of the generated sequences (Figures 2A and 2B). This effect could be mitigated by reducing θ and β when reweighting Equation (2) and optimizing the loss function (Equation (3)) or other parameters. Furthermore, extending the flexibility of RfamGen to simultaneously train all RNA families and transfer information between RNA families, for example through large-scale language models, could improve overall performance and broaden the scope of RNA design beyond existing RNA families.

[0074] Using gel electrophoresis and massively parallel assays, we demonstrated RfamGen's high design success rate for diverse RNA families. Surprisingly, we generated sequences with high activity that is rare among natural sequences (Figure 5C). To further explore this observation in other RNA families, we compared the predicted activity of natural and generated sequences across 628 RNA families. The computational results also suggested slightly lower loss values ​​for the generated sequences, indicating their high activity (data not shown). This is likely due to the noise reduction effect of VAE, which absorbs biological sequence variation as statistical noise. This finding indicates that naturally occurring RNAs may leave evolved sequence space to further enhance specific functions, and that RfamGen can be used to promote or optimize the activity of functional RNAs. As previously discussed for protein design, RfamGen's ability to predict variant effects and its semantically rich latent space are useful for functional RNA design. Taken together, our computational and experimental results can provide a foundation for deep generative RNA engineering.

[0075] The applicability of RfamGen extends beyond the natural functional RNAs in the Rfam database. In principle, RfamGen can be easily applied to any type of functional RNA, provided there is a consensus secondary structure for constructing the MSA and CM. As a demonstration, we constructed CMs from the consensus secondary structure of the MSA and ZWD pistol ribozyme, subsequently training RfamGen and verifying the functionality of the generated sequences. Interesting applications of RfamGen to genetic engineering and RNA medicine could include the design of tRNAs, ribosomal RNAs, or synthetic RNAs obtained through in vitro selection or experimental evolution. Indeed, we demonstrated that exploring the potential space identified sequences similar to the protein-responsive riboswitches experimentally discovered by Auslander, S. et al., A general design strategy for protein-responsive riboswitches in mammalian cells. Nat. Methods 11, 1154-1160 (2014). This is one of the most studied types of artificial RNAs. These results suggest the utility of RfamGen in RNA engineering to discover functional RNAs with additional properties. In addition to functional RNAs, another promising application of RfamGen is the discovery and optimization of RNA aptamers that specifically interact with various ligands, and their potential for medical applications has attracted considerable attention. RfamGen's high data efficiency and high design success rate make it a powerful and versatile platform for RNA engineering in biotechnology and medicine.

[0076] [Supplement 1] Details of Pre- and Post-Processing: RfamGen data processing is described in more detail. RfamGen can be considered an extension of grammar VAE to covariance modeling (CM). Within the framework of grammar VAE (Wilson, DS & Szostak, JW. In vitro selection of functional nucleic acids. Annu. Rev. Biochem. 68, 611-647 (1999)), Kusner et al.'s parse tree was converted to a one-hot representation. Similarly, the traceback information of alignments to CM was converted to a one-hot representation and further decomposed into three matrices X = (tr, ss, bp). tr, ss, and bp contain information on the shape of secondary structures, single-stranded nucleotides, and base-paired nucleotides, respectively. Conversely, CM was reconstructed from X = (tr, ss, bp), and an argmax sequence was generated from the reconstructed CM. RfamGen's pre- and post-processing are described in the following three subsections. (1) One-hot representation of alignment to CM (2) Splitting the one-hot representation into three matrices (3) Decoding the alignment and generating the sequence

[0077] (1) One-hot representation of alignment to CM The covariance model is a kind of context-free grammar. As shown in Figure 6(1), all transition / output rules are called a "CM grammar." The CM grammar contains 76 grammars consisting of the following 56 transition rules, 4 single-strand output rules, and 16 base pair output rules.

[0078] In this grammar, the base pair output n is used so that RfamGen can learn base pair information directly. l , n r The alignment to the CM was performed by the "cmalign" command in the Infernal software. As shown in Figure 6 (2), the alignment on the CM is the path from the S state to the E state. Next, as shown in Figure 6 (3) and (4), the jth non-terminal state jThe i-th rule r of the CM grammar in i We created a binary matrix traceback(i, j) by ordering the use of varies depending on the consensus structure of the CM.

[0079] (2) Splitting the One-Hot Representation into Three Matrices In this section, we explain how to compress the traceback information into three matrices. As shown in Figure 7, we found that this traceback(i, j) contains many zero or missing values, which cause training instability. For example, there can be no single-strand output in the state column {S,MP,D,B}. To mitigate training instability and ignore these meaningless values, we defined reindexing to move values ​​to the left side of the matrix (Figure 7, left).

[0080] For single-strand and base pair outputs, the problem was solved by defining the reindexing functions σ and β as follows: where I A are the indicator functions of set A. The role of these functions is to ignore and remove information on single strands or states that do not output base pairs. Next, to reduce missing values, traceback(i, j) was converted into a triplet of three matrices (tr, ss, bp) as follows:

[0081] Here, trans, ss, and bp correspond to the transition, single-strand output, and base pair output. Since trans was not processed, there were still many missing values. Therefore, we further utilized the node information of the covariance model (Figure 8). A node is a coarse-grained block of states given internally by the Infernal software. Since a single node does not contain states of the same type, the intra-node transition rules are not used between different states. Therefore, the transition information of states within a node is summarized into a single vector, generating a column of transition information. Formally, a node is a set of states, and every node has a λ We assume that the nodes are indexed by s.t. We used the order of the state index s.t to index the nodes. Using this node indexing, we were able to define a reindexing function from the trans indices to the reformatted matrix, as follows: Using reindexing, I was able to convert trans to tr without losing any information, like this:

[0082] Columns containing only missing values ​​are dropped, and the remaining missing values ​​are filled with zeros. These three matrices were used in training. To reconstruct the logits from the decoder output (tr', ss', bp'), the above preprocessing steps were completely reversed. This is achieved by "memorizing" the jth non-terminal list.

[0083] (3) Decoding the alignment and generating a sequence The output from the decoder is (tr', ss', bp'), which is reformatted into logits as described above (Fig. 9(1) and (2)). If CM has no value, we mask it (Fig. 9(3)), and replace the CM value with the remaining value of the logits (Fig. 9(4)). To generate a sequence from the reconstructed alignment, we use the ith rule r i from argmax i+1 p(r i+1 ; r i ) to collect the most probable rules of the CM grammar (Fig. 9(5)) and collect the terminal symbols (Fig. 9(6)).

[0084] [Supplement 2] Description of benchmark models To evaluate the effectiveness of incorporating alignment and structural information, we constructed benchmark models using ablation information, character VAE (CVAE), grammar VAE (GVAE), and gapped-character VAE (GCVAE). The information that each model can consider is shown in Table 4 below.

[0085] Character VAE (CVAE) The CVAE was trained using one-hot representations of sequences. Exceptionally short or long sequences whose number was less than 1% of the total sequence size were excluded. The "AUCG" in the RNA sequence was converted to a one-hot representation. The architecture of the CVAE is similar to that of RfamGen. The encoder converts the one-hot representation X to a latent variable Z using three one-dimensional convolutional layers followed by two fully connected layers. The decoder recovers triplets X' from the latent variable Z using two fully connected layers followed by three one-dimensional deconvolutional layers. The convolutional / deconvolutional parameters, hidden dimension, and latent dimension are the same as those of RfamGen. To sample sequences, we sampled the latent variables, decoded them into one-hot logits, and obtained an argmax array from the logits.

[0086] Grammar VAE (GVAE) GVAE is a VAE using a context-free grammar (CFG) proposed by Kusner et al. GVAE has previously been applied to small molecule compounds but has not yet been used for RNA design. However, because it was necessary for the validation of RfamGen, we performed the first RNA design incorporating structure using GVAE. The preprocessing procedure was essentially the same as that of Kusner et al., except that the grammar was specialized for RNA structure. More specifically, the sequences were preprocessed as follows: (1) Align the sequences to the CM and assign a consensus structure to each sequence; (2) Remove gaps between the corresponding structures; (3) Analyze the corresponding structures using the G6 grammar; and (4) Create a one-hot representation from the analysis tree. This preprocessing method was also the first to be investigated and implemented by the inventors for the purpose of RNA generation. We chose the G6 grammar because it can analyze a larger number of sequences than other grammars. Sequences that could not be analyzed using the G6 grammar were excluded. This time, we performed structural analysis using G6 grammar, but any grammar that can analyze bracket notation, including G1 and G2, can also be used.

[0087] The architecture of GVAE is similar to that of RfamGen. The encoder converts the one-hot representation X into a latent variable Z using three one-dimensional convolutional layers followed by two fully connected layers, and the decoder recovers the triplet X. The latent variable Z is then generated using two fully connected layers followed by three one-dimensional deconvolutional layers. The convolutional / deconvolutional parameters, hidden dimensions, and latent dimensions are the same as those of RfamGen. To sample the array, we sampled the latent variables and decoded them into logits of the one-hot representation of the parse tree. To generate the array, we stepwise selected the G6 grammar rule with the highest probability from among the rules that matched the last previously selected rule. Details are disclosed in Kusner et al.

[0088] Gapped-Character VAE (GCVAE) GCVAE is a modified CVAE. The sequences were preprocessed as follows. This preprocessing method was also the first we have considered and implemented for the purpose of RNA generation. (1) The sequences were aligned to CM and gap characters in the insertion state were removed. (2) "AUCG-" in the aligned RNA sequence was converted to a one-hot representation. The architecture of GCVAE is similar to that of RfamGen. The encoder converts the one-hot representation X into a latent variable Z using three one-dimensional convolutional layers followed by two fully connected layers, and the decoder recovers the triplet X'. The latent variable X is generated using two fully connected layers followed by three one-dimensional deconvolutional layers. The convolutional / deconvolutional parameters, hidden dimension, and latent dimension are the same as those of RfamGen. To sample the sequences, the latent variables were sampled and decoded into logits of the one-hot representation of the parse tree, and an argmax array was obtained from the logits.

[0089] Comparison of Benchmark Models: In Figure 2 of the main text, we first carefully selected 18 RNA families with complete alignments consisting of at least 10,000 sequences. From these RNA families, we excluded RNA families with too short base pairs (less than 8 bp). To compare RNA Family RfamGen with the benchmark models, we prepared sequences suitable for all models, especially CVAE and GVAE. That is, we prepared sequences that did not contain ambiguous nucleotides such as N / Y / R, and were neither too long nor too short (cutoff length <1% frequency for CVAE), and their structures could be analyzed using the G6 grammar (for GVAE).

Claims

1. A method for generating functional RNA sequences, comprising: a) obtaining input data X that digitizes secondary structure information stored in a training RNA sequence; b) inputting the input data X for a plurality of training RNA sequences into a training model to obtain generated output data X'; and c) presenting the generated plurality of RNA sequences based on the output data X'.

2. The method for generating functional RNA sequences according to claim 1, wherein the step a) is a step of obtaining input data X in which alignment information of the learning RNA sequence is quantified in addition to secondary structure information of the learning RNA sequence.

3. A method for generating a functional RNA sequence, wherein the step a) is a step of obtaining the digitized input data X based on a covariance model (CM).

4. The method according to claim 3, wherein the step a) includes the steps of: aligning the training RNA sequence to a CM and converting a path on the CM into a numerical matrix representation; and separating non-terminal symbol parts, base pair parts, and non-base pair parts from the matrix of the numerical matrix representation, thereby obtaining the digitized input data X, wherein the numerical matrix representation is a matrix with CM state indexes as columns and transition / output rules used as rows, in which 1 is assigned to the rule that was used and 0 is assigned to the others, and in the case where multiple insertion transition rules are used, the number of times the insertion transition rule was used is assigned.

5. The method of claim 3, wherein the input data X is represented as X=(tr, ss, bp), where tr represents non-terminal transition information, ss represents single strand output information, and bp represents base pair output information.

6. The method according to claim 3, wherein the step b) includes a step of transforming the plurality of input data X into a latent variable Z using a first nonlinear function, and obtaining output data X' from the latent variable Z using a second nonlinear function, wherein the first and second nonlinear functions are parameterized, and the parameters are learned so as to minimize a first term that brings the input data X and the output data X' closer together, plus a second term that brings the distribution of the latent variable Z obtained by transforming the input data X closer to a prior distribution.

7. The method of claim 3, wherein the learning model performs machine learning of the input data X using a variational autoencoder.

8. The method of claim 6, wherein the first term represents the cross-entropy error between the input data X and the output data X', and the second term represents a weighted term of the Kullback-Leibler distance between the distribution of Z and a normal distribution with mean 0 and variance 1.

9. The method of claim 6, wherein the first nonlinear function is a first neural network and the second nonlinear function is a second neural network.

10. The method according to claim 3, wherein the step c) comprises: creating an output CM by replacing parameters of the CM in the step a) based on the output data X'; and obtaining an RNA sequence having a maximum probability from the output CM.

11. The method of any one of claims 1 to 10, wherein the produced functional RNA comprises linear or circular, mRNA or non-coding RNA.

12. A method according to any one of claims 1 to 10, wherein the generated functional RNA exhibits activity equivalent to that of the learning RNA sequence in vitro, in vivo, ex vivo, and / or in cellulo (e.g., in a test tube including in an aqueous solution or on a paper surface, in a liposome, in a cell, in an organoid, in a tissue, and / or in a living organism).

13. The method according to any one of claims 1 to 10, wherein the functionality includes an RNA activity selected from an RNA activity corresponding to an RNA family registered in the Rfam Database, an RNA activity corresponding to all RNA families known from literature, etc., and any other RNA activity.

14. The method of any one of claims 1 to 10, wherein the functionality comprises riboswitch activity, RNA interference activity, aptamer activity, tRNA activity, signal recognition activity, snRNA activity, snoRNA activity, guide RNA activity, Y RNA activity, long non-coding activity, eRNA activity, microRNA activity, piRNA activity, siRNA activity, rRNA activity, genome editing activity, RNA editing activity, transposon activity, ribosome activity, spliceosome activity, and / or ribozyme activity (including, for example, self-cleaving ribozyme activity and RNA methylation ribozyme activity).

15. The method of claim 13 or 14, wherein step a) comprises obtaining input parameters X for a plurality of training RNA sequences having a common sequence or structure.

16. The method according to claim 6, wherein in step b), in a coordinate space visualized using t-SNE of the latent space represented by the latent variable Z, the relative distance of variables into which the multiple pieces of training RNA sequence information are transformed represents the similarity in characteristics between each of the training RNA sequences.

17. A method for generating a functional RNA sequence, comprising: A) obtaining a new latent variable Z' for a trained latent variable Z, converting the latent variable Z', and obtaining output data X'; and B) presenting a plurality of generated RNA sequences based on the output data X', wherein the latent variable Z has information on the sequence, structure, biological species, and / or activity intensity of the input data X based on the trained RNA sequence, and the latent variable Z and the new latent variable Z' have similar information.

18. A computer program for causing an information processing device having at least one processor to execute each step of the method according to any one of claims 1 to 17.

19. A non-transitory computer readable medium comprising the computer program of claim 18.

20. A method according to any one of claims 1 to 17, executed by an information processing device having at least one processor.

21. An information processing device comprising at least one processor, the processor executing a method according to any one of claims 1 to 17.

Images