SpeciesLM training<\/p>\n
For SpeciesLM metazoa, we obtained metazoan genomes comprising 494 different species from the Ensembl 110 database36<\/a>. For each annotated protein-coding gene, we extracted 2,000 bases 5\u2032 to the start codon and trained a species-aware masked language model on this region. We followed the training and tokenization procedure outlined in Species-aware gLMs4<\/a>, but kept the batch size at 2,304, despite increasing the input sequence length, resulting in approximately twice as many tokens seen during training as in SpeciesLM fungi 5\u2032. We used rotary positional encoding to inject positional information into the Transformer blocks.<\/p>\n For SpeciesLM fungi, we deviated from the above recipe by tokenizing each base of the sequences discussed in ref. 4<\/a> separately (single nucleotide, 1-mer tokenization) and using learned absolute positional encodings. To stabilize training, we increased dropout in the multilayer perceptron layers of the transformer to 0.2 and set it to 0.1 for attention dropout.<\/p>\n Overall, we improved the training efficiency by fusing biases of the linear layers, the multilayer perceptron in the transformer and the optimizer using Nvidia Apex. We used FlashAttention2 (ref. 51<\/a>) to train all models.<\/p>\n Nucleotide dependencies and variant influence score<\/p>\n We define the dependency between a variant nucleotide kalt at position i and a target position j as<\/p>\n $${e}_{i,j,{k}_{\\mathrm{alt}}}=\\max {\\left\\{\\left|{\\log }_{2}\\left(\\frac{{\\rm{o}}\\hat{{\\rm{d}}}\\mathrm{ds}\\left({n}_{\\!j}={k|}{n}_{1},\\ldots ,{n}_{i}={k}_{\\mathrm{alt}},\\ldots ,{n}_{N}\\right)}{{\\rm{o}}\\hat{{\\rm{d}}}\\mathrm{ds}\\left({n}_{\\!j}={k|}{n}_{1},\\ldots ,{n}_{i}={k}_{\\mathrm{ref}},\\ldots ,{n}_{N}\\right)}\\right)\\right|\\right\\}}_{k{\\rm{\\in }}\\left\\{{\\rm{A}},{\\rm{C}},{\\rm{G}},{\\rm{T}}\\right\\}}$$<\/p>\n where k is one of the four possible nucleotides A, C, G or T; ni and nj are the nucleotides at position i and j, respectively; kref is the nucleotide in the reference, nonaltered input sequence, and kalt is the nucleotide in the alternative sequence. The odds estimates are computed from the predictions of a gLM under consideration. For this computation, none of the nucleotides (including the target nucleotide) is masked.<\/p>\n The variant influence score \\({{e}_{i,k}}_{\\mathrm{alt}}\\), for a sequence of N nucleotides, is defined by averaging the dependencies on a variant nucleotide at position i across all positions j = 1, \u2026, N such that \\(j\\ne i\\).<\/p>\n A nucleotide dependency ei,j between a query position i and a target position j on a sequence of N nucleotides is given by:<\/p>\n $${e}_{i,j}=\\max {\\left\\{\\left|{\\log }_{2}\\left(\\frac{{\\rm{o}}\\hat{{\\rm{d}}}\\mathrm{ds}\\left({n}_{\\!j}={k|}{n}_{1},\\ldots ,{n}_{i}{\\rm{\\ne }}{k}_{\\mathrm{ref}},\\ldots ,{n}_{N}\\right)}{{\\rm{o}}\\hat{{\\rm{d}}}\\mathrm{ds}\\left({n}_{\\!j}={k|}{n}_{1},\\ldots ,{n}_{i}={k}_{\\mathrm{ref}},\\ldots ,{n}_{N}\\right)}\\right)\\right|\\right\\}}_{k{\\rm{\\in }}\\left\\{{\\rm{A}},{\\rm{C}},{\\rm{G}},{\\rm{T}}\\right\\}}$$<\/p>\n We compute dependencies for all i,j pairs such that \\(i\\ne j,\\) that is we do not consider self-dependencies.<\/p>\n In autoregressive models, a query variant cannot directly affect the prediction of a target position located 5\u2032 of the query. Thus, to obtain the lower triangular matrix of the dependency map, we also run the model on the reverse strand.<\/p>\n In the SpeciesLM metazoa, which predicts nucleotides as overlapping 6-mers, the procedure needs to be adapted to yield one prediction for each target nucleotide. This is achieved by first computing for each of the six 6-mers that overlap the target nucleotide of interest, which probability it implies for this target nucleotide, as previously described4<\/a>. We then average these six probabilities to obtain a single probability.<\/p>\n For the Nucleotide Transformer models, which predict only nonoverlapping 6-mers, we use a similar approach. Consider the case of predicting the probability of observing nucleotide n at position i of the sequence. In the tokenized sequence, this nucleotide has position p in the kth 6-mer where:<\/p>\n $$k=\\left\\lfloor \\frac{i}{6}\\right\\rfloor$$<\/p>\n $$\\begin{array}{cc}p=i & \\mathrm{mod}6\\end{array}$$<\/p>\n The model predicts a distribution over all 46 possible 6-mers at position k. We first discard all predictions corresponding to 6-mers that contain a nucleotide that differs from the reference sequence at any location other than p\u2014which leaves only four 6-mers. We renormalize so that the predicted probability of these remaining 6-mers sums to one. We then record the (renormalized) probability of the 6-mer that has the desired nucleotide n at position p.<\/p>\n Apart from extracting nucleotide-level probabilities with the above-mentioned method, we have also experimented with computing the probability for a nucleotide at position i as the sum of all k-mers containing that nucleotide at that position. Evaluation of nucleotide dependencies within tRNAs revealed a worse performance with this method.<\/p>\n Variant impact benchmarks<\/p>\n As our metric of variant impact, we used the variant influence score. This average is computed over the full receptive field of the model for the SpeciesLM. For Nucleotide Transformer models, we only average over the central 2\u2009kb, so as to facilitate comparisons. Nevertheless, we provide the full sequence context for which this model has been trained.<\/p>\n For comparison, we also calculated a variant effect score based on the gLM reconstruction at the query variant. Specifically, this score is the log ratio between the predicted probability of the variant nucleotide and the predicted probability of the reference nucleotide5<\/a>,6<\/a>.<\/p>\n Finally, we downloaded conservation scores (PhyloP and PhastCons) for human and S. cerevisiae from the University of California, Santa Cruz genome browser database21<\/a>,22<\/a>,23<\/a>,24<\/a>,52<\/a>. For humans, these include the conservation scores based on the 100-way, 447-way and 470-way alignment.<\/p>\n Promoter saturation mutagenesis<\/p>\n Promoter saturation mutagenesis (ref. 20<\/a>) data mapped to hg38 were provided by V. Agarwal (mRNA Center of Excellence, Sanofi, Waltham, MA, USA). As discussed in ref. 29<\/a>, we excluded the FOXE1 promoter due to the low replicability of the measurements, leaving nine promoters and comprising 8,635 variants. Variants were then intersected with the human gene 5\u2032 regions (that is, the regions 2-kb 5\u2032 of annotated start codons). Then, the variant influence score was calculated for each variant measured in the assay from the LM dependencies for these regions. The variant influence score was then correlated with the absolute value of the measured log2 fold change in expression. This correlation was computed for each promoter and then averaged across promoters.<\/p>\n To determine confidence intervals, we performed 100 bootstrap samples per promoter and recomputed the correlation for each bootstrap sample. The confidence interval was defined by adding\/subtracting two standard deviations of the average correlation.<\/p>\n eQTL variants<\/p>\n For human eQTL, we downloaded SUSIE26<\/a> fine-mapped GTEx eQTL data from EBI. We then intersected these data with the human gene 5\u2032 regions. This procedure, by design, enriches for promoter eQTL. Similar to the details in ref. 29<\/a>, we considered every eQTL variant with a posterior inclusion probability higher than 0.9 as putative causal and we considered any eQTL variant with posterior inclusion probability lower than 0.01 as putative noncausal. We only considered putative noncausal eQTL intersecting regions, which also include at least one causal eQTL. This procedure gave 2,958 eQTL variants, of which 1,631 were classified as putative causal. Then, the influence score for each variant was computed based on the nucleotide dependencies in these regions. We ranked variants according to the influence score. Confidence intervals were computed using bootstrapping as before.<\/p>\n For yeast eQTL, we downloaded the results of an MPRA study assessing candidate cis-eQTL variants27<\/a>. After this study, we classify any eQTL variant with false discovery rate\u2009<\u20090.05 in the MPRA assay as causal and we classify any eQTL with (unadjusted) P value of >0.2 as noncausal. This yielded 3,056 eQTL variants, of which 379 were classified as causal. These eQTL variants were then intersected with yeast gene 5\u2032 regions and influence scores were computed from the SpeciesLM fungi dependency maps. Confidence intervals were computed using bootstrapping as before.<\/p>\n Clinvar<\/p>\n We used ClinVar version 2023_07_17 (ref. 19<\/a>), previously downloaded from https:\/\/ftp.ncbi.nlm.nih.gov\/pub\/clinvar\/vcf_GRCh38\/<\/a>. We considered noncoding any variant in the categories \u2018intron_variant\u2019, \u20185_prime_UTR_variant\u2019, \u2018splice_acceptor_variant\u2019, \u2018splice_donor_variant\u2019, \u20183_prime_UTR_variant\u2019, \u2018non_coding_transcript_variant\u2019, \u2018genic_upstream_transcript_variant\u2019 and \u2018genic_downstream_transcript_variant\u2019. As discussed in ref. 53<\/a>, we considered as pathogenic any variant classified as pathogenic or likely pathogenic and as benign any variant classified as benign or likely benign. We excluded variants with fewer than one review star. This resulted in 385,572 variants, of which 22,313 were classified as pathogenic.<\/p>\n As most ClinVar variants fall outside the 5\u2032 regions of genes, we chose not to intersect with these regions. Instead, we computed the dependency map centered on the variant of interest. Confidence intervals were computed using bootstrapping as before.<\/p>\n Borzoi<\/p>\n We ran Borzoi in mixed precision to reduce computational overhead using the PyTorch Borzoi package. Replicate zero of Borzoi was used for all analyses. For the eQTL analysis, we computed the L2 score as discussed in ref. 25<\/a>. We used the tissues of borzoi predictions matching the eQTLs. If several Borzoi tracks matched the tissue, we averaged the scores across these tracks. For ClinVar, we followed a similar approach, except that we collected Borzoi predictions for all tissues and assays. We then computed the L2 score across tracks to give a tissue-agnostic and mechanism-agnostic variant-effect score.<\/p>\n For the Kircher saturation mutagenesis dataset, we computed the logSED score as discussed in ref. 25<\/a>. We mapped the cell types used in the assay to Borzoi tracks as follows: for the GP1BB, HBB, NHBG1 and PKLR promoters, we used \u2018RNA:K562\u2019; for the F9 and LDLR promoters, we used \u2018RNA:HepG2\u2019; for HNF4A and MSMB, we used \u2018RNA:kidney\u2019 (as HEK293 is originally a kidney cell) and for TERT, we used \u2018RNA:astrocyte\u2019 (as glioblastoma are cancerous astrocytes).<\/p>\n Integrative model using Borzoi and the influence score<\/p>\n We integrated Borzoi and the influence score using logistic regression\u2014for the eQTL and ClinVar predicitions\u2014and using linear regression for the mutagenesis data using fivefold cross-validation scheme for all benchmarks. Notably, for the Kircher saturation mutagenesis task, model fitting and cross-validation were performed separately for each promoter, and performance was then averaged across folds and promoters.<\/p>\n Alternative dependency metrics<\/p>\n All benchmarks on alternative dependency metrics were performed on the SpeciesLM fungi.<\/p>\n Gradient-based<\/p>\n We computed the gradient of the prediction for each nucleotide at position i with respect to each nucleotide at position j yielding a 4\u2009\u00d7\u20094 matrix. To achieve this, we first replaced the tokenization layer with a one-hot encoding and a linear layer, which map the one-hot encoded nucleotides to their respective token embeddings. We then propagated gradients from each target nucleotide prediction to each one-hot encoded input nucleotide. As a metric of nucleotide dependency, we then used the maximum absolute value across the 4\u2009\u00d7\u20094 matrix of each i,j position.<\/p>\n Mask-based<\/p>\n Masked-based dependencies are computed as:<\/p>\n $${e}_{i,j}=\\max {\\left\\{\\left|{\\log }_{2}\\left(\\frac{{\\rm{o}}\\hat{{\\rm{d}}}\\mathrm{ds}\\left({n}_{\\!j}={k|}{n}_{1},\\ldots ,{n}_{i}=\\left[\\mathrm{MASK}\\right],\\ldots ,{n}_{N}\\right)}{{\\rm{o}}\\hat{{\\rm{d}}}\\mathrm{ds}\\left({n}_{\\!j}={k|}{n}_{1},\\ldots ,{n}_{i}={k}_{\\mathrm{ref}},\\ldots ,{n}_{N}\\right)}\\right)\\right|\\right\\}}_{k{\\rm{\\in }}\\left\\{{\\rm{A}},{\\rm{C}},{\\rm{G}},{\\rm{T}}\\right\\}}$$<\/p>\n where \u2018[MASK]\u2019 stands for the mask token, k belongs to one of the four possible nucleotides A, C, G or T; ni and nj are the nucleotides at position i and j, respectively; kref is the nucleotide in the reference, nonaltered input sequence.<\/p>\n S. cerevisiae tRNA structure benchmark<\/p>\n S. cerevisiae genome assembly version R64-1-1 and annotation version R64-1-1.53 were downloaded from EnsemblFungi36<\/a>. The S. cerevisiae tRNA secondary structures were downloaded from GtRNAdb54<\/a>. We considered only the tRNAs overlapping the 1\u2009kb 5\u2032 regions to any yeast start codon, yielding 172 tRNA sequences. Subsequently, dependency maps on tRNAs were processed by taking the maximum between ei,j and ej,i. This symmetrizes the dependency map and achieves one unique score per pair of positions in the tRNA sequence. We then used this score to predict whether a pair of nucleotides belonged to a secondary structure contact.<\/p>\n Assessment of donor\u2013acceptor dependencies in S. cerevisiae<\/p>\n We extracted intron sequences by selecting the regions within annotated gene intervals that lie between exon annotations. This resulted in 380 sequences. We then retained only introns bounded by canonical splice site dinucleotides GT and AG, yielding 272 sequences. We then computed the average dependency between every donor and acceptor nucleotide within the intron as a measure of dependency between the donor and acceptor sites. We designed two negative sets for a given intron. For the negative set \u2018Decoy acceptor\u2019, we compute the average dependency between donor nucleotides and each AG dinucleotide within the intron that does not include the acceptor site. For the negative set \u2018Matched distance\u2019, we sampled four random dependencies between nucleotides that were as distant from each other as the donor was from the acceptor, without including the donor or acceptor themselves.<\/p>\n TF motif mapping<\/p>\n We downloaded FIMO PWM scan results from http:\/\/www.yeastss.org<\/a> (ref. 55<\/a>) and Chip-Exo TF binding peaks from http:\/\/www.yeastepigenome.org<\/a> (ref. 31<\/a>). We then extracted all Chip-Exo peaks for the available PWMs. We excluded PWM matches for which no Chip-Exo data were available for the corresponding factor. This procedure yielded data for 68 TFs. We annotated every nucleotide within 1\u2009kb 5\u2032 of a start codon as part of a binding TF motif if it is (1) part of a PWM match with P value of <0.01 and (2) this PWM match is within ten bases of a Chip-Exo peak of the corresponding TF. We defined the positive class in this way to ensure that we capture nucleotides relevant for determining binding (that is, motif) rather than all nucleotides close to a Chip-Exo peak, regardless of their role in binding. This resulted in 92,117 binding nucleotides out of a total of 6,538,427. We designated a nucleotide as repeat if it was masked by RepeatMasker. We extracted this information from the soft-masked GTF provided by Ensembl36<\/a>.<\/p>\n The 69-way alignment<\/p>\n We used progressive cactus56<\/a> to align 69 budding yeast species using default parameters and specifying S. cerevisiae as reference quality genome. We then extracted fourfold degenerate sites and used phyloFit with the EM algorithm to estimate a neutral model. Using this neutral model and the alignment, we ran phastCons with the parameters –rho 0.3 –estimate-rho –target-coverage 0.4 –expected-length 23, which correspond to the parameters used in the seven-way alignment21<\/a>. We also ran phyloP22<\/a>, with –method LRT –mode CONACC.<\/p>\n Dependencies in rare-variant-associated aberrant splicing<\/p>\n We computed dependency maps for all rare SNVs associated with splicing outliers in GTEx28<\/a> as described earlier32<\/a>. Because the input length of SpliceBert14<\/a> is limited to 1,024\u2009bp, the complete set of variant outlier pairs (n\u2009=\u200918,371) was filtered such that the variant and associated outlier junction were located within an 800-bp window (n\u2009=\u20091,811) and 100\u2009bp of sequence was added from the maximum and minimum positions of the variant and outlier junction splice sites. For each variant location, we extracted the average value of the dependency map at the intersection of either variant and outlier donor dinucleotide or variant and outlier acceptor dinucleotide. This variant effect score was compared against a background score. This background score was computed as the mean over all dependencies that were as distant from each other as the variant was from the outlier donor (matched distance) or the outlier acceptor. The scores were filtered for a minimum distance of 5\u2009bp between the variant and splicing dinucleotide to filter values near the diagonal corresponding to self-interactions. Variant categories were annotated with the Ensembl variant effect predictor (VEP)57<\/a>. For each variant, the most severe VEP annotation was considered. For the \u2018exon\u2019 category, the following VEP categories were grouped together: synonymous_variant, missense_variant, stop_lost, stop_gained.<\/p>\n Genome-wide search for parallel and antiparallel dependencies<\/p>\n We scanned dependency maps for parallel and antiparallel dependencies using 5\u2009\u00d7\u20095 convolutional filters. We constructed the antiparallel filter by populating the antidiagonal of a zero-filled 5\u2009\u00d7\u20095 matrix with ones, and for the parallel filter, by populating the diagonal with ones. We then centered each filter by subtracting the mean value from each position to ensure that a convolution on a uniform 5\u2009\u00d7\u20095 region yields a result of zero. We applied these filters to dependency maps from SpeciesLM fungi (both filters) and RiNALMo (antiparallel filter only)4<\/a>,7<\/a>.<\/p>\n Search for parallel and antiparallel dependencies in fungi using the SpeciesLM fungi<\/p>\n For the SpeciesLM fungi, we have computed dependency maps spanning 1\u2009kb 5\u2032 of each annotated start codon on a set of representative fungi species, including Agaricus bisporus, Candida albicans, Debaryomyces hansenii, Kluyveromyces lactis, Neurospora crassa, S. cerevisiae, S. pombe and Yarrowia lipolytica. The genomes and annotation files for each species were downloaded from EnsemblFungi release 53 with accessions GCA_000300555.1<\/a>, GCA000182965v3<\/a>, GCA_000006445.2<\/a>, GCA000002515.1<\/a>, GCA_000182925.2<\/a>, GCA_003046715.1<\/a>, GCA_000002945.2<\/a> and GCA_000002525.1<\/a>, respectively.<\/p>\n All regions annotated as \u2018five_prime_utr\u2019, \u2018three_prime_utr\u2019, \u2018intron\u2019, \u2018CDS\u2019, \u2018pseudogene_with_CDS\u2019 and other regions (for example, nonannotated introns) inside an annotated gene interval were categorized as protein-coding gene. All regions annotated as \u2018tRNA\u2019, \u2018tRNA_pseudogene\u2019, \u2018rRNA\u2019, \u2018snRNA\u2019, \u2018ribozyme\u2019, \u2018SRP_RNA\u2019, \u2018snoRNA\u2019, \u2018RNase_P_RNA\u2019 and \u2018RNase_MRP_RNA\u2019 were categorized as structured RNA. Finally, all regions annotated as \u2018transposable_element\u2019, \u2018pseudogene\u2019 and regions without any annotation were considered as intergenic.<\/p>\n Search for antiparallel dependencies and RNA structure in E. coli using RiNALMo<\/p>\n For RiNALMo, we computed dependency maps for regions 100, 200 and 500\u2009bp before each annotated start codon in E. coli str. K-12 substr. MG1655, whose genome and annotation were downloaded from GenBank58<\/a> with accession U00096.3<\/a>.<\/p>\n As candidates for a new RNA structure, we first filtered positions whose convolution value is greater or equal to 25 to select only high-value antiparallel dependencies, resulting in a filtered convolved dependency map. Next, we counted the unique number of antidiagonals potentially belonging to one stem by extracting the unique i + j nonzero positions supported by at least three nonzero values.<\/p>\n As candidates for a new structure, we selected maps suggesting the existence of at least two potential stems.<\/p>\n RNA secondary structure benchmarking<\/p>\n We downloaded the database of secondary structures Archive II37<\/a>, which includes 3,865 curated RNA structures across nine families (5S rRNA, SRP RNA, tRNA, tmRNA, RNase P RNA, group I intron, 16S rRNA, telomerase RNA and 23S rRNA). For each structure, we generated the dependency map with the pretrained RiNALMo and retained the largest of the two dependency map entries for each pair of nucleotides (maximum of i,j and j,i). The AUROC curve was computed for each structure against the Archive II secondary structure annotations.<\/p>\n Benchmarking of canonical and noncanonical RNA contacts<\/p>\n We downloaded the database of RNA structures CompaRNA38<\/a>, which is a compilation of RNA contacts based on 201 available RNA structures in the Protein Data Bank by RNAView59<\/a>. Contacts are classified either as \u2018standard\u2019 or as \u2018extended\u2019. While the first includes only canonical AU, GC and wobble GU pairs in the cis-Watson\u2013Crick\/Watson\u2013Crick conformation60<\/a>, the latter calls all interacting bases regardless of their conformation, including noncanonical or tertiary contacts. Of the 201 structures, 196 had a length below the maximum input length of RiNALMo (1,022 nt). For each structure, we generated the dependency map using the pretrained RiNALMo and retained the largest entry from the two dependency maps for each pair of nucleotides. Similarly, the same structures were also evaluated with the fine-tuned RiNALMo model version rinalmo_giga_ss_bprna_ft, resulting in a predicted value for each pair of nucleotides. To evaluate their performance in predicting noncanonical contacts, we excluded all canonical contacts and computed the AUROC curve for all remaining positions across all structures. Significance between ROC AUCs was determined by bootstrapping over 10,000 permutations.<\/p>\n Comparison with RNAalifold<\/p>\n We evaluated the performance of the dependency maps against RNAalifold61<\/a>, a standard alignment-based method for predicting a consensus RNA structure by incorporating sequence covariation from a set of aligned RNA sequences as input. For this, we use the 201 PDB entries in CompaRNA that had at least one Rfam match and consider two subsets. The first subset consisted of the 33 PDB sequences that contained an exact sequence match between the PDB entry and at least one Rfam seed alignment. In case of multiple matching Rfam seed alignments (for example, ribosomal RNA), we considered an arbitrarily chosen single Rfam seed alignment to avoid confounding the evaluations by duplicates. The second subset consisted of the remaining 168 sequences. For this, we used nhmmer (v3.1b2)62<\/a> to find homologous sequences within a database of 220,478 bacterial and archaeal genomes and plasmids downloaded from NCBI. After removing sequences longer than 1,022 nt (the maximum context length for which the gLM RiNALMo has been trained), this resulted in 67 sequences with hits in the database.<\/p>\n On the first subset, we use the Rfam seed alignments as input to RNAalifold. To assess the robustness of the analyses to the alignment procedure, we additionally realigned the sequences in the seed alignments using Clustal-Omega (v1.2.4)63<\/a> and MAFFT (v7.525)64<\/a>. For the second subset, we performed sequence alignments using both Clustal-Omega and MAFFT, limiting the alignments to a maximum of 1,000 sequences (by aligning the PDB sequence to the top 999 nhmmer hits) to reduce computation time. On both subsets and from each alignment, a base-pair probability matrix corresponding to the predicted RNA structure was generated using RNAalifold available through the ViennaRNA (v2.6.4) package65<\/a>. RNAalifold was run in the following two modes: using the default energy model (command: RNAalifold -p) and with RIBOSUM scoring (command: RNAalifold -p -r).<\/p>\n Pseudoknot benchmark<\/p>\n We downloaded the compendium dataset bpRNA-1m(90) that contains 28,370 annotated RNA structures with less than 90% sequence similarity obtained from the databases CRW, tmRNA, SRP, tRNAdb2009, RNP, RFAM and PDB39<\/a>. From these, we extracted all structures that contain pseudoknot contacts and are no longer than 1,022 nt (the maximum context length for which the gLM RiNALMo has been trained). These resulted in 2,530 structures of varying lengths and sources. We then extracted the pseudoknot contacts from the dot-bracket notation provided by bpRNA that takes into account non-nested pairs39<\/a>. Finally, we computed the dependency maps for each one of these structures and evaluated their ability to predict whether a pair of nucleotides belongs to a pseudoknot contact (positive set) or does not belong to a structure contact (pseudoknot or canonical structure contact\u2014negative set).<\/p>\n DMS-MaPseq analysis of E. coli cells<\/p>\n E. coli TOP10 cells were grown in LB broth at 37\u2009\u00b0C with shaking until OD600\u2009=\u20090.5, after which dimethyl sulfate (DMS; Sigma-Aldrich, D186309), prediluted 1:4 in ethanol, was added to a final concentration of 200\u2009mM. Bacteria were incubated for 2\u2009min at 37\u2009\u00b0C, and reaction was quenched by addition of 0.5\u2009M final DTT. Bacteria were pelleted by centrifugation at 17,000g for 1\u2009min at 4\u2009\u00b0C, after which they were resuspended in cell pellets in 12.5-\u03bcl resuspension buffer (20\u2009mM Tris\u2013HCl pH 8.0; 80\u2009mM NaCl; 10\u2009mM EDTA pH 8.0), supplemented with 100\u2009\u03bcg\u2009ml\u22121 final lysozyme (L6876, Merck) and 20 U SUPERase\u00b7In RNase Inhibitor (Thermo Fisher Scientific, A2696), by vortexing. After 1\u2009min, 12.5-\u03bcl lysis buffer (0.5% Tween-20; 0.4% sodium deoxycholate; 2\u2009M NaCl; 10\u2009mM EDTA) were added, and samples were incubated at room temperature for 2 additional min. Then 1\u2009ml TRIzol Reagent (Thermo Fisher Scientific, 15596018) was added, and RNA extracted as per the manufacturer\u2019s instructions. rRNA depletion was performed on 1\u2009\u03bcg total RNA using the RiboCop for Bacteria kit (Lexogen, 126). DMS-MaPseq library preparation was performed as previously described41<\/a>. After sequencing, reads were aligned to the E. coli str. K-12 substr. MG1655 genome (GenBank, U00096.3<\/a>), using the rf-map module of the RNA framework66<\/a> and Bowtie2 (ref. 67<\/a>). Count of DMS-induced mutations and coverage and reactivity normalization were performed using the rf-count-genome and rf-norm modules of the RNA framework. Experimentally informed structure modeling was performed using the rf-fold module of the RNA framework and ViennaRNA (v2.5.1)67<\/a>.<\/p>\n RNA structure covariation analysis<\/p>\n Covariation analysis was performed using the cm-builder pipeline (https:\/\/github.com\/dincarnato\/labtools<\/a>) and a nonredundant database of 7,598 representative archaeal and bacterial genomes (and associated plasmids, when present) from RefSeq68<\/a>.<\/p>\n Evaluation of artificial forward and inverted duplications<\/p>\n We generated random sequences of 100 nucleotides by sampling from regions 1\u2009kb 5\u2032 of the start codon in S. cerevisiae to ensure a representative GC content and shuffling the sequences to destroy potential functional elements. Additionally, we created 100 unique duplicated sequences, ranging from 2 to 20 nucleotides in length, by randomly sampling each nucleotide with equal probability. Each duplicated sequence was then inserted into a uniquely generated 100-nucleotide sequence at a random distance from each other, ensuring no overlaps occurred. We used the SpeciesLM fungi to generate dependency maps for each sequence. We then computed average dependencies by taking the mean of the dependencies between nucleotides and their duplicates. This involved averaging across a parallel diagonal for forward duplications and an antiparallel diagonal for inverted duplications.<\/p>\n For tRNA-sized sequences, we followed a similar method but generated each sequence by shuffling each unique tRNA sequence in S. cerevisiae once. We computed the average number of inverted duplications by averaging the occurrences of duplicated sequences of specific lengths across 10,000 shuffled versions of each tRNA sequence.<\/p>\n Genome-wide analysis of dependency distribution<\/p>\n Using the SpeciesLM fungi, we computed dependency maps across the genomes of S. cerevisiae and S. pombe. Because the SpeciesLM fungi was pretrained on sequences of 1,003 nucleotides, including the start codon at the end, we discarded dependencies involving the last three nucleotides of each sequence, yielding dependencies for 1,000 nucleotides. Genome-wide dependency maps of 1-kb span were obtained with a tiling approach. Along each chromosome, we computed 1-kb square dependency maps every 500\u2009bp and averaged overlapping entries.<\/p>\n To ensure that the same number of targets is computed before and after a specific query nucleotide, we considered dependencies involving nucleotides at most 500 positions away from each other. For each map, we sampled 1,000 dependencies and averaged dependencies mapping to the same genomic positions but computed from different overlapping maps. Due to limitations in numerical precision, we considered only dependencies larger than 0.001.<\/p>\n To compute the power\u2013law coefficients, a linear regression was fitted to predict the logarithm of the dependency from the logarithm of its corresponding distance in nucleotides. The scaling coefficient was then obtained by exponentiating the fitted intercept of the linear regression, and the decay rate was obtained directly from the fitted slope. The scaling coefficient and decay rate were computed for different regions in the genome which are as follows: (1) nuclear\u2014involving all dependencies belonging to nuclear DNA; (2) mitochondria\u2014involving all dependencies within mitochondrial DNA; (3) structured RNA\u2014belonging to the annotations \u2018tRNA\u2019, \u2018tRNA_pseudogene\u2019, \u2018rRNA\u2019, \u2018snRNA\u2019, \u2018ribozyme\u2019, \u2018SRP_RNA\u2019, \u2018snoRNA\u2019, \u2018RNase_P_RNA\u2019 or \u2018RNase_MRP_RNA\u2019; (4) protein-coding gene\u2014belonging to the annotations \u2018five_prime_utr\u2019, \u2018three_prime_utr\u2019, \u2018CDS\u2019 or \u2018pseudogene_with_CDS\u2019; (5) intron\u2014belonging to the regions inside an annotated gene interval but not to exons and (6) intergenic\u2014belonging to all regions annotated as \u2018transposable_element\u2019, \u2018pseudogene\u2019, as well as regions without any annotation.<\/p>\n Model comparison<\/p>\n All other models used were downloaded from Huggingface or from their publicly available repositories. Human tRNA sequences were downloaded from GtRNAdb54<\/a>. Exact duplicate sequences were removed, leaving 266 tRNAs.<\/p>\n Reporting summary<\/p>\n