RNAStructuromeDB: A genome-wide database for RNA structural inference

RNA plays important roles in almost every aspect of biology, and every aspect of RNA biology is influenced by its folding. This is a particularly important consideration in the era of high-throughput sequencing, when the discovery of novel transcripts far outpaces our knowledge of their functions. To gain a comprehensive picture of biology requires a structural framework for making functional inferences on RNA. To this end we have developed the RNA Structurome Database (https://structurome.bb.iastate.edu), a comprehensive repository of RNA secondary structural information that spans the entire human genome. Here, we compile folding information for every base pair of the genome that may be transcribed: coding, noncoding, and intergenic regions, as well as repetitive elements, telomeres, etc. This was done by fragmenting the GRCh38 reference genome into 154,414,320 overlapping sequence fragments and, for each fragment, calculating a set of metrics based on the sequence’s folding properties. These data will facilitate a wide array of investigations: e.g. discovery of structured regulatory elements in differential gene expression data or noncoding RNA discovery, as well as allow genome-scale analyses of RNA folding.

families (Rfam) database, for example, contains entries for families of ncRNAs linked by homology. Rfam entries contains information describing ncRNA biosynthesis, localization, phylogenetic distribution and functional roles, as well as evolutionary conservation of primary sequence and, importantly, secondary structure. Conservation of secondary structure is a defining feature of ncRNAs and is used in essentially all ncRNA prediction programs 36 .
A powerful, and popular, program for genome-wide ncRNA prediction is RNAz [37][38][39][40] . This program is based on a support vector machine (SVM) that is trained on data from known ncRNAs contained within Rfam. Two primary training parameters are used for ncRNA classification: a structure conservation index (SCI), which measures conservation of secondary structure and a thermodynamic z-score, which measures the propensity of a particular sequence to form a defined and energetically stable structure. Typical ncRNAs have structures with high conservation and propensity to form structure. RNAz was previously used to scan through whole genome alignments comparing human vs. animal genomes (ranging in similarity down to zebrafish) to identify putative ncRNAs. This yielded >30,000 high-confidence predictions with ~1,000 that were conserved throughout all vertebrates 41 . The identification of so many deeply conserved structured RNAs highlights their likely ubiquity and importance.
In both coding and noncoding RNAs, secondary structure plays key roles throughout their functions. The diversity of RNA functions and potential for physiological impact (much like with proteins) is made possible by the ability of RNA to fold into unique functional structures. Functional RNA structures form thermodynamically stable base pairs that have been selected for by evolution 42 . This is the key principle behind the thermodynamic z-score implemented in RNAz: functional RNAs have a more stable folding energy than randomized sequences. Specific folds can be recognized by regulatory proteins 43 , occlude/present functional motifs 44 , or alter the distance between functional sites 45 .The impact of RNA structure on alternative splicing, for example, has been particularly well-studied 46,47 . Additionally, awareness of the importance of non-specific RNA structure (e.g. regions that do not form particular folds or adopt dynamic structures) is growing. For example, thermodynamically stable regions within open reading frames are proposed to modulate the speed of translation and thus affect protein folding 48,49 . RNA intramolecular thermodynamic stability is important in mediating its accessibility for intermolecular interactions: e.g. stable regions in UTRs are less accessible to miRNA binding, thus affecting miRNA-mediated gene silencing 50 . Dynamic RNA structure also has significance to disease: single nucleotide polymorphisms (SNPs) can affect RNA folding in ways that impede healthy function by disrupting specific motifs or altering conformational equilibria 51,52 .
Advances in sequencing technology will continue to massively expand the list of interesting RNA sequences and, excitingly, also provide information on secondary structure. Several in-vivo RNA structure probing techniques have been developed to acquire transcriptome-wide folding information 53 . Snapshots of the human "RNA structurome" (ranging across tissue/cell types, disease states, and treatments with drugs) will become more common in the near future; informing our knowledge of human biology and advancing our understanding of pathogenesis. For this reason, and the reasons discussed above, it is critical to have a knowledge framework in place to understand the roles of RNA structure in human biology. This is the motivation for the creation of the RNAStructuromeDB. Here, we have compiled computed RNA folding information across the entire human genome, irrespective of whether or not it is known to be transcribed. The RNAStructuromeDB is a web-accessible (https://structurome.bb.iastate.edu) repository for investigators to obtain structural metrics for any RNA sequence originating from the human genome. To further aid investigators we have put the data into context by incorporating comprehensive Gencode annotations 54 using the biological database schema Chado 55 . This allows the rapid comparison of differential gene expression data (e.g. to identify regulatory RNA structures) or transcriptome-wide RNA biochemical probing data against the RNAStructuromeDB facilitating these, and other types of studies.

Results and Discussion
The RNAStructuromeDB holds the results of a genome-wide computational analysis in which we folded the entire human genome. The results of this analysis are comprised of folding metrics which indicate every region of the genome's propensity to generate structured RNA molecules. Tools have been built by which users can quickly download these metrics or use them to search for novel RNA structural elements or to assess the global folding properties of a transcript. Here we present the metrics that will aid in the investigation of RNA structure/function, and the tools which have been developed to view these metrics effectively.

Data types.
To analyze a genome as large as the human genome requires fragmenting the sequence data into smaller pieces that are both computationally tractable, as well as biologically meaningful. The fragmentation approach taken in the construction of the RNAStructuromeDB is similar to those that were successfully used for the analysis of RNA structure in the human and other genomes [39][40][41][56][57][58] as well as in specific lncRNAs 59 . Here we consider discrete, overlapping windows whose size was selected based on several considerations. Due to the kinetics of RNA folding, RNAs typically fold into structures composed of smaller motifs; for example, lncR-NAs like HOTAIR (2,421 nt) are comprised of locally folding domains (4 domains), each containing many local substructures (e.g. RNA hairpins 60 ). Additionally, the majority of known RNA sequences/structures archived in Rfam are short (<200 nt long). Incidentally, a window size spanning 100 to 150 nt was shown to be most accurate at predicting base pairs in known cis-regulatory structures of large mRNA molecules (>85,000 nts long) 61 . Finally, folding algorithms scale O(N 3 ) in time and O(N 2 ) in memory, where N is the sequence length 62 ; thus, shorter sequences, in addition to yielding more accurate predictions, can be folded more quickly and efficiently. With these factors in mind, a window size of 120 nt was selected to maximize the chance of encapsulating structured elements, while also optimizing prediction accuracy and computational resources. A step size of 40 nt was selected to ensure the entire genome is canvassed with multiple frames, increasing the likelihood of capturing all Scientific REPORtS | 7: 17269 | DOI:10.1038/s41598-017-17510-y structured elements. This window and step size resulted in the generation of 154,414,320 windows which were analyzed in both strand orientations.
For each window, five folding metrics were calculated and archived. The Gibb's minimum free energy (MFE) of folding (ΔG), which estimates the thermodynamic stability of the most stable 2D conformation of an RNA given Turner nearest neighbor energy parameters 63,64 (a set of experimentally measured values taken from the analysis of many small RNA motifs). The ΔG was calculated using the program RNAfold, which is a component of the ViennaRNA package 63 ; predictions were made at 37 °C (human body temperature) and values are reported in kcal/mol. 2D structures associated with window MFE predictions were captured and archived in "dot-bracket" notation; here, paired nt are represented by matched brackets "()" and single stranded nt by dots ". ". To determine if the MFE ΔG depends on the nt order or composition, we compared the native sequence ΔG vs. in silico randomized sequences with the same nt composition. In general, structured ncRNA molecules have lower (more stable) MFE values than random sequences of the same nucleotide composition 65 ; a property which can be exploited for structured RNA motif discovery 66 . For each window, we report the thermodynamic z-score. The z-score is calculated by taking the difference between native MFE (ΔG native ) and the mean MFE of random sequences (ΔG ) random , then normalizing by the standard deviation, as described in equation (1) in the Materials & Methods section. The z-score sign indicates if the native MFE is either lower (negative) or higher (positive) than ΔG random and its magnitude indicates the standard deviations ΔG native is from ΔG random . A z-score less than −1, for example, indicates a ΔG native which is one standard deviation lower than random. Negative z-score suggests that a sequence's order was selected by evolution to fold into a stable structure 42 : e.g. if the order is important, shuffling nt will lead to less stable folds by disrupting native (evolved) pairing contacts that give the RNA thermodynamic stability. Positive z-scores are harder to interpret, however, they could suggest regions of RNAs that are evolved to be more accessible/unstructured 59 . The z-score can also be qualified using the p-value, which is calculated as the fraction of ΔG random values with more thermodynamic stability than ΔG native . This provides a measure of the z-score quality, where p-values close to zero indicate higher prediction confidence: few or no randomized sequences are more stable than the native sequence; this also indicates if the randomization number is high enough to yield reliable z-scores.
Two values are recorded which measure statistical properties of the RNA 2D folding ensemble: the ensemble diversity (ED) and the frequency of the MFE (fMFE) metrics. Structured RNAs have rough folding landscapes (with many "suboptimal" folds that are near in energy to the native MFE fold) and, in some cases, are capable of occupying several conformations. The probably of finding any of these structures within the Boltzmann ensemble is mathematically described by the partition function 67 , which was calculated using RNAfold 63 . From this partition the ED is calculated by taking the average base-pair "distance" between all structures in the ensemble, where distance is defined as the number of base pairs different between structures 65 . The ED then, is a metric to mathematically signify the variety of folding structures within the ensemble (a low ED implies a small number of similar structures are present, while a high ED suggests several alternative folds or a lack of overall structure in the ensemble 68 ). The fMFE metric is the probability of finding the MFE within the Boltzmann distribution of structures, where a high probability suggests the MFE structure is more likely to be the dominant fold.
The five archived metrics each suggest regions of the human genome that may generate RNAs with interesting folding properties. The MFE ΔG defines the thermodynamic stability, while the z-score suggests whether or not that stability is unusual (given the sequence composition). A region can have a very stable ΔG, but mediocre (or positive) z-scores because the order is unimportant vs., for example, the GC content. The p-value estimates the quality of the z-score. If almost every randomized sequence is less stable than native, then the likelihood of the order being significant is higher. The ED and fMFE indicate the diversity of the RNA folding ensemble in a region and how well-represented the MFE prediction is in the 2D structural ensemble, respectively. These can suggest which regions have well-defined folds (e.g. with low ED and high fMFE) or where folding may be diverse (e.g. dynamic regions or regulatory structural switches that have higher ED and low fMFE). Taken together, these metrics can also suggest which regions may be harboring functional structures. Overlapping windows with interesting folding metrics can be concatenated to define larger regions of interest, which can then be used for comparative sequence/structure modeling 69 or ncRNA prediction (e.g. using approaches such as RNAz [37][38][39][40] . Examples of how the metrics contained in the RNAStructuromeDB can be used are discussed below in the "Examples" subsections. Data accessibility. There are three main methods for accessing the data: flat files, data tables, and via an interactive genome browser. Users can freely download raw data as flat files in generic feature format version three (GFF3) directly from the download page (https://structurome.bb.iastate.edu/downloads). GFF3 is a standardized file format developed by The Sequence Ontology 70 for the sharing of genomic information in a unified manner and is compatible with several bioinformatics programs. This file contains each window's sequence, coordinates, strand orientation (forward or reverse), MFE, z-score, p-value, ED, fMFE and dot-bracket 2D structure. This data format will allow interested end-users to parse out trends using their preferred methods or protocols; however, users may want to filter the metrics, downloading only information relevant to their needs/interests. Therefore, we have also created an interface where a user can filter and search for data based on any combination of genomic coordinates, Ensembl IDs, gene symbols, or folding metrics-facilitating varied inquiries. Investigators interested in multiple genes can quickly search for, populate, and download folding metrics coinciding with each gene as a single consolidated CSV (comma-separated values) file. Conversely, a user may have discovered an RNA molecule lacking any established annotations; as long as coordinates are supplied, this user will find folding metrics for their novel transcript as well. This genome-wide analysis provides the freedom to obtain data from any known or novel transcript.
Scientific REPORtS | 7: 17269 | DOI:10.1038/s41598-017-17510-y The JBrowse genome browser provides significant insights into the structural landscape of each chromosome via a customizable visual interface. Within the JBrowse interface, folding metrics have been adapted to allow for rapid structural inferences to be made: they are displayed as bar graphs canvasing each chromosome sequence. Each folding metric (MFE, z-score, p-value, ED, and fMFE) has been split into separate graphical tracks that can be displayed parallel to genomic annotations. A user interested in a specific transcript will be able to observe the folding metrics of its constitutive features (exons, introns, and UTRs are all displayed with unique visualizations), while at the same time be able to see the folding metrics of the entire gene within greater genomic context. The human eye is very good at identifying patterns (e.g. within the structure of the data represented in the JBrowse tracks) and visualizing the genome, transcriptome and predicted RNA structurome simultaneously can facilitate discovery. Examples are given below.
Example 1: The MALAT1 lncRNA. Currently, thousands of lncRNAs have been discovered in the human transcriptome and a growing list has been assigned functions and are implicated in diseases such as cancer [71][72][73] ; however, most lncRNAs and their functions remain uncharacterized 73 . One aspect of lncRNA function that remains controversial, is the role played by RNA secondary structure. Many reports describe extensive local or global folding for lncRNAs 59,60,74 ; yet, statistical evidence of structure conservation remains tenuous 75 . It is worth noting that, even in the absence of high structural conservation, RNA folding can play functional roles 76 . Analyzing lncRNA folding landscapes and building 2D structural models of these molecules will aid in understanding their mechanisms of action (e.g. in detecting regulatory structural motifs) and, possibly, in developing therapeutic strategies 77 to modulate function. lncRNAs tend be modular in structure, consisting of multiple structured domains 78 . Sliding window approaches for RNA structure detection are able to roughly define the extent of these domains 59 . The RNAStructuromeDB suggests these domains for all potential human lncRNAs with the folding metrics and local structural models necessary to begin to decipher the structure and function of these transcripts.
For example, MALAT1 (metastasis-associated lung adenocarcinoma transcript 1) is a highly conserved (throughout 20 mammalian species, including mouse and human 79 ) lncRNA, which is involved in numerous cellular processes (e.g. transcriptional regulation 80 , alternative splicing 81 , and cellular localization 82 ) and implicated in disease states: such as cancer 83 and diabetes 84 . Here, the RNAStructuromeDB data table interface is used to define the structured domains of MALAT1 (summarized in Fig. 1). Upon inputting the MALAT1 target (specified by genomic coordinates, gene symbol, or Ensembl ID), all overlapping folding metric windows are extracted from the DB. These can be browsed on the website or downloaded as a single CSV file. With these data, investigators can use any method to define structured regions; in this example, regions with overlapping windows having z-scores 1 σ more negative (< −1.70) than the average MALAT1 z-score (−0.51) were concatenated into domains (similar to previous work on the Xist lncRNA 59 ). This resulted in seven domains (labeled I-VII in Table 1) likely to generate structured RNAs and 11 individual windows with z-scores less than −1.70. Concatenated domains range from 160 to 280 nt (comprising two to six overlapping windows). Interestingly, in addition to their low z-scores, the defined domains also have lower than average ΔG native (in all but domain II) and ED scores (Table 1), suggesting stable folding with one (or few) dominant conformations in the structural ensemble. Structure models for individual windows can be analyzed directly from the RNAStructuromeDB. For example, domain VII contains a window (positions chr11:65,506,081-65,506,200) that contains a known structured element important to MALAT1 maturation and that also has independent functions in the cell: the MALAT1-associated small cytoplasmic (masc)RNA 85 . The terminal window in this domain contains the predicted model of the mascRNA (highlighted in Fig. 1d), which correctly predicts the mascRNA tRNA-like folding this is essential to its function. Sequences corresponding to each longer domain can be used for structure modeling: e.g. using RNAfold locally or through the RNAfold web server 86 , be used for BLAST 87 searches to identify homologs, aligned to related sequences, and used for consensus folding, comparative sequence/structure analyses and ncRNA prediction (e.g. using the RNAz server [38][39][40]. Links to all of these tools can be found on the RNAStructuromeDB website. Example 2: The VEGFA gene. In 2008, the first human riboswitch was discovered in the mRNA of vascular endothelial growth factor-A (VEGFA) 44 . Within this mRNA, is a region of the 3′ UTR which can adopt two structural conformations, each of which leads to different translational levels of VEGFA. The adoption of either conformation is mediated by protein binding which occurs as a result of hypoxic signaling. When expressed in high abundance (under hypoxic conditions) the regulatory protein hnRNP L binds to the VEGFA riboswitch causing two sequences, known as the GAIT (gamma interferon inhibitor of translation) element and the stem stability sequence (Fig. 2d), to anneal to each other: inhibiting association with the GAIT complex and stimulating VEGFA expression. When hnRNP L is not highly expressed (under normoxic conditions), these elements form their own discrete hairpins (Fig. 2d) and the GAIT element is able to bind the GAIT complex repressing VEGFA expression. Figure 2 shows the basic pipeline involved in using the genome browser. Upon opening the JBrowse tool, a user is able to input genomic coordinates or Ensembl ID and navigate to their region of interest. In the case of VEGFA, the Ensembl ID (ENSG00000112715.21) was entered to navigate to the genome coordinates. Graphical representations of folding metrics are displayed as tracks, which flank the VEGFA transcript model (Fig. 2a). Transcript models display exons (yellow), introns (black line), and UTRs (blue) as separate entities, allowing for quick visual analyses. VEGFA has a multitude of low z-score windows (197 windows below average), the majority of which (85%) lie within introns and UTRs. The lowest z-score window in the VEGFA gene resides in the 3′ UTR; excitingly, this window fully encompasses the VEGFA riboswitch, highlighting the utility of the metrics archived in the RNAStructuromeDB as well as the ease of picking out regions of interest visually using JBrowse. Users can highlight this region to "zoom in" to investigate that particular window (Fig. 2b).
Scientific REPORtS | 7: 17269 | DOI:10.1038/s41598-017-17510-y Individual windows can be selected to view a "pop-up" window ( Fig. 2c) containing: folding metrics, RNA sequence, and a dot-bracket structure model (dot-bracket structures can be conveniently represented in 2D using the VARNA java applet 88 ). Using VARNA, we annotated the 2D image with the key riboswitch elements (Fig. 2d). The model MFE structure archived in the RNAStructuromeDB recapitulates the translation permissive (TP) conformation that dominates under hypoxia, where the GAIT element and stem stability sequence are annealed to each other. The differences in the MFE model and the TP conformation described in the literature occur primarily at the hnRNP L binding site. The literature model is based on enzymatic probing of an in vitro generated construct that terminates 5 nt upstream of the hnRNP L interaction site. The MFE model, however, includes 13 additional upstream nt because of the window size used. These additional nt allow for the formation of a short hairpin stem that sequesters part of the hnRNP L binding site in a short hairpin loop (Fig. 2d). It is worth noting that the in vitro model for the TP conformer is poorly defined in this region: high reactivity spans the 5′ end and a strong  92 . Inosine is recognized by most cellular machinery (as well as the enzymes used in RNA sequencing) as guanosine; indeed it is the A-to-G transitions observed as variants in sequencing reads that allows one to quantify the percent A-to-I editing at a nt 93 . Efforts towards creating an "inosinome Atlas" (by conducting whole genome and transcriptome sequencing of human cells to detect A-to-I editing events) resulted in the identification of ~3 million A-to-I editing sites 90 . Given their quantity, likely roles in regulating gene expression, and their implication in disease states it is important to gain insight into the structural features of A-to-I "hot-spots".
ADAR enzymes have a strong preference for double stranded (ds)RNA regions 94 . ADAR editing was first described in helical regions of very long hairpins found in 3′ UTRs 95,96 Subsequently, many editing sites were found; for example, within regions containing Alu elements 97 whose inverted repeat structure facilitates hairpin formation. This structure-preference for ADAR makes the folding metrics within the RNAStructuromeDB particularly useful for interpreting patterns of A-to-I editing data. We utilized the customizability of JBrowse to visualize A-to-I editing sites alongside folding metrics in order to examine the structural landscape corresponding to edited regions. We prepared tracks using data from a study which looked at RNA editing in human B cells 91 where several transcripts were discovered to be hyperedited. For example, the formin binding protein 1 (FNBP1) transcript was found to have the most editing events (291) of any gene, and when viewed alongside folding metrics it was clear that predicted structured regions aligned well to editing sites (Fig. 3). A particularly striking example is the window with metrics shown in Fig. 3b (as well as the genome browser tracks shown in Fig. 3c and d); this window has not only a highly-negative z-score (almost five standard deviations more stable than random), but low (favorable) MFE, ED and high fMFE values as well. This window overlaps a cluster of five highly-edited sites, which are annotated on the predicted structure shown in Fig. 3e. Interestingly, in addition to the highly-stable hairpin in the individual window, concatenating adjacent windows which overlap editing sites (or are within 40 nt and have a less than average z-score) can define a larger hairpin structure (677 nt long) where this particular window forms the terminal hairpin stem loop ( Figure S1). Nine out of the 15 windows overlapping this region had z-scores lower than the transcript average (Supplementary Table 1) showing how the sliding window approach can still be used to define domains > the window size used (120 nt in this case). In this particular region there have been two inverted Alu element insertions (Fig. 3d), which provide the complementarity for forming such a large stem structure. As this whole region is transcribed as part of an intron, it is possible for it to fold as predicted  Table 1. Low z-score windows and regions of MALAT1. Each region was defined as a low z-score region if it contained two or more consecutive windows with z-scores 1 σ lower (< −1.71) than the MALAT1 average (−0.50). Also reported is each region's minimum and average ΔG native (kcal/mol) and ensemble diversity (ED) value. Single windows with low z-scores (< −1.71) are also shown. Region names (I-VII) correspond to those shown in Fig. 1.
into the long hairpin, which is structurally similar and similar in length to the long 3′ UTR hairpins described in earlier studies of ADAR editing 95,96 .
In addition to helping define the extent of structured domains recognized by ADAR (which could facilitate the design of assays to study editing, by removing/adding domains to reporters, etc.) the information in the RNAStructuromeDB can also offer insights into the detailed structural contexts of A-to-I editing sites. The primary sequence context of these editing sites is being investigated, with insights being made into flanking nucleotide preferences 90 and tools have been created that predict editing sites based on these primary sequence rules 98 ; the RNAStructuromeDB may enhance these investigations by providing secondary structure context. For example, we can compare the percent editing at an individual site to its structural context. In Fig. 3e, three editing sites occur in AU pairs (two flanked by Watson-Crick pairs, and the other at the end of a helix upstream of a 3 nt bulge-loop) and two occur in internal loops opposite to cytosine residues. In all cases the inosine substitution would be expected to minimally disrupt folding energy (AU to IU pairs) or enhance stability (AC mismatch to IC pair). Additionally, the flanking nearest-neighbor nt can also affect stability. Thus, the folding information in the RNAStructuromeDB might be able to help discern sequence/structural features of editing hot spots that can help predict sites of editing as well as the effect of editing on RNA folding.

Conclusion
The RNAStructuromeDB is a repository of useful RNA folding metrics and a powerful vehicle for exploring the human genome via RNA structure. It allows users to browse, access, and retrieve the data quickly and flexibly, which will facilitate a wide array of researches. We presented three examples of how this database can be used: to generate a map of RNA folding throughout the MALAT1 lncRNA, VEGFA gene, and FNBP1 transcript. In each example, global properties of folding (e.g. the enrichment for stable folding in intronic or UTR regions), as well as the identification of functional motifs (e.g. the identification and modeling of the VEGFA riboswitch and stable hairpin in FNBP1) could be quickly deduced. We anticipate that this tool will have many applications both in basic research and in the therapeutic targeting of disease-associated human RNAs.

Materials and Methods
Overview. Our   and genomic data stored within the Chado schema. Folding metrics were calculated on Iowa State University's High Powered Compute cluster using RNAfold (version 2.2.10) and Perl (version 5). The website pages were built using Iowa State University's Luggage platform (http://luggagedocs.info/), which is constructed on an underlying Drupal 7 framework. Tripal 99,100 (version 2.1) was used to upload all data into the Chado schema and populate page "views" by later pulling relevant data from the Chado schema. JBrowse API (https://github.com/isubit/tri-pal_jbrowse_api) was used to generate JBrowse tracks directly from the Chado schema. An overview of this procedure is indicated in Fig. 4.
Genomic sequence and annotations files. This analytical approach utilizes FASTA chromosomal DNA sequence files as input. For the human genome, the standardized reference sequence is generated by the Human Genome Reference Consortium. The latest major release of this human genome reference sequence (GRCh38/ hg38) was acquired from Ensembl FTP server (ftp://ftp.ensembl.org/pub/release-90/fasta/homo_sapiens/dna/) as 24 separate chromosomal FASTA files (consisting of unmasked DNA sequence). Regions of sequence ambiguity and/or difficult to sequence nucleotides are depicted as the character "N" and were not considered in any calculations, but remain as placeholders within corresponding windows. Comprehensive Gencode genome annotations 54 (release 26) were acquired from Gencode directly in the form of GFF3 files (ftp://ftp.sanger.ac.uk/pub/ gencode/Gencode_human/release_26/gencode.v26.annotation.gff3.gz).
RNA structure, partition function and z-score calculations. Each chromosomal FASTA file supplied the underlying sequence for metric calculations, and was analyzed every 40 nt using 120 nt windows. Each window fragment is run through RNAfold to calculate its MFE (value and fold) and partition function at 37 °C. To account for the reverse strand, each fragment is also converted to its reverse complement and run through the same process in a parallel script. Partition function calculations output an ensemble diversity score and an fMFE in the ensemble value for each window. Each "native" window sequence is then scrambled to produce 30 randomized versions which are run through RNAfold to calculate their MFE values as well. The MFE values of the native sequences are then compared to the MFE values of random sequences to generate a thermodynamic z-score, calculated with a method adapted from Clote et al. 66 as illustrated in the following equation: The randomization number was optimized to yield converged z-score values using the minimal randomization number. As a measure of z-score quality, we also calculate the fraction of MFE values that were lower than the native (the p-value) using Perl operations. Database web development. Web-based interfaces were prepared to put folding metric output into genomic context. An example of the pipeline used to generate the RNAStructuromeDB is shown in Fig. 4. Data tables were constructed by linking genome annotations to folding metrics based on genomic coordinates (whenever such relationships existed) and made available as searchable and downloadable CSV files; a preview of the data can be found alongside download links on the "Data Search" and "Gene Search" web pages, where users can view the data from individual windows and find links to view the corresponding MFE structures in forna 101 . Folding metric windows and Gencode annotations were also prepared as feature tracks for the JBrowse genome browsing interface; each folding metric window is displayed as a feature with directionality and each gene was visualized as "processed transcripts" (in order to visualize directionality, underlying exons, and UTRs when appropriate). Additionally, each folding metric (MFE, z-score, p-value, ED, and fMFE) was extracted and prepared as separate big-Wig tracks for both forward and reverse strands -with each folding metric window trimmed to span only its initial 40 nt, allowing for proper visualization of each metric as a bar graph. We have left JBrowse customizable: users are given the option to upload tracks (in BAM, BED, bigWig, VCF, or FASTA formats) alongside RNAStructuromeDB metrics. Users are also able to download data from every track (whole sets or selected regions).

MALAT1
Example. MALAT1 metrics were obtained from https://structurome.bb.iastate.edu/folding-metricfeature-search by filtering using Ensembl ID (ENSG00000251562.7) and downloading via the CSV link in the webpage. Statistical analyses were performed and graphs/tables were generated in Excel.
VEGFA Example. VEGFA was navigated to at https://structurome.bb.iastate.edu/jbrowse/ by inputting its Ensembl ID (ENSG00000112715.21). All data was gathered from the JBrowse interface directly. All structures were generated using VARNA and annotated as described in Ray, et al.

FNBP1 Example.
A-to-I editing sites were obtained from ref. 88 supplemental table five. Genomic coordinates were converted from human genome reference version hg18 to version hg38 using the UCSC reference conversion tool: https://genome.ucsc.edu/cgi-bin/hgLiftOver. Converted coordinates were then used to create bigWig tracks for each nucleotide site. These tracks were uploaded directly to JBrowse from the user interface for visualization along with folding metrics.
Data availability. The datasets generated during the current study are available in the RNAStructuromeDB downloads repository, https://structurome.bb.iastate.edu/downloads or from corresponding author on reasonable request.